chapter2.Rmd

---
title: "Chapter 2"
author: "Nate"
date: "5/9/2020"
output: html_document
---

```{r setup, message = F}
knitr::opts_chunk$set(comment = NA, message=F, warning=F)
library(faraway)
library(skimr)
library(magrittr)
library(tidyverse)
```

# 1. 
```{r teengamb}
skim(teengamb)
m <- lm(gamble ~ ., teengamb)
```

## a.
```{r teengamb_a}
sumary(m)
```

53%, reported by **the coefficient of determination**, aka 'R-squared'.

## b.
```{r teengamb_b}
res <- resid(m)

which.max(res)
```

Case (row-number) 24 gambles more than the model expects.

## c.
```{r teengamb_c}
mean(res)
median(res)
```

The mean of the residuals should always be zero, we get this very-small (practically zero) value because of how computer math works. The median is a product of have only 47 cases, as more observations are collected the median will trend towards zero too.

## d.

```{r teengamb_d}
predict(m) %>% cor(res)
```

Again a practically zero value, indicating no correlation. This is means our linear model assumptions were not violated.

## e.

```{r teengamb_e}
teengamb$income %>% cor(res)
```

## f.
```{r teengamb_f}
?teengamb
# 0 = male
# gamb in pounds per year

sumary(m)
```

An average female would spend 22 pounds less per year that an average male, if we consider status, income and verbal as nuissance variables.

# 2.
```{r}
data("uswages")
?uswages
skim(uswages)

m <- lm(wage ~ educ + exper, uswages)
sumary(m)
```

For every year of education the model estimates a \$51 per week increase in wage. For every year of experience the model esimates only \$10.

```{r log_wage}
update(m, log(wage, 10) ~ .) %>% sumary()
```

Log transformation of the response, means we interpret the coefficients as magitude (multiplicative) changes. So for each year of education the model expects a 4% increase in wage. For each year of experience < 1%.

# 3.

```{r hand_calc}
x <- 1:20
y <- x + rnorm(20)

m <- lm(y ~ I(x^2))
sumary(m)

direct_calc <- function(x, y, degree = 2) {
  x_mat <- model.matrix(~I(x^degree))
  solve(crossprod(x_mat), crossprod(x_mat, y))
}

direct_calc(x, y)

map(set_names(3:7), possibly(~ direct_calc(x, y, degree = .), "Error"))
```

# 4.

```{r prostate}
data("prostate")
skim(prostate)

m <- lm(lpsa ~ lcavol, prostate)

get_resid_error <- function(m) {
  s <- invisible(summary(m))
  data.frame(
    "mod" = formula(m) %>% deparse(),
    "sigma" = s$sigma,
    "r_squared" = s$r.squared
    )
} 

resid_error_dat <- get_resid_error(m)

predictors_to_add <- c("lweight", "svi", "lbph", "age", "lcp", "pgg45", "gleason")

for (p in predictors_to_add) {
  f <- formula(m) %>% deparse() %>% paste("+", p)
  m <- update(m, f)
  
  resid_error_dat <- bind_rows(resid_error_dat, get_resid_error(m))
}

resid_error_dat %>% 
  pivot_longer(-mod, "key") %>% 
  ggplot(aes(mod, value, color = key, group = key)) +
  geom_path() +
  theme(axis.text.x = element_text(angle= 25, hjust = 1))
```

As more predictors are added, the estimate of population variance decreases and the coefficient of determinaiton (ammount of variance explained by the model) increases. Unadjuted R-squared will always increase as predictors are added to the model, the [adjusted R-squared](https://en.wikipedia.org/wiki/Coefficient_of_determination#Adjusted_R2) adds a penalty for n-predictors to compensate for that. 

Note the rapid improvment in fit up to the model `lpsa ~ lcavol + lweight + svi`, perhaps that's a good place to stop adding more predictors.

# 5.

```{r prostate2}
m <- lm(lcavol ~ lpsa, prostate)
m2 <- lm(lpsa ~ lcavol, prostate)

ggplot(prostate, aes(lcavol, lpsa)) +
  geom_point(alpha = .5) +
  geom_line(aes(x = predict(m), color = "lcavol ~ lpsa")) +
  geom_line(aes(y = predict(m2), color = "lpsa ~ lcavol")) +
  geom_point(aes(y = 2.48, x = 1.36), shape = 1, size = 10)
```

Algebra not showm but you can calculate the intesection if you solve for either `lcavol` or `lpsa` in the system of equations represented by the model coefficients.

# 6.

```{r cheddar}
data("cheddar")
skim(cheddar)
```

## a.

```{r cheddar_a}
m <- lm(taste ~ ., cheddar)
sumary(m)
```

## b.

```{r cheddar_b}
cor(m$fitted.values, cheddar$taste)^2
```

R-Squared, the coeffecient of determination or percentage of variance explained

## c. 

```{r cheddar_c}
m2 <- update(m, . ~ -1 + .)
cor(m2$fitted.values, cheddar$taste)^2
```

## d.

```{r cheddar_d}
m_mat <- model.matrix(m)

qr_decomp <- qr(m_mat)

backsolve(
  qr.R(qr_decomp), # upper-right
  t(qr.Q(qr_decomp)) %*% cheddar$taste
)
```

# 7.

```{r wafer}
data("wafer")
skim(wafer)
m <- lm(resist ~ ., wafer)
```

## a.

```{r wafer_a}
x <- model.matrix(m)
x
```

As [dummy variables](https://en.wikipedia.org/wiki/Dummy_variable_(statistics)), 1s and 0s.

## b.

```{r wafer_b}
cor(x)
```

The intercept is constant and has no variance, so `cor()` throws a `NA` back.

```{r}
sd(x[,1])
cor(x[,1], x[,1])
```

## c.

```{r wafer_c}
sumary(m)
```

Resistivity will increase by 26 unknown units. The SI unit of resistivity are ohm-meter.[Wikipedia](https://en.wikipedia.org/wiki/Electrical_resistivity_and_conductivity)

## d. 

```{r wafer_d}
m2 <- update(m, . ~ . - x4)
sumary(m)
sumary(m2)
```

The coefficient estimates are idential but everything else is differnt. There is more estimated variance without `x4` in the model so the `t.values` get closer to zero and the `p.values` increase slightly. None of these changes impact the inferences from the model.

## e.

The corrlations matrix shows that none of the predictors are correlated, so removing `x4` has no impact on the value estimates for the remaining predictors.


# 8. 

```{r truck}
data("truck")
skim(truck)

truck %<>% mutate_if(is.factor, ~ ifelse(. == "-", -1, 1))
```


## a.

```{r truck_a}
m <- lm(height ~ ., truck)
sumary(m)
```

## b.

```{r truck_b}
m2 <- update(m, . ~ . -O)
sumary(m2)
```

The coefficients are identical, because there is no correlation between predictors, `cor(model.matrix(m))` shows us that.

## c. 

```{r truck_c}

truck %<>% 
  rowwise() %>% 
  mutate(A = sum(B, C, D, O))

m3 <- update(m, . ~ . + A)
sumary(m3)
```

`A` is not present becuse it is a linear combination the other predictors. Because the model realizes it can't identify all of the predictors it warns us and drops the last predictor listed in the formula, until the model becomes identifiable. In this case `A` gets cut, but `update(m, . ~ A + .) %>% sumary()` would drop `O`.

## d.

```{r truck_d, error=TRUE}
x <- model.matrix(m3)
solve(t(x) %*% x)
```

Singularities are a real problem.

## e.

```{r truck_e}
qr_decomp <- qr(x)
backsolve(
  qr.R(qr_decomp), # upper-right
  t(qr.Q(qr_decomp)) %*% truck$height
)
```

No, most of the estimates are insanely large, which indicate a problem.

## f. 

```{r truck_f}
qr.coef(qr_decomp, truck$height)
```