QK Box 10.8

Wamelink et al. (2014) compared the growth of 14 species of plant grown in soil collected from near the Rhine River with simulant regolith representing the soil conditions on the moon and on Mars. We will focus on total biomass of Sedum reflexum as the response variable. Twenty spatial blocks were established in a glasshouse and for this species, there were 3 small pots (each with five seeds) in each block, representing the three soils types. Soil type is considered a fixed factor and block is a random factor. The full layout of the experiment is shown clearly in Figures 1 and 2 of their paper.

[]

Sedum. Mick Keough,

The paper is here

Wamelink, G. W. W., Frissel, J. Y., Krijnen, W. H. J., Verwoert, M. R. & Goedhart, P. W. (2014). Can plants grow on Mars and the Moon: A growth experiment on Mars and Moon soil simulants. PLoS One, 9, e103138.

Preliminaries

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

Import wamelink data file (wamelink1.csv)

wamelink <- read.csv("../data/wamelink1.csv")
head(wamelink,5)

NA

Set contrasts from afex & make block a factor

set_sum_contrasts()

wamelink$block <- factor(wamelink$block)

Check assumptions with boxplots

boxplot(totbiomass~soil, data=wamelink)

Fit OLS model with default aov SS

wamelink.aov <- aov(totbiomass~soil+block, data=wamelink)

Check residuals - some evidence for interaction

plot(wamelink.aov)

Do interaction plot

interaction.plot(wamelink$block,wamelink$soil,wamelink$totbiomass)

Order of treatments consistent but differences vary greatly between blocks

boxplot(log10(totbiomass)~soil, data=wamelink)

Transform to logs due to variance heterogeneity and to minimise interaction

wamelink$ltotbiomass <- log10(wamelink$totbiomass)

Recheck diagnostics - much better with less interaction

boxplot(ltotbiomass~soil, data=wamelink)

interaction.plot(wamelink$block,wamelink$soil,wamelink$ltotbiomass)

Fit OLS model with default aov SS

wamelink.aov1 <- aov(ltotbiomass~soil+block, data=wamelink)

Check residuals - look OK

plot(wamelink.aov1)

Examine results

options(digits = 3)
summary(wamelink.aov1)

            Df Sum Sq Mean Sq F value Pr(>F)    
soil         2  19.39    9.69  196.02 <2e-16 ***
block       19   0.66    0.03    0.71   0.79    
Residuals   38   1.88    0.05                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Use VCA package to get anova var comps

Note that VCs can be -ve when based on OLS model fits

wamelink.vca <- anovaMM(ltotbiomass~soil+(block), wamelink)

Convert variable soil from "character" to "factor"!

wamelink.vca

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

    [Fixed Effects]

      int soilEarth  soilMars  soilMoon 
    0.438     0.753     1.391     0.000 


    [Variance Components]

  Name  DF SS       MS       VC       %Total SD       CV[%]    
1 total 57                   0.049458 100    0.222391 19.298591
2 block 19 0.664948 0.034997 0*       0*     0*       0*       
3 error 38 1.879389 0.049458 0.049458 100    0.222391 19.298591

Mean: 1.15 (N = 60) 

Experimental Design: balanced  |  Method: ANOVA | * VC set to 0 | adapted MS used for total DF

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

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

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

    [Fixed Effects]

      int soilEarth  soilMars  soilMoon 
    0.438     0.753     1.391     0.000 


    [Variance Components]

  Name  DF SS     MS     VC     %Total SD     CV[%]   Var(VC)
1 total 57               0.0495 100    0.2224 19.2986        
2 block 19 0.6649 0.035  0*     0*     0*     0*      0      
3 error 38 1.8794 0.0495 0.0495 100    0.2224 19.2986 1e-04  

Mean: 1.15 (N = 60) 

Experimental Design: balanced  |  Method: ANOVA | * VC set to 0 | adapted MS used for total DF


> VC:
-----
      Estimate CI LCL CI UCL One-Sided LCL One-Sided UCL
total 0.0495   0.0353 0.0741 0.0373        0.0694       
block 0        0*     0.0105 0*            0.0088       
error 0.0495   0.033  0.0821 0.0352        0.0755       

> SD:
-----
      Estimate CI LCL CI UCL One-Sided LCL One-Sided UCL
total 0.2224   0.188  0.2723 0.1931        0.2634       
block 0        0*     0.1024 0*            0.0938       
error 0.2224   0.1817 0.2866 0.1876        0.2748       

> CV[%]:
--------
      Estimate CI LCL  CI UCL  One-Sided LCL One-Sided UCL
total 19.2986  16.3152 23.6275 16.7546       22.8536      
block 0        0*      8.8866  0*            8.141        
error 19.2986  15.7717 24.8716 16.2822       23.8483      


95% Confidence Level  | * CI-limits constrained to be >= 0
SAS PROC MIXED method used for computing CIs 

Fit random intercept model using lme4 and REML

wamelink.lmer1 <- lmer(ltotbiomass~soil + (1|block), REML=TRUE, wamelink)

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

summary(wamelink.lmer1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: ltotbiomass ~ soil + (1 | block)
   Data: wamelink

REML criterion at convergence: -4.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6076 -0.6731 -0.0887  0.5306  2.7223 

Random effects:
 Groups   Name        Variance Std.Dev.
 block    (Intercept) 9.83e-26 3.14e-13
 Residual             4.46e-02 2.11e-01
Number of obs: 60, groups:  block, 20

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   1.1524     0.0273 57.0000    42.2   <2e-16 ***
soil1         0.0385     0.0386 57.0000     1.0     0.32    
soil2         0.6762     0.0386 57.0000    17.5   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
      (Intr) soil1 
soil1  0.000       
soil2  0.000 -0.500
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

AICc(wamelink.lmer1)

[1] 6.83

F tests of fixed effects (Type III SS) - matches t-test from summary command

anova(wamelink.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)    
soil   19.4    9.69     2    38     217 <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

CI on variance components (remembering to square CIs from lmer with are in SD units)

wamelink.ci3 <- confint.merMod(wamelink.lmer1)

Computing profile confidence intervals ...
Warning: slightly lower deviances (diff=-3.19744e-14) detectedWarning: Last two rows have identical or NA .zeta values: using minstepWarning: Last two rows have identical or NA .zeta values: using minstepWarning: slightly lower deviances (diff=-3.19744e-14) detectedWarning: non-monotonic profile for .sig01Warning: bad spline fit for .sig01: falling back to linear interpolationWarning: collapsing to unique 'x' values

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

              2.5 % 97.5 %
.sig01      0.00000 0.0093
.sigma      0.03025 0.0620
(Intercept) 1.20872 1.4528
soil1       0.00132 0.0129
soil2       0.36152 0.5641

Fit model that omits blocks

wamelink.gls2 <- gls(ltotbiomass~soil,wamelink)
summary(wamelink.gls2)

Generalized least squares fit by REML
  Model: ltotbiomass ~ soil 
  Data: wamelink 

Coefficients:

 Correlation: 
      (Intr) soil1
soil1  0.0        
soil2  0.0   -0.5 

Standardized residuals:
    Min      Q1     Med      Q3     Max 
-1.6076 -0.6731 -0.0887  0.5306  2.7223 

Residual standard error: 0.211 
Degrees of freedom: 60 total; 57 residual

AICc(wamelink.gls2)

[1] 4.45