Massey University — School of Mathematical and Computational Sciences

161.251 Regression Modelling

Computer Laboratory Exercise 2B Solutions

In this exercise you will:

• add the fitted line from a model to a scatter plot
• take the model fitted on the transformed scale back to the original scale and present it on a graph

Continuing the 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 km²; 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)

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. Reproduce a scatterplot of Population against Area using:

library(ggplot2)
Kittiwake.scatter <- ggplot(Kittiwake, aes(x=Area, y=Population)) + geom_point()
Kittiwake.scatter

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.

2. Create a new variable lPop by

library(tidyverse)
Kittiwake |> mutate(lPop = log(Population)) -> Kittiwake
glimpse(Kittiwake)

Now produce a scatterplot of the log of Population against Area by modifying the ggplot() command used previously. Store it using the name Kittiwake.scatter2.

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

3. Reproduce the simple linear regression model to the transformed data by

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

4. Extract the parameter estimates for the fitted model equation using:

coef(Kittiwake.lm)

 (Intercept)         Area 
6.0611805028 0.0009102822 

The fitted model is \(\mathbb{E}[lPop] = 6.0611805 + 9.1028218\times 10^{-4} Area\)

5. Look at the model specified slightly differently, using:

Kittiwake.lm1 <- lm(log(Population) ~ Area, data=Kittiwake)

and compare its coefficients.

coef(Kittiwake.lm1)

 (Intercept)         Area 
6.0611805028 0.0009102822 

6. We can update our graph using:

Kittiwake.scatter2 + geom_smooth(method="lm", se=FALSE)

`geom_smooth()` using formula = 'y ~ x'

Note that we “added” the line fitted by lm() on the fly here. The geom_smooth() didn’t use the model we fitted, but fitted its own one — the same one as it happens.

7. No one thinks of populations on a log scale.

\(\log{y} = a + bx\) (same form as our model) is equivalent to y = e^a e^bx where e is the exponential function that is the inverse of the natural logarithm. This is often re-written as y = Ae^bx where A=e^a.

This means our fitted line is equivalent to \[ \mathbb{E}[\mbox{Population}] = 428.9 e^{0.00091\mbox{Area}} \]