---
title: "Multiple explanatory variables: R code for Chapter 18 examples"
author: "Michael Whitlock and Dolph Schluter"
output: 
  html_document:
  toc: yes
toc_depth: 3
---

_Note: This document was converted to R-Markdown from [this page](http://whitlockschluter.zoology.ubc.ca/r-code/rcode18) by M. Drew LaMar.  You can download the R-Markdown [here](https://qubeshub.org/collections/post/1400/download/chap18.Rmd)._

Download the R code on this page as a single file [here](http://whitlockschluter.zoology.ubc.ca/wp-content/rcode/chap18.r)

## New methods

Hover over a function argument for a short description of its meaning.  The variable names are plucked from the examples further below.

**Fit a linear model for two fixed factors**:

> algaeFullModel <- lm(<u><a title="A model formula. height * herbivores is shorthand for height + herbivores + height:herbivores, where the first two terms are the main effects and height:herbivores is the interaction.">sqrtArea ~ height * herbivores</a></u>, <u><a title="The name of the data frame containing the variables named in the formula.">data = algae</a></u>)

**Residual plot**:

> plot(<u><a title="The first part of the formula refers to the residuals for each of the data points.">residuals(moleRatNoInteractModel)</a></u> ~ <u><a title="The second part of the formula refers to the fitted (predicted) values for each of the data points.">fitted(moleRatNoInteractModel)</a></u>)

**Other new methods**:

 * Linear model for randomized block design.
 * Linear model for one categorical and one numeric variable (ANCOVA).
 * Interaction plot.
 * Visualize fit of null and full models to the data.
 * Type I and Type III sums of squares.

## Figure 18.1-1 [Modeling with linear regression](http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter17/chap17e3PlantDiversityAndStability.csv)

*__Compare the fits of the null and univariate regression models__ to data on the relationship between stability of plant biomass production and the initial number of plant species assigned to plots. <br \>The data are from Example 17.3.*

Read and inspect data.

```{r}
prairie <- read.csv(url("http://www.zoology.ubc.ca/~schluter/WhitlockSchluter/wp-content/data/chapter17/chap17e3PlantDiversityAndStability.csv"))
head(prairie)
```

**Take the log-transformation** of stability.

```{r}
prairie$logStability <- log(prairie$biomassStability)
head(prairie)
```

**Fit the null model to the data**, which in simple linear regression is a line of 0 slope.

```{r}
prairieNullModel <- lm(logStability ~ 1, data = prairie)
```

**Fit the full model**, which includes the treatment variable.

```{r}
prairieRegression <- lm(logStability ~ nSpecies, data = prairie)
```

**Scatter plot to compare models visually**.

```{r}
plot(logStability ~ nSpecies, data = prairie, bty = "l", col="firebrick", pch = 1, las = 1, cex = 1.5, xlim = c(0,16), xaxp = c(0,16,8), xlab = "Species number treatment", ylab = "Log-transformed ecosystem stability")
abline(prairieNullModel, lty = 2)
abline(prairieRegression)
```

**The F-test of improvement in fit** of the full (regression) model.

```{r}
anova(prairieRegression)
```

-----

## Figure 18.1-2. [Generalizing linear regression](http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter15/chap15e1KneesWhoSayNight.csv)

*__Compare the fits of the null and single-factor ANOVA model__ to data on phase shift in the circadian rhythm of melatonin production in participants given alternative light treatments.<br \>The data are from Example 15.1.*

Read and inspect data.

```{r}
circadian <- read.csv(url("http://www.zoology.ubc.ca/~schluter/WhitlockSchluter/wp-content/data/chapter15/chap15e1KneesWhoSayNight.csv"))
```

**Set the order of groups** for tables and graphs.

```{r}
circadian$treatment <- factor(circadian$treatment, levels = c("control", "knee","eyes")) 
```

**Fit the null model to the data**. In single-factor ANOVA, this fits a constant (grand mean) to all groups.

```{r}
circadianNullModel <- lm(shift ~ 1, data = circadian)
```

**Fit the full model to the data**, which includes the treatment effect.

```{r}
circadianAnova <- lm(shift ~ treatment, data = circadian)
```

**Strip chart to compare models visually**. The dashed line is the null model. The value `adjustAmount` is used to control the width of the horizontal line segments.

```{r}
par(bty="l")
adjustAmount <- 0.15
stripchart(shift ~ treatment, data = circadian, method = "jitter", vertical = TRUE, las = 1, pch = 1, xlab = "Light treatment", ylab = "Shift in circadian rhythm (h)", col = "firebrick", cex = 1.2)

xpts <- as.numeric(circadian$treatment)
ypts <- predict(circadianNullModel)
segments(xpts - adjustAmount, ypts, xpts + adjustAmount, ypts, lty = 2)

xpts <- as.numeric(circadian$treatment)
ypts <- predict(circadianAnova)
segments(xpts - adjustAmount, ypts, xpts + adjustAmount, ypts, col = "firebrick")
```

**$F$-test of improvement in fit** of the full model.

```{r}
anova(circadianAnova)
```

-----

## Example 18.2. [Zooplankton depredation](http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter18/chap18e2ZooplanktonDepredation.csv)

*__Analyze data from a randomized block experiment__, which measured the effects of fish abundance (the factor of interest) on zooplankton diversity. Treatments were repeated at 5 locations in a lake (blocks).*

Read and inspect data.

```{r}
zooplankton <- read.csv(url("http://www.zoology.ubc.ca/~schluter/WhitlockSchluter/wp-content/data/chapter18/chap18e2ZooplanktonDepredation.csv"))
head(zooplankton)
```

**Set the order of groups** in tables and plots.

```{r}
zooplankton$treatment <- factor(zooplankton$treatment, levels = c("control","low","high"))
```

**Table of the data** (Table 18.2-1). One measurement is available from each combination of treatment and block

```{r}
tapply(zooplankton$diversity, list(Treatment = zooplankton$treatment,  Location = zooplankton$block), unique)
```

A **blocking variable is typically analyzed as a random effect**. We need to load the `nlme` package.

```{r}
library(nlme)
```

**Fit the null model**, which includes block but leaves out the treatment variable.

```{r}
zoopNullModel <- lme(diversity ~ 1, random = ~ 1| block, data = zooplankton)
```

**Fit the full model**, which includes treatment.

```{r}
zoopRBModel <- lme(diversity ~ treatment, random = ~ 1| block, data = zooplankton)
```

**Visualize model fits, beginning with the null model**. The result here differs from that in the book, which shows block means. Lines are coincident when we analyze the data using `lme`, which estimates hardly any variance among the block means.

```{r}
interaction.plot(zooplankton$treatment, zooplankton$block, response = predict(zoopNullModel), ylim = range(zooplankton$diversity), trace.label = "Block", las = 1)

points(zooplankton$treatment, zooplankton$diversity, pch = as.numeric(zooplankton$block))
```

**Visualize model fits, continuing with the full model**. With the treatment term included, `lme` estimates that there is indeed variation among block means.

```{r}
interaction.plot(zooplankton$treatment, zooplankton$block, response = predict(zoopRBModel), ylim = range(zooplankton$diversity), trace.label = "Block", las = 1)

points(zooplankton$treatment, zooplankton$diversity, pch = as.numeric(zooplankton$block))
```

**$F$-test of improvement in fit** of the full model. This represents the test of treatment effect. Notice that R will not test the random (block) effect of an `lme` model.

```{r}
anova(zoopRBModel)
```

-----

## Example 18.3 [Intertidal interaction zone](http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter18/chap18e3IntertidalAlgae.csv)

*__Analyze data from a factorial experiment__ investigating the effects of herbivore presence, height above low tide, and the interaction between these factors, on abundance of a red intertidal alga using __two-factor ANOVA__.*

Read and inspect the data.

```{r}
algae <- read.csv(url("http://www.zoology.ubc.ca/~schluter/WhitlockSchluter/wp-content/data/chapter18/chap18e3IntertidalAlgae.csv"))
head(algae)
```

**Fit the null model having both main effects but no interaction term**. Note that we use `lm` because both factors are fixed effects.

```{r}
algaeNoInteractModel <- lm(sqrtArea ~ height + herbivores, data = algae)
```

**Fit the full model**, with interaction term included.

```{r}
algaeFullModel <- lm(sqrtArea ~ height * herbivores, data = algae)
```

**Visualize the model fits, beginning with the no-interaction model.**

```{r}
interaction.plot(algae$herbivores, algae$height, response = predict(algaeNoInteractModel), ylim = range(algae$sqrtArea), trace.label = "Height", las = 1, ylab = "Square root surface area (cm)", xlab = "Herbivore treatment")

adjustAmount = 0.05

points(sqrtArea ~ c(jitter(as.numeric(herbivores), factor = 0.2) + adjustAmount), data = subset(algae, height == "low"))

points(sqrtArea ~ c(jitter(as.numeric(herbivores), factor = 0.2) - adjustAmount), data = subset(algae, height == "mid"), pch = 16)
```

**Visualize the model fits, continuing with the full model including the interaction term.**

```{r}
interaction.plot(algae$herbivores, algae$height, response = predict(algaeFullModel), ylim = range(algae$sqrtArea), trace.label = "Height", las = 1, ylab = "Square root surface area (cm)", xlab = "Herbivore treatment")

adjustAmount = 0.05

points(sqrtArea ~ c(jitter(as.numeric(herbivores), factor = 0.2) + adjustAmount), data = subset(algae, height == "low"))

points(sqrtArea ~ c(jitter(as.numeric(herbivores), factor = 0.2) - adjustAmount), data = subset(algae, height == "mid"), pch = 16)
```

**Test the improvement in fit** of the model including the interaction term. 

```{r}
anova(algaeNoInteractModel, algaeFullModel)
```

**Test all terms in the model in a single ANOVA table**. Most commonly, this is done using either "Type III" sums of squares (see footnote 5 on p 618 of the book) or "Type I" sums of squares (which is the default in R). In the present example the two methods give the same answer, because the design is completely balanced, but this will not generally happen when the design is not balanced.

Here is how to test all model terms using **"Type III" sums of squares**. We need to include a `contrasts` argument for the two categorical variables in the `lm` command. Then we need to load the `car` package and use its `Anova` command. Note that "A" is in upper case in `Anova` - a very subtle difference.

```{r, warning=FALSE}
algaeFullModelTypeIII <- lm(sqrtArea ~ height * herbivores, data = algae, contrasts = list(height = contr.sum, herbivores = contr.sum))

library(car)

Anova(algaeFullModelTypeIII, type = "III") # note "A" in Anova is capitalized
```

Here is how we test all model terms using **"Type I" (sequential) sums of squares**. Note that "a" is in lower case in `anova`.

```{r}
anova(algaeFullModel)
```

**A residual plot** from the full model.

```{r}
plot( residuals(algaeFullModel) ~ fitted(algaeFullModel) )
abline(0,0)
```

**A normal quantile plot of the residuals**.

```{r}
qqnorm(residuals(algaeNoInteractModel), pch = 16, col = "firebrick", las = 1, ylab = "Residuals", xlab = "Normal quantile", main = "")
```

-----

## Example 18.4 [Mole-rat layabouts](http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter18/chap18e4MoleRatLayabouts.csv)

*__Analyze a factor while adjusting for a covariate__, comparing energy expenditure of two castes of naked mole-rat while adjusting for differences in body mass using __analysis of covariance ANCOVA__.*

Read and inspect the data.

```{r}
moleRat <- read.csv(url("http://www.zoology.ubc.ca/~schluter/WhitlockSchluter/wp-content/data/chapter18/chap18e4MoleRatLayabouts.csv"))
head(moleRat)
```

We are going to **sort the data** according to the value of the ln of body mass. This simplifies graphing of the model fits. The graph commands below assume that the data are sorted in this way.

```{r}
moleRatSorted <- moleRat[ order(moleRat$lnMass), ]
```

**Scatter plot** of the data.

```{r}
plot(lnEnergy ~ lnMass, data = moleRat, type = "n", las = 1, bty = "l")

points(lnEnergy ~ lnMass, data = subset(moleRatSorted, caste == "worker"), pch = 1, col = "firebrick")
points(lnEnergy ~ lnMass, data = subset(moleRatSorted, caste == "lazy"), pch = 16, col = "firebrick")
```

Fit models to the data, **beginning with the model lacking an interaction term**. Use `lm` because caste and mass are fixed effects. Save the predicted values in the data frame.

```{r}
moleRatNoInteractModel <- lm(lnEnergy ~ lnMass + caste, data = moleRatSorted)
moleRatSorted$predictedNoInteract <- predict(moleRatNoInteractModel)
```

**Fit the full model, which includes the interaction term**. Again, save the predicted values in the data frame.

```{r}
moleRatFullModel <- lm(lnEnergy ~ lnMass * caste, data = moleRatSorted)
moleRatSorted$predictedInteract <- predict(moleRatFullModel)
```

**Visualize the model fits**, beginning with the fit of the no-interaction model. Redraw the scatter plot (see commands above), if necessary, before issuing the following commands.

```{r}
# BEGIN: Redrawing scatter plot
plot(lnEnergy ~ lnMass, data = moleRat, type = "n", las = 1, bty = "l")

points(lnEnergy ~ lnMass, data = subset(moleRatSorted, caste == "worker"), pch = 1, col = "firebrick")
points(lnEnergy ~ lnMass, data = subset(moleRatSorted, caste == "lazy"), pch = 16, col = "firebrick")
# END: Redrawing scatter plot

lines(predictedNoInteract ~ lnMass, data = subset(moleRatSorted, caste == "worker"), lwd = 1.5)
lines(predictedNoInteract ~ lnMass, data = subset(moleRatSorted, caste == "lazy"), lwd = 1.5)
```

**Visualize the fit of the full model**, including the interaction term. Redraw the scatter plot, if necessary, before issuing the following commands.

```{r}
# BEGIN: Redrawing scatter plot
plot(lnEnergy ~ lnMass, data = moleRat, type = "n", las = 1, bty = "l")

points(lnEnergy ~ lnMass, data = subset(moleRatSorted, caste == "worker"), pch = 1, col = "firebrick")
points(lnEnergy ~ lnMass, data = subset(moleRatSorted, caste == "lazy"), pch = 16, col = "firebrick")
# END: Redrawing scatter plot

lines(predictedInteract ~ lnMass, data = subset(moleRatSorted, caste == "worker"), lwd = 1.5, lty = 2)
lines(predictedInteract ~ lnMass, data = subset(moleRatSorted, caste == "lazy"), lwd = 1.5, lty = 2)
```

**Test the improvement in fit of the full model**, including the interaction term. This is a test of the interaction term only.

```{r}
anova(moleRatNoInteractModel, moleRatFullModel)
```

**Test for differences** in ln body mass between castes, **assuming that no interaction term** is present in the mole rat population (i.e., assuming that the two castes hve equal slopes). Most commonly this is done using either "Type III" sums of squares (see footnote 5 on p 618 of the book) or "Type I" sums of squares (the default in R). The two methods do not give identical answers here because the design is not balanced (in a balanced design, each value of the x-variable would have the same number of y-observations from each group).

Test using **"Type III" sums of squares**. We need to include a `contrasts` argument for the categorical variable in the `lm` command. Then we need to load the `car` package and use its `Anova` command. Note that "A" is in upper case in `Anova()`.

```{r}
moleRatNoInteractModelTypeIII <- lm(lnEnergy ~ lnMass + caste, data = moleRat, contrasts = list(caste = contr.sum))

library(car)
Anova(moleRatNoInteractModelTypeIII, type = "III") # note "A" in Anova is capitalized
```

Test using **"Type I" (sequential) sums of squares**. Make sure that the covariate (`lnMass`) comes before the factor (`caste`) in the `lm` formula, as shown. Note that "a" is in lower case in `anova`.

```{r}
moleRatNoInteractModel <- lm(lnEnergy ~ lnMass + caste, data = moleRat)
anova(moleRatNoInteractModel)
```

**Residual plot** from the linear model. 

```{r}
plot( residuals(moleRatNoInteractModel) ~ fitted(moleRatNoInteractModel) )
abline(0,0)
```

**Normal quantile plot of residuals.**

```{r}
qqnorm(residuals(moleRatNoInteractModel), pch = 16, col = "firebrick", las = 1, ylab = "Residuals", xlab = "Normal quantile", main = "")
```