Caballes et al. (2016) examined the effects of maternal nutrition (three treatments: starved or fed one of two coral genera: Acropora or Porites) on the larval biology of crown-of-thorns seastars. There were three female seastars nested within each treatment, 50 larvae reared from each female were placed into each of three glass culture jars and the lengths of ten larvae from each jar after four days were measured after 4 days. This fully balanced design has maternal nutrition as a fixed factor with three random levels of nesting: females within nutrition treatment, jars within females and individual larvae within jars.
The paper is here
Caballes, C. F., Pratchett, M. S., Kerr, A. M. & Rivera-Posada, J. A. (2016). The role of maternal nutrition on oocyte size and quality, with respect to early larval development in the coral-eating starfish, Acanthaster planci. PLoS One, 11, e0158007.
First, load the required packages (afex, car, lattice, lme4, lmerTest, nlme, VCA, Rmisc)
Import caballes data file (caballes_length.csv)
caballes_length <- read.csv("../data/caballes_length.csv")
head(caballes_length,10)
set_sum_contrasts()
setting contr.sum globally: options(contrasts=c('contr.sum', 'contr.poly'))
Make female and jar factors
caballes_length$female <- factor(caballes_length$female)
caballes_length$jar <- factor(caballes_length$jar)
Reorder nutrition treatments so starved is first for default lm contrasts
caballes_length$diet <- factor(caballes_length$diet, levels = c("starved","aa","pr"))
Note: can’t get the aov commands to work for 3 level nested design
caballes.aov <- aov(length~diet+Error(female/jar), caballes_length)
Warning: Error() model is singularError in qr.qty(qr.e, resp) : NA/NaN/Inf in foreign function call (arg 1)
caballes.lm <- lm(length ~ diet/female/jar, caballes_length)
plot(caballes.lm)
summary(caballes.lm)
Call:
lm(formula = length ~ diet/female/jar, data = caballes_length)
Residuals:
Min 1Q Median 3Q Max
-0.3405 -0.0840 0.0065 0.0801 0.3631
Coefficients: (702 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.6899 0.0242 28.52 <2e-16 ***
diet1 -0.1025 0.0342 -3.00 0.003 **
diet2 0.0799 0.0342 2.34 0.020 *
dietstarved:female1 0.0244 0.0592 0.41 0.681
dietaa:female1 NA NA NA NA
dietpr:female1 0.0129 0.1026 0.13 0.900
dietstarved:female2 -0.1375 0.0592 -2.32 0.021 *
dietaa:female2 NA NA NA NA
dietpr:female2 NA NA NA NA
dietstarved:female3 NA NA NA NA
dietaa:female3 NA NA NA NA
dietpr:female3 NA NA NA NA
dietstarved:female4 NA NA NA NA
dietaa:female4 -0.0557 0.0592 -0.94 0.348
dietpr:female4 NA NA NA NA
dietstarved:female5 NA NA NA NA
dietaa:female5 -0.0326 0.0592 -0.55 0.583
dietpr:female5 NA NA NA NA
dietstarved:female6 NA NA NA NA
dietaa:female6 NA NA NA NA
dietpr:female6 NA NA NA NA
dietstarved:female7 NA NA NA NA
dietaa:female7 NA NA NA NA
dietpr:female7 0.0373 0.0592 0.63 0.530
dietstarved:female8 NA NA NA NA
dietaa:female8 NA NA NA NA
dietpr:female8 NA NA NA NA
dietstarved:female1:jar1 -0.0649 0.0592 -1.10 0.274
dietaa:female1:jar1 NA NA NA NA
dietpr:female1:jar1 NA NA NA NA
dietstarved:female2:jar1 NA NA NA NA
dietaa:female2:jar1 NA NA NA NA
dietpr:female2:jar1 NA NA NA NA
dietstarved:female3:jar1 NA NA NA NA
dietaa:female3:jar1 NA NA NA NA
dietpr:female3:jar1 NA NA NA NA
dietstarved:female4:jar1 NA NA NA NA
dietaa:female4:jar1 NA NA NA NA
dietpr:female4:jar1 NA NA NA NA
dietstarved:female5:jar1 NA NA NA NA
dietaa:female5:jar1 NA NA NA NA
dietpr:female5:jar1 NA NA NA NA
dietstarved:female6:jar1 NA NA NA NA
dietaa:female6:jar1 NA NA NA NA
dietpr:female6:jar1 NA NA NA NA
dietstarved:female7:jar1 NA NA NA NA
dietaa:female7:jar1 NA NA NA NA
dietpr:female7:jar1 NA NA NA NA
dietstarved:female8:jar1 NA NA NA NA
dietaa:female8:jar1 NA NA NA NA
dietpr:female8:jar1 NA NA NA NA
dietstarved:female9:jar1 NA NA NA NA
dietaa:female9:jar1 NA NA NA NA
dietpr:female9:jar1 -0.1506 0.1026 -1.47 0.144
dietstarved:female1:jar2 -0.1039 0.0592 -1.75 0.081 .
dietaa:female1:jar2 NA NA NA NA
dietpr:female1:jar2 NA NA NA NA
dietstarved:female2:jar2 NA NA NA NA
dietaa:female2:jar2 NA NA NA NA
dietpr:female2:jar2 NA NA NA NA
dietstarved:female3:jar2 NA NA NA NA
dietaa:female3:jar2 NA NA NA NA
dietpr:female3:jar2 NA NA NA NA
dietstarved:female4:jar2 NA NA NA NA
dietaa:female4:jar2 NA NA NA NA
dietpr:female4:jar2 NA NA NA NA
dietstarved:female5:jar2 NA NA NA NA
dietaa:female5:jar2 NA NA NA NA
dietpr:female5:jar2 NA NA NA NA
dietstarved:female6:jar2 NA NA NA NA
dietaa:female6:jar2 NA NA NA NA
dietpr:female6:jar2 NA NA NA NA
dietstarved:female7:jar2 NA NA NA NA
dietaa:female7:jar2 NA NA NA NA
dietpr:female7:jar2 NA NA NA NA
dietstarved:female8:jar2 NA NA NA NA
dietaa:female8:jar2 NA NA NA NA
dietpr:female8:jar2 NA NA NA NA
dietstarved:female9:jar2 NA NA NA NA
dietaa:female9:jar2 NA NA NA NA
dietpr:female9:jar2 NA NA NA NA
dietstarved:female1:jar3 NA NA NA NA
dietaa:female1:jar3 NA NA NA NA
dietpr:female1:jar3 NA NA NA NA
dietstarved:female2:jar3 NA NA NA NA
dietaa:female2:jar3 NA NA NA NA
dietpr:female2:jar3 NA NA NA NA
dietstarved:female3:jar3 NA NA NA NA
dietaa:female3:jar3 NA NA NA NA
dietpr:female3:jar3 NA NA NA NA
dietstarved:female4:jar3 NA NA NA NA
dietaa:female4:jar3 NA NA NA NA
dietpr:female4:jar3 NA NA NA NA
dietstarved:female5:jar3 NA NA NA NA
dietaa:female5:jar3 NA NA NA NA
dietpr:female5:jar3 NA NA NA NA
dietstarved:female6:jar3 NA NA NA NA
dietaa:female6:jar3 NA NA NA NA
dietpr:female6:jar3 NA NA NA NA
dietstarved:female7:jar3 NA NA NA NA
dietaa:female7:jar3 NA NA NA NA
dietpr:female7:jar3 NA NA NA NA
dietstarved:female8:jar3 NA NA NA NA
dietaa:female8:jar3 NA NA NA NA
dietpr:female8:jar3 NA NA NA NA
dietstarved:female9:jar3 NA NA NA NA
dietaa:female9:jar3 NA NA NA NA
dietpr:female9:jar3 NA NA NA NA
dietstarved:female1:jar4 NA NA NA NA
dietaa:female1:jar4 NA NA NA NA
dietpr:female1:jar4 NA NA NA NA
dietstarved:female2:jar4 0.0089 0.0592 0.15 0.881
dietaa:female2:jar4 NA NA NA NA
dietpr:female2:jar4 NA NA NA NA
dietstarved:female3:jar4 NA NA NA NA
dietaa:female3:jar4 NA NA NA NA
dietpr:female3:jar4 NA NA NA NA
dietstarved:female4:jar4 NA NA NA NA
dietaa:female4:jar4 NA NA NA NA
dietpr:female4:jar4 NA NA NA NA
dietstarved:female5:jar4 NA NA NA NA
dietaa:female5:jar4 NA NA NA NA
dietpr:female5:jar4 NA NA NA NA
dietstarved:female6:jar4 NA NA NA NA
dietaa:female6:jar4 NA NA NA NA
dietpr:female6:jar4 NA NA NA NA
dietstarved:female7:jar4 NA NA NA NA
dietaa:female7:jar4 NA NA NA NA
dietpr:female7:jar4 NA NA NA NA
dietstarved:female8:jar4 NA NA NA NA
dietaa:female8:jar4 NA NA NA NA
dietpr:female8:jar4 NA NA NA NA
dietstarved:female9:jar4 NA NA NA NA
dietaa:female9:jar4 NA NA NA NA
dietpr:female9:jar4 NA NA NA NA
dietstarved:female1:jar5 NA NA NA NA
dietaa:female1:jar5 NA NA NA NA
dietpr:female1:jar5 NA NA NA NA
dietstarved:female2:jar5 -0.0302 0.0592 -0.51 0.611
dietaa:female2:jar5 NA NA NA NA
dietpr:female2:jar5 NA NA NA NA
dietstarved:female3:jar5 NA NA NA NA
dietaa:female3:jar5 NA NA NA NA
dietpr:female3:jar5 NA NA NA NA
dietstarved:female4:jar5 NA NA NA NA
dietaa:female4:jar5 NA NA NA NA
dietpr:female4:jar5 NA NA NA NA
dietstarved:female5:jar5 NA NA NA NA
dietaa:female5:jar5 NA NA NA NA
dietpr:female5:jar5 NA NA NA NA
dietstarved:female6:jar5 NA NA NA NA
dietaa:female6:jar5 NA NA NA NA
dietpr:female6:jar5 NA NA NA NA
dietstarved:female7:jar5 NA NA NA NA
dietaa:female7:jar5 NA NA NA NA
dietpr:female7:jar5 NA NA NA NA
dietstarved:female8:jar5 NA NA NA NA
dietaa:female8:jar5 NA NA NA NA
dietpr:female8:jar5 NA NA NA NA
dietstarved:female9:jar5 NA NA NA NA
dietaa:female9:jar5 NA NA NA NA
dietpr:female9:jar5 NA NA NA NA
dietstarved:female1:jar6 NA NA NA NA
dietaa:female1:jar6 NA NA NA NA
dietpr:female1:jar6 NA NA NA NA
dietstarved:female2:jar6 NA NA NA NA
dietaa:female2:jar6 NA NA NA NA
dietpr:female2:jar6 NA NA NA NA
dietstarved:female3:jar6 NA NA NA NA
dietaa:female3:jar6 NA NA NA NA
dietpr:female3:jar6 NA NA NA NA
dietstarved:female4:jar6 NA NA NA NA
dietaa:female4:jar6 NA NA NA NA
dietpr:female4:jar6 NA NA NA NA
dietstarved:female5:jar6 NA NA NA NA
dietaa:female5:jar6 NA NA NA NA
dietpr:female5:jar6 NA NA NA NA
dietstarved:female6:jar6 NA NA NA NA
dietaa:female6:jar6 NA NA NA NA
dietpr:female6:jar6 NA NA NA NA
dietstarved:female7:jar6 NA NA NA NA
dietaa:female7:jar6 NA NA NA NA
dietpr:female7:jar6 NA NA NA NA
dietstarved:female8:jar6 NA NA NA NA
dietaa:female8:jar6 NA NA NA NA
dietpr:female8:jar6 NA NA NA NA
dietstarved:female9:jar6 NA NA NA NA
dietaa:female9:jar6 NA NA NA NA
dietpr:female9:jar6 NA NA NA NA
dietstarved:female1:jar7 NA NA NA NA
dietaa:female1:jar7 NA NA NA NA
dietpr:female1:jar7 NA NA NA NA
dietstarved:female2:jar7 NA NA NA NA
dietaa:female2:jar7 NA NA NA NA
dietpr:female2:jar7 NA NA NA NA
dietstarved:female3:jar7 -0.0364 0.0592 -0.61 0.540
dietaa:female3:jar7 NA NA NA NA
dietpr:female3:jar7 NA NA NA NA
dietstarved:female4:jar7 NA NA NA NA
dietaa:female4:jar7 NA NA NA NA
[ reached getOption("max.print") -- omitted 529 rows ]
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.132 on 243 degrees of freedom
Multiple R-squared: 0.446, Adjusted R-squared: 0.387
F-statistic: 7.52 on 26 and 243 DF, p-value: <2e-16
anova(caballes.lm)
Analysis of Variance Table
Response: length
Df Sum Sq Mean Sq F value Pr(>F)
diet 2 2.53 1.267 72.16 <2e-16 ***
diet:female 6 0.37 0.062 3.55 0.0022 **
diet:female:jar 18 0.53 0.029 1.67 0.0460 *
Residuals 243 4.27 0.018
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Get F and P values using correct denominators
#Diet F
f <- 1.26668/0.06224
pf(f, df1 = 2, df2 = 6, lower.tail = FALSE)
[1] 0.00212
#Females F
f <- 0.06224/0.02925
pf(f, df1 = 6, df2 = 18, lower.tail = FALSE)
[1] 0.1
#Jars F
f <- 0.02925/0.01755
pf(f, df1 = 18, df2 = 243, lower.tail = FALSE)
[1] 0.0459
Variance components from VCA
caballes.vca <- anovaMM(length ~ diet/(female)/(jar), caballes_length)
caballes.vca
ANOVA-Type Estimation of Mixed Model:
--------------------------------------
[Fixed Effects]
int dietaa dietpr dietstarved
0.527 0.219 0.189 0.000
[Variance Components]
Name DF SS MS VC %Total SD CV[%]
1 total 200.914152 0.019823 100 0.140793 21.242557
2 diet:female 6 0.37346 0.062243 0.0011 5.547827 0.033162 5.003435
3 diet:female:jar 18 0.526523 0.029251 0.00117 5.90134 0.034202 5.160385
4 error 243 4.265431 0.017553 0.017553 88.550832 0.132489 19.989554
Mean: 0.663 (N = 270)
Experimental Design: balanced | Method: ANOVA
VCAinference(caballes.vca, alpha = 0.05, VarVC = TRUE, ci.method = "satterthwaite")
Inference from Mixed Model Fit
------------------------------
> VCA Result:
-------------
[Fixed Effects]
int dietaa dietpr dietstarved
0.527 0.219 0.189 0.000
[Variance Components]
Name DF SS MS VC %Total SD CV[%] Var(VC)
1 total 200.9142 0.0198 100 0.1408 21.2426
2 diet:female 6 0.3735 0.0622 0.0011 5.5478 0.0332 5.0034 0
3 diet:female:jar 18 0.5265 0.0293 0.0012 5.9013 0.0342 5.1604 0
4 error 243 4.2654 0.0176 0.0176 88.5508 0.1325 19.9896 0
Mean: 0.663 (N = 270)
Experimental Design: balanced | Method: ANOVA
> VC:
-----
Estimate DF CI LCL CI UCL One-Sided LCL One-Sided UCL
total 0.0198 200.91 0.0165 0.0244 0.0169 0.0235
diet:female 0.0011 1.57 0.0003 0.1039 0.0003 0.0427
diet:female:jar 0.0012 2.80 0.0004 0.0188 0.0004 0.0112
error 0.0176 243.00 0.0148 0.0211 0.0152 0.0205
> SD:
-----
Estimate DF CI LCL CI UCL One-Sided LCL One-Sided UCL
total 0.1408 200.91 0.1283 0.156 0.1302 0.153
diet:female 0.0332 1.57 0.0164 0.322 0.0184 0.207
diet:female:jar 0.0342 2.80 0.0191 0.137 0.0210 0.106
error 0.1325 243.00 0.1217 0.145 0.1233 0.143
> CV[%]:
--------
Estimate DF CI LCL CI UCL One-Sided LCL One-Sided UCL
total 21.24 200.91 19.35 23.5 19.64 23.2
diet:female 5.00 1.57 2.47 48.6 2.77 31.2
diet:female:jar 5.16 2.80 2.88 20.7 3.16 16.0
error 19.99 243.00 18.36 21.9 18.61 21.6
95% Confidence Level
Satterthwaite methodology used for computing CIs
caballes.lmer <- lmer(length ~ diet + (1|female/jar), caballes_length)
summary(caballes.lmer)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: length ~ diet + (1 | female/jar)
Data: caballes_length
REML criterion at convergence: -289
Scaled residuals:
Min 1Q Median 3Q Max
-2.6901 -0.6966 0.0139 0.6414 2.5937
Random effects:
Groups Name Variance Std.Dev.
jar:female (Intercept) 0.00117 0.0342
female (Intercept) 0.00110 0.0332
Residual 0.01755 0.1325
Number of obs: 270, groups: jar:female, 27; female, 9
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.6628 0.0152 5.9999 43.65 9.7e-09 ***
diet1 -0.1359 0.0215 5.9999 -6.33 0.00073 ***
diet2 0.0830 0.0215 5.9999 3.87 0.00831 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) diet1
diet1 0.000
diet2 0.000 -0.500
Get F-ratio for diet test using lmerTest
anova(caballes.lmer, ddf = "Kenward-Roger")
Type III Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
diet 0.714 0.357 2 6 20.4 0.0021 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
CI on variance components (remembering to square CIs from lmer which are in SD units)
caballes.ci <- confint.merMod(caballes.lmer)
Computing profile confidence intervals ...
caballes.vc <- (caballes.ci)^2
print(caballes.vc)
2.5 % 97.5 %
.sig01 0.00000 0.00405
.sig02 0.00000 0.00318
.sigma 0.01477 0.02108
(Intercept) 0.40405 0.47600
diet1 0.03036 0.00951
diet2 0.00199 0.01474