lec07.Rmd

---
title: "Practice Lecture 7 MATH 342W Queens College"
author: "Professor Adam Kapelner"
---

## Simple Linear Regression (p = 1)

To understand what the algorithm is doing - best linear fit by minimizing the squared errors, we can draw a picture. First let's make up some very simple training data $\mathbb{D}$.

```{r}
set.seed(1984)
n = 20
x = runif(n)
beta_0 = 3
beta_1 = -2
h_star = beta_0 + beta_1 * x
epsilons = rnorm(n, mean = 0, sd = 0.33)
y = h_star + epsilons 
```

And let's plot the data:


```{r}
pacman::p_load(ggplot2)
simple_df = data.frame(x = x, y = y)
simple_viz_obj = ggplot(simple_df, aes(x, y)) + 
  geom_point(size = 2)
simple_viz_obj
```

And its true $h^*$ line:

```{r}
true_hstar_line = geom_abline(intercept = beta_0, slope = beta_1, color = "green")
simple_viz_obj + true_hstar_line
```

Now let's calculate the simple least squares coefficients:

```{r}
r = cor(x, y)
s_x = sd(x)
s_y = sd(y)
ybar = mean(y)
xbar = mean(x)

b_1 = r * s_y / s_x
b_0 = ybar - b_1 * xbar
b_0
b_1
```

Note how $b_0$ and $b_1$ are not exactly the same as $\beta_0$ and $\beta_1$. Why?

And we can plot it:


```{r}
simple_ls_regression_line = geom_abline(intercept = b_0, slope = b_1, color = "red")
simple_viz_obj + simple_ls_regression_line + true_hstar_line
```

Review of the modeling framework:

The difference between the green line and red line is the "estimation error". The difference between the green line and the points is a combination of error due to ignorance and error due to misspecification of $f$ as a straight line. In most real-world applications, estimation error is usually small relative to the other two. In the era of "big data", $n$ is usually big so estimation error is pretty small.

Recall that the noise (epsilons) are the difference between the data and the green line:

```{r}
simple_df$hstar = beta_0 + beta_1 * simple_df$x
simple_viz_obj = ggplot(simple_df, aes(x, y)) + 
  geom_point(size = 2)
epsilon_line_segments = geom_segment(aes(xend = x, yend = hstar), position = position_nudge(x = 0.002))
simple_viz_obj + epsilon_line_segments + true_hstar_line
```

And that the residuals (e's) are the difference between the measurements of the response in the actual data and the green line:

```{r}
simple_df$gs = b_0 + b_1 * simple_df$x
simple_viz_obj = ggplot(simple_df, aes(x, y)) + 
  geom_point(size = 2)
e_line_segments = geom_segment(aes(xend = x, yend = gs), color = "purple")
simple_viz_obj + simple_ls_regression_line + e_line_segments
```

Examining both at the same time:

```{r}
simple_viz_obj + simple_ls_regression_line + true_hstar_line + e_line_segments + epsilon_line_segments
```


## Assessing quality of a simple linear regression

Regenerate the data from last week:


```{r}
set.seed(1984)
n = 20
x = runif(n)
beta_0 = 3
beta_1 = -2
h_star = beta_0 + beta_1 * x
epsilons = rnorm(n, mean = 0, sd = 0.33)
y = h_star + epsilons 

pacman::p_load(ggplot2)
simple_df = data.frame(x = x, y = y)
simple_viz_obj = ggplot(simple_df, aes(x, y)) + 
  geom_point(size = 2)
simple_viz_obj
```

Note: our $\mathcal{A}$ was ordinary least squares. What follows below is a method of assessing model fit quality not only for the least squares line, or any linear fit, but any regression fit.

```{r}
simple_df$yhat = b_0 + b_1 * simple_df$x
simple_df$e = y - simple_df$yhat

sse = sum(simple_df$e^2)
mse = sse / (n - 2)
rmse = sqrt(mse)
sse
mse
rmse

s_sq_y = var(y)
s_sq_e = var(simple_df$e)
rsq = (s_sq_y - s_sq_e) / s_sq_y
rsq

#calculated in a different, but equivalent way
sse_0 = (n - 1) * s_sq_y
(sse_0 - sse) / sse_0
```

Let's take a look at $R^2$ visually. We compute null residuals (the $e_0$'s) and model residuals (the $e$'s) and plot a them.

```{r}
simple_df$e_0 = y - mean(y)

ggplot(simple_df) + 
  geom_histogram(aes(x = e), fill = "darkgreen", alpha = 0.3) + 
  geom_histogram(aes(x = e_0, fill = "red", alpha = 0.3)) +
  theme(legend.position = "none")

ggplot(simple_df) + 
  stat_density(aes(x = e), fill = "darkgreen", alpha = 0.3) + 
  stat_density(aes(x = e_0, fill = "red", alpha = 0.3)) +
  theme(legend.position = "none")
```

Note residuals always have sample average = 0 (modulo numeric error):

```{r}
mean(simple_df$e_0)
mean(simple_df$e)
```

We will prove this fact later in class.

Let's take a look at predictions by truth:


```{r}
ggplot(simple_df, aes(x = yhat, y = y)) + 
  geom_point() + 
  xlim(0, max(simple_df$yhat, y)) + 
  ylim(0, max(simple_df$yhat, y)) + 
  xlab("yhat") + 
  coord_fixed() + 
  geom_abline(intercept = 0, slope = 1, color = "orange")
```

Linear regression is pretty popular so there's obviously support for this in R. Before we talk about this, we need to discuss another object type in R. It is called the "formula" object. Here's an example:

```{r}
simple_model_formula = as.formula("y ~ x")
simple_model_formula
```

You can use a convenience:

```{r}
simple_model_formula = y ~ x
simple_model_formula
```

How did this work? R interprets this as a formula because it sees the tilde inside (you have to dig pretty deep into the R language to understand the tilde operator but feel free to ignore it now).

By default the formula object when executed prints out the string you supplied. But obviously it is not just a string. This object contains instructions to model `y` with `x`. This may seem opaque now but you'll see how this works soon.

Getting back to support for the default linear model. The popular function for that implements least squares linear regression is loaded into R automatically in the package `stats`. Here is a list of packages that are loaded into R by default (save the `ggplot2` and RStudio's addition):

```{r}
search()
```

The function `lm` runs least squares. Let's see it:

```{r}
simple_linear_model = lm(simple_model_formula)
simple_linear_model
```

You can skip a step by putting the formula in-line (this is how it's usually done):


```{r}
simple_linear_model = lm(y ~ x)
simple_linear_model
class(simple_linear_model)
```

What did this do? By specifying the formula that $y$ should be modeled with $x$ sing the simple linear model $y = w_0 + w_1 x$
By default it prints out $b_0$ and $b_1$. You can store the vector via:

```{r}
b = coef(simple_linear_model)
b
names(b)
class(b) #i.e. a vector of numbers dimension 2 where each entry is named
```

You can query the linear model about its fit as well:

```{r}
names(summary(simple_linear_model))
summary(simple_linear_model)$r.squared #the R^2
summary(simple_linear_model)$sigma #the RMSE
```

Cleanup...

```{r}
rm(list = ls())
```


## Simple Linear Regression with an example data set

Load up the famous Boston Housing data

```{r}
?MASS::Boston
Xy = MASS::Boston
head(Xy)
```

We would like to see how each feature relates to the response, `medv`. This is a quick and dirty way to do it:

```{r}
for (feature in setdiff(colnames(Xy), "medv")){
  plot(ggplot(Xy, aes(x = Xy[, feature], y = medv)) + geom_point() + xlab(feature))
}
```

Let's try to explain `medv` using the feature `rm` in a simple linear regression (least squares) model.

```{r}
x = Xy$rm
y = Xy$medv

r = cor(x, y)
s_x = sd(x)
s_y = sd(y)
ybar = mean(y)
xbar = mean(x)

b_1 = r * s_y / s_x
b_0 = ybar - b_1 * xbar
b_0
b_1
```

and we can plot this line atop the data:

```{r}
simple_viz_obj = ggplot(Xy, aes(x = rm, y = medv)) + geom_point()
simple_ls_regression_line = geom_abline(intercept = b_0, slope = b_1, color = "red")
simple_viz_obj + simple_ls_regression_line
```

And how well did we do?

```{r}
yhat = b_0 + b_1 * x #this is the g(x^*) function!
e = y - yhat
sse = sum(e^2)
mse = sse / length(y)
rmse = sqrt(mse)
sse
mse
rmse

s_sq_y = var(y)
s_sq_e = var(e)
rsq = (s_sq_y - s_sq_e) / s_sq_y
rsq
```

SSE is not a super useful number alone. MSE is not super useful alone. RMSE is... what does it mean? What does $R^2$ mean?


```{r}
Xy$null_residuals = y - mean(y)
Xy$residuals = e

ggplot(Xy) + 
  stat_density(aes(x = residuals), fill = "darkgreen", alpha = 0.6, adjust = 0.5) + 
  stat_density(aes(x = null_residuals, fill = "red", alpha = 0.6, adjust = 0.5)) +
  theme(legend.position = "none")
```

This is not a great model. Why? Three sources of error... what do you think are the biggest sources of error?

```{r}
rm(list = ls())
```

Let's do this again using R's `lm` function. Here we can leverage the data frame object by using the names of the variables in our model formula:

```{r}
mod = lm(medv ~ rm, data = MASS::Boston)
```

I read this as "build a linear model where we explain median household value using the average number of rooms in the Boston housing dataset". One line! We can of course ping the model for everything we've been talking about:

```{r}
coef(mod)
summary(mod)$r.squared
summary(mod)$sigma
```

And ggplot has amazing integration. Here it runs the model internally and gives smoothed confidence bands (we did not discuss this):

```{r}
ggplot(MASS::Boston, aes(rm, medv)) + 
  geom_point() + 
  geom_smooth(method = 'lm')
```

```{r}
rm(list = ls())
```


Let's take a look at another dataset.

```{r}
?MASS::Cars93
cars = MASS::Cars93
cars
```

Usually, we are trying to build a model for `Price`. Let's see how `Horsepower` is related to price:

```{r}
pacman::p_load(ggplot2)
ggplot(cars, aes(Horsepower, Price)) + 
  geom_point() + 
  geom_smooth(method = 'lm')
```


```{r}
simple_linear_model = lm(Price ~ Horsepower, data = cars)
coef(simple_linear_model)
summary(simple_linear_model)$r.squared
summary(simple_linear_model)$sigma
```

62\% is pretty good $R^2$! But the RMSE is about \$6,000. Using the empirical rule heuristic, that means you can only predict within $\pm \$12,000$ around 95\% of the time. Not so good!

# Predictions with linear models in R

After the model is fit, you may want to predict with it using the $g$ function. Of course R can do this:

```{r}
predict(simple_linear_model, data.frame(Horsepower = 200)) 
#i.e. yhat = g(400)
predict(simple_linear_model, data.frame(Horsepower = c(200, 300, 500))) 
#i.e. the yhat vector = [g(200), g(300), g(500)]
```