[R-sig-ME] How do I interpret linear mixed model contrast estimates from multcomp::glht()?
Lenth, Russell V
russell-lenth at uiowa.edu
Wed Dec 17 19:29:11 CET 2014
What concerns me about this question is that I don't see it expressed what contrasts you WANT. Rather, the focus is on discovering what contrasts are being produced. Surely you have some particular objectives in mind.
I don't know much about the contrast package, but this sort of thing is easy to do in lsmeans. With the model you fitted (by the way there is a warning message when fitting it), suppose you want to see predictions for each 'female' (I'd have named the variable 'sex') for each of the years, when size is 380:
R> library(lsmeans)
R> female.year <- lsmeans(math.lmm, ~ female | year,
+ at = list(year = c(0.5, 1.5, 2.5), size = 380))
Loading required namespace: pbkrtest
R> female.year
year = 0.5:
female lsmean SE df lower.CL upper.CL
Female -0.3177517 0.08298837 72.58 -0.4831633 -0.1523400
Male -0.3224042 0.08328416 73.30 -0.4883778 -0.1564306
year = 1.5:
female lsmean SE df lower.CL upper.CL
Female 0.4508660 0.08363388 74.80 0.2842516 0.6174805
Male 0.4462135 0.08390500 75.43 0.2790818 0.6133452
year = 2.5:
female lsmean SE df lower.CL upper.CL
Female 1.2194837 0.08510503 80.13 1.0501235 1.3888440
Male 1.2148312 0.08534946 80.67 1.0450019 1.3846605
Confidence level used: 0.95
### Then if you want to compare the sexes in each year...
R> pairs(female.year)
year = 0.5:
contrast estimate SE df t.ratio p.value
Female - Male 0.004652517 0.04245527 1675.57 0.11 0.9128
year = 1.5:
contrast estimate SE df t.ratio p.value
Female - Male 0.004652517 0.04245527 1675.57 0.11 0.9128
year = 2.5:
contrast estimate SE df t.ratio p.value
Female - Male 0.004652517 0.04245527 1675.57 0.11 0.9128
(same result each year since year doesn't interact with female)
### Or if you prefer using glht to do the testing:
R> summary(as.glht(female.year))
Note: df set to 76
$`year = 0.5`
Simultaneous Tests for General Linear Hypotheses
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
Female == 0 -0.31775 0.08299 -3.829 0.000440
Male == 0 -0.32240 0.08328 -3.871 0.000381
(Adjusted p values reported -- single-step method)
$`year = 1.5`
Simultaneous Tests for General Linear Hypotheses
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
Female == 0 0.45087 0.08363 5.391 1.36e-06
Male == 0 0.44621 0.08390 5.318 1.81e-06
(Adjusted p values reported -- single-step method)
$`year = 2.5`
Simultaneous Tests for General Linear Hypotheses
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
Female == 0 1.21948 0.08511 14.33 <1e-10
Male == 0 1.21483 0.08535 14.23 <1e-10
(Adjusted p values reported -- single-step method)
### Or if you want to see the slope of the 'size' trend for each of those years:
R> lstrends(math.lmm, ~ year, var = "size",
+ at = list(year = c(0.5,1.5,2.5), size = 380))
year size.trend SE df lower.CL upper.CL
0.5 -0.0003041466 0.0001910213 57.49 -0.0006865890 7.829584e-05
1.5 -0.0003648159 0.0001921332 58.82 -0.0007492979 1.966618e-05
2.5 -0.0004254852 0.0001947970 62.12 -0.0008148634 -3.610690e-05
Results are averaged over the levels of: female
Confidence level used: 0.95
At any rate, the main thing is to figure out what contrasts you want first.
Russ
Russell V. Lenth - Professor Emeritus
Department of Statistics and Actuarial Science
The University of Iowa - Iowa City, IA 52242 USA
Voice (319)335-0712 (Dept. office) - FAX (319)335-3017
Message: 2
Date: Tue, 16 Dec 2014 22:24:57 -0700
From: Matthew Van Scoyoc <scoyoc at gmail.com>
To: Ken Beath <ken.beath at mq.edu.au>
Cc: "r-sig-mixed-models at r-project.org"
<r-sig-mixed-models at r-project.org>
Subject: Re: [R-sig-ME] How do I interpret linear mixed model contrast
estimates from multcomp::glht()?
Message-ID:
<CALx9ERU6pvc79a-aZbE8B=rOzn6h650HyftE989Q=CDPq2DUFw at mail.gmail.com>
Content-Type: text/plain; charset="UTF-8"
Thanks Ken and Thierry,
The rows in summary.glht() correspond to the row names of the contrast matrix given in the linfct argument as Thierry stated in my post on r-sig-ecology <http://r-sig-ecology.471788.n2.nabble.com/How-do-I-interpret-linear-mixed-model-contrast-estimates-from-multcomp-glht-td7579236.html#a7579237>
(thanks
again Thierry). I'm using the contrast function in the contrast package to help compute the contrast matrix. The contrast function from the contrast package produces this...
> # Calculate the contrasts
> cc<-contrast(math.lm, a = list(year = c(.5, 1.5, 2.5), size = 380,
> female
= levels(egsingle$female)),
+ b = list(year = c(.5, 1.5, 2.5), size = 800,
+ female
= levels(egsingle$female)))
> cc
lm model parameter contrast
Contrast S.E. Lower Upper
t df Pr(>|t|)
0.1671791 0.01739038 0.1330889 0.2012694 9.61 7225 0
0.2148280 0.02235506 0.1710055 0.2586504 9.61 7225 0
0.2624768 0.03167964 0.2003754 0.3245781 8.29 7225 0
0.1671791 0.01739038 0.1330889 0.2012694 9.61 7225 0
0.2148280 0.02235506 0.1710055 0.2586504 9.61 7225 0
0.2624768 0.03167964 0.2003754 0.3245781 8.29 7225 0
So, how do the rows of the contrast object relate to the coefficients of the mixed or linear model? The summary of *math.lm* produces 5 rows of coefficients...
> summary(math.lm)
Call:
lm(formula = math ~ year * size + female, data = egsingle)
Residuals:
Min 1Q Median 3Q Max
-3.5631 -0.7775 -0.0523 0.7183 5.0195
Coefficients:
Estimate Std. Error t value
Pr(>|t|)
(Intercept) -5.590e-01 3.687e-02 -15.163 < 2e-16 ***
year 8.521e-01 2.391e-02 35.644 < 2e-16 ***
size -3.413e-04 4.251e-05 -8.030 1.13e-15 ***
femaleFemale -2.232e-02 2.576e-02 -0.867 0.38624
year:size -1.135e-04 2.948e-05 -3.849 0.00012 ***
I keep looking at cc$X (what I'm passing to glht() as the contrast matrix in the linfct argument), but my brain starts to melt at this point...
> cc$X
(Intercept) year size femaleFemale year:size
1 0 0 -420 0 -210
2 0 0 -420 0 -630
3 0 0 -420 0 -1050
4 0 0 -420 0 -210
5 0 0 -420 0 -630
6 0 0 -420 0 -1050
attr(,"assign")
[1] 0 1 2 3 4
attr(,"contrasts")
attr(,"contrasts")$female
[1] "contr.SAS"
This is my first try at mixed models, and it's been several years since my linear regression course that was taught in SAS, so any help is appreciated.
Cheers,
MVS
=====
Matthew Van Scoyoc
<mvanscoyoc at aggiemail.usu.edu>https://sites.google.com/site/scoyoc/
=====
Think SNOW!
More information about the R-sig-mixed-models
mailing list