title: 'Linear Regression Simulation Project'

author: "Nov-27-2020, Oran Chan"

date: ''

output:

github_document:

pandoc_args: --webtex

# html_document:

# keep_md: true

# theme: readable

# toc: yes

pdf_document: default

urlcolor: cyan

```{r setup, echo = FALSE, message = FALSE, warning = FALSE}

In this simulation study we will investigate the significance of regression test. We will simulate from two different models:

1. The **"significant"** model

Y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + \epsilon_i

where $\epsilon_i \sim N(0, \sigma^2)$ and

- $\beta_0 = 3$,

- $\beta_1 = 1$,

- $\beta_2 = 1$,

- $\beta_3 = 1$.

2. The **"non-significant"** model

Y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + \epsilon_i

where $\epsilon_i \sim N(0, \sigma^2)$ and

- $\beta_0 = 3$,

- $\beta_1 = 0$,

- $\beta_2 = 0$,

- $\beta_3 = 0$.

For both, we will consider a sample size of $25$ and three possible levels of noise. That is, three values of $\sigma$.

- $n = 25$

- $\sigma \in (1, 5, 10)$

Simulation will be used to obtain an empirical distribution for each of the following values, for each of the three values of $\sigma$, for both models.

- The **$F$ statistic** for the significance of regression test.

- The **p-value** for the significance of regression test

- **$R^2$**

Data found in [`study_1.csv`](study_1.csv) will be the values of the predictors. These will be kept constant for the entirety of this study. The `y` values in this data are a blank placeholder.

For each model and $\sigma$ combination, there will be $2000$ simulations. For each simulation, we will fit a regression model of the same form used to perform the simulation.

The `y` vector $2 (models)×3 (sigmas)×2000 (sims)=12000$ will be simulated times.

- The null and alternative hypotheses

- $H_0: \beta_1 = \beta_2 = \beta_3 = 0$

- $H_1:$ At Least one of $\beta_j \neq 0,j=1,2,3$

- Simulation results will be plotted by below functions

```{r fig.height=5, fig.width=10}

```{r fig.height=5, fig.width=10}

```{r fig.height=5, fig.width=10}

- $H_0: \beta_1 = \beta_2 = \beta_3 = 0$

- $H_1:$ At Least one of $\beta_j \neq 0,j=1,2,3$

- Since know the true distribution of F-statistics with two of the degress of freedom, it has been plotted with the empirical value in the F-statistics grid of plots in coming discussion.

- F-statistics has a distribution under null hypothesis. It means that assuming null hypothesis is true, the test statistics follows an F-distribution.

Since F-statistics has the above assumption, the empirical value of **"significant"** model does not follow the true distribution.

- We can see in the plot on the left below, most of the F-value fall in the range from 20 to 50, which are considered as some very large F values. Such large F value corresponds to a very small p-value, that is to reject the null hypothesis.

- When we look at the plot of $R^2$ for $\sigma$=1, majority of value falls above 0.7. It means the althernative hypothesis model is explaining well by the predictors on the relationship with the response. It is also a indicator for the decision to reject $H_0$.

- The average F-statistics, $R^2$ and p-value are listed for each $\sigma$.

```{r fig.height=5, fig.width=10}

- In the above plots of F-value for the other two $\sigma$(s), the empirical values lean to the left of the graph, which represents majority of the F-value are small. It also means the corresponding p-value becomes larger.

- In the plots of $R^2$, majority of values fall in the range of 0.1 to 0.3 for $\sigma$=5, 0.05 to 0.15 for $\sigma$=10. The model with these $\sigma$ do not performe well in explaining the relatinoship.

- With sufficiently larger p-value, the chance of FTR $H_0$ will be much higher than that with $\sigma$=1.

- They partially align with the characteristic of F-distribution. Nonetheless, even the empirical values becomes closer to the true value, it is still far from it. We will see how it goes with **"Non-significant "** model in the next section.

P-value distribution of 3 $\sigma$ varies.

- In the first plots with $\sigma$ = 1, as discussed above, it comes with a critically small p-value so the decision is to reject the null hypothesis. The peak close to zero as shown in the plot is where the alternative hypothesis located.

```{r fig.height=3, fig.width=10}

- In the second plot with $\sigma$ = 5. More p-value is falling into larger values, which contribute to null hypothesis. The frequency of the value close to zero is lower comparing to the 1st plot.

- We should be reminded that null hypothesis contains low p-value as well.

- In the thrid plot with $\sigma$ = 10. Even p-value is falling into larger values, we can see a close to uniform distribution when $p-value > 0.2$. For p-value < 0.1, frequency from 0 to 40 might contain null hypothesis as well.

- The interpretation of p-value aligns with our previous explaination that the chance of FTR $H_0$ will be much higher.

```{r fig.height=5, fig.width=10}

- As shown from the plots above, the F-value distribution of all $\sigma$ are critically close to the true F-distribution.

- The average F values are rougly between 1.1 and 1.2, which is considered very small comparning to the significant model.

- The average $R^2$ of all $\sigma$(s) are very small. It implies the **"non-significant"** model is not performing well in explaining the relationship between the predictors and reponse.

- The average P-value of all $\sigma$(s) are considered very large as they are close to 0.5.It alignes with the small F-values. Therefore, the decision is very likely to be FTR the null hypothesis.

- Hence, the very similar behaviour between the true F-distrbution and the **"non-significant"** model could be explained by the above findings.

- We can also interpret the result by looking at setting of $\beta_{\texttt{1-3}}$. As they are all zero, so we expect the $H_0$ is very likely to be true.

P-value distribution of 3 $\sigma$ are very similar.

- The P-value follow a uniform distribution in all three plots.

```{r fig.height=3, fig.width=10}

- It means that 5% of the p-value will be above the 9th percentile, another 5% will be between 90th and 95th percentile, and so forth. 5% of the value will be from 0 to 0.05, and another will be from 0.05 to 0.1.

- The uniformly distributed p-value is one of the definition of a p-value under null hypothesis.

- The above observation aligns with our discussed devision to be FTR $H_0$.

We will look into the relatinoship between $\sigma$, and $F$ statistic, the p-value, $R^2$.

- When we perform a significant test on regression, the relationship varies differently according to the settings of the model. Specifically, the relationship is affected by the value of $\beta$.

In **"Significant"** model, we have positive signal strength for $\beta_{\texttt{1-3}}$.

- When $\sigma$ becomes larger, the F-statistic distribution becomes closer to the true F-distribution. The value of $R^2$ drops significantly and p-value increases significantly.

- We can interpret that, when $\sigma$ becomes larger in **"Signigficant"** model, it is more likely to be FTR $H_0$.

In **"Non-Significant"** model, we have zero signal strength for $\beta_{\texttt{1-3}}$ and $\beta_0$ = 3.

- When $\sigma$ varies, the effect on $F$ statistic, the p-value and $R^2$ are considerably small.

- Hence, given $\beta_0$ = 3, under the settings of **"Non-Significant"**, we can interpret that, $\sigma$ has very slight influcence on $F$ statistic, the p-value and $R^2$. The decision of FTR $H_0$ would remain unchanged.

In this simulation study we will investigate how well this procedure works. Since splitting the data is random, we don’t expect it to work correctly each time. We could get unlucky. But averaged over many attempts, we should expect it to select the appropriate model.

We will simulate from the model

Y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + \beta_4 x_{i4} + \beta_5 x_{i5} + \beta_6 x_{i6} + \epsilon_i

where $\epsilon_i \sim N(0, \sigma^2)$ and

- $\beta_0 = 0$,

- $\beta_1 = 3$,

- $\beta_2 = -4$,

- $\beta_3 = 1.6$,

- $\beta_4 = -1.1$,

- $\beta_5 = 0.7$,

- $\beta_6 = 0.5$.

Data in [`study_2.csv`](study_2.csv) will be the values of the predictors. These will be kept constant for the entirety of this study. The `y` values in this data are a blank placeholder.

We will consider a sample size of $500$ and three possible levels of noise. That is, three values of $\sigma$.

- $n = 500$

- $\sigma \in (1, 2, 4)$

Each time the data is simulated, they will be randomly split into train and test sets of equal sizes (250 observations for training, 250 observations for testing).

For each, we will fit **nine** models, with forms:

- `y ~ x1`

- `y ~ x1 + x2`

- `y ~ x1 + x2 + x3`

- `y ~ x1 + x2 + x3 + x4`

- `y ~ x1 + x2 + x3 + x4 + x5`

- `y ~ x1 + x2 + x3 + x4 + x5 + x6`, the correct form of the model as noted above

- `y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7`

- `y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8`

- `y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9`

For each model, we will calculate Train and Test RMSE.

\text{RMSE}(\text{model, data}) = \sqrt{\frac{1}{n} \sum_{i = 1}^{n}(y_i - \hat{y}_i)^2}

- In each simluation, one out of nine models will be selected for each sigma.

- The one with lowest RMSE value in test set will be selected.

This above process will be repeated with $1000$ simulations for each of the $3$ values of $\sigma$. For each value of $\sigma$, there will be a plot that shows how average Train RMSE and average Test RMSE changes as a function of model size. We can also tell number of times the model of each size was chosen for each value of $\sigma$.

We will simulate the $y$ vector $3×1000=3000$ times and fit $9×3×1000=27000$ models.

- Simulation results will be plotted by below functions

```{r fig.height=4, fig.width=10}

We will discuss whether the RMSE method **always** select the correct model.

- RMSE can be interpreted as the smaller the better.

- From the above 3 plots, we can see a considerably high RMSE value when model size is 1 and 2. It is a case of underfitting for both train and test set.

- RMSE value starts to diminish slowely after sitting at the model with size 3.

- Meanwhile, the test and train set RMSE start to diverge.

- Train set RMSE remains in a state of good performance and it is further improved at model with size 6. We don't see any issue of underfitting as train set RMSE converges to the minimum at model with size 9.

- RMSE value starts to diminish slowely after model with size 3, and even slower than Train set RMSE.

- Test RMSE set converges to minimum point at model with size 6.

- However, it starts to rebound from model size 7 to 9. It is a case of overfitting.

```{r fig.height=5, fig.width=10}

- Hence, we conclude that model with size 6 is the correct model for explaining the dataset. Although train set RMSE is the lowest for model 9, test set RMSE is higher, which indicates the model does not work well with new data and suffers from overfitting.

- In the above plots, the red line indicates the correct model size.

- The model with size 6, is the most selected model for all 3 $\sigma$(s).

- But the portion of the selection is not ideal. The highest percentage rougly **54%** of $\sigma$ = 1. The lowest is only **34%** of $\sigma$ = 4.

```{r fig.width=5, fig.align='center'}

- The RMSE method **does not** always select the correct model. It is not a reliable method for model selection in current study.

We will discuss whether the level of noice would affect our selection result.

```{r fig.height=5, fig.width=10}

- RMSE

- When $\sigma$ becomes larger, the average RMSE value becomes larger, but the variability of RMSE value becomes lower. There is a positive relationship between $\sigma$ and average RMSE value.

- When $\sigma$ becomes larger, the performace of test set will start to decrease earlier. Test and train set will diverge more than that of smaller $\sigma$.

- Model selection

- When $\sigma$ becomes larger, the selection becomes broader.

- $\sigma$ = 1, 4 out of 9 models were selected

- $\sigma$ = 2, 7 out of 9 models were selected

- $\sigma$ = 4, 8 out of 9 models were selected

We can interpret as, with higher variance, there is higher probability other models will perform better, but generally over 1000 simulations, they perform poorly.

- When $\sigma$ becomes larger, the lower chance the correct model is selected.

We can interpret that some other models would perform better due to high variance. Therefore, we prefer a lower $\sigma$ value during model selection process.

In this simulation study we will investigate the **power** of the significance of regression test for simple linear regression.

H_0: \beta_{1} = 0 \ \text{vs} \ H_1: \beta_{1} \neq 0

Recall, we had defined the *significance* level, $\alpha$, to be the probability of a Type I error.

\alpha = P[\text{Reject } H_0 \mid H_0 \text{ True}] = P[\text{Type I Error}]

Similarly, the probability of a Type II error is often denoted using $\beta$; however, this should not be confused with a regression parameter.

\beta = P[\text{Fail to Reject } H_0 \mid H_1 \text{ True}] = P[\text{Type II Error}]

*Power* is the probability of rejecting the null hypothesis when the null is not true, that is, the alternative is true and $\beta_{1}$ is non-zero.

\text{Power} = 1 - \beta = P[\text{Reject } H_0 \mid H_1 \text{ True}]

Essentially, power is the probability that a signal of a particular strength will be detected. Many things affect the power of a test. In this case, some of those are:

- Sample Size, $n$

- Signal Strength, $\beta_1$

- Noise Level, $\sigma$

- Significance Level, $\alpha$

We'll investigate the first three.

For simplicity, we will let $\beta_0 = 0$, thus $\beta_1$ is essentially controlling the amount of "signal." We will then consider different signals, noises, and sample sizes:

- $\beta_1 \in (-2, -1.9, -1.8, \ldots, -0.1, 0, 0.1, 0.2, 0.3, \ldots 1.9, 2)$

- $\sigma \in (1, 2, 4)$

- $n \in (10, 20, 30)$

We will hold the significance level constant at $\alpha = 0.05$.

We will simulate from the model

Y_i = \beta_0 + \beta_1 x_i + \epsilon_i

where $\epsilon_i \sim N(0, \sigma^2)$.

For each possible $\beta_1$ and $\sigma$ combination, simulate from the true model at least $1000$ times. Each time, perform the significance of the regression test.

The following sample code is the method to generate the predictor values, `x`: values for different sample sizes.

```{r eval=FALSE}

In the functions below, before plotting the result, we will calculate the estimate of Power by this formula,

\hat{\text{Power}} = \hat{P}[\text{Reject } H_0 \mid H_1 \text{ True}] = \frac{\text{# Tests Rejected}}{\text{# Simulations}}