\hrule

## Instructions

* This file serves as a template for your assignment. \textcolor{red}{\textbf{Remember to set eval = TRUE below when compiling the report.}}
```{r setup}
knitr::opts_chunk$set(echo = TRUE, eval = FALSE,
                      fig.height = 4, fig.width = 7, fig.align = "center",
                      message=FALSE)
```
* Use comments to annotate and explain your code whenever possible.
* This assignment covers **Training models** (bias-variance, resampling) and **Regression I** (subset selection, ridge, lasso). It uses one data set: `Boston Housing`, saved as `assign1_regData.csv`
\newline
\hrule

## Theory (10 points)

### T1. Mathematical basis of variance-bias tradeoff (10 points)

Assume $Y = f(X) + \epsilon$, where $E(\epsilon) = 0, \text{Var}(\epsilon) = \sigma^2$ and $\epsilon$ and $f(X)$ are uncorrelated. Prove the key formula for bias-variance tradeoff:
$$ \text{MSE} = (\text{bias})^2 + \text{variance} + \text{irreducible error} $$
where $\text{MSE} = E\left[(y - \hat{f}(x))^2\right]$, $\text{bias} = E[\hat{f}(x)] - f(x)$, $\text{variance} = \text{Var}[\hat{f}(x)]$ and $\text{irreducible error} = \sigma^2$.

<!-- Hint: expand $E\left[(y - \hat{f}(x))^2\right] = E\left[(y - f(x) + f(x) - \hat{f}(x))^2\right]$ and expand again $E\left[(f(x) - \hat{f}(x))^2\right] = E\left[(f(x) - E(\hat{f}(x)) + E(\hat{f}(x)) - \hat{f}(x))^2\right]$. -->

## Application (80 points)
Dataset: `Boston Housing`, saved as `assign1_regData.csv`.

[Harrison and Rubinfeld (1978)](https://doi.org/10.1016/0095-0696(78)90006-2)\footnote{David Harrison, Jr and Daniel L Rubinfeld. Hedonic housing prices and the demand for clean air. Journal of Environmental Economics and Management, 5(1):81–102, 1978.} studied families' willingness to pay for clean air using
Boston area housing market data in 1978. The cleaned data consists of 394 rows with **complete** data on census tracts in Boston, with variables collected from a mixture of surveys, administrative records, and other research papers. The variables are described below.

- `CRIM`:	per-capita crime rate by town (units and crimes included not specified)
- `ZN`:	percentage of residential land zoned for lots over 25,000 sq.ft.
- `INDUS`:	proportion of non-retail business acres per town (nonretail businesses include industrial and commercial businesses)
- `CHAS`:	Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
- `NOX`:	nitric oxides concentration (parts per 10 million)
- `RM`:	average number of rooms per dwelling
- `AGE`:	percentage of owner-occupied units built prior to 1940
- `DIS`:	weighted log distances to five Boston employment centers
- `RAD`:	index of accessibility to radial highways, where a higher number indicates higher accessibility
- `TAX`:	full-value property-tax rate per $10,000
- `PTRATIO`:	pupil-teacher ratio by school district
- `B`:	$1000(B_k - 0.63)^2$ where $B_k$ is the proportion of Black people by town
- `LSTAT`:	% lower status of the population, defined as half of the sum of proportion of adults without some high school education and proportion of male workers classified as laborers
- `MEDV`: Median value of owner-occupied homes, in $1000's

We will load the dataset and set the outcome variable to be `MEDV`.
```{r, echo = TRUE, eval = TRUE, message=FALSE}
# Load required libraries and dataset
pacman::p_load(tidyverse, caret, boot, glmnet, pls, leaps)
set.seed(42)
data <- read_csv("assign1_regData.csv") |> drop_na()
data$CHAS <- as.factor(data$CHAS)
```

### A1. Data Exploration & Baseline Model (5 + 5 + 8 = 18 points)

1. Use the `summary` function to provide descriptive statistics for the dataset. Randomly split your data into a training set (80%) and a test set (20%). Make sure that the distribution of the outcome is similar in both sets (please read the documentation for `caret::createDataPartition` to know how). Use a fixed seed (use `seed = 42`) for reproducibility.
```{r}
# your code here
```

2. Perform a brief exploratory data analysis: comment on how the outcome varies with the `RM` predictor - is the behavior expected?
```{r}
# your code here
```

3. Fit a multiple linear regression model using all predictors on the training data. Report the test set Root Mean Squared Error (RMSE) and interpret the coefficient on `RM`.
```{r}
# your code here
```


### A2. Resampling and Model Assessment (6 + 10 + 4 = 20 points)

1. Using the training data, perform 10-fold cross-validation on the full multiple linear regression model from A1. What is the average cross-validated RMSE?
```{r}
# your code here
```

2. Use the bootstrap method (with $B=1000$ resamples, set `seed = 42`) on the training data to estimate the $95\%$ confidence interval (CI) for the `RM` coefficient. Compare the bootstrap-based $95\%$ CI with the $95\%$ CI using the standard formula from our linear regression model.
```{r}
# your code here
```

3. Interpretation: *Write your interpretation of the bootstrap confidence interval here. How does the theoretically-derived interval compare to the one you just estimated with the bootstrap method? If they differ, what could be some potential causes?*


### A3. Feature Selection & Regularization (12 + 12 + 12 + 6 = 42 points)

1. Stepwise Selection: On the training data, perform forward stepwise selection to find the best model according to BIC. List the variables selected by the procedure, and report the test set RMSE using this smaller model.
```{r}
# your code here
```


2. Ridge Regression: Fit a ridge regression model, using 10-fold CV to determine the optimal lambda for the training set. Report this lambda and the test set RMSE.
```{r}
# your code here
```

3. Lasso Regression: Fit a lasso regression model. Use 10-fold CV to determine the optimal lambda. Report optimal lambda, list variables shrunk to zero, and calculate test set RMSE. Finally, plot coefficient paths for varying lambda.
```{r}
# your code here
```

4. Comparison: Make a small table comparing the test set RMSE of the full linear model (A1.3), the stepwise-selected model, ridge, and lasso. Which model would you recommend, and why?
```{r}
# your code here
```

## Optional / Bonus (unscored, will still provide feedback if you answer)

1. LOO-CV: Prove formula (5.2) from ISLR, Section 5.1.2, page 202.

2. Basis of $0.632$ bootstrap: Problem 5.4.2 from ISLR, page 219.

3. Subgroup fairness (Boston housing): The `B` variable in this data set is defined as $1000(B_k - 0.63)^2$ where $B_k$ is the proportion of Black people by town — due to the squared term this treats towns with very low and very high proportions of Black residents similarly, which is a problematic construction. Split the data into two groups by the median of `B`, fit separate linear models on each group's `RM` and `LSTAT`, and compare in-group vs. cross-group test RMSE. What does this suggest about the dangers of training and deploying models on demographically specific or unrepresentative data?
