---
title: "Lab 4: Data Wrangling and Linear Models"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(pdxTrees)
library(tidyverse)
library(moderndive)
```

# Due: 12am night before lab

# Instructions
Work through the problems below, and submit this document as a knitted .pdf file with answers to the Lab 4 assignment on your lab section gradescope. 

**In some problems, an incomplete R chunk has been provided. To evaluate the chunk, follow the given instructions, and then change `eval=FALSE` to `eval=TRUE`.

**For each problem, put your solution between the bars of red stars.**


# Background
For this lab set, we will first run through an example of a simple linear regression, getting some practice with data wrangling along the way. Then we will work through a regression analysis independently.


Run the code chunk below to create the `near_reed` dataset from the `pdxTrees` package.

```{r}
# Grab trees near Reed

pdxTrees_parks <- get_pdxTrees_parks()
near_reed <- filter(pdxTrees_parks, Park %in% c("Crystal Springs Rhododendron Garden",
                                          "Kenilworth Park", 
                                          "Eastmoreland Garden",
                                          "Berkeley Park"))
#Remove trees with no functional type
near_reed <- drop_na(near_reed, Functional_Type)
```



# An Example/Demo
## Visualization
We will start by investigating the relationship between the diameter of a tree's trunk and its height.

We will do this by creating a scatterplot with `Tree_Height` as the outcome variable on the y-axis and `DBH` as the explanatory variable on the x-axis.

Note the use of the `geom_smooth()` function, that tells R to add a regression line. While the points do scatter/vary around the blue regression line, of all possible lines we can draw in this point of clouds, the blue line is the “best-fitting” line in that in minimizes the sum of the squared residuals (see moderndive 6.3.3 for a fuller explanation).


```{r}
ggplot(data = near_reed, aes(x = DBH, y = Tree_Height)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE) +
  labs(x = "Diameter at Breast Height (inches)", 
       y = "Height (feet)", 
       title = "Tree Height and Trunk Diameter") 
```

### Exercise 1

Comment on the direction, strength and linearity of the relationship shown in the preceding plot.

********************************



********************************


## The correlation coefficient (R)
We can numerically quantify the strength of the linear relationship between the two variables with the correlation coefficient. The following tells R to `summarize()` the correlation coefficient between the numerical variables `DBH` and `Tree_Height`. Note that the correlation coefficient only exists for pairs of numerical variables.


```{r}
near_reed %>%
  summarise(r = cor(DBH, Tree_Height))
```


## Running a linear regression model

In R we can fit a linear regression model (a regression line), like so:

```{r}
tree_mod <- lm(Tree_Height ~ DBH, data = near_reed)
```

Note that:

- the function `lm()` is short for “linear model”
- the first argument is a formula in the form `y ~ x` or in other words `outcome variable ~ explantory variable`.
- the second argument is the data frame in which the outcome and explanatory variables can be found.
- we SAVED THE MODEL RESULTS as an object called `tree_mod`

This object `tree_mod` contains all of the information we need about the linear model that was just fit and we'll be accessing this information again later.

## Get the regression table
The `get_regression_table()` function from the moderndive package will output a regression table. Let's focus on the value in the second column: an estimate for 1) an intercept, and 2) a slope for the `Tree_Height` variable. We'll revisit what the other columns mean in a future problem set.

```{r}
get_regression_table(tree_mod)
```


We can interpret the `intercept` and `poverty` slope like so:

- When DBH is 0, the predicted average tree height 7.938 feet
- For every increase DBH by 1 inch, there is an **associated increase** in tree height of 2.633 feet.

 

### Exercise 2

We found a positive correlation coefficient. Is it reasonable for us to conclude that increasing a tree's diameter will **cause** an increase in tree height? Explain why or why not?

***********************************************



***********************************************

## Making predictions
Based on the R output of our model, the following is our least squares regression line for the linear model. Let $Y$ denote the `Tree_Height` variable and $X$ denote the `DBH` variable:

$$\hat{Y}=7.938 + 2.633 X$$

We can use the line from our graph of the `Tree_Height`, `DBH` relationship to visually make predictions. For example, when DBH is 15 inches, the line shows a value of just over 47 feet.

To get a more accurate prediction, we could actually plug 15 into the regression equation like so:

```{r}
y_hat = 7.938 + 2.633 * 15
y_hat
```


### Exercise 3

What height would you expect for a tree with trunk diameter of 20 inches?

***********************************************



***********************************************


### Exercise 4

Do you think it is a good idea to predict Tree height at 80 inch diameter, based on this regression equation? Explain your reasoning.

***********************************************


***********************************************

## Residuals

Recall that model residuals $(E)$ are the difference between the observed values in your data set and the values predicted by the line:

$$E = Y - \hat{Y}$$ 

For instance, below, one data point is highlighted in red...the residual is the difference between the y value of the data point (here 125), and the y value predicted by the line (roughly 165). Here the residual is roughly -40 (125-165=-40). The regression equation has over-estimated tree height support, compared to this data point.

```{r}
ggplot(data = near_reed, aes(x=DBH, y = Tree_Height)) + 
  geom_point() + 
  geom_point(x = 60, y = 125, color = "red", size = 3)+
  geom_segment(x = 14, y = 59.5, xend = 14, yend = 69, color = "black", size =.35, linetype = "dotted")+
  geom_smooth(method = "lm", se = FALSE) +
  labs(x = "Diameter at Breast Height (inches)", 
       y = "Height (feet)", 
       title = "Tree Height and Trunk Diameter") 
```


The following function `get_regression_points()` gives you the *fitted* (also known as *predicted*) value for every data point, and the residual for every data point. The first row in the output is the first data point, you see that Tree Height was 110, DBH was 39.3, the regression equation predicted 111.396 tree height, and the residual was -1.396 (110 - 111.396).

```{r}
get_regression_points(tree_mod)
```

## Residual Plots

It is often easier to determine the size of residuals by way of a residual plot, which displays points with the same horizontal value as in the scatterplot, but with vertical values equal to the correspond residual. A residual plot can be created in ``ggplot` using the data from the `get_regression_points` function.

```{r}
get_regression_points(tree_mod) %>% 
  ggplot( aes(x =DBH, y = residual) )+
    geom_point() +
    geom_smooth(method = "lm", se = FALSE)
```

To assess the distribution of residuals, we can use a histogram:

```{r}
get_regression_points(tree_mod) %>% 
ggplot(aes(x = residual)) +
  geom_histogram(bins = 30, color = "white")
```

# Independent Practice


### Exercise 5

The crown of a tree refers to its total aboveground parts, including stems, branches, and leaves. In the `near_reed` dataframe, `Crown_Width_EW` represents the greatest width of the crown in the east-west direction, while `Crown_Width_NS` represents the greatest width of the crown in the north-south direction. Add a new variable to the `near_reed` dataframe called `Crown_Width` which is the average of `Crown_Width_EW` and `Crown_Width_NS`.

********************************




********************************

### Exercise 6

Generate a scatterplot with `Crown_Width` as the explanatory variable and `Tree_Height` as the response. Be sure to include an informative title and axis labels to your plot. This will help contextualize it.

********************************




********************************


### Exercise 7

D you expect the correlation coefficient (for `Crown_Width` and `Tree_Height`) to be positive or negative? Write code chunk testing your prediction.

***********************************************




***********************************************

### Exercise 8

Using the `near_reed` data frame, run a linear regression using `Crown_Width` as the explanatory variable and `Tree_Height` as the outcome variable. Interpret the Intercept and slope estimates.

***********************************************





***********************************************

### Exercise 9

Make a numerical prediction for the `Tree_Height` for a tree that has 40 foot crown width, using the formula for the linear model you created in the previous exercise.

***********************************************




***********************************************


### Exercise 10

Create a residual plot with `Tree_Height` modeled as a function of `Crown_Width`. Do the residuals appear to have relatively constant variability?

***********************************************




***********************************************


### Exercise 11

Create a histogram of the residual values. Does the distribution appear to be bell-shaped, unimodal, symmetric and centered at $0$?

***********************************************



***********************************************


### Exercise 12

Based on the evidence you have so far (scatterplots and correlation coefficients), which of the explanatory variables we have considered (`DBH` or `Crown_Width`) seems to be a better predictor `Tree_Height`? Explain.

********************************



********************************



# Additional `dplyr` Exercises

### Exercise 13

Functional type describes what type of leaves a tree has.  Unfortunately, the categories of the `Functional_Type` variable are not the most informative.  Change "BE" to "Broadleaf Evergreen", "BD" to "Broadleaf Deciduous.", "CE" to "Coniferous Evergreen" and "CD" to "Coniferous Deciduous.".  Save the results as a dataframe named `near_reed`.

********************************



********************************



## Exercise 14

Create a scatterplot of `Tree_Height` as a function of `Crown_Width` and color by `Functional_Type`. What trends are evident in this plot that were not evident in the plot in Exercise 7?

********************************



********************************


### Exercise 15

How many, and what proportion of trees are there for each functional type?

********************************





********************************


### Exercise 16

Create a variable `Functional_Type_2` that classifies a tree either as `Deciduous` or `Evergreen`. Save this data frame as `near_reed`

********************************




********************************

### Exercise 17

Create two new linear models for `Tree_Height` as a function of `Crown_Width`, one that just uses `Deciduous` trees, and one that just uses `Evergreen` trees. Be sure to include the regression table for each model.

********************************


********************************


### Exercise 18

Compute the correlation coefficients for `Tree_Height` and `Crown_Width` for each model. How do the slopes, intercepts and correlation coefficients for the two linear models compare (a) to each other, and (b) to those from Exercises 7 and 8.

********************************



********************************
