---
title: "Linear Models"
author: "Mathieu Ribatet"
format: html
editor: visual
---

![](Varsovie.jpg){fig-align="center" width="703"}

## 🏘️ Apartment prices in Warsaw

In this Lab we are going to play with a data set which contains the prices of the apartments which were sold in Warsaw, Poland, during the calendar years 2007 to 2009. Overall we have 409 observations on the following 6 variables :

-   **surface** area of the apartment in square meters.

-   **district** a factor corresponding to the district of Warsaw with levels `Mokotow`, `Srodmiescie`, `Wola` and `Zoliborz`

-   **n.rooms** number of rooms in the apartment

-   **floor** floor on which the apartment is located

-   **construction.date** year that the apartment was constructed

-   **areaPerMzloty** area in square meters per million zloty

The dataset can be retrieved from the following piece of code:

```{r}
rm(list=ls())##clear memory before working
url <- "https://mribatet.perso.math.cnrs.fr/CentraleNantes/Data/Warsaw.csv"
data <- read.csv(url, row.names = 1)
```

*Hint: It may be tedious to write code like `data$areaPerMzloty` so we might want to "attach" the data set to your R session so that we can just use variable names, e.g., `areaPerMzoloty`* . This is done using the following code

```{r}
attach(data)
```

Of course you can "detach" the data set by invoking

```{r eval=FALSE}
detach(data)
```

1.  Familiarized yourself with the dataset. What are the type of the variables? Which ones are the explanatory variables, i.e., $x_j$ in the lecture ? Which one is the response variable, i.e., $Y$ in the lecture ?

2.  Plot the evolution of the price w.r.t. the district and comment.

    ```{r}
    ##Write code here
    ```

3.  Plot the evolution of the price w.r.t. construction date and comment.

    ```{r}
    ##Write code here
    ```

4.  Plot the evolution of the price w.r.t. the number of rooms and comment.

    ```{r}
    ##Write code here
    ```

5.  Plot the evolution of the price w.r.t. the surface and comment.

    ```{r}
    ##Write code here
    ```

6.  Read the documentation of the `lm` function and fit your first linear model where you model the price as a function of all the covariates (you will store it in a object called `fit`)

    ```{r}
    ##Write code here
    ```

7.  Run the following code and say what it does. Compute by yourself the hat matrix and compute its trace. Was it expected? Compute the sum of each row? Again was it expected?

    ```{r eval=FALSE}
    X <- model.matrix(fit)

    ##Write code here
    ```

8.  Run the following code and comment. Make sure you understand every single output.

    ```{r eval=FALSE}
    summary(fit)
    ```

    What are the meaning of the "\*"? Is the `surface` statistically significant? Why do we have 401 degrees of freedom for the residual standard error? Look at the bottom line where one can read `p-value: < 2.2e-16`, what does it say?

9.  Why are we having lines named `districtSrodmiescie`, `districtWola`and `districtZoliborz`? Why there is no `districtMokotow`? Is the variable `district` statistically significant?

10. What can you say about your fitted model? Produce diagnostic plots using the following code

    ```{r eval=FALSE}
    plot(fit)
    ```

11. We see that our model is not very accurate. We may need to work with a more complex model by considering "new variables" such as `log(surface)`. Nowadays, this stage is known as *feature engineering*–a few years ago we didn't use a name for such a basic stage ;-)

    Using the `update` function, try to get a better fitted model by incorporating new variables; you will call this new model `new_fit`.

    For instance, you may want to consider the log, square or even product between two variables such as the `surface * n.rooms` (in statistics it is called an **interaction**). You may also want to "merge" some levels of a categorical variable, e.g., merging levels `Wola`, `Zoliborz` and `Mokotow` into a single level.

    ```{r}
    ## Write code here
    ```

12. Perform model selection between models `fit` and `new_fit` using an information criterion and, if relevant, analysis of variance.

    ```{r}
    ##Write code here
    ```

13. Read the documentation of the `step` function and next perform stepwise model selection. From this output plot the price of the apartments to predictions. Interpret the effect of the `floor` variable on the price.

    ```{r}
    ##Write code here
    ```

14. Perform residual analysis for your best model. Comment and do a scatterplot of prices w.r.t. predictions (the function `predict` may be useful for this purpose).

    ```{r}
    ## Write code here
    ```

# 🎁 Bonus: Additive models

It is a bit unfortunate that your best linear model is not extremely powerful. The linear model is widely used and highly interpretable but often lacks some flexibility. More flexible models exist but the price to pay is a decrease in interpretability. Neural networks are the typical example of powerful models but not interpretable at all. A good compromise are **additive models**:

$$
Y = f_1(X_1) + \ldots + f_p(X_p) + \varepsilon, \qquad \varepsilon \sim N(0, \sigma^2),
$$

where $f_j$ are unknown functions to be estimated. We will not cover the theory beyond this model but I just want you to know that those functions are defined in a given basis and, as so, parameters of the above models are the coefficient in that basis for each function $f_j$ (and consequently we get a linear smoother).

Fitting such model is a piece of cake using the `mgcv` library (that comes with R but is not loaded in memory):

```{r}
library(mgcv)
add_mod <- gam(areaPerMzloty ~ s(surface) + s(construction.date) + district + s(floor), data = data)
summary(add_mod)
```

Without knowing the theory (which is unfortunate of course), we can guess what's going on with the above output. For instance one can see that the (adjusted) $R^2$ is around 45% and that the `floor` variable, even within a additive form, is not statistically significant. The variables `surface` and `construction.date` are strongly significant. So we can update our model by removing the `s(floor)` term and next have a look the effect of the two aforementioned variables.

```{r}
add_mod <- update(add_mod, . ~ . - s(floor))
par(mfrow = c(1, 2), mar = c(4, 5, 0.5, 0.5))
plot(add_mod)
```

From the above (right) graph we can see that the cheapest apartments are those built around 1980 while much older or recent ones are more expensive.
