Computer Laboratory Exercise 2A Solutions

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)
            Area Population
W. Unst    207.7        311
Hermaness 1570.0       3872
N.E. Unst 1588.0        495
W. Yell    125.7        134
Buravoe    353.4        485
Fetlar     931.0        372
tail(Kittiwake)
             Area Population
Papa Stour  808.9       1036
Foula      2927.0       5570
Eshaness   1069.0       2430
Uyea        898.2        731
Gruney      564.8       1364
Fair Isle  3957.0      17000
  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)
Rows: 22
Columns: 3
$ Area       <dbl> 207.70, 1570.00, 1588.00, 125.70, 353.40, 931.00, 1616.00, …
$ Population <int> 311, 3872, 495, 134, 485, 372, 284, 10767, 1975, 970, 3243,…
$ lPop       <dbl> 5.739793, 8.261526, 6.204558, 4.897840, 6.184149, 5.918894,…

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.

Kittiwake.scatter2 <- ggplot(Kittiwake, aes(x=Area, y=lPop)) + geom_point() + labs(ylab="log(Population)")

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) 

Call:
lm(formula = lPop ~ 1 + Area, data = Kittiwake)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.88322 -0.60450 -0.01131  0.74978  2.02422 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.0611805  0.2902016  20.886 4.71e-15 ***
Area        0.0009103  0.0002151   4.233 0.000408 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9511 on 20 degrees of freedom
Multiple R-squared:  0.4725,    Adjusted R-squared:  0.4461 
F-statistic: 17.92 on 1 and 20 DF,  p-value: 0.0004082

to take a look at your results.

  1. Write down the fitted model equation.

Hence the fitted model is \(\mathbb{E}[lPop] = 6.061 + 0.001Area\)

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

The coefficient of determination is R2 = 47.3%. This is very respectable for the sort of observational study that we are looking at. Nonetheless, if you look at the scatterplot again (or at the model diagnostics later on) you should be able to see that the apparent strength of the linear relationship is driven quite heavily by two data points with large x-values (and hence high leverage).

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

If we denote the regression slope by \(\beta_1\), then a formal test of \(H_0:,\beta_1 = 0\) versus \(H_1:,\beta_1 \ne 0\) has the t-statistic 4.233 and corresponding P-value of P=0.000408, providing overwhelming evidence of a relationship between (log) population and colony area

  1. 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)
Rows: 21
Columns: 4
$ rowname    <chr> "W. Unst", "Hermaness", "N.E. Unst", "W. Yell", "Buravoe", …
$ Area       <dbl> 207.70, 1570.00, 1588.00, 125.70, 353.40, 931.00, 1616.00, …
$ Population <int> 311, 3872, 495, 134, 485, 372, 284, 1975, 970, 3243, 500, 2…
$ lPop       <dbl> 5.739793, 8.261526, 6.204558, 4.897840, 6.184149, 5.918894,…

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

Kittiwake2.lm <- lm(lPop ~ 1 + Area, data=Kittiwake2)
summary(Kittiwake2.lm)

Call:
lm(formula = lPop ~ 1 + Area, data = Kittiwake2)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7612 -0.5718  0.0719  0.7792  1.0518 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.0011972  0.2609201  23.000 2.47e-15 ***
Area        0.0008719  0.0001931   4.514 0.000237 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8513 on 19 degrees of freedom
Multiple R-squared:  0.5175,    Adjusted R-squared:  0.4921 
F-statistic: 20.38 on 1 and 19 DF,  p-value: 0.0002373

Hence the new fitted model is \({\mathbb{E}}[{\tt lPop}] = 6.001 + 0.001 {\tt Area}\)

This is little different from the original model, which should be of little surprise, since data point 8 has modest leverage (and Cook’s distance).