I have used lme to generate a mixed effect model of the response of cells to a certain stimulus. The stimulus is applied 3 times in a row (coded by the Exposure factor that can be 1, 2, or 3) and the Genotype factor has 3 levels (wt, A, and B)
My model looks like:
model <- lme(AUC ~ Exposure + Sex + Genotype + Exposure:Sex + Exposure:Genotype + Basal, random = ~ 1 | Date / Experiment / Cell, data = mydata) anova returns significant interactions between the exposure number and both sex and genotype
numDF denDF F-value p-value (Intercept) 1 168 55.1581 <.0001 Exposure 2 168 10.1711 0.0001 Sex 1 22 4.1581 0.0536 Genotype 2 22 3.7246 0.0404 Basal 1 41 1142.2234 <.0001 Exposure:Sex 2 168 5.5801 0.0045 Exposure:Genotype 4 168 2.0960 0.0835 For the fixed effects, summary returns:
Value Std.Error DF t-value p-value (Intercept) 204417.8 109088.33 168 1.87387 0.0627 Exposure2 -192542.9 58653.05 168 -3.28274 0.0013 Exposure3 -332725.9 58653.05 168 -5.67278 0.0000 SexM -232599.7 104210.57 22 -2.23202 0.0361 GenotypeA -184772.3 125723.74 22 -1.46967 0.1558 GenotypeB -40073.2 128715.16 22 -0.31133 0.7585 Basal 1650.4 48.83 41 33.79680 0.0000 Exposure2:SexM 135000.3 58133.66 168 2.32224 0.0214 Exposure3:SexM 203637.4 58133.66 168 3.50292 0.0006 Exposure2:GenotypeA 106377.5 71472.21 168 1.48838 0.1385 Exposure3:GenotypeA 159548.1 71472.21 168 2.23231 0.0269 Exposure2:GenotypeB 101246.7 70397.58 168 1.43821 0.1522 Exposure3:GenotypeB 191111.2 70397.58 168 2.71474 0.0073 Am I correct in interpreting the interaction terms of this output as:
There is a statistically significant difference in response of M vs F both during exposure 2 and 3.
There is a statistically significant difference between genotype A and the reference genotype (wt) for exposure 3 (p=0.02) but not for exposure 2 (p=0.13)?
And similarly for genotype B?
What if I wanted to compare exposure 2 to exposure 3?
I read about using estimable in the gmodels package, but I struggle to understand its syntax.
Best Answer
If you look at the summary of your fixed effects portion of the model, you can label each row as follows:
Value Std.Error DF t-value p-value beta0 (Intercept) 204417.8 109088.33 168 1.87387 0.0627 beta1 Exposure2 -192542.9 58653.05 168 -3.28274 0.0013 beta2 Exposure3 -332725.9 58653.05 168 -5.67278 0.0000 beta3 SexM -232599.7 104210.57 22 -2.23202 0.0361 beta4 GenotypeA -184772.3 125723.74 22 -1.46967 0.1558 beta5 GenotypeB -40073.2 128715.16 22 -0.31133 0.7585 beta6 Basal 1650.4 48.83 41 33.79680 0.0000 beta7 Exposure2:SexM 135000.3 58133.66 168 2.32224 0.0214 beta8 Exposure3:SexM 203637.4 58133.66 168 3.50292 0.0006 beta9 Exposure2:GenotypeA 106377.5 71472.21 168 1.48838 0.1385 beta10 Exposure3:GenotypeA 159548.1 71472.21 168 2.23231 0.0269 beta11 Exposure2:GenotypeB 101246.7 70397.58 168 1.43821 0.1522 beta12 Exposure3:GenotypeB 191111.2 70397.58 168 2.71474 0.0073 So the fixed effects portion of your model looks like this:
beta0 + beta1*Exposure2 + beta2*Exposure3 + beta3*SexM + beta4*GenotypeA + beta5*GenotypeB + beta6*Basal + beta7*Exposure2*SexM + beta8*Exposure3*SexM + beta9*Exposure2*GenotypeA + beta10*Exposure3*GenotypeA + beta11*Exposure2*GenotypeB + beta12*Exposure3*GenotypeB where the symbol * denotes multiplication and all the beta's are true fixed effects (i.e., unknown but estimable from the data via the values listed in the Value column of your fixed effects summary).
If you are interested in describing the effect of Sex, for instance, all you have to do is to find all the terms in the model which include the dummy variable SexM (which is equal to 1 for Males and 0 for Females) and group them together. The coefficient of SexM obtained after this grouping denotes the effect of Sex:
(beta3 + beta7*Exp2 + beta8*Exp3)*SexM (1) From the above, we can see that the effect of Sex depends on the value of Exp. Recall that Exp2 and Exp3 are dummy variables defined as follows: Exp2 = 1 if Exp = 2 and 0 otherwise; Exp3 = 1 if Exp = 3 and 0 otherwise. We can exploit this to spell out the effect of SexM for:
Exp = 1 (that is, for Exp2 = 0 and Exp3 = 0); Exp = 2 (that is, for Exp2 = 1 and Exp3 = 0); Exp = 3 (that is, for Exp2 = 0 and Exp3 = 1). When Exp = 1, substituting Exp2 = 0 and Exp3 = 0 in expression (1) above yields that the coefficient of SexM is beta3. This is the effect of Sex on your outcome variable when Exp = 1, all else being equal.
When Exp = 2, substituting Exp2 = 1 and Exp3 = 0 in expression (1) above yields that the coefficient of SexM is beta3 + beta7. This is the effect of Sex on your outcome variable when Exp = 2, all else being equal.
When Exp = 3, substituting Exp2 = 0 and Exp3 = 1 in expression (1) above yields that the coefficient of SexM is beta3 + beta8. This is the effect of Sex on your outcome variable when Exp = 3, all else being equal.
So if you wanted to set up linear combinations of parameters which encapsulate the effect of SexM for each value of Exp, all you need to do is to specify these combinations so they reflect the parameters you are interested in: beta3, beta3 + beta7 and beta3 + beta8. (Again, the betas are true fixed effects, not estimated fixed effects.)
Since the fixed effects portion of your model includes the parameters beta0 through beta12, you are going to set up each combination via a row vector which includes 13 components. The first component of this vector corresponds to beta0, the second component corresponds to beta1, …, the last component corresponds to beta12.
In R, beta3 is a combination of all the fixed effects model parameters with weights given by the components of the row vector c1:
c1 <- rep(0, 13) names(c1) <- paste0("beta",0:12) c1[names(c1)=="beta3"] <- 1 c1 In other words, beta3 = 0*beta0 + 0*beta1 + 0*beta2 + 1*beta3 + 0*beta4 + 0*beta5 + 0*beta6 + 0*beta7 + 0*beta8 + 0*beta9 + 0*beta10 + 0*beta11 + 0*beta12.
The linear combination of model parameters beta0 through beta12 that will encapsulate the parameter beta3 + beta7 is given by:
c2 <- rep(0, 13) names(c2) <- paste0("beta",0:12) c2[names(c2)=="beta3"] <- 1 c2[names(c2)=="beta7"] <- 1 c2 In other words, beta3 + beta7 = 0*beta0 + 0*beta1 + 0*beta2 + 1*beta3 + 0*beta4 + 0*beta5 + 0*beta6 + 1*beta7 + 0*beta8 + 0*beta9 + 0*beta10 + 0*beta11 + 0*beta12.
The linear combination of model parameters beta0 through beta12 that will encapsulate the parameter beta3 + beta8 is given by:
c3 <- rep(0, 13) names(c3) <- paste0("beta",0:12) c3[names(c3)=="beta3"] <- 1 c3[names(c3)=="beta8"] <- 1 c3 such that beta3 + beta8 = 0*beta0 + 0*beta1 + 0*beta2 + 1*beta3 + 0*beta4 + 0*beta5 + 0*beta6 + 0*beta7 + 1*beta8 + 0*beta9 + 0*beta10 + 0*beta11 + 0*beta12.
If you now want to simultaneously test hypotheses such as:
H0: beta3 = 0 versus Ha: beta3 != 0 H0: beta3 + beta7 = 0 versus Ha: beta3 + beta7 != 0 H0: beta3 + beta8 = 0 versus Ha: beta3 + beta8 != 0 you can achieve this by using the multcomp package in R. (Here, != 0 stands for "not equal to zero".) These hypotheses will enable you to determine whether Sex has an effect for any of the specific levels of Exp (i.e., if males differ from females with respect to the average response, all else being the same). This package can also be used to compute simultaneous confidence intervals for beta3, beta3 + beta7 and beta3 + beta8.
All you need to do to use multcomp in conjunction with your lme model is something like this:
library(multcomp) c <- rbind(c1, c2, c3) g <- glht(model, linfct = c) s <- summary(g, test=adjusted("holm")) s ci <- confint(summary(g, test=adjusted("holm"))) ci If you don't want to adjust p-values or confidence levels for multiplicity, just use adjusted("none") in the above.
Of course, there are R packages which will do all of the above for you with a minimal number of commands. But it helps to know how to test your own simple effects (e.g., beta3, beta3 + beta7, beta3 + beta8) when interested in probing a significant interaction.
In the approach I presented here, contrasts are specified via linear combinations of all model parameters. There are other ways to specify contrasts, which I will leave to others on this forum to address.
Similar Posts:
- Solved – In what sense is the interpretation of coefficients in a GLMM subject-specific
- Solved – In what sense is the interpretation of coefficients in a GLMM subject-specific
- Solved – In what sense is the interpretation of coefficients in a GLMM subject-specific
- Solved – In what sense is the interpretation of coefficients in a GLMM subject-specific
- Solved – Using R and plm to estimate fixed-effects models that include interactions with time