# <b>Základy regrese | NMFM 334</b>

<b>Letný semester 2024 | Cvičenie 4 | 19.03.2024</b><br>
HTML Markdown -- <a href = "cvicenie4.Rmd">Rmd source file</a> (UTF encoding)


**************
**************

<b>Outline of the lab session (individual work):</b><br><br>
<ul>
<li>Multiple regression model for the contnuous response vector</li> 
<li>Multiple regression model for a categorical response variable</li> 
</ul>

<br><br>


The dataset `mtcars` (available under the standard R installation) will be used 
for the purposes of this lab session. Some insight about the data can be taken from 
the help session (command `?mtcars`) or from some outputs provided below. 

```{R}
head(mtcars)
str(mtcars)
summary(mtcars)
```

<br><br>

### <b>1. Multiple linear regression model</b> 


For simplicity, we will only consider three explanatory (continuous) variables from the whole `mtcars` dataset, i.e. $\boldsymbol{x} \in \mathbb{R}^3$. Any generalization for general $\boldsymbol{X} \in \mathbb{R}^p$
for some $p \in \mathbb{N}$ is straightforward. 
In particular, the (continuous) dependent variable $Y \in \mathbb{R}$ is the car consumption (variable in the dataset denoted as `mpg` -- miles per gallon) 
and three independent variables will be the car's displacement (`disp`), horse powers (`hp`), and the rear axle ratio information (`drat`). 

<br><br>

Three simple scatterplots showing the dependent variable $Y$ agains one independent variable (based on the random sample $\{(Y_1, \boldsymbol{X}_i^\top)^\top;~i = 1, \dots, 32\}$) are shown below. 
Only a limited information about the whole multivariate structure of the underlying data can be, howver, taken from the series of two-dimensional scatterplots. 

```{r, fig.width = 10}
par(mfrow = c(1,3))
plot(mpg ~ disp, data = mtcars, pch = 22, cex = 0.8, bg = "lightblue", ylab = "Consumption [miles per gallon]", xlab = "Displacement [cubic inches]")
lines(lowess(mtcars$mpg ~ mtcars$disp), col = "red", lwd = 2)
plot(mpg ~ hp, data = mtcars, pch = 22, cex = 0.8, bg = "lightblue", ylab = "", xlab = "Horse Power")
lines(lowess(mtcars$mpg ~ mtcars$hp), col = "red", lwd = 2)
plot(mpg ~ drat, data = mtcars, pch = 22, cex = 0.8, bg = "lightblue", ylab = "", xlab = "Rear axle ratio")
lines(lowess(mtcars$mpg ~ mtcars$drat), col = "red", lwd = 2)
```

The red curves are non-parametric estimates of the conditional mean---these estimates offer a usefull insight 
that the dependence of the consumptions is very likely negative (regressive) wrt. the displacement and horsepower while it seems to positive (progresive) for the rear axle ratio. 

A simple multiple regression model of the form 
$$
Y = a + \beta_1 \textrm{`disp`} + \beta_2\textrm{`hp`} + \beta_3 \textrm{`drat`} + \varepsilon
$$
or mathematically (the model for the empirical data, for $i = 1, \dots, 32$)
$$
Y_i = \boldsymbol{X}_i^\top\boldsymbol{\beta} + \varepsilon_i
$$
where $\boldsymbol{X} = (1, disp, hp, drat)^\top$ and $\boldsymbol{\beta} \in \mathbb{R}^4$ will be considered. 
In the R software, the model can be fitted using the following R code: 

```{r}
fit1 <- lm(mpg ~ disp + hp + drat, data = mtcars)
print(fit1)
```

The information in the output above gives values for four estimated parameters -- elements of the parameter vector $\boldsymbol{\beta} \in \mathbb{R}^4$. The estimated parameters 
are not the same as estimated three independent simple linear regression model. Compare the output above with the following outputs: 

```{r}
fit_disp <- lm(mpg ~ disp, data = mtcars)
fit_hp <- lm(mpg ~ hp, data = mtcars)
fit_drat <- lm(mpg ~ drat, data = mtcars)

print(fit_disp)
print(fit_hp)
print(fit_drat)
```

Indeed, while the sigs of the estimated parameters correspond, the magnitues are substantially different. The ordinary regressions can be directly 
ploted into the corresponding scatterplots: 

```{r, fig.width = 10}
par(mfrow = c(1,3))
plot(mpg ~ disp, data = mtcars, pch = 22, cex = 0.8, bg = "lightblue", ylab = "Consumption [miles per gallon]", xlab = "Displacement [cubic inches]")
lines(lowess(mtcars$mpg ~ mtcars$disp), col = "red", lwd = 2)
abline(fit_disp, col = "blue")
plot(mpg ~ hp, data = mtcars, pch = 22, cex = 0.8, bg = "lightblue", ylab = "", xlab = "Horse Power")
lines(lowess(mtcars$mpg ~ mtcars$hp), col = "red", lwd = 2)
abline(fit_hp, col = "blue")
plot(mpg ~ drat, data = mtcars, pch = 22, cex = 0.8, bg = "lightblue", ylab = "", xlab = "Rear axle ratio")
lines(lowess(mtcars$mpg ~ mtcars$drat), col = "red", lwd = 2)
abline(fit_drat, col = "blue")
```

However, the multiple regression model does not estimate a line but a three dimensional subspace. For instance, the intercept value $\widehat{a} = 29.59$ from the simple regression 
model `fit_disp` stands for the expected consuption (miles driven per one gallon) if the displacement of the engine is equal to zezo cubic inches (i.e., an unrealistic scenario).

Nevertheless, the intercept estimate $\widehat{\beta}_0 = 19.34$ stands for the expected car consumption (again in  miles driven per one gallon) but if all three regressors -- the information about 
the displacement, the horse powers, and, also, the rear to axel ratio -- are all equal to ze. Therefore, to draw simple line for two-dimensional scatterplot, it is needed to pre-select some reasonable 
information (realistic case) for all other covariates in the model. 

For instance, consider the multiple model from above and, in particular, the scatterplot of the consumption against the displacement. Various regression lines can be considered: 

```{r, fig.width = 10}
plot(mpg ~ disp, data = mtcars, pch = 22, cex = 0.8, bg = "lightblue", ylab = "Consumption [miles per gallon]", xlab = "Displacement [cubic inches]")

xGrid <- seq(min(mtcars$disp), max(mtcars$disp, 100))
mCoeff <- fit1$coeff

yMeanMean <- mCoeff[1] + mCoeff[2] * xGrid + mCoeff[3] * mean(mtcars$hp) + mCoeff[4] * mean(mtcars$drat) 
lines(yMeanMean ~ xGrid, col = "red", lwd = 2)

yMinMin <- mCoeff[1] + mCoeff[2] * xGrid + mCoeff[3] * min(mtcars$hp) + mCoeff[4] * min(mtcars$drat) 
yMaxMax<- mCoeff[1] + mCoeff[2] * xGrid + mCoeff[3] * max(mtcars$hp) + mCoeff[4] * max(mtcars$drat) 
lines(yMinMin ~ xGrid, col = "blue", lwd = 2, lty = 2)
lines(yMaxMax ~ xGrid, col = "darkgreen", lwd = 2, lty = 2)

legend("topright", legend = c("Average HP & average rear axle ratio", "Min HP & min rear axle ratio", "Max HP & max rear axle ratio"), col =c("red", "blue", "darkgreen"), lwd = c(3,3,3), lty = c(1,2,2))
```

<br><br>

### <font color = "#880000">Indidividual work</font> 
The whole set of different regression lines can be plotted depending on particular choices for the remaining covariates that are not used 
in the scatterplot (in this case, the information about the horse powers and the rear axle ratio). 
<ul>
<li>For all three considered covariates draw a separate scatterplot with a set of regression lines that will reasonably explain the underlying 
dependence of the car consumption.</li>
<li>In practical applications, it is reasonable to use the exploratory analysis to pre-select some reasonable values for the covariates 
that are used in the model formula, the relationship of the dependent covariate is not directly illustrated in the plot. For instance, 
one can use a series of nice curves for the estimated deciles of, let's say the horse power, while always using the median value of the rear axle ratio.</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>


### <b>2. Basic inference on $\boldsymbol{\beta} \in \mathbb{R}^p$ </b>

The main purpose of statistic is to do inference -- i.e., to say something relevant about the whole (unknwon) population (from which we have 
a small random sample $\{(Y_1, \boldsymbol{X}_i^\top)^\top;~i = 1, \dots, 32\}$) using only the infomarion (data) in the sample and the techniques 
of statistical inference. There are two main tools for the statistical inference -- condidence intervals for the uknown paramters and the hypothesis tests. 

While the confidence intervals must be calculated manually -- however, using just the standard output from the regression model -- some statistical 
tests are provided automaticaly. A statistical output of the multiple linear regression model fitted above can be obtained by the command `summary()`: 


```{r}
summary(fit1)
```

Four (the number of estimated parameters) statistical tests -- $p$-values of the tests are provided in the output (the fourth column of the table) by default.
In each case, the corresponding null hypothesis stantes that the estimated parameter (not the estimate, but the unknown true value) is equal to zero, i.e., 

$$
H_0: \beta_j = 0
$$
against a general alternative
$$
H_A: \beta_j \neq 0
$$
where $j \in \{1, 2, 3, 4\}$ corresponds to the line in the output table, respectively the element of the estimated paramter vector $\boldsymbol{\beta} \in \mathbb{R}^4$. 

<br><br>

### <font color = "#880000">Indidividual work</font> 
<ul>
<li>Recall, how a standard statistical test is provided. What is the corresponding test statistic, the distribution of the test statistic 
under the null hypothesis, or the critical region for the test? </li>
<li>Can you use some of the information provided in the output to reconstruct the quantities needed for the test (e.g., the test statistic?). </li>
<li>What should be some reasonable assumptions to be considered together with the test? </li>
<li>What is the interpretation of the statistical test given in the first line of the summary output (the line denoted as `Intercept`)?</li>
<li>What is the interpretation of the statistical tests for the remaining three lines of the summary output (the lines denoted as `disp`, `hp`, and `drat`)?</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

### <b>3. Model selection strategies</b> 
Considering the critical level of $\alpha = 0.05$, we can not reject the null hypothesis  $H_0; \beta_4 = 0$  which (in words) means that 
the consumption of the car is independent of the rear axle ratio. Thus, instead of the model with three predictor variablees we can consider 
a smaller (restricted) model with only the information on the displacement (`disp`) and the horse power (`hp`). 

```{r}
fit2 <- lm(mpg ~ disp + hp, data = mtcars)
summary(fit2)
```

In the restricted model, using again the critical value of $\alpha = 0.05$, we can again see that the null phypothesis $H_0: \beta_3 = 0$
is not rejected -- the corresponding $p$-value $0.0737$ is slightly higher that $0.05$. Thus, we can again consider a submodel, where the car consumption 
is only explained by the information about the engine displacement: 

```{r}
fit3 <- lm(mpg ~ disp, data = mtcars)
summary(fit3)
```

In the model above, both statistical tests return significant $p$-value and the model can not be further reduced. 
Thus, <b>we performed the model selection</b> procedure. From all considered models we selected one final model that is (in some 
well-defined sense) the most appropriate one. This particular selection strategy is called <b>the stepwise backward</b> selection strategy. 
<br><br>

#### <b>Stepwise backward selection</b> 
The idea is to fit a full model firstly and, consecutively, to reduce the model by taking out the less significant predictors. 
The final model is the one in which no other predictor can be added. This is the procedure that we just applied above. 
<br><br>

#### <b>Stepwise forward selection</b> 
The idea is reverse. We start with all possible (three in this case) simple regresion models (as we fitted them above) and we conclude that 
the most significant dependence is for the model with the dependence of the consumption on the displacement -- compare the $p$-values below: 

```{r}
summary(fit_disp)
summary(fit_hp)
summary(fit_drat)
```

Thus, the model that we use to proceed further is the model where the consumption is explained by the displacement. In the second step of 
the stepwise forward selection procedure we look for another covariate (out of two remaining ones, `hp`, or `drat`) which can be added into the model (meaning 
that the corresponding $p$-value of the test is significant). We consider two models:

```{r}
fit4 <- lm(mpg ~ disp + hp, data = mtcars)
summary(fit4)
```

```{r}
fit5 <- lm(mpg ~ disp + drat, data = mtcars)
summary(fit5)
```

In both cases we observe the $p$-value for the added covariate which is larger than $0.05$. Thus, the final model is the model denoted as `fit3` (which is, of course, equivalent to the model denoted as `fit_disp`).

### <font color = "#880000">Indidividual work</font> 
<ul>
<li>Note, that both model selection strategies (backward and forward) retured the same final model. This is, however, not always the case. </li>
<li>Discuss some obvious advantages and disadvantages of both selection strategies.</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>


#### <b>Model selection criteria</b> 
There are also other selection strategies based on various selection criteria. Typically used criteria are, for instance, 
the Akaike's Information Criterion (AIC) or Bayesian Information Criterion (BIC). But there are also others. 

For the AIC criterion, there is a standard R function `AIC()`: 
```{r}
AIC(fit3)

### compare also with the AIC for the remaining two (simple) regression models
AIC(fit_hp)
AIC(fit_drat)
```


The BIC criterion is implemented, for instance, in the library `flexmix` (installation with the command `install.packages("flexmix")` is needed)

```{r}
library(flexmix)
BIC(fit3)

### compare also with the BIC for the remaining two (simple) regression models
BIC(fit_hp)
BIC(fit_drat)
```

<br><br>

### <font color = "#880000">Indidividual work</font> 
<ul>
<li>Consider both, AIC and BIC criterion and use them to perform the model selection approach. The smaller the value the better the fit.</li>
<li>Use some liteature or references and find out how AIC and BIC are calculated. </li>
<li>What are some discussed advantages/disadvantages of all four model selection approaches?</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>


### <b>4. Categorical predictor variable</b> 

Finally, we will consider a regression dependence of a continuous variable $Y$ on some categorical covariate $X$ (taking more than 2 levels). 
We will consider the number of cylinders (there are four, six, and eight cylinder cars in the data set). Note, that the categorical covariate 
can be either assumed to be nominal or, also, ordinal. 

```{r, fig.width = 10}
boxplot(mpg ~ cyl, data = mtcars, col = "lightblue")
```

Compare and explain the following models: 

```{r}
fit6 <- lm(mpg ~ cyl, data = mtcars)
print(fit6)

fit7 <- lm(mpg ~ as.factor(cyl), data = mtcars)
print(fit7)

fit8 <- lm(mpg ~ cyl + I(cyl^2), data = mtcars)
summary(fit8)
```


<br><br>


### <font color = "#880000">Indidividual work</font> 
<ul>
<li>Try to use a diffent parametrization for the number of cylinders. For instance, consider values 4, 6 and 10
and refit models `fit6`, `fit7`, and `fit8`. </li>
<li>What are other possible parametrizations that can be used? Compare the following model with the model denoted as `fit7`. 
````{r}
mtcars$dummy1 <- as.numeric(mtcars$cyl == 4)
mtcars$dummy1[mtcars$cyl == 8] <- -1

mtcars$dummy2 <- as.numeric(mtcars$cyl == 6)
mtcars$dummy2[mtcars$cyl == 8] <- -1

fit9 <- lm(mpg ~ dummy1 + dummy2, data = mtcars)
print(fit9)
```
</li>

<li>Can you reconstruct model `fit7` from the model above (`fit9`)? Look at the output of the `summary()` function 
and explain the statistical tests that are provided for each column of the summary table (in model `fit7` and model `fit9`).</li>

<li>Think of some other parametrization that could be suitable and try to implement them in the R software. </li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

