# Analysing the residuals from regression analysis on individual pines from the Arne 2019 data

```{r setup_arne2, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE)
```

## Introduction

In today's class you will learn how to run a regression analysis in R using the data collected from Arne as an example.

There are several handouts in both the online course book for this unit and AQM that cover the theory of regression in more depth. The purpose of this class is to provide the tools in R to run a regression and interpret the output.

### Load R libraries 

The libraries that are almost always used in any R analysis include ggplot2 (for graphics) dplyr (for data manipulation) and readr/readxl for data import.

```{r}
library(ggplot2)
library(aqm)
library(dplyr)
library(mgcv)

library(mapview)
library(sf)

```

### Load the data

When the data is loaded from the aqm package it also contains the gps way points that you collected.

```{r}
## This loads data from the aqm package
data(arne_2019)


```


```{r}
mapview(gps)
```

You can look at the points as a data frame




```{r}
gps %>% st_set_geometry(NULL) %>% dt()
```

Notice that the waypoints alone do not contain a variable that indicates the site. 


One way to obtain this (there are a lot of possibilities) is to cluster the points by proxmity.

```{r}
clust<-hclust(dist(st_coordinates(gps)))  ## Default option for a distance matrix uses euclidaen distance which is what we want.
plot(clust)
```

There are clearly two main clusters of points.

The cutree function in R cuts this tree into a given number of groups and returns an identifying number for the group. So we can now find the two sites easily and map the results out.


```{r}
gps$site<-as.factor(cutree(clust,2))
mapview(gps,zcol="site", burst=TRUE)
```

## Adding the GPS points to the pines

When we have two different data tables in R of any type we can always join them together if we can find a unique identifying feature that the two data sets have in common. 

Let's look at the names of the variables in the gps data


```{r}
str(gps)



```

The element that uniquely identifies the point is the name of the waypoint combined with the group that collected the data.

The pine_hts data has similar fields.


```{r}
str(pine_hts)
```
Notice that the Site column is currently blank.
We can make an id column for each set of data, then merge the results. 


```{r}
gps$id<-paste(gps$group,as.numeric(gps$name), sep="_")
pine_hts$id<-paste(pine_hts$Group,pine_hts$WP, sep="_")
pine_hts_gps<-merge(gps,pine_hts)
```



```{r}
mapview(pine_hts_gps,zcol="site")
```

## Looking at the residual

Let's make a working data frame to simplify typing.



```{r}
d<-pine_hts_gps

```

These data are the same as we used previously, so the scatterplot should look the same.

```{r}
fig1<-ggplot(d,aes(x=DiMM, y=HeCM)) +geom_point() + geom_smooth()
fig1
```

### Plotting the data by site

If we add a colour aesthetic to our ggplot then the data will be split into two labeled groups which correspond with the site.

```{r, message=FALSE}
fig2<-ggplot(d,aes(x=DiMM, y=HeCM, colour=site)) +geom_point() + geom_smooth()
fig2
```

What can you see from this figure?

### Obtaining residuals

Last time we tried fitting both lines and splines to the height diameter relationship. Regression lines were found to be useful for summarising the reslationship, but splines provided a closer fit to the actual data. So we can take the residuals from the spline.  



```{r}
d$residuals<-residuals(gam(data=d,HeCM~s(DiMM)))
```

This is the equivalent of measuring the distance between the fitted line shown in figure one and each data point. Some will be positive and others negative. Negative values imply that the pine height is lower than expected for a given diameter.


### Plotting the residuals by site

```{r}
fig3<-ggplot(d,aes(x=site,y=residuals)) 
fig3 + geom_boxplot()
```

### Confidence intervals

```{r}
aqm::ci(fig3)
```

### Simple unpaired t-test

```{r}
t.test(d$residuals~d$site)
```

