---
title: Massey University --- School of Mathematical and Computational Sciences  
subtitle: 161.251 Regression Modelling
author: "Staff member responsible for this workshop: Jonathan Godfrey a.j.godfrey@massey.ac.nz"
date: "last updated 11 July 2025"
output:
    html_document
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(comment="", dev="png")
library(tidyverse)
```


N.B. Please note the following:

- Read the full workshop before using this template. External images are not included in your template.
- The above chunk is setting some global options for the rest of this document. It isn't strictly necessary.
- The template worked on my computer, but I have access to all external files, including the data and images.
- Before you proceed, test this template works for you.


# Computer Laboratory  Exercise 2A 

An [R markdown template](https://R-Resources.massey.ac.nz/161251/labs/Lab2A.Rmd) 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 km^2^; 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

```{r}
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  `<a href="data:text/csv;base64,IkNvbG9ueSIsIkFyZWEiLCJQb3B1bGF0aW9uIg0KIlcuIFVuc3QiLDIwNy43LDMxMQ0KIkhlcm1hbmVzcyIsMTU3MCwzODcyDQoiTi5FLiBVbnN0IiwxNTg4LDQ5NQ0KIlcuIFllbGwiLDEyNS43LDEzNA0KIkJ1cmF2b2UiLDM1My40LDQ4NQ0KIkZldGxhciIsOTMxLDM3Mg0KIk91dCBTa2VycmllcyIsMTYxNiwyODQNCiJOb3NzIiwxMzE3LDEwNzY3DQoiTW91c3NhIiw2MTQsMTk3NQ0KIkRhbHNldHRlciIsNjAuMTIsOTcwDQoiU3VtYnVyZ2giLDEyNzMsMzI0Mw0KIkZpdGZ1bCBIZWFkIiw1OTUuNyw1MDANCiJTdCBOaW5pYW5zIiwxMDUuNiwyNTANCiJTLiBIYXZyYSIsMjQxLjksOTI1DQoiUmVhd2ljayIsMTExLjEsOTcwDQoiVmFpbGEiLDMwMi40LDI3OA0KIlBhcGEgU3RvdXIiLDgwOC45LDEwMzYNCiJGb3VsYSIsMjkyNyw1NTcwDQoiRXNoYW5lc3MiLDEwNjksMjQzMA0KIlV5ZWEiLDg5OC4yLDczMQ0KIkdydW5leSIsNTY0LjgsMTM2NA0KIkZhaXIgSXNsZSIsMzk1NywxNzAwMA0K" download="kittiwake.csv">Download kittiwake.csv</a>`{r} then save it on your
computer, use a command like: 


``` r
## 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:

``` r
head(Kittiwake)
tail(Kittiwake)
```

2. Produce a scatterplot of `Population` against `Area` using:


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

3. Create a new variable `lPop` by


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

4. Fit a simple linear regression model to the transformed data by

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

Then use

``` r
summary(Kittiwake.lm) 
```
to take a look at your results. 

5. Write down the fitted model equation.


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


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


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


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


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

or by using another package:


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


``` r
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.](https://R-Resources.massey.ac.nz/161251/labs/Lab2A-sols.html)