---
title: "Linear Models"
author: "Julian Faraway"
date: "`r format(Sys.time(), '%d %B %Y')`"
output:
  github_document
---

```{r echo=FALSE}
library(knitr)
opts_chunk$set(comment=NA, warning=FALSE, fig.path="/tmp/Figs/")
```

We start by loading the `tidyverse` and the data.

```{r}
library(tidyverse)
load("hospdata.rda")
hospdata
```

The data for this example was generated using:

- air pollution data from station 8 in 2018
- one hospital located near this station where we take only the last two quarters of the year
- the two datasets were joined so that we know the air pollution readings on the date of admission for each patient
- a more extensive analysis would use more data - this is just for an example

The first few lines in the data:

```{r}
hospdata
```

and here is a summary:

```{r}
summary(hospdata)
```

A linear model expresses a linear relationship between a *response* and one or more *predictors*. For
example, we might model the `Cost` as a linear function of `Age` and `TreatDays` as follows:

```
Cost = beta0 + beta1 * Age + beta2 * TreatDays + error
```

The *parameters* are `beta0`, `beta1` and `beta2`. Usually, `beta0` is called the intercept. We do
not expect that the data will obey this relation exactly so we include an error term.

**Before fitting a model, it is important to investigate the data in an exploratory way with plots and numerical
summaries. We have already done much of this in the previous modules**

# Fitting the model

We can fit this model with:

```{r}
lmod = lm(Cost ~ Age + TreatDays, hospdata)
```

The model formula is `Cost ~ Age + TreatDays`. We can get this from the mathematical formula 
above by just deleting the parameters and the error. The model formula is an easier way to represent the same model.

We can extract estimates of the parameters:

```{r}
coef(lmod)
```

or we can view a more detailed summary of the model fit:


```{r}
summary(lmod)
```

## Exercises

Add `NO2` to the model as a predictor. Compute the parameter estimates.

# Checking the model

Before using a model to explain a relationship or predict future outcomes, we should check the model
is an adequate fit to the data. The best way to do this is to make diagnostic plots.

We will find the `broom` package useful for this purpose. 

```{r}
library(broom)
```

Now we compute the diagnostic quantities for the model.

```{r}
dmod = augment(lmod)
dmod
```

You may recognize some of these quantities if you have done linear modelling in the past.

The first and most important plot shows the relationship between the fitted values (the response
predicted by the model) and the residuals (the difference between the observed response and
the fitted values). This plot checks two assumptions of the linear model:

1. The linear relation between the response and the predictors. If this is true, the plot
will show no relationship between the fitted values and the predictors.

2. The errors should vary equally. If this is true, the variation of the residuals should
be constant as the fitted values vary.

Here is the so-called residuals-fitted plot for this model:

```{r}
ggplot(dmod,aes(.fitted,.resid)) + geom_point() + geom_hline(yintercept = 0) 
```

We see that the second assumption is broken because the variation in the residuals increases as the fitted values increase. We also suspect that the first assumption is also broken.

The second diagnostic plots you should make checks for the normality of the errors. We do this using a Q-Q plot which is constructed as follows:

```{r}
ggplot(dmod,aes(sample=.resid)) + stat_qq() + stat_qq_line()
```

If the residuals are normally distributed, they will stay close to the line drawn on the plot. In this case, they do not. Large positive residuals are much bigger than expected, indicating that the model does not account well for incomes that are much larger than expected.

Our model is inadequate. From our previous plots of the data, we know that the `Cost` is skewed and that a log-transformation is helpful. We have added one to avoid the problem of log(0). Let's do that here and repeat the plots:

```{r}
lmod = lm(log(Cost+1) ~ Age + TreatDays, hospdata)
dmod = augment(lmod)
ggplot(dmod,aes(.fitted,.resid)) + geom_point() + geom_hline(yintercept = 0) 
```

Why are there lines of residuals on the plot? The transformation has helped.

Let's try the Q-Q plot

```{r}
ggplot(dmod,aes(sample=.resid)) + stat_qq() + stat_qq_line()
```

This is improved. We do not require the points to follow the line exactly. Approximate normality is sufficient for our purposes.

We can try plotting the residuals and the predictors. First we try `Age`

```{r}
ggplot(dmod,aes(Age,.resid)) + geom_point() + geom_hline(yintercept = 0) 
```

This is an example of a perfect plot - no pattern and equal variance. So no problem with `Age`. Let's try `TreatDays`:

```{r}
ggplot(dmod,aes(TreatDays,.resid)) + geom_point() + geom_hline(yintercept = 0) 
```

No problem here.


We should also check for unusual points. Large residuals are an example of this. We would see these
in the Q-Q plot. We did not seen any such points.

Leverage measures how extreme a point is in the space of predictors. We can plot these with
the standardized residuals. Leverages are also called hat-values.

```{r}
ggplot(dmod, aes(.hat,.std.resid)) + geom_point()
```

We can see that there are two points with large leverage. This means they have unusual predictor values. We used the standardized residuals since this makes it easier to judge the size. 

*Cook statistics* measure influential points. A large value indicates a combination of large leverage and large residual.
We can plot these as:


```{r}
ggplot(dmod,aes(seq_along(.cooksd),.cooksd)) + geom_col() + xlab("Case Number") + ylab("Cook Statistic")
```

We see that there is several large Cook statistics. Let us investigate:

```{r}
dmod %>% arrange(desc(.hat))
```

We have listed the diagnostic information in descending order of the leverage. Can you see what is unusual about the first (and second) case? If you identify influential points, it is best to see how much the model would change if you eliminated these points:

Here are the parameters estimates with all the data:

```{r}
coef(lmod)
```

and here are the parameter estimates with only `TreatDays <40`.

```{r}
lmod2 = lm(log(Cost+1) ~ Age + TreatDays, hospdata, subset = TreatDays < 40)
coef(lmod2)
```

Naturally, they are different but are they different enough to concern us? It is hard to judge now. Bear in mind that dropping cases from your data is a serious choice that you will need to justify.

## Exercises

Fit a linear model with the formula: `TreatDays ~ Age + NO2`. Recreate all the plots above and
interpret them.

# Interpreting the model

Our current model is

```{r}
lmod = lm(log(Cost+1) ~ Age + TreatDays, hospdata)
```

and we can obtain a summary of the output as

```{r}
summary(lmod)
```

First we might question whether any of the predictors have a relationship
to the response. For this we look at the *p-values* - Pr(>|t|) in the output.

- The p-value for `TreatDays` is 1e-04. This is quite small so we conclude
that `TreatDays` has a statistically significant relationship to the response. 

- The p-value for `Age` is 0.9. Since this is larger (greater than 0.05) we conclude
that this predictor is not helpful in predicting the response *if* we already know
`TreatDays`.

- The p-value for the intercept is of no interest. We usually do not care about
the statistical significance of the intercept since we plan to include it in the model
regardless.

Now consider the meaning of the parameter estimates starting with `Age`. The
coefficient is -0.00016.

- If `Age` increases by 1, log(Cost) is predicted to decrease by 0.00016.

- Decreasing log(ingreso) by 0.00016 means that ingreso is multiplied by exp(-0.00016) 

- exp(-0.00016)=0.99984 so the cost is reduced by 0.016%

- We see that the effect of age is very small

The effect of treatment days is

```{r}
exp(0.02814)
```

Each extra day raises the cost by 3%. Also not a large effect.


Overall, there is evidence that this model does not predict very
accurately:

- The R-squared of 0.0117 is very small

- The residual standard error of 0.832 corresponds to a multiplier of exp(0.832) which
is about 2.3. In other words, our predictions will tend to be wrong by a factor of 2. 

## Exercises

Continue with linear model with the formula: `TreatDays ~ Age + NO2`. Interpret
the parameter estimates. Describe the fit of the model.


# Prediction

The fitted values are obtained by using the predictor values in the data
to predict the response using the fitted model. These are useful in assessing
the model but we want to predict the future and not the past (which we already know).

Consider a new household with `Age`=40 and `TreatDays`=5. We can make a prediction
of `log(Cost+1)` with:

```{r}
predict(lmod, tibble(Age=40,TreatDays=5))
```

If we want this on the original scale of income, we need to exponentiate:

```{r}
exp(predict(lmod, tibble(Age=40,TreatDays=5)))-1
```

We can construct a 95% prediction interval for a single patient:

```{r}
exp(predict(lmod, tibble(Age=40,TreatDays=5), interval="prediction"))-1
```

We see that this interval is wide - we cannot predict with much certainty.

We might also ask what is mean cost for patients aged 40 staying for 5 days:

```{r}
exp(predict(lmod, tibble(Age=40,TreatDays=5), interval="confidence"))-1
```

The point prediction is the same but the interval is much narrower because
we are asking about a mean rather than a single observation. Usually, the prediction
interval is what is really wanted.

Often we want to make more than one prediction. Suppose we want to predict the costs
for a range of pa where edad=40. First, we construct the predictor values:

```{r}
x = tibble(Age = 40, TreatDays=1:20)
x
```

Now we make the predictions with confidence. We bind on the predictor values for convenience:

```{r}
predv = as.tibble(exp(predict(lmod, x, interval="confidence"))-1)
predv = bind_cols(predv,x)
predv
```

We can make a plot of this:

```{r}
ggplot(predv, aes(x=TreatDays,y=fit, ymin=lwr, ymax=upr)) + geom_line() + geom_ribbon(alpha=0.25)
```

The `geom_ribbon` geometry adds the shaded area - we need to make it transparent using `alpha=0.25` or the area will be black. We can see how the cost is related to the treatment days. Notice how the bands become wider for larger longer stays where there is not much data.

It is interesting to make a similar plot for 'Age':

```{r}
x = tibble(Age = 0:80, TreatDays=5)
predv = as.tibble(exp(predict(lmod, x, interval="confidence"))-1)
predv = bind_cols(predv,x)
ggplot(predv, aes(x=Age,y=fit, ymin=lwr, ymax=upr)) + geom_line() + geom_ribbon(alpha=0.25)
```

We see that the response does not vary much. A constant response would fit within the bands
- this is a sign that the predictor is not statistically significant.

## Exercises

Continue with linear model with the formula: `TreatDays ~ Age + NO2`. Make a prediction for a patient with age=40 and NO2=10. Construct a plot that shows how the predicted response changes as `Age` varies with NO2 fixed at 10.

# Categorical Variables

Let's add a categorical predictor to the model - `Gender`:


```{r}
lmod = lm(log(Cost+1) ~ Age + TreatDays + Gender, hospdata)
coef(lmod)
```

Why is there only Male and not Female? The coefficient of -0.124 represents
the difference between men and women. We predict that for the same level of `Age`
and `TreatDays`, the predicted `log(Cost)` of men is -0.124 relative to women.

Let's see how that works - we predict the outcome for two households which
are the same except one is headed by a man and one by a woman:

```{r}
x = tibble(Age=40,TreatDays=5,Gender=c("Female","Male"))
p = predict(lmod, x)
p
```

The difference between the two predictions is

```{r}
p[2] - p[1]
```

which is precisely the coefficient seen in the model.

The level which does not appear in the coefficients is called the *reference* level. R
chooses the first alphabetical level by default. We can change this as seen below. The
function `fct_relevel` comes from the `forcats` package which is another
member of the `tidyverse`.

```{r}
hospdata %>% mutate(Gender2 = fct_relevel(Gender,"Male")) -> hospdata
lmod2 = lm(log(Cost+1) ~ Age + TreatDays + Gender2, hospdata)
coef(lmod2)
```

Now `Male` is the reference level and we see that women are predicted
to be 0.124 less than men. The outcome is essentially the same -
you can choose either one.

Since exp(0.124) = 1.132, we can say males cost 13% more than females.

Let's see how this works for categorical predictors with more
than two levels using a four level version of the illness as created
in the merging and tidying module:

```{r}
hospdata %>% mutate(shortill = str_sub(Illness,1,3)) %>% 
  mutate(shortill = ifelse(shortill %in% c("J18","J34","J44"),shortill,"other")) -> hospdata
table(hospdata$shortill)
```

Now add this variable to the model:

```{r}
lmod = lm(log(Cost+1) ~ Age + TreatDays + Gender + shortill, hospdata)
summary(lmod)
```

The illness with code J18 is the reference level. We can see the relative costs to
this level with:

```{r}
exp(coef(lmod)[5:7])
```

We see that illness J34 results in 69% higher costs.

Are the differences between the languages statistically significant? We might look
at the summary but this is not helpful in answering the question. The test are all pairwise
with respect to the reference level. We need to make an F-test:

```{r}
drop1(lmod, test="F")
```

We see that the illness is very significant.

## Exercises

We can analyse all the illness codes rather than grouping the less common codes with:

```{r}
hospdata %>% mutate(shortill = str_sub(Illness,1,3)) -> hospdata2
```

What difference does this make to the model?

# Model Selection

Now let's put all the variables into the model (except `fex` because this is a weight not a predictor).

```{r}
lmod = lm(log(Cost+1) ~ Age + TreatDays + Gender + shortill + NO2 + PM10, na.omit(hospdata))
summary(lmod)
```

We see from the R-squared of 0.388 that the predictive ability of the model is much improved over
the models used previously. Even so, it is still not very accurate.

We might ask whether all these variables are necessary. Can we make the same or better
predictions with fewer variables?

The *Akaike Information Criterion* (AIC) is commonly used to select a good model. The AIC
is a measure that balances the fit of the model against complexity. We want a good fit without
using too many predictors. We may find a better choice of model by optimising AIC which can be
done like this:

```{r}
lmod2 = step(lmod)
summary(lmod2)
```

The process removed age and NO2. Note that although Gender and PM10 are
not statistically signifcant but are still retained.

In analyses where explanation is important, model selection is best done with user input and not
using an automated method such as AIC. We include this here to show what to do when the purpose
the analysis is prediction.

## Exercises

We previously identified two very long staying patients as having an influence
on the fit. Can you check how much influence they have on the model selection?






