Computer Laboratory Exercise 2A

An R markdown template is available for this workshop. Copy and paste its content into a new R markdown file in RStudio.

In this exercise you will:

  • practice implementing linear regression modelling in R;
  • gain experience in the use of model diagnostics;
  • apply these techniques to the analysis of data on the sizes of kittiwake colonies.

Regression Analysis of Kittiwake Colony Data

This exercise is concerned with the relation between population size and foraging area for seabird colonies. Data are available on 22 black-legged kittiwake (a northern gull) colonies on Scotland’s Shetland and Orkney Islands. The variables recorded are the colony name; Area, in km2; and Population, the number of breeding pairs. (The data source is Cairns, D. K. (1988). “The regulation of seabird colony size: a hinterland model.” The American Naturalist, 134, 141-146.)

You should be able to read the data directly in R, and stored as the data frame Kittiwake, by

Kittiwake <- read.csv(file="https://r-resources.massey.ac.nz/data/161251/kittiwake.csv", header=TRUE, row.names=1)

Notice here that we are specifying that the first column of the file, which contains the colony names, should be used by R to provide names for the observations (rather than as a normal column of data).

If you Download kittiwake.csv then save it on your computer, use a command like:

## Kittiwake <- read.csv(file="<file path>/kittiwake.csv", header=TRUE, row.names=1)

where you should replace <file path> with the appropriate address corresponding to the location where you stored the file on your computer.

  1. Look at the data frame by simply typing its name, but also try:
head(Kittiwake)
tail(Kittiwake)
  1. Produce a scatterplot of Population against Area using:
library(ggplot2)
Kittiwake.scatter <- ggplot(Kittiwake, aes(x=Area, y=Population)) + geom_point()
Kittiwake.scatter

What are your initial impressions?

You should see that the scatterplot suggests a non-linear relationship between the variables, with the curve increasing rapidly for values of Area. We will therefore need to transform the data if we are to apply a simple linear regression model.

  1. Create a new variable lPop by
library(tidyverse)
Kittiwake |> mutate(lPop = log(Population)) -> Kittiwake
glimpse(Kittiwake)

This new variable contains the natural logarithms of Population. Now produce a scatterplot of the log of Population against Area by modifying the ggplot() command used previously.

You should see that the relationship Area and lPop looks reasonably linear.

  1. Fit a simple linear regression model to the transformed data by
Kittiwake.lm <- lm(lPop ~ 1 + Area, data=Kittiwake)

Then use

summary(Kittiwake.lm)

to take a look at your results.

  1. Write down the fitted model equation.

  2. What percentage of the variation in lPop is explained by Area?

  3. Is there statistically significant evidence that lPop is (linearly) related to Area? (Hint: think of an appropriate hypothesis test regarding the regression slope parameter.)

  4. Look at the standard diagnostic plots for your model by

plot(Kittiwake.lm)

Some people prefer to look at all four diagnostic plots simultaneously. This can be achieved by splitting up the graph window before issuing the plot() command:

par(mfrow=c(2,2))
plot(Kittiwake.lm)

or by using another package:

library(ggfortify)
autoplot(Kittiwake.lm)

N.B. You may need to reset the graphics window to a single plot by par(mfrow=c(1,1)).

Do you think there are any problems with your model?

You should note from the residual plots that observations 7 and 8 have the most extreme residuals. However, their standardized residual values are only a little larger than 2, which provides no particular cause for concern. It is worth noting that there are two data points with relatively high leverage. These data points have high Area values (take a look back at the scatterplot).

The model diagnostics do not give any real causes for concern regarding the validity of the assumptions underlying our model. (Nonetheless, you should appreciate that with only 22 data points, the residual plots are relatively crude diagnostic tools.) Nonetheless, as an exercise we will see what happens if we remove the data point with largest residual - data point 8, the Noss colony.

Create new Area and lPop variables with the Noss colony omitted by

Kittiwake2 = Kittiwake |> rownames_to_column() |> filter(rowname !="Noss") |> mutate(lPop = log(Population)) 
glimpse(Kittiwake2)

Fit a linear regression model using the new data. How much does the new fitted model differ from the original one?

Solutions

You should compare your work with the solutions for this workshop.