5
$\begingroup$

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.

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ May I suggest that if you would like a "canonical answer" to your question, you should consider sharing your data (or a modified version of the data to protect your research findings)? $\endgroup$
    – Isabella Ghement
    May 31, 2018 at 14:35

3 Answers 3

Reset to default
10
+150
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Wow Isabella, thank you for taking the time to write such a clear explanation, I really appreciate this. $\endgroup$
    – nico
    May 31, 2018 at 15:28
  • $\begingroup$ @nico: Great! Note that I edited the ci command in my answer. $\endgroup$
    – Isabella Ghement
    May 31, 2018 at 17:51
2
$\begingroup$

You could also reparameterise your model and look at the parameter estimates Something like model2 <- lme(AUC ~ Exposure -1 + Exposure:Sex + Exposure:Genotype + Basal, random = ~ 1 | Date / Experiment / Cell, data = mydata) then look differences between sexes with Exposure categories and same for Genotype. And try again with exposure nested within Sex and Genotype BUT be cautious of rising Type I error rates.

model3 <- lme(AUC ~ Sex/Exposure -1 + Genotype/Exposure + Basal, random = ~ 1 | Date / Experiment / Cell, data = mydata)

Improve this answer
$\endgroup$
2
$\begingroup$

I think for the kind of questions you are asking, you want to use post-hoc comparisons, such as with emmeans, rather than trying to interpret the summary output.

I have sample code below. I created some data and used a simpler model, so obviously my results will be different than yours.

Because your model is complex, I recommend reading the documentation and vignettes for the emmeans package before relying on the results.

For example, for the "What if I wanted to compare exposure 2 to exposure 3?" question, you can do pairwise comparisons of all levels of Exposure.

For the "There is a statistically significant difference in response of M vs F both during exposure 2 and 3.", you can compare levels of Sex within each level of Exposure. Or you can do pairwise comparisons of each combination of Sex and Exposure. These results could be summarized in a compact letter display (cld).

if(!require(nlme)){install.packages("nlme")}
if(!require(car)){install.packages("car")}
if(!require(emmeans)){install.packages("emmeans")}

set.seed(1234)

Exposure = factor(rep(c("1", "2", "3"), 1, each=6))
Sex      = factor(rep(c("F", "M"), 9))
Subject  = factor(rep(letters[1:6],3))

AUC      = as.numeric(Exposure) * 2 +
           as.numeric(Sex) +
           as.numeric(Exposure) * as.numeric(Sex) * 1.2
AUC      = AUC + rnorm(length(Exposure), 0, 1)

Data = data.frame(Exposure, Sex, Subject, AUC) 

str(Data) 

library(nlme) 

model = lme(AUC ~ Exposure * Sex, random=~1|Subject, data=Data)

library(car)

Anova(model)


### What if I wanted to compare exposure 2 to exposure 3?

library(emmeans) 

marginal = emmeans(model, ~ Exposure)

pairs(marginal)

   ### NOTE: Results may be misleading due to involvement in interactions

   ###  contrast  estimate        SE df t.ratio p.value
   ###  1 - 2    -3.334045 0.4925335  8  -6.769  0.0004
   ###  1 - 3    -7.595154 0.4925335  8 -15.421  <.0001
   ###  2 - 3    -4.261108 0.4925335  8  -8.651  0.0001


### There is a statistically significant difference in response of M vs F both during exposure 2 and 3.

marginal = emmeans(model, ~ Sex | Exposure)

pairs(marginal) 

   ### Exposure = 1:
   ###  contrast  estimate        SE df t.ratio p.value
   ###  F - M    -1.577096 0.7533083  4  -2.094  0.1044
   ### 
   ### Exposure = 2:
   ###  contrast  estimate        SE df t.ratio p.value
   ###  F - M    -3.127110 0.7533083  4  -4.151  0.0142
   ### 
   ### Exposure = 3:
   ###  contrast  estimate        SE df t.ratio p.value
   ###  F - M    -4.390249 0.7533083  4  -5.828  0.0043


### Or you could compare all levels of the interaction

marginal = emmeans(model, ~ Sex + Exposure)

pairs(marginal)

cld(marginal, Letters=letters)

   ###  Sex Exposure    emmean        SE df  lower.CL  upper.CL .group
   ###  F   1         4.302167 0.5326694  5  2.932896  5.671437  a    
   ###  M   1         5.879262 0.5326694  4  4.400335  7.358190  ab   
   ###  F   2         6.861205 0.5326694  5  5.491935  8.230475   bc  
   ###  M   2         9.988315 0.5326694  4  8.509387 11.467242    cd 
   ###  F   3        10.490744 0.5326694  5  9.121473 11.860014     d 
   ###  M   3        14.880993 0.5326694  4 13.402065 16.359920      e
   ### 
   ### Confidence level used: 0.95 
   ### P value adjustment: tukey method for comparing a family of 6 estimates 
   ### significance level used: alpha = 0.05
Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you Sal, emmeans seems indeed quite a handy package. Something more to read about! $\endgroup$
    – nico
    Jun 4, 2018 at 10:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.