inconsistent results in two-way repeated measure analysis using mixed models

Question

I have two factors: interval with two levels, and spd_des with three levels. Each of the subjects (grouping variable ratID) repeats the experiment for each level of spd_des. In each of these experiments, I observe the variable cc for the two levels of interval (one observation for each level). I would like to study the effect of interval and spd_des (and their possible interaction) on the dependent variable cc. Basically, this is a two-way repeated measure analysis. The data look like in the image below.

I would like to use mixed effect models with the package glmmTMB. However, to start with something simpler, I used paired t-tests for each level of spd_des and I corrected the p-values with Bonferroni.

> s1    <- t.test(x1, y1, paired = TRUE)
> s2    <- t.test(x2, y2, paired = TRUE)
> s3    <- t.test(x3, y3, paired = TRUE)

> p_all <- c(s1$p.value,s2$p.value,s3$p.value)
[1] 0.17689345 0.02041172 0.02103782

> p.adjust(p_all,method="bonferroni")
[1] 0.53068035 0.06123517 0.06311345

It seems that there is no difference between the two levels of interval for the first level s1 of spd_des, and there is a tendency (albeit not significant after Bonferroni) for the levels s2 and s3. This makes sense if one looks at the image above.

Unfortunately, I get a completely different picture when I use the more correct mixed-effect models, and I wonder why. Given the design described above (in which the observations obtained for the two levels of interval are obtained within each level of spd_des within each subject), I nest spd_des within ratID

> linM <- glmmTMB(cc ~ interval*spd_des + (1|ratID/spd_des), data=dat_trf, na.action=na.omit,  control = glmmTMBControl(optCtrl=list(iter.max=1e3,eval.max=1e3)))

Running Anova on this model (with 'sum' contrast and Type III sum of square), I obtain:

Anova.glmmTMB(linM, type = 3) Analysis of Deviance Table (Type III Wald chisquare tests)

Response: cc
                    Chisq Df Pr(>Chisq)    
(Intercept)      434.5436  1  < 2.2e-16 ***
interval          23.6361  1  1.164e-06 ***
spd_des            0.7465  2     0.6885    
interval:spd_des   3.6512  2     0.1611    

Now it seems that there is an effect of interval independently on spd_des. This is not consistent with the previous result, and it does not convince me much by looking at the image above (I have also tried other random effect structures: using random slope on spd_des or just random intercept, and I get equivalent results). Question 1: Is there anything wrong with this model? Here is the summary:

summary(linM)

Family: gaussian  ( identity )
Formula:          cc ~ interval * spd_des + (1 | ratID/spd_des)
Data: dat_trf

     AIC      BIC   logLik deviance df.resid 
    22.2     40.1     -2.1      4.2       45 

Random effects:

Conditional model:
 Groups        Name        Variance  Std.Dev. 
 spd_des:ratID (Intercept) 5.092e-12 2.256e-06
 ratID         (Intercept) 1.849e-02 1.360e-01
 Residual                  5.191e-02 2.278e-01
Number of obs: 54, groups:  spd_des:ratID, 27; ratID, 10

Dispersion estimate for gaussian family (sigma^2): 0.0519 

Conditional model:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)         1.11193    0.05334  20.846  < 2e-16 ***
interval1           0.15137    0.03114   4.862 1.16e-06 ***
spd_des1            0.03030    0.04314   0.702   0.4824    
spd_des2           -0.03407    0.04481  -0.761   0.4469    
interval1:spd_des1 -0.07871    0.04283  -1.838   0.0661 .  
interval1:spd_des2  0.01630    0.04394   0.371   0.7108

[UPDATE] Turns out that if I run post-hoc tests on this model, even if there is no significant interaction according to Anova (see above), I obtain that there is significant difference between the two levels of interval only for some of the levels of spd_des. These results seem to make sense based on the picture and the simple paired ttests above. Question 2: Should I run post-hoc on the full models even if there is no significant interaction according to Anova? If so, how do I interpret the two different results?

> ph_conditional <- c("interval1+interval1:spd_des1=0",
                      "interval1+interval1:spd_des2=0",
                      "interval1-interval1:spd_des1-interval1:spd_des2=0")


> linM.ph <- glht(linM, linfct = ph_conditional);
> summary(linM.ph,test = adjusted("bonferroni"))
Simultaneous Tests for General Linear Hypotheses

Fit: glmmTMB(formula = cc_marg ~ interval * spd_des + (1 | ratID/spd_des), 
    data = dat_trf, na.action = na.omit, control = glmmTMBControl(optCtrl = list(iter.max = 1000, 
        eval.max = 1000)), ziformula = ~0, dispformula = ~1)

Linear Hypotheses:
                                                         Estimate Std. Error z value Pr(>|z|)    
interval1 + interval1:spd_des1 == 0                       0.07266    0.05095   1.426 0.461521    
interval1 + interval1:spd_des2 == 0                       0.16766    0.05370   3.122 0.005389 ** 
interval1 - interval1:spd_des1 - interval1:spd_des2 == 0  0.21379    0.05696   3.753 0.000524 ***

[END OF UPDATE]

Then, I tried to fit a full model where I defined random effects of intercept*spd_des, and I obtained something more similar to the picture above.

> linM2 <- glmmTMB(cc ~ interval*spd_des + (interval*spd_des|ratID), data=dat_trf, na.action=na.omit,  control = glmmTMBControl(optCtrl=list(iter.max=1e3,eval.max=1e3)))
> Anova.glmmTMB(linM2, type = 3)
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: cc
                    Chisq Df Pr(>Chisq)    
(Intercept)      448.8546  1  < 2.2e-16 ***
interval          35.5090  1  2.539e-09 ***
spd_des            1.2995  2   0.522170    
interval:spd_des  11.3635  2   0.003408 ** 

There seems to be a significant interaction between the two factors. When I run the post-hoc tests to compare the two levels of interval for each level of spd_des (adjusted with Bonferroni):

> ph_conditional <- c("interval1+interval1:spd_des1=0",
                      "interval1+interval1:spd_des2=0",
                      "interval1-interval1:spd_des1-interval1:spd_des2=0")


> linM.ph <- glht(linM, linfct = ph_conditional);
> summary(linM.ph,test = adjusted("bonferroni"))
Simultaneous Tests for General Linear Hypotheses

Fit: glmmTMB(formula = cc ~ interval * spd_des + (interval * 
    spd_des | ratID), data = dat_trf, na.action = na.omit, control = glmmTMBControl(optCtrl = list(iter.max = 1000, 
    eval.max = 1000)), ziformula = ~0, dispformula = ~1)

Linear Hypotheses:
                                                         Estimate Std. Error z value Pr(>|z|)   
interval1 + interval1:spd_des1 == 0                       0.07266    0.04704   1.545  0.36740   
interval1 + interval1:spd_des2 == 0                       0.19857    0.06155   3.226  0.00377 **
interval1 - interval1:spd_des1 - interval1:spd_des2 == 0  0.25354    0.08649   2.931  0.01012 *

This now tells that for the first level of spd_des there is no significant difference between the two levels of interval, but this difference becomes significant for the other two levels of spd_des. I had somewhat assumed that this full model was not correct and could not be fitted because I had a single observation for each interval:spd_des but the results seem to be more reasonable than model linM above. Question 3: Is this model (linM2) correct? and why?

Note that I have checked all the assumptions for both models and they are all met: constant variance of residuals, normal residuals and normal random errors. Also, I know I could use nlme or lme4, but the former has issues at estimating the denominator DF for random slope models, and the latter provides less flexibility than glmmTMB in defining the structure of the var-cov matrix for the random effect (in case I needed to change it). I would use to glmmTMB, if possible. Thanks!

Answer 1

I don't see why the results are inconsistent, especially when you take the full RE structure, which corresponds most closely to the paired t-tests. The bits I notice are:

1. You're using Type-III tests. Why? It's very easy to test nonsense hypotheses with Type-III tests. Check out the classic Exegesis on Linear Models from Venables.

2. The paired t-tests each have only a piece of the data, while the mixed model has all of the data. As such, there's more data to estimate the effect of interval and more power, at least for the marginal effect of interval across all spd_des.

3. You don't need post hoc tests for the mixed model, if you're clever with your contrasts. Or if you want something that looks more like your p-values from your paired t-test, try the emmeans package.

Answer 2

I don't know if it's a way toward the answer, but I do not understand the logic behind the first model. If RatID is the random effect, and assuming each rat really has a different ID, I do not see what the sp_des:Rat additional random effect term really means. A syntax like (1|RatID) would have been expected, for random effects on intercepts only. Especially if each rat experiences the three sp_des levels... Indeed, if you look at the results of the model, you can see two random effect terms for the intercept: RatID and spd_des:ratID; I guess the second somehow catches up all interaction effects. In these aspects, the second, « full », model does not have this strange random effect syntax...