# rcode18

## Multiple explanatory variables: R code for Chapter 18 examples

Click on 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(, )`

Residual plot:
`plot( ~ )`

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

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. The data are from Example 17.3.

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

Take the log-transformation of stability.

```prairie\$logStability <- log(prairie\$biomassStability)
```

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

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

Fit the full model, which includes the treatment variable.

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

Scatter plot to compare models visually.

```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.

```anova(prairieRegression)
```

### Figure 18.1-2. Generalizing linear regression

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. The data are from Example 15.1.

```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.

```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.

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

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

```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.

```par(bty="l")
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)

```

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

```anova(circadianAnova)
```

### Example 18.2. Zooplankton depredation

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).

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

Set the order of groups in tables and plots.

```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

```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.

```library(nlme)
```

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

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

Fit the full model, which includes treatment.

```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.

```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.

```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.

```anova(zoopRBModel)
```

### Example 18.3 Intertidal interaction zone

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.

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

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

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

Fit the full model, with interaction term included.

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

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

```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")
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.

```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")
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.

```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.

```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`.

```anova(algaeFullModel)
```

A residual plot from the full model.

```plot( residuals(algaeFullModel) ~ fitted(algaeFullModel) )
abline(0,0)
```

A normal quantile plot of the residuals.

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

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.

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

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.

```moleRatSorted <- moleRat[ order(moleRat\$lnMass), ]
```

Scatter plot of the data.

```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.

```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.

```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.

```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.

```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.

```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()`.

```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`.

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

Residual plot from the linear model.

```plot( residuals(moleRatNoInteractModel) ~ fitted(moleRatNoInteractModel) )
abline(0,0)
```

Normal quantile plot of residuals.

```qqnorm(residuals(moleRatNoInteractModel), pch = 16, col = "firebrick",
las = 1, ylab = "Residuals", xlab = "Normal quantile", main = "")
```
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.`
The name of the data frame containing the variables named in the formula.
The first part of the formula refers to the residuals for each of the data points.
The second part of the formula refers to the fitted (predicted) values for each of the data points.