Back to Repository Home

Setting up ANCOVA in R

ANCOVA can involved mixed effects or not, but they are always nested designs. We simply collapse the design into an analysis that involves the two highest order levels.

d <- read.csv(file="Coral_interaction_experiment.csv", header=T, sep=",")
library(DT)
datatable(d)

PROTIP: You should only use aov for balanced designs. For unbalanced designs you should use lm, lme, or nlme procedures. As an alternative, you can also use the anova function and specify that you do NOT want to run a Type I (i.e., balanced design) ANOVA.


Conducting ANCOVA

Mixed effects specification using nlme

#Specifying and using mixed effects:

#use the nlme package for an unbalanced design:
library(nlme)

#this specifies that size is the outcome, predicted by site, with a random effect of transect)
#model_EXAMPLE <- lme(fixed=size~site, random=~1|tran, data=d)

#test these data using a similar approach
model2 <- lme(fixed=growth~treatment, 
              random=~1|tank, data=d)
summary(model2)
## Linear mixed-effects model fit by REML
##  Data: d 
##        AIC      BIC    logLik
##   1175.817 1188.067 -583.9085
## 
## Random effects:
##  Formula: ~1 | tank
##         (Intercept) Residual
## StdDev:    1.007869 9.434722
## 
## Fixed effects: growth ~ treatment 
##                Value Std.Error  DF   t-value p-value
## (Intercept) 149.7438  1.113396 151 134.49285       0
## treatment   -48.7800  1.491761 151 -32.69962       0
##  Correlation: 
##           (Intr)
## treatment -0.67 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.15108347 -0.63854855  0.01520533  0.60215427  2.78297564 
## 
## Number of Observations: 160
## Number of Groups: 8


Test for group normality

We can also test the normality of groups assumption using the shapiro test. You should run this test at the lowest level of grouping of your data. You can use the tapply function to run a test based on some grouping fairly quickly. 

#conduct a shapiro test on growth, grouped by the list of tank values.

tapply(d$growth, list(d$tank), shapiro.test)
## $`1`
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.94179, p-value = 0.2592
## 
## 
## $`2`
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.90415, p-value = 0.04936
## 
## 
## $`3`
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.8632, p-value = 0.008942
## 
## 
## $`4`
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.83855, p-value = 0.003447
## 
## 
## $`5`
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.89386, p-value = 0.03168
## 
## 
## $`6`
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.89912, p-value = 0.03971
## 
## 
## $`7`
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.86328, p-value = 0.008972
## 
## 
## $`8`
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.82189, p-value = 0.001867