fact2.crd <- read.csv("twofactorial.crd.csv")
head(fact2.crd)## Pesticide Variety Yield
## 1 P1 V1 49
## 2 P1 V1 39
## 3 P1 V2 55
## 4 P1 V2 41
## 5 P1 V3 66
## 6 P1 V3 68out1 <- aov(Yield ~ Pesticide*Variety,
data = fact2.crd) #Implements ANOVA
summary(out1) # Displays ANOVA table## Df Sum Sq Mean Sq F value Pr(>F)
## Pesticide 3 2227 742.5 17.556 0.00011 ***
## Variety 2 3996 1998.0 47.244 2.05e-06 ***
## Pesticide:Variety 6 457 76.2 1.801 0.18168
## Residuals 12 508 42.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(out1) # Alternative function to display ANOVA table## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Pesticide 3 2227.5 742.49 17.5563 0.0001098 ***
## Variety 2 3996.1 1998.04 47.2443 2.048e-06 ***
## Pesticide:Variety 6 456.9 76.15 1.8007 0.1816844
## Residuals 12 507.5 42.29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
To implement ANOVA for experiments with random effects, one can still use the aov() function but he/she has to recompute the F ratios according to the Expected Mean Squares.
anova(out1)## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Pesticide 3 2227.5 742.49 17.5563 0.0001098 ***
## Variety 2 3996.1 1998.04 47.2443 2.048e-06 ***
## Pesticide:Variety 6 456.9 76.15 1.8007 0.1816844
## Residuals 12 507.5 42.29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1From the above ANOVA table, if both Pesticides (Factor A) and Variety (Factor B) are random factors, then the correct F ratios for the main effects are: \[ F_A = \frac{742.49}{76.15} = 9.75 \\ F_B = \frac{1998.04}{76.15} = 26.24 \]
with associated p-values of:
\[ p_A = 0.01007 \\ p_B = 0.00108 \]
Alternatively, we can use the lmer() function in the lme4 package to analyze random effects models. By default, it will display the estimates of the variance components.
library(lme4)## Loading required package: Matrixout2 <- lmer(Yield ~ (1 | Pesticide) + (1 | Variety) + (1 | Pesticide:Variety),
data = fact2.crd)
summary(out2)## Linear mixed model fit by REML ['lmerMod']
## Formula: Yield ~ (1 | Pesticide) + (1 | Variety) + (1 | Pesticide:Variety)
## Data: fact2.crd
##
## REML criterion at convergence: 174.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3961 -0.5260 -0.1458 0.5919 1.6536
##
## Random effects:
## Groups Name Variance Std.Dev.
## Pesticide:Variety (Intercept) 16.93 4.115
## Pesticide (Intercept) 111.06 10.538
## Variety (Intercept) 240.24 15.500
## Residual 42.29 6.503
## Number of obs: 24, groups: Pesticide:Variety, 12; Pesticide, 4; Variety, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 61.46 10.54 5.833The lmer() function will not display the ANOVA table. But if we are really interested in testing the significance of the variance components,we can generate, say, the 95% confidence intervals and use these in testing hypotheses.
confint(out2, level = 0.95, oldNames = FALSE)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## sd_(Intercept)|Pesticide:Variety 0.000000 11.69133
## sd_(Intercept)|Pesticide 3.615636 26.93050
## sd_(Intercept)|Variety 6.268525 40.34079
## sigma 4.568343 10.14869
## (Intercept) 38.031705 84.88503In the above output, the 95% CI for the variance of the interaction term includes 0, hence, not it is not significant.
We can still use the aov() function to generate ANOVA for for experiments with mixed effects. However, just like for random effects model, we have to recompute some of the F ratios based on the expected mean squares. This is illustrated below.
anova(out1)## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Pesticide 3 2227.5 742.49 17.5563 0.0001098 ***
## Variety 2 3996.1 1998.04 47.2443 2.048e-06 ***
## Pesticide:Variety 6 456.9 76.15 1.8007 0.1816844
## Residuals 12 507.5 42.29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Recall that if A (Pesticide) is fixed and B (Variety) is random, the expected mean squares are: \[ \begin{aligned} E[MSE] &= \sigma_{\epsilon}^2 \\ E[MSAB] &= r\sigma_{\alpha\beta}^2 + \sigma_{\epsilon}^2 \\ E[MSB] &= ar \sigma_{\beta}^2 + \sigma_{\epsilon}^2 \\ E[MSA] &= br \frac{\sum_{i=1}^a\alpha_i^2}{a-1} + r\sigma_{\alpha\beta}^2 + \sigma_{\epsilon}^2 \end{aligned} \]
Consequently, the correct F ratios are:
\[ \begin{aligned} F_{AB} &= \frac{MSAB}{MSE} \\ F_B &= \frac{MSB}{MSE} \\ F_A &= \frac{MSA}{MSAB} \end{aligned} \]
Therefore, from the ANOVA table generated using aov() function, the F ratios for the interaction effect (AB) and the Variety effect (B) are already correct. We only need to recompute the F ratio for the fixed effect (A=Pesticide). It is given by
\[ F_A = \frac{MSA}{MSAB} = \frac{742.49}{76.15} = 9.7504 \]
with associated p-value given by
\[ p_A = 0.0101 \]
Therefore,
There is no significant interaction effect (p>0.05).
There is significant Variety effect (p<0.01).
There is significant Pesticide effect (p<0.01).
Alternatively, we can also use the lmer() function to analyze mixed effects models. Suppose Pesticide is fixed and Variety is random.
out3 <- lmer(Yield ~ Pesticide + (1 | Variety) + (1 | Pesticide:Variety),
data = fact2.crd)
summary(out3)## Linear mixed model fit by REML ['lmerMod']
## Formula: Yield ~ Pesticide + (1 | Variety) + (1 | Pesticide:Variety)
## Data: fact2.crd
##
## REML criterion at convergence: 150.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3220 -0.5453 -0.1250 0.6821 1.5452
##
## Random effects:
## Groups Name Variance Std.Dev.
## Pesticide:Variety (Intercept) 16.93 4.115
## Variety (Intercept) 240.24 15.500
## Residual 42.29 6.503
## Number of obs: 24, groups: Pesticide:Variety, 12; Variety, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 53.000 9.632 5.503
## PesticideP2 14.833 5.038 2.944
## PesticideP3 -1.833 5.038 -0.364
## PesticideP4 20.833 5.038 4.135
##
## Correlation of Fixed Effects:
## (Intr) PstcP2 PstcP3
## PesticideP2 -0.262
## PesticideP3 -0.262 0.500
## PesticideP4 -0.262 0.500 0.500
We can display the ANOVA table for the fixed effects using the anova() function, but there is no p-value.
anova(out3)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## Pesticide 3 1237 412.34 9.75The p-value can be generated using the lmerTest package.
library(lmerTest)
out4 <- lmer(Yield ~ Pesticide + (1 | Variety) + (1 | Pesticide:Variety),
data = fact2.crd)
anova(out4)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Pesticide 1237 412.34 3 6 9.75 0.01007 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1