QK Box 12.2

Chadha et al. (2019) set up a completely randomized experiment studying the effects of different soil water holding capacities (WHCs) on various characteristics of the agronomic weed Lactuca serriola. Plants were grown individually in pots allocated to one of four WHCs (100%, 75%, 50% and 25%). There were seven plants for each WHC, although we just use a subset of the data from the 100% treatment group so there was no between-subjects factor. The number of leaves on each of the seven plants in 100% WHC soil was recorded weekly for nine weeks (counts at the start of the experiment, week 0, were omitted). Time was the within-subjects (repeated measures) fixed factor and individual plants were the random subjects. With nine weeks and reasonably linear trends through time for each plant (Figure 12.2), it made sense to treat time as a continuous covariate for analysis. Although week 0 was not included in the analysis, we did not centre time for analysis, so intercepts represent the number of leaves for week 0.

Prickly lettuce or Milk thistle. Mick Keough, CC SA-BY 4.0

Chadha, A., Florentine, S., Chauhan, B. S., Long, B. & Jayasundera, M. (2019). Influence of soil moisture regimes on growth, photosynthetic capacity, leaf biochemistry and reproductive capabilities of the invasive agronomic weed; Lactuca serriola. PLoS One, 14, e0218191.

Link to paper

Preliminaries

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

Import chadha data file

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

set contrasts from afex

set_sum_contrasts()

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

select 100%WHC

chadha1 <- subset(chadha, treat=="100WHC")

make individual a factor

chadha1$plant <- factor(chadha1$plant)

plot slopes - very consistent

xyplot(noleaves~week|plant, type=c("p","r"), auto.key=T, chadha1)

Fit OLS “ancova” model

chadha1.aov <- aov(noleaves~week*plant, data=chadha1)
plot(chadha1.aov)

summary(chadha1.aov)

            Df Sum Sq Mean Sq F value  Pr(>F)    
week         1   4999    4999  310.55 < 2e-16 ***
plant        6    710     118    7.35 1.2e-05 ***
week:plant   6    110      18    1.14    0.35    
Residuals   49    789      16                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get Type III SS)

Anova(lm(chadha1.aov), type='III')

Anova Table (Type III tests)

Response: noleaves
            Sum Sq Df F value Pr(>F)    
(Intercept)   3568  1  221.63 <2e-16 ***
week          4999  1  310.55 <2e-16 ***
plant          131  6    1.35   0.25    
week:plant     110  6    1.14   0.35    
Residuals      789 49                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get correct F-ratio and P value for week (tested against week by plant)

4999.050/18.400

[1] 272

1-pf(4999.050/18.400, 1, 6, lower.tail = TRUE, log.p = FALSE)

[1] 3.18e-06

no GG and HF adjustments as week is continuous with 1 df

Get variance components

chadha1.vca <- anovaMM(noleaves~week+(plant)+(week*plant), chadha1)
chadha1.vca

ANOVA-Type Estimation of Mixed Model:
--------------------------------------

    [Fixed Effects]

  int  week 
16.40  3.45 


    [Variance Components]

  Name       DF        SS         MS         VC        %Total    SD       CV[%]    
1 total      15.101789                       26.537662 100       5.151472 15.308619
2 plant      6         710.095238 118.349206 10.401893 39.196719 3.225197 9.584311 
3 week:plant 6         110.4      18.4       0.038377  0.144613  0.1959   0.582156 
4 error      49        788.772222 16.097392  16.097392 60.658668 4.012156 11.922915

Mean: 33.7 (N = 63) 

Experimental Design: balanced  |  Method: ANOVA

VCAinference(chadha1.vca, alpha=0.05, VarVC=TRUE, excludeNeg=FALSE, constrainCI=FALSE)

Inference from Mixed Model Fit
------------------------------

> VCA Result:
-------------

    [Fixed Effects]

  int  week 
16.40  3.45 


    [Variance Components]

  Name       DF      SS       MS       VC      %Total  SD     CV[%]   Var(VC)
1 total      15.1018                   26.5377 100     5.1515 15.3086        
2 plant      6       710.0952 118.3492 10.4019 39.1967 3.2252 9.5843  77.6065
3 week:plant 6       110.4    18.4     0.0384  0.1446  0.1959 0.5822  0.0343 
4 error      49      788.7722 16.0974  16.0974 60.6587 4.0122 11.9229 10.5766

Mean: 33.7 (N = 63) 

Experimental Design: balanced  |  Method: ANOVA


> VC:
-----
           Estimate CI LCL CI UCL One-Sided LCL One-Sided UCL
total       26.5377 14.506 63.339        15.949        54.661
plant       10.4019 -6.864 27.668        -4.088        24.892
week:plant   0.0384 -0.325  0.401        -0.266         0.343
error       16.0974 11.232 24.997        11.890        23.247

> SD:
-----
           Estimate CI LCL CI UCL One-Sided LCL One-Sided UCL
total         5.152   3.81  7.959         3.994         7.393
plant         3.225  -2.62  5.260        -2.022         4.989
week:plant    0.196  -0.57  0.633        -0.516         0.586
error         4.012   3.35  5.000         3.448         4.822

> CV[%]:
--------
           Estimate CI LCL CI UCL One-Sided LCL One-Sided UCL
total        15.309  11.32  23.65         11.87         21.97
plant         9.584  -7.79  15.63         -6.01         14.83
week:plant    0.582  -1.69   1.88         -1.53          1.74
error        11.923   9.96  14.86         10.25         14.33


95% Confidence Level  
SAS PROC MIXED method used for computing CIs 

Drop interaction

chadha2.aov <- aov(noleaves~week+plant, data=chadha1)
Anova(lm(chadha2.aov), type='III')

Anova Table (Type III tests)

Response: noleaves
            Sum Sq Df F value Pr(>F)    
(Intercept)   3568  1  218.22 <2e-16 ***
week          4999  1  305.78 <2e-16 ***
plant          710  6    7.24  1e-05 ***
Residuals      899 55                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

get new variance components

chadha2.vca <- anovaMM(noleaves~week+(plant), chadha1)
chadha2.vca

ANOVA-Type Estimation of Mixed Model:
--------------------------------------

    [Fixed Effects]

  int  week 
16.40  3.45 


    [Variance Components]

  Name  DF        SS         MS         VC        %Total    SD       CV[%]    
1 total 23.462944                       27.681988 100       5.261368 15.635196
2 plant 6         710.095238 118.349206 11.333402 40.941432 3.366512 10.004257
3 error 55        899.172222 16.348586  16.348586 59.058568 4.043338 12.015581

Mean: 33.7 (N = 63) 

Experimental Design: balanced  |  Method: ANOVA

VCAinference(chadha2.vca, alpha=0.05, VarVC=TRUE, excludeNeg=FALSE, constrainCI=FALSE)

Inference from Mixed Model Fit
------------------------------

> VCA Result:
-------------

    [Fixed Effects]

  int  week 
16.40  3.45 


    [Variance Components]

  Name  DF      SS       MS       VC      %Total  SD     CV[%]   Var(VC)
1 total 23.4629                   27.682  100     5.2614 15.6352        
2 plant 6       710.0952 118.3492 11.3334 40.9414 3.3665 10.0043 57.76  
3 error 55      899.1722 16.3486  16.3486 59.0586 4.0433 12.0156 9.7191 

Mean: 33.7 (N = 63) 

Experimental Design: balanced  |  Method: ANOVA


> VC:
-----
      Estimate CI LCL CI UCL One-Sided LCL One-Sided UCL
total     27.7  16.79   54.0         18.17          48.3
plant     11.3  -3.56   26.2         -1.17          23.8
error     16.3  11.62   24.7         12.27          23.1

> SD:
-----
      Estimate CI LCL CI UCL One-Sided LCL One-Sided UCL
total     5.26   4.10   7.35          4.26          6.95
plant     3.37  -1.89   5.12         -1.08          4.88
error     4.04   3.41   4.97          3.50          4.80

> CV[%]:
--------
      Estimate CI LCL CI UCL One-Sided LCL One-Sided UCL
total     15.6  12.18   21.8         12.67          20.7
plant     10.0  -5.61   15.2         -3.21          14.5
error     12.0  10.13   14.8         10.41          14.3


95% Confidence Level  
SAS PROC MIXED method used for computing CIs 

Now mixed effects modelcomparing random slopes and random intercepts using ML

chadha.lmer1 <- lmer(noleaves~week + (week|plant), REML=FALSE, chadha1)

boundary (singular) fit: see help('isSingular')

 # singular fit so set correlation between random effects to zero
chadha.lmer1 <- lmer(noleaves~week + (week||plant), REML=FALSE, chadha1)
chadha.lmer2 <- lmer(noleaves~week + (1|plant), REML=FALSE, chadha1)

Compare models

anova(chadha.lmer1, chadha.lmer2)

Data: chadha1
Models:
chadha.lmer2: noleaves ~ week + (1 | plant)
chadha.lmer1: noleaves ~ week + (week || plant)
             npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
chadha.lmer2    4 375 383   -183      367                    
chadha.lmer1    5 376 386   -183      366  0.83  1       0.36

AICc(chadha.lmer1, chadha.lmer2)

Focus on random intercept model - refit with REML

chadha.lmer3 <- lmer(noleaves~week + (1|plant), REML=TRUE, chadha1)
summary(chadha.lmer3, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: noleaves ~ week + (1 | plant)
   Data: chadha1

REML criterion at convergence: 366

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.923 -0.509  0.086  0.320  3.927 

Random effects:
 Groups   Name        Variance Std.Dev.
 plant    (Intercept) 11.3     3.37    
 Residual             16.3     4.04    
Number of obs: 63, groups:  plant, 7

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)   16.401      1.689 13.433    9.71  1.9e-07 ***
week           3.450      0.197 55.000   17.49  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
week -0.584

anova(chadha.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)    
week   4999    4999     1    55     306 <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get variance components

chadha.ci3 <- confint.merMod(chadha.lmer3, oldNames=FALSE)

Computing profile confidence intervals ...

chadha.vc3 <- (chadha.ci3)^2
print(chadha.vc3)

                      2.5 % 97.5 %
sd_(Intercept)|plant   2.77   38.2
sigma                 11.33   23.8
(Intercept)          170.12  390.4
week                   9.36   14.7

Random slopes output for comparison - refit using REML

chadha.lmer4 <- lmer(noleaves~week + (week||plant), REML=TRUE, chadha1)
summary(chadha.lmer4, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: noleaves ~ week + (week || plant)
   Data: chadha1

REML criterion at convergence: 365

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.017 -0.521  0.077  0.375  3.925 

Random effects:
 Groups   Name        Variance Std.Dev.
 plant    (Intercept)  6.797   2.607   
 plant.1  week         0.141   0.375   
 Residual             15.644   3.955   
Number of obs: 63, groups:  plant, 7

Fixed effects:
            Estimate Std. Error    df t value Pr(>|t|)    
(Intercept)    16.40       1.47 10.28    11.2  4.4e-07 ***
week            3.45       0.24 10.28    14.4  3.8e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
week -0.530

anova(chadha.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)    
week   3245    3245     1  10.3     207 3.8e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Check whether AR(1) covariance structure would improve the fit using nlme

chadha.lme1 <- lme(noleaves~week, random=~1|plant, method="ML", chadha1)
chadha.lme2 <- lme(noleaves~week, random=~1|plant, method="ML", correlation=corAR1(form=~1|plant), chadha1)
anova(chadha.lme1, chadha.lme2)

AICc(chadha.lme1, chadha.lme2)

no improvement with AR(1)

chadha <- lm(noleaves~week, chadha1)
summary(chadha)

Call:
lm(formula = noleaves ~ week, data = chadha1)

Residuals:
   Min     1Q Median     3Q    Max 
-9.101 -3.151  0.249  2.949 14.549 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   16.401      1.410    11.6   <2e-16 ***
week           3.450      0.251    13.8   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.14 on 61 degrees of freedom
Multiple R-squared:  0.756, Adjusted R-squared:  0.752 
F-statistic:  189 on 1 and 61 DF,  p-value: <2e-16