QK Box 11.3

The rat pup dataset from Pinheiro and Bates (2000) was introduced in Chapter 10. The between-plots factor was drug treatment (fixed with three doses: control, low dose, high dose) and the plots were ten female rats, and their subsequent litters (level 2), assigned to each treatment. The pups from a litter produced by each female were sexed and weighed, so pups were the sub-plots (level 1) with sex (fixed: male vs females) as the within-plots factor and pup weight as the response variable. Litter size (number of pups) was recorded as a continuous covariate for each litter (level 2) as we might reasonably expect larger litters to have smaller pups on average. The design was unbalanced for both plots (only seven out of ten rats survived in the high dose group) and subplots (litter size ranged between two and 18 pups). Additionally, one litter (number 12) only had female pups so this litter was omitted from the analysis to avoid a missing cell.

This analysis uses a subset of the data - it’s here

Preliminaries

First, load the required packages (afex, car, lattice, lme4, lmerTest, nlme, VCA, ez, emmeans, Rmisc, MuMIn)

Import ratpupno12 data file (ratpupno12.csv)

ratpupno12 <- read.csv("../data/ratpupno12.csv")
ratpupno12

Set contrasts from afex Make litter a factor

set_sum_contrasts()

setting contr.sum globally: options(contrasts=c('contr.sum', 'contr.poly'))

ratpupno12$litter <- factor(ratpupno12$litter)

Check boxplots - some variance difference but not related to mean

boxplot(weight~treatment*sex, ratpupno12)

Not suited to OLS fit because of imbalance

Without covariate

First, model using lmer without litter size as covariate - this model includes litter(treatment) x sex interaction

ratpup.lmer1 <- lmer(weight~treatment+sex+treatment*sex+(1|litter)+(1|litter:sex), REML=TRUE, ratpupno12)
plot(ratpup.lmer1)

summary(ratpup.lmer1, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: weight ~ treatment + sex + treatment * sex + (1 | litter) + (1 |      litter:sex)
   Data: ratpupno12

REML criterion at convergence: 420.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.8869 -0.5151  0.0131  0.5374  2.8411 

Random effects:
 Groups     Name        Variance Std.Dev.
 litter:sex (Intercept) 0.01564  0.1251  
 litter     (Intercept) 0.25296  0.5030  
 Residual               0.15727  0.3966  
Number of obs: 320, groups:  litter:sex, 52; litter, 26

Fixed effects:
                Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)      6.13527    0.10450 23.14353  58.709  < 2e-16 ***
treatment1       0.27970    0.14181 22.99973   1.972   0.0607 .  
treatment2      -0.06747    0.15614 23.56229  -0.432   0.6696    
sex1            -0.18000    0.03096 25.15990  -5.815  4.5e-06 ***
treatment1:sex1 -0.03385    0.04090 22.29970  -0.828   0.4167    
treatment2:sex1  0.02127    0.04845 30.21454   0.439   0.6637    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmn1 trtmn2 sex1   trt1:1
treatment1  -0.117                            
treatment2   0.155 -0.529                     
sex1         0.002  0.017 -0.010              
trtmnt1:sx1  0.017  0.015 -0.005 -0.194       
trtmnt2:sx1 -0.010 -0.004 -0.005  0.287 -0.578

anova (ratpup.lmer1, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
              Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
treatment     0.6938  0.3469     2 23.091  2.2057    0.1329    
sex           5.3173  5.3173     1 25.160 33.8089 4.501e-06 ***
treatment:sex 0.1081  0.0540     2 23.749  0.3434    0.7128    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get variance components with profile CIs

ratpup.ci <- confint.merMod(ratpup.lmer1,oldNames=FALSE)

Computing profile confidence intervals ...

ratpup.vc <- (ratpup.ci)^2
print(ratpup.vc)

                                 2.5 %       97.5 %
sd_(Intercept)|litter:sex 0.000000e+00  0.046258112
sd_(Intercept)|litter     1.234467e-01  0.429932353
sigma                     1.335858e-01  0.188435253
(Intercept)               3.522724e+01 40.136473884
treatment1                9.926651e-05  0.305713065
treatment2                1.339955e-01  0.053504195
sex1                      5.675929e-02  0.014566084
treatment1:sex1           1.254337e-02  0.001852688
treatment2:sex1           4.848962e-03  0.012933502

Now simplify model by testing if litter(treatment) x sex term should be retained

ratpup.lmer1a <- lmer(weight~treatment+sex+treatment*sex+(1|litter)+(1|litter:sex), REML=FALSE, ratpupno12)
ratpup.lmer2 <- lmer(weight~treatment+sex+treatment*sex+(1|litter), REML=FALSE, ratpupno12)
anova(ratpup.lmer1a, ratpup.lmer2)

Data: ratpupno12
Models:
ratpup.lmer2: weight ~ treatment + sex + treatment * sex + (1 | litter)
ratpup.lmer1a: weight ~ treatment + sex + treatment * sex + (1 | litter) + (1 | litter:sex)
              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
ratpup.lmer2     8 415.72 445.86 -199.86   399.72                     
ratpup.lmer1a    9 417.09 451.00 -199.54   399.09 0.6278  1     0.4282

AICc(ratpup.lmer1a, ratpup.lmer2)

Use simpler model

ratpup.lmer2a <- lmer(weight~treatment+sex+treatment*sex+(1|litter), REML=TRUE, ratpupno12)
summary(ratpup.lmer2a, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: weight ~ treatment + sex + treatment * sex + (1 | litter)
   Data: ratpupno12

REML criterion at convergence: 422.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-7.3671 -0.4894  0.0064  0.5477  3.1089 

Random effects:
 Groups   Name        Variance Std.Dev.
 litter   (Intercept) 0.2641   0.5139  
 Residual             0.1634   0.4042  
Number of obs: 320, groups:  litter, 26

Fixed effects:
                 Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)       6.13428    0.10509  23.12458  58.373  < 2e-16 ***
treatment1        0.28459    0.14266  23.01002   1.995    0.058 .  
treatment2       -0.06812    0.15696  23.52001  -0.434    0.668    
sex1             -0.17674    0.02555 295.61185  -6.917 2.86e-11 ***
treatment1:sex1  -0.03098    0.03317 294.19956  -0.934    0.351    
treatment2:sex1   0.02024    0.04091 296.97605   0.495    0.621    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmn1 trtmn2 sex1   trt1:1
treatment1  -0.116                            
treatment2   0.154 -0.529                     
sex1         0.005  0.018 -0.010              
trtmnt1:sx1  0.019  0.019 -0.007 -0.242       
trtmnt2:sx1 -0.009 -0.006 -0.003  0.351 -0.600

anova(ratpup.lmer2a, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
              Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
treatment     0.7382  0.3691     2  23.076  2.2590    0.1271    
sex           7.8167  7.8167     1 295.612 47.8417 2.862e-11 ***
treatment:sex 0.1437  0.0718     2 295.032  0.4397    0.6446    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get var comps

ratpup.ci <- confint.merMod(ratpup.lmer2,oldNames=FALSE)

Computing profile confidence intervals ...

ratpup.vc <- (ratpup.ci)^2
print(ratpup.vc)

                             2.5 %       97.5 %
sd_(Intercept)|litter 1.310156e-01  0.437332351
sigma                 1.381952e-01  0.191155575
(Intercept)           3.521029e+01 40.138977305
treatment1            1.448129e-04  0.310379990
treatment2            1.351348e-01  0.053876450
sex1                  5.113174e-02  0.015893474
treatment1:sex1       9.272754e-03  0.001119387
treatment2:sex1       3.483295e-03  0.010190477

Fit ratpup.lmer2 model including covariate with simple effect

ratpup.lmer4 <- lmer(weight~treatment+sex+treatment*sex+lsize+(1|litter), REML=TRUE, ratpupno12)
summary(ratpup.lmer4, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: weight ~ treatment + sex + treatment * sex + lsize + (1 | litter)
   Data: ratpupno12

REML criterion at convergence: 406.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-7.4805 -0.5008  0.0275  0.5704  3.0167 

Random effects:
 Groups   Name        Variance Std.Dev.
 litter   (Intercept) 0.1001   0.3164  
 Residual             0.1630   0.4037  
Number of obs: 320, groups:  litter, 26

Fixed effects:
                 Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)       7.67662    0.27842  27.04136  27.572  < 2e-16 ***
treatment1        0.42523    0.09489  22.62806   4.482 0.000175 ***
treatment2       -0.42228    0.11751  23.18499  -3.594 0.001520 ** 
sex1             -0.17434    0.02541 300.85249  -6.861 3.90e-11 ***
lsize            -0.12722    0.02206  25.89312  -5.766 4.58e-06 ***
treatment1:sex1  -0.03163    0.03303 298.05341  -0.958 0.339005    
treatment2:sex1   0.02164    0.04062 303.45006   0.533 0.594629    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmn1 trtmn2 sex1   lsize  trt1:1
treatment1   0.221                                   
treatment2  -0.451 -0.577                            
sex1        -0.023  0.020  0.001                     
lsize       -0.970 -0.257  0.502  0.026              
trtmnt1:sx1  0.025  0.034 -0.020 -0.237 -0.018       
trtmnt2:sx1 -0.025 -0.017  0.010  0.346  0.023 -0.599

anova(ratpup.lmer4, type="3", ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
              Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
treatment     3.5205  1.7603     2  22.200 10.7987 0.0005302 ***
sex           7.6728  7.6728     1 300.852 47.0720 3.903e-11 ***
lsize         5.4189  5.4189     1  25.893 33.2448 4.583e-06 ***
treatment:sex 0.1499  0.0750     2 299.778  0.4599 0.6317951    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Fit ratpup.lmer2 model including covariate (match West analysis)

ratpup.lmer3 <- lmer(weight~treatment+sex+treatment*sex+lsize+treatment*lsize+(1|litter), REML=TRUE, ratpupno12)
summary(ratpup.lmer3, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: weight ~ treatment + sex + treatment * sex + lsize + treatment *      lsize + (1 | litter)
   Data: ratpupno12

REML criterion at convergence: 411.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-7.4770 -0.5069  0.0096  0.5761  3.0265 

Random effects:
 Groups   Name        Variance Std.Dev.
 litter   (Intercept) 0.08843  0.2974  
 Residual             0.16287  0.4036  
Number of obs: 320, groups:  litter, 26

Fixed effects:
                   Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)        7.429896   0.323326  22.059200  22.980  < 2e-16 ***
treatment1         0.595882   0.389028  23.029597   1.532 0.139215    
treatment2         0.405952   0.402879  24.376536   1.008 0.323532    
sex1              -0.174655   0.025395 301.020191  -6.878 3.52e-11 ***
lsize             -0.115517   0.025105  22.152562  -4.601 0.000137 ***
treatment1:sex1   -0.030903   0.033007 298.289775  -0.936 0.349909    
treatment2:sex1    0.018630   0.040602 303.337744   0.459 0.646664    
treatment1:lsize  -0.005981   0.029544  22.655861  -0.202 0.841387    
treatment2:lsize  -0.072522   0.034670  24.334329  -2.092 0.047065 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmn1 trtmn2 sex1   lsize  trt1:1 trt2:1 trtm1:
treatment1  -0.460                                                 
treatment2  -0.360  0.001                                          
sex1        -0.026  0.029 -0.023                                   
lsize       -0.969  0.474  0.237  0.033                            
trtmnt1:sx1  0.027 -0.011  0.013 -0.237 -0.025                     
trtmnt2:sx1 -0.018  0.011 -0.041  0.346  0.023 -0.599              
trtmnt1:lsz  0.484 -0.966  0.071 -0.027 -0.523  0.022 -0.020       
trtmnt2:lsz  0.214  0.063 -0.956  0.027 -0.069 -0.021  0.047 -0.179

anova(ratpup.lmer3, type="3", ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
                Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
treatment       0.5476  0.2738     2  23.614  1.6807 0.2077737    
sex             7.7038  7.7038     1 301.020 47.3013 3.525e-11 ***
lsize           3.4484  3.4484     1  22.153 21.1732 0.0001367 ***
treatment:sex   0.1454  0.0727     2 299.889  0.4464 0.6403272    
treatment:lsize 0.7689  0.3845     2  22.907  2.3602 0.1169492    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Contrast control to other 2 groups

ratpup3.emm <- emmeans(ratpup.lmer3, ~treatment)

NOTE: Results may be misleading due to involvement in interactions

ratpup3.con <- contrast(ratpup3.emm, "trt.vs.ctrl", ref="Control")
summary(ratpup3.con, adjust="none")

 contrast       estimate    SE   df t.ratio p.value
 High - Control   -1.082 0.225 20.9  -4.801  0.0001
 Low - Control    -0.466 0.148 19.7  -3.138  0.0052

Results are averaged over the levels of: sex 
Degrees-of-freedom method: kenward-roger 

Get var comps

ratpup.ci <- confint.merMod(ratpup.lmer3,oldNames=FALSE)

Computing profile confidence intervals ...

ratpup.vc <- (ratpup.ci)^2
print(ratpup.vc)

                             2.5 %       97.5 %
sd_(Intercept)|litter  0.035086693 1.308154e-01
sigma                  0.137587610 1.900248e-01
(Intercept)           46.899112019 6.419454e+01
treatment1             0.010371330 1.690570e+00
treatment2             0.103954079 1.278899e+00
sex1                   0.049471899 1.516133e-02
lsize                  0.025826883 4.948339e-03
treatment1:sex1        0.009221154 1.067045e-03
treatment2:sex1        0.003524194 9.702075e-03
treatment1:lsize       0.003531796 2.208381e-03
treatment2:lsize       0.018197201 9.723359e-05