---
title: "Ozone dataset"
output: html_notebook
---

We will (briefly) analyze the *airquality* dataset. Here we will try to predict the ozone concentration level w.r.t. some features.

## Descriptive analysis

## Linear modelling

Here we will use a simple model
$$
\mathbb{E}(Y \mid X) = \beta_0 + \beta_1 \times \text{solar radiation} + \beta_2 \times \text{wind speed} + \beta_3 \times \text{temperature}
$$
```{r}
fit <- lm(log(Ozone) ~ Solar.R + Wind + Temp,
          data = airquality)
summary(fit)


```
From the output we can see that the ozone level depends on the the solar radiation (since p-value is very very small).

Now do we need wind speed or temperature? For this we need to perform analysis of variance
```{r}
fit1 <- lm(log(Ozone) ~ Solar.R + Wind + Temp,
          data = airquality)
fit2 <- lm(log(Ozone) ~ Solar.R, data = airquality)
anova(fit2, fit1)
```

From the above output we should rather use the most complex model, i.e., *fit1*.

**Remark: When doing anova, models need to be NESTED**

## Model checking

Here we will perform residual analysis. Remember here we have to check 3 assumptions:
  - mean of error is 0
  - constant variance (homoscedasticity)
  - normality of error
  
In *R* you just invoke
```{r}
plot(fit1)
```

Observation number 21 appears to have a large error. Unfortunately Cook's distance is smaller than 0.5 indicating that, from a statistical standpoint, there is no leverage issue with it. But still I will try to see what happens when we drop it from our dataset
```{r}
fit1_sub <- update(fit1, subset = -21)
plot(fit1_sub)
```

An attempt to fix the heteroscedasticity issue, i.e., not working on the log of ozone levels.
```{r}
attempt <- lm(Ozone ~ Solar.R + Wind + Temp, data = airquality)
plot(attempt)
```

We can see a quatratic trend on the first plot, suggesting that maybe squaring some feature may be relevant.
```{r}
fit_square <- lm(Ozone ~ Solar.R + Wind + Temp + 
                   I(Solar.R^2) + I(Wind^2) + I(Temp^2),
          data = airquality)
summary(fit_square)
```
Interestingly, the use of quadratic term seems very relevant here and we finally work with the following model
```{r}
final_model <- update(fit_square, . ~ . - Solar.R - I(Solar.R^2))
summary(final_model)
```

Let's do model checking
```{r}
plot(final_model)
```
Observation 117 might have a leverage impact. Discard it and stop please the analysis.
```{r}
final_model_I_swear <- update(final_model, subset = -117)
plot(final_model_I_swear)
```

## Predictions and interpretations

To do predictions from my best model, I just use the *predict* function. Invoking
```{r}
predict(final_model_I_swear)
```
will give me predictions for all my dataset. But if I want to do prediction for a new observation I need to do
```{r}
new_obs <- data.frame(Wind = 100, Temp = 60)
predict(final_model_I_swear, newdata = new_obs)
```

From the outpout below
```{r}
final_model_I_swear
```

**All things being equal**, a one unit increase in wind speed decrease (in mean) the amount of ozone level by around 10. 
