Below are the solutions to these exercises on Two way ANOVA.

if (!require(car)){install.packages(car, dep=T)}

## Warning: package 'car' was built under R version 3.3.2

library(car) if (!require(ggplot2)){install.packages(ggplot2, dep=T)}

## Warning: package 'ggplot2' was built under R version 3.3.3

library(ggplot2) if (!require(dplyr)){install.packages(dplyr, dep=T)}

## Warning: package 'dplyr' was built under R version 3.3.3

library(dplyr) if (!require(lattice)){install.packages(lattice, dep=T)} library(lattice) ############### # # # Exercise 1 # # # ############### # Load data Chick <- read.csv(file.choose()) str(Chick)

## 'data.frame': 24 obs. of 3 variables: ## $ Food : Factor w/ 2 levels "Even","Patchy": 1 1 1 1 1 1 1 1 1 1 ... ## $ Pen : Factor w/ 4 levels "P1","P2","P3",..: 1 1 1 2 2 2 3 3 3 4 ... ## $ Cortisol: num 15.5 16.2 15 18.7 15 20.8 12 15 13.2 12.4 ...

```
Chick
```

## Food Pen Cortisol ## 1 Even P1 15.5 ## 2 Even P1 16.2 ## 3 Even P1 15.0 ## 4 Even P2 18.7 ## 5 Even P2 15.0 ## 6 Even P2 20.8 ## 7 Even P3 12.0 ## 8 Even P3 15.0 ## 9 Even P3 13.2 ## 10 Even P4 12.4 ## 11 Even P4 10.4 ## 12 Even P4 9.5 ## 13 Patchy P1 16.1 ## 14 Patchy P1 7.3 ## 15 Patchy P1 10.5 ## 16 Patchy P2 7.5 ## 17 Patchy P2 7.2 ## 18 Patchy P2 10.8 ## 19 Patchy P3 5.2 ## 20 Patchy P3 9.4 ## 21 Patchy P3 5.5 ## 22 Patchy P4 8.8 ## 23 Patchy P4 10.3 ## 24 Patchy P4 9.9

# see balance...... table(Chick$Food, Chick$Pen)

## ## P1 P2 P3 P4 ## Even 3 3 3 3 ## Patchy 3 3 3 3

############### # # # Exercise 2 # # # ############### #boxplot par(mfrow = c(2, 1)) boxplot(Cortisol ~ Food, data = Chick, xlab = "Food", ylab = "Cortisol") boxplot(Cortisol ~ Pen, data = Chick, xlab = "Pen", ylab = "Cortisol")

#outlier presence #histogram P1<- Chick[Chick$Pen == "P1", ] P2<- Chick[Chick$Pen == "P2", ] P3<- Chick[Chick$Pen == "P3", ] P4<- Chick[Chick$Pen == "P4", ] par(mfrow = c(2, 2)) hist(P1$Cortisol, xlim = c(0, 20), main = "", col = "Darkblue", xlab = "P1") hist(P2$Cortisol, xlim = c(0, 20), main = "", col = "Darkblue", xlab = "P2") hist(P3$Cortisol, xlim = c(0, 20), main = "", col = "Darkblue", xlab = "P3") hist(P4$Cortisol, xlim = c(0, 20), main = "", col = "Darkblue", xlab = "P4")

#coplot par(mfrow = c(1, 1)) coplot(Cortisol ~ Pen | Food, panel = panel.smooth, xlab = "Pen", ylab = "Cortisol", data = Chick)

############### # # # Exercise 3 # # # ############### #normality qqnorm(Chick$Cortisol) qqline(Chick$Cortisol)

#shapiro test shapiro.test(Chick$Cortisol)

## ## Shapiro-Wilk normality test ## ## data: Chick$Cortisol ## W = 0.9669, p-value = 0.5915

#OK #levene test leveneTest(Cortisol ~ Pen, data = Chick)

## Levene's Test for Homogeneity of Variance (center = median) ## Df F value Pr(>F) ## group 3 3.4458 0.03624 * ## 20 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

leveneTest(Cortisol ~ Food, data = Chick) #not OK?

## Levene's Test for Homogeneity of Variance (center = median) ## Df F value Pr(>F) ## group 1 0.104 0.7502 ## 22

############### # # # Exercise 4 # # # ############### #Interaction between explanatory variables interaction.plot(Chick$Food, Chick$Pen, Chick$Cortisol, col=c(2,3), xlab = "Pen", ylab = "Cortisol", trace.label = "Food")

#all factors must be significant #check explanatory variables interaction effect using two-way ANOVA One <- aov(Cortisol ~ Food * Pen, data = Chick) # Look at the results summary(One)

## Df Sum Sq Mean Sq F value Pr(>F) ## Food 1 177.13 177.13 32.655 3.19e-05 *** ## Pen 3 63.49 21.16 3.902 0.0288 * ## Food:Pen 3 59.50 19.83 3.656 0.0352 * ## Residuals 16 86.79 5.42 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# The results suggest an highly signficant food effect, a significant pen effect and a significant interaction effect. ############### # # # Exercise 5 # # # ############### #Using a nested analysis Two <- aov(Cortisol ~ Food + Error(Food / Pen), data = Chick) # Look at the results summary(Two)

## ## Error: Food ## Df Sum Sq Mean Sq ## Food 1 177.1 177.1 ## ## Error: Food:Pen ## Df Sum Sq Mean Sq F value Pr(>F) ## Residuals 6 123 20.5 ## ## Error: Within ## Df Sum Sq Mean Sq F value Pr(>F) ## Residuals 16 86.79 5.424

# There are no p values and F ratios # To find those you need to compare the nested models (between pen and food) and we'll do it manually ############### # # # Exercise 6 # # # ############### # STEP ONE - the food effect comparison # 1. You need the F statistic for the food effect so divide the mean sq on the top level # (Food) by the mean sq for the level below (Food : Pen) 177.1/20.5 # = 8.639024

## [1] 8.639024

# 2. Now compute the p value using the properties of the F distribution # Use the formula 1 - pf(F value, Df top level test, DF second level test) 1 - pf(8.63902, 1, 6) # gives a food effect of (p=0.02597392 or 0.026 if we round it)

## [1] 0.02597392

############### # # # Exercise 7 # # # ############### # STEP TWO - the pen effect. Divide the food:pen residuals by the error within residuals 20.5/5.424 # this is 3.779499

## [1] 3.779499

# Compute the p value 1 - pf(3.779499, 6, 16) # p=0.01549547 or p=0.015 if we round it)

## [1] 0.01549547

# Do the same thing using code forcing R to generate the residual comparisons # The food effect FoodEffect <- aov(Cortisol ~ Food + Error(Food : Pen), data = Chick)

## Warning in aov(Cortisol ~ Food + Error(Food:Pen), data = Chick): Error() ## model is singular

# look at results and ignore the warning summary(FoodEffect)

## ## Error: Food:Pen ## Df Sum Sq Mean Sq F value Pr(>F) ## Food 1 177.1 177.1 8.641 0.026 * ## Residuals 6 123.0 20.5 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Error: Within ## Df Sum Sq Mean Sq F value Pr(>F) ## Residuals 16 86.79 5.424

# Food effect is significant at p=0.026 ############### # # # Exercise 8 # # # ############### # now consider pen effect by comparing the un-nested and nested models using an ANOVA command PenEffect <- anova(aov(Cortisol ~ Food, data = Chick), aov(Cortisol ~ Food : Pen, data = Chick), test = "F") # Look at the results PenEffect # Notice we don't need the summary command because we used an ANOVA call to compare the models

## Analysis of Variance Table ## ## Model 1: Cortisol ~ Food ## Model 2: Cortisol ~ Food:Pen ## Res.Df RSS Df Sum of Sq F Pr(>F) ## 1 22 209.772 ## 2 16 86.787 6 122.98 3.7789 0.0155 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# So we can see that mean cortisol differed among by pens by treatment. Over and above this effect at this spatial scale # mean cortisol levels were effected by levels of food predictability. # Cortisol levels were lower on average in pens with patchy food than even distribution

**What's next:**

- Become a Top R Programmer Fast with our Individual Coaching Program
- Explore all our (>4000) R exercises
- Find an R course using our R Course Finder directory
- Subscribe to receive weekly updates and bonus sets by email
- Share with your friends and colleagues using the buttons below

Jeremy Kraft says

As per:

https://www.r-bloggers.com/raccoon-ch-2-5-unbalanced-and-nested-anova/

and

http://www.flutterbys.com.au/stats/tut/tut9.2a.html#h2_23

the correct aov call for nested model should be:

summary(aov(Cortisol~Food+Food:Pen, data=chick)) or

summary(aov(Cortisol~Food/Pen, data=chick))

and the results should not be much different between nested / non-nested model (ie summary(aov(Cortisol~Food*Pen, data=chick)) ).