Setup

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
library(ggfortify)
library(dplyr)

Part 1 - Two-Way ANOVA

A. Import and Inspect the Data

growth.moo <- read.csv("growth.csv")

glimpse(growth.moo)
## Rows: 48
## Columns: 3
## $ supplement <chr> "supergain", "supergain", "supergain", "supergain", "contro…
## $ diet       <chr> "wheat", "wheat", "wheat", "wheat", "wheat", "wheat", "whea…
## $ gain       <dbl> 17.37125, 16.81489, 18.08184, 15.78175, 17.70656, 18.22717,…
levels(growth.moo$diet)
## NULL
levels(growth.moo$supplement)
## NULL

The response variable is cow and weight gain, while diet and supplement are categorical explanatory variables. Meaning, two-way ANOVA is appropriate.

B. Set Reference Level

growth.moo <- growth.moo %>%
  mutate(
    supplement = factor(supplement),
    supplement = relevel(supplement, ref = "control")
  )

C. Calculate Summary Statistics

sumMoo <- growth.moo %>%
  group_by(diet, supplement) %>%
  summarise(
    meanGrow = mean(gain),
    seGrow = sd(gain) / sqrt(n()),
    .groups = "drop"
  )

D. Make an interaction plot

ggplot(sumMoo,
       aes(x = supplement,
           y = meanGrow,
           colour = diet,
           group = diet)) +
  geom_point() +
  geom_line() +
  theme_bw()

If the lines are parallel, there is little evidence of an interaction. If they differ in shape or cross, this suggests that the effect of supplement may depend on diet.

E. Build Two-Way ANOVA

model_cow <- lm(gain ~ diet * supplement,
                data = growth.moo)

diet * supplement = diet + supplement + diet:supplement

F. Check Assumptions

autoplot(model_cow, smooth.colour = NA)
## Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 48 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 48 rows containing missing values or values outside the scale range
## (`geom_line()`).

G. Examine ANOVA table

anova(model_cow)
## Analysis of Variance Table
## 
## Response: gain
##                 Df  Sum Sq Mean Sq F value    Pr(>F)    
## diet             2 287.171 143.586 83.5201 2.999e-14 ***
## supplement       3  91.881  30.627 17.8150 2.952e-07 ***
## diet:supplement  6   3.406   0.568  0.3302    0.9166    
## Residuals       36  61.890   1.719                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Part 2 - ANCOVA

A. Import data and inspect

limp <- read.csv("limpet.csv")

glimpse(limp)
## Rows: 24
## Columns: 3
## $ DENSITY <int> 8, 8, 8, 8, 8, 8, 15, 15, 15, 15, 15, 15, 30, 30, 30, 30, 30, …
## $ SEASON  <chr> "spring", "spring", "spring", "summer", "summer", "summer", "s…
## $ EGGS    <dbl> 2.875, 2.625, 1.750, 2.125, 1.500, 1.875, 2.600, 1.866, 2.066,…

Egg production is the response variable, density is the continous explanatory variable, and season is a categorical explanatory variable.

B. Plot the data first

ggplot(limp,
       aes(x = DENSITY,
           y = EGGS,
           colour = SEASON)) +
  geom_point() +
  theme_bw()

C. Fit the ANCOVA

limp.mod <- lm(EGGS ~ DENSITY * SEASON,
               data = limp)

Density * Season = Density + Season + Density:Season

D. Check Assumptions

autoplot(limp.mod, smooth.colour = NA)
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).

E. ANOVA table

anova(limp.mod)
## Analysis of Variance Table
## 
## Response: EGGS
##                Df Sum Sq Mean Sq F value    Pr(>F)    
## DENSITY         1 5.0241  5.0241 30.1971 2.226e-05 ***
## SEASON          1 3.2502  3.2502 19.5350 0.0002637 ***
## DENSITY:SEASON  1 0.0118  0.0118  0.0711 0.7925333    
## Residuals      20 3.3275  0.1664                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Egg production changes significantly with density and differs significantly between seasons. There is no significant density x season interaction, meaning the relationship between density and egg production does not alter significantly between seasons.

Part 3 - Final Plots

ggplot(sumMoo,
       aes(x = supplement,
           y = meanGrow,
           colour = diet,
           group = diet)) +
  geom_point() +
  geom_line() +
  geom_errorbar(
    aes(ymin = meanGrow - seGrow,
        ymax = meanGrow + seGrow),
    width = 0.1
  ) +
  theme_bw()