1
$\begingroup$

When modelling QPCR data using LME4 I am getting a result that tells me my treatment effect is insignificant. When I plot the data this looks wrong and if I use JMPpro the p value for Treatment is tiny, as expected.

I have this data:

Strain <- c("A", "B", "C")
Strain <- rep(Strain, each = 18)


Treatment <- c("Water", "dsRNA")
Treatment <- rep(Treatment, each=3, times=9)

Time <- c("24", "48", "72")
Time <- rep(Time, each=6, times=3)

Cohort <- c("1", "2", "3")
Cohort <- rep(Cohort, times = 18)

deltaCT <- c(-2.772848129, -2.382880211, -2.802474976, -2.545639038, -2.363836288, -2.280316353, -2.418494225, -2.333683968, -2.999289513, -0.828756332, -2.032594681, -1.340758324, -2.595740318,
-2.231285095, -2.265996933, -1.700395584, -1.744121552, -1.32080555, -0.910932541, -2.313327789, -2.346955299, -1.538050652, -0.690214157, -1.678500175, -2.221671104, -3.087368965, -3.176133156,
-1.53543663, -2.560225487, -2.625535011, -3.341409683, -3.042736053, -3.546187401, -0.139810562, -2.315425873, -2.392004967, -2.916553497, -2.504823685, -3.253599167, -2.005877495, -2.682225227,
-3.007395744, -3.070969582, -3.156041145,-3.030291557,-0.097309113, -2.800593376, -2.803215981, -2.930361748, -2.752057076, -1.931899071, -1.238015175, -1.467027664,-1.421812057)

mdat <- data.frame(Strain, Treatment, Time, Cohort, deltaCT)

mdat$Strain <- factor(mdat$Strain)
mdat$Treatment <- factor(mdat$Treatment)
mdat$Time <- factor(mdat$Time)
mdat$Cohort <- factor(mdat$Cohort)

Using lme4 to run a mixed model, calculating the significance of the three fixed effects with interactions (Strain, Treatment and Time) and one random effect (Cohort) on deltaCT. My code looks like this:

m1 <- lmer(deltaCT ~ Treatment * Strain * Time + (1|Cohort), data = mdat)

Anova(m1, test.statistic = 'F',  type =3)

My results show that the Treatment effect is not significant but Strain is and the interaction between Strain and time is too:

Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)

Treatment p = 0.59228

Strain p = 0.02364 *

Strain:Time p = 0.04189 *

If I plot the data I can see the treatment effect is clearly very strong and if I parse the same data with JMPpro, the p value for Treatment is <0.001. The Strain and Strain:Time interaction may very well have significant effects on deltaCT but the Treatment effect should be stronger.

What am I doing wrong here ? Thanks in advance.

$\endgroup$
2
  • 1
    $\begingroup$ You didn't compute the type 3 anova correctly: you need to set the sum contrast for the three factors. Also -- are you certain you want the type 3 anova in the first place? $\endgroup$
    – dipetkov
    Commented Jun 14 at 17:08
  • $\begingroup$ I think type 3 is necessary when looking at interactions. How would I set the sum contrast ? $\endgroup$
    – Mikeed
    Commented Jun 14 at 17:29

1 Answer 1

2
$\begingroup$

Your experimental design is balanced: there are two treatments, three strains and three time points, with 3 observations for every (treatment, strain, time) combination. In this situation, type 2 and type 3 ANOVAs produce the same result. Instead here's what we get with your code:

car::Anova(m1, test.statistic = "F", type = 2)
#> Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
#> 
#> Response: deltaCT
#>                             F Df Df.res    Pr(>F)    
#> Treatment             29.5543  1     34 4.656e-06 ***
#> Strain                 0.8296  2     34  0.444877    
#> Time                   0.6083  2     34  0.550118    
#> Treatment:Strain       0.1545  2     34  0.857424    
#> Treatment:Time         2.4734  2     34  0.099341 .  
#> Strain:Time            4.4768  4     34  0.005163 ** 
#> Treatment:Strain:Time  0.7148  4     34  0.587654

car::Anova(m1, test.statistic = "F", type = 3)
#> Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
#> 
#> Response: deltaCT
#>                             F Df Df.res  Pr(>F)    
#> (Intercept)           43.5907  1 26.205 5.1e-07 ***
#> Treatment              0.2923  1 34.000 0.59228    
#> Strain                 4.1894  2 34.000 0.02364 *  
#> Time                   2.4952  2 34.000 0.09747 .  
#> Treatment:Strain       0.1086  2 34.000 0.89737    
#> Treatment:Time         0.9618  2 34.000 0.39236    
#> Strain:Time            2.7870  4 34.000 0.04189 *  
#> Treatment:Strain:Time  0.7148  4 34.000 0.58765

We can compute type 3 ANOVA correctly if we first "center" the predictors. We center continuous variables by subtracting the mean: x.c = x - mean(x); we center categorical variables by using the sum contrast (aka deviation coding): contrasts(A) = contr.sum(nlevels(A)). The default in R is the treatment contrast (aka dummy coding); it compares factor levels to one level chosen as reference rather than to the mean of all the levels.

So a correct way to fit the model if we want to use type 3 ANOVA is:

m1.sum.contr <- lmer(
  deltaCT ~ Treatment * Strain * Time + (1 | Cohort),
  data = mdat,
  contrasts = list(
    Treatment = "contr.sum",
    Strain = "contr.sum",
    Time = "contr.sum"
  )
)

Now we get the same table for type 2 and type 3 ANOVAs, as expected.

Note also that we get the correct table with type 2 whether we center the predictors or not. This is one reason to prefer type 2 over type 3. See this great answer by @gung to learn more: How to interpret type I, type II, and type III ANOVA and MANOVA? I also like this blog post: Everything You Always Wanted to Know About ANOVA (But Were Afraid to Ask). Though I would say that you should know all that about ANOVA before using it to analyze your data.

car::Anova(m1.sum.contr, test.statistic = "F", type = 2)
#> Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
#> 
#> Response: deltaCT
#>                             F Df Df.res    Pr(>F)    
#> Treatment             29.5543  1     34 4.656e-06 ***
#> Strain                 0.8296  2     34  0.444877    
#> Time                   0.6083  2     34  0.550118    
#> Treatment:Strain       0.1545  2     34  0.857424    
#> Treatment:Time         2.4734  2     34  0.099341 .  
#> Strain:Time            4.4768  4     34  0.005163 ** 
#> Treatment:Strain:Time  0.7148  4     34  0.587654

car::Anova(m1.sum.contr, test.statistic = "F", type = 3)
#> Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
#> 
#> Response: deltaCT
#>                             F Df Df.res    Pr(>F)    
#> (Intercept)            0.0000  1  2.014  1.000000    
#> Treatment             29.5543  1 34.000 4.656e-06 ***
#> Strain                 0.8296  2 34.000  0.444877    
#> Time                   0.6083  2 34.000  0.550118    
#> Treatment:Strain       0.1545  2 34.000  0.857424    
#> Treatment:Time         2.4734  2 34.000  0.099341 .  
#> Strain:Time            4.4768  4 34.000  0.005163 ** 
#> Treatment:Strain:Time  0.7148  4 34.000  0.587654

PS: In R, one way to check the contrasts used in fitting a model is with the model.matrix function. For example, model.matrix(m1) reports "contr.treatment" for the three predictors while model.matrix(m1.sum.contr) reports "sum.contr".

PPS: In a comment the OP writes that "type 3 is necessary when looking at interactions". This is a misunderstanding. See this discussion by @EdM: How does type II vs type III ANOVA work when there are more than two categorical predictors? and references therein.

$\endgroup$
1
  • 1
    $\begingroup$ That's brilliant, I'll look into those. Thanks very much $\endgroup$
    – Mikeed
    Commented Jun 14 at 17:49

Your Answer

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

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