# Instructions

+ Upload a PDF file, named with your UC Davis email ID and homework number (e.g., `xjw18_hw7.pdf`), to Gradescope (accessible through Canvas). You will give the commands to answer each question in its own code block, which will also produce output that will be automatically embedded in the output file. When asked, answer must be supported by written statements as well as any code used. 

+ All code used to produce your results must be shown in your PDF file (e.g., do not use `echo = FALSE` or `include = FALSE` as options anywhere). `Rmd` files do not need to be submitted, but may be requested by the TA and must be available when the assignment is submitted.

+ Students may choose to collaborate with each other on the homework, but must clearly indicate with whom they collaborated. Every student must upload their own submission.

+ Start to work on it as early as possible. Finishing this homework can help prepare midterm 1.

+ When you want to show your result as a vector that is too long, slice the first 10 objects. When you want to show your result as a data frame, use `head()` on it. Failure to do so may lead to point deduction.

+ Directly knit the Rmd file will give you an html file. Open that file in your browser and then you can print it into a PDF file.

```{r setup, include=FALSE}
# Don't change anything here
knitr::opts_chunk$set(
  comment = "", prompt = F, message=F, warning = F
)
```


```{r}
# Load necessary libraries
library(tidyverse)
library(ggplot2)
library(MASS)
```

# Problem 1: the Boston data set

There are 14 variables (columns) and 506 rows in the Boston house values data frame in the `MASS` package. To look at those variables you may run `names(Boston)` however we are only going to look at 2 of those variables:

* `lstat` lower status of the population (percent).

* `medv` median value of owner-occupied homes in $1000

We want to build a model use the lower status of the population (percent) to predict the median value of owner-occupied homes in $1000.

1. Find the estimated regression line. Use `summary` to show the result.


2. Interpret the slope of the estimated regression line in terms of the problem.


3. Interpret the intercept of the estimated regression line in terms of the problem (if appropriate). 


4. Predict the median value of owner-occupied homes in $1000 when lower status of the population is 25%.


5. Interpret the value of R squared (`Multiple R-squared` in `summary`) in terms of the problem.


6. What is the $s_e$ ($\hat \sigma$) from the regression result?



7. Find the 99% confidence interval for the slope and intercept using `confint`.


# Problem 2: Sampling distribution of regression

In this problem, we will use simulation to study the distribution of the regression estimates.

Suppose our true model is:

$$Y = \beta_0 +\beta_1 X +\varepsilon$$

Where $X \sim Uniform(-3,3)$, $\varepsilon \sim N(0,\sigma^2)$.

1. When $\beta_0 = 1$,$\beta_1 = 2$, $\sigma = 1$, write R code to simulate 100 pairs of $X$ and $Y$, and store them into a data frame. Use scatter plot to show the data.


2. Wrap the code to generate the dataframe as a function, with name `simulate_1`. The function should take inputs as: beta_0, beta_1, sigma and n. Use the values from part 1 as the default argument value. [About default argument value](https://www.javatpoint.com/r-function-default-arguments)


3. Now use the `simulate_1()` simulate 1 dataset (blank in the `()` to use default argument value). Then first use the `lm` function fit the model, and use the [slide formula](https://xjw1001001.github.io/lecture/Lecture%2023/lecture23.html#22) to fit the model, check the result are the same. 


4. Here I provide the code to fit a linear model and obtain its estimate of slope, intercept and the estimation of the [sample SD of residual](https://xjw1001001.github.io/lecture/Lecture%2024/lecture24.html#17) (sigma_hat in code):

```
#suppose the dataset name is df:
set.seed(123)
df = simulate_1()

fit = lm(Y~X, df)
slope = fit$coefficients[2]
intercept = fit$coefficients[1]
sigma_hat = sigma(fit)
```

Now calculate the `sigma_hat` by the definition in the lecture note [sample SD of residual](https://xjw1001001.github.io/lecture/Lecture%2024/lecture24.html#17) and check it's the same with the result `sigma(fit)`. In your calculation, the dataset `df` and the fitted value of slope and intercept just use the ones provided in the previous code chunk.

Hint: $e_i$ are the difference of the actual $y_i$ with the predicted value $\hat y_i$. The predicted value can be obtained by plug in $x_i$ into the fitted regression line. 


5. Each time we simulate 100 pairs of the points, we can calculate the estimated slope and intercept, plus the estimated error standard deviation sigma_hat for these 100 pairs of the points. Now use for loop to simulate 500 sample of the 100 pairs of the points, calculate the estimated slope and intercept, sigma_hat for each of the sample. Store the result into a dataframe named `regression_simulation`.

```{r}
regression_simulation = data.frame(
  index = 1:500, slope = numeric(500),
  intercept = numeric(500), sigma_hat = numeric(500)
)
#you need to fill in the corresponding result for slope,
# intercept, sigma_hat in the for loop.
```


6. Now we have the 500 samples of regression line estimation: slope, intercept and sigma_hat. Now summary their mean and standard deviation using `summarise` function.


7. Call `summary(fit)`, compare the mean from part 6 with the actual value of $\beta_0,\beta_1,\sigma$. What do you find? Compare the sd from part 6 with the `Std. Error` of intercept and X, what do you find?

8. Understand the p value in summary. In the lecture we know in the summary we test the null hypothesis about the slope coefficient is 0 against the alternative hypothesis the slope coefficient is not 0. The test statistic is:

$$T = \frac{\hat \beta_1 - \beta_1}{SE(\hat \beta_1)} = \frac{\hat \beta_1 - 0 }{SE(\hat \beta_1)} \sim t_{n-2} \text{   (Under null hypothesis)}$$

Now we will use our simulation of the distribution of $\beta_1$ to check the p value of $H_0: \beta_1 = 1.9$. 

Now draw the histogram of the `regression_simulation$slope`, and add a vertical line at 1.9. 

Also calculate the $T$ statistics when $H_0: \beta_1 = 1.9$ by plug in $\beta_1 = 1.9$ (Use Problem 7 summary result), and calculate the corresponding p value by using `pt(t_statistic, df = n-2, lower.tail = F)`. What should the `df` be beased on the summary? Is the p value explain the histogram and the vertical line?

# Problem 3: violations of the regression assumption: linearity.

Now in this problem, we will look at how the violations of the regression assumption: linearity will affect the linear regression model.

Suppose our true model is:

$$Y = \beta_0 +\beta_1 (X-1.5)^2 +\varepsilon$$

Where $X \sim Uniform(-3,3)$, $\varepsilon \sim N(0,\sigma^2)$, $\beta_0 = 1$,$\beta_1 = 2$, $\sigma = 4$

1. write R code to simulate 100 pairs of $X$ and $Y$, and store them into a data frame. Use scatter plot to show the data. Add layer 1: `geom_smooth(color = "blue")` to add the true fitted line, and add layer 2: `geom_smooth(method = lm, color = "red")` to add the linear regression fitted line. Remember to `set.seed(123)`


2. Is the linear regression still a good model for this dataset? Why?

3. Now fit the linear model on the data you generated. Save the fitted model as `fit`. Then use `plot(fit, 1)` to generate the [residual VS fitted plot](https://online.stat.psu.edu/stat462/node/117/). Interpret the residual vs fitted plot. 


4. Now fit the polynomial model on the data you generated use: `fit2 = lm(Y ~ X + I(X^2), df)`. Then use `plot(fit2, 1)` to generate the residual VS fitted plot. Check `summary(fit2)` in your console (don't print it in rmd) and then write down the model in mathematical form $Y = \beta_0 + \beta_1 X + \beta_2 X^2$. Also check the residual vs fitted plot. 


5. Draw the $Y = \beta_0 + \beta_1 X + \beta_2 X^2$ and the dataset using ggplot.


# Problem 4: violations of the regression assumption: equal residual variance.

Now in this problem, we will look at how the violations of the regression assumption: equal residual variance will affect the linear regression model.

Suppose our true model is:

$$Y = \beta_0 +\beta_1 X +\varepsilon$$

Where $X \sim Uniform(-3,3)$, $\varepsilon \sim N(0,(\sigma*|X|)^2)$, $\beta_0 = 1$,$\beta_1 = 2$, $\sigma = 4$

1. write R code to simulate 100 pairs of $X$ and $Y$, and store them into a data frame. Use scatter plot to show the data. Add layer 1: `geom_smooth(method = lm, color = "red")` to add the true fitted line. Remember to `set.seed(123)`.

Hint: `rnorm(n, mean_vector, sd_vector)` will generate n normal RV, the ith mean and ith sd is the corresponding in the `mean_vector` and `sd_vector`.


2. Is the linear regression still a good model for this dataset? Why?

3. Now fit the linear model on the data you generated. Save the fitted model as `fit`. Then use `plot(fit, 1)` to generate the [residual VS fitted plot](https://online.stat.psu.edu/stat462/node/117/). Interpret the residual vs fitted plot. 



## Include the person you coop with:

Names:

## Appendix

```{r eval = TRUE}
sessionInfo()
```











