<!--- rmarkdown::render("/Volumes/Files/courses/cruk/LinearModelAndExtensions/git_linear-models-r/multiple_regression+.Rmd") --->

```{r message = FALSE, warning = FALSE, echo = FALSE} 
# change working directory: should be the directory containg the Markdown files:
#setwd("/Volumes/Files/courses/cruk/LinearModelAndExtensions/git_linear-models-r/")

```

# Section 1: Multiple Regression

The in-built dataset `trees` contains data pertaining to the `Volume`, `Girth` and `Height` of 31 felled black cherry trees. In the Simple Regression session, we constructed a simple linear model for `Volume` using `Girth` as the independent variable. Now we will expand this by considering `Height` as another predictor.

Start by plotting the dataset:
```{r}
plot(trees)
```

This plots all variables against each other, enabling visual information about correlations within the dataset.

Re-create the original model of `Volume` against `Girth`:
```{r}
m1 = lm(Volume~Girth,data=trees)
summary(m1)
```

Now include `Height` as an additional variable:
```{r}
m2 = lm(Volume~Girth+Height,data=trees)
summary(m2)
```

Note that the R^2 has improved, yet the `Height` term is less significant than the other two parameters.

Try including the interaction term between `Girth` and `Height`:
```{r}
m3 = lm(Volume~Girth*Height,data=trees)
summary(m3)
```

All terms are highly significant. Note that the `Height` is more significant than in the previous model, despite the introduction of an additional parameter.

We'll now try a different functional form - rather than looking for an additive model, we can explore a multiplicative model by applying a log-log transformation (leaving out the interaction term for now).
```{r}
m4 = lm(log(Volume)~log(Girth)+log(Height),data=trees)
summary(m4)
```

All terms are significant. Note that the residual standard error is much lower than for the previous models. However, this value cannot be compared with the previous models due to transforming the response variable. The R^2 value has increased further, despite reducing the number of parameters from four to three.
```{r}
confint(m4)
```

Looking at the confidence intervals for the parameters reveals that the estimated power of `Girth` is around 2, and `Height` around 1. This makes a lot of sense, given the well-known dimensional relationship between `Volume`, `Girth` and `Height`!

For completeness, we'll now add the interaction term.
```{r}
m5 = lm(log(Volume)~log(Girth)*log(Height),data=trees)
summary(m5)
```

The R^2 value has increased (of course, as all we've done is add an additional parameter), but interestingly none of the four terms are significant. This means that none of the individual terms alone are vital for the model - there is duplication of information between the variables. So we will revert back to the previous model.

Given that it would be reasonable to expect the power of `Girth` to be 2, and Height to be 1, we will now fix those parameters, and instead just estimate the one remaining parameter.
```{r}
m6 = lm(log(Volume)-log((Girth^2)*Height)~1,data=trees)
summary(m6)
```

Note that there is no R^2 (as only the intercept was included in the model), and that the Residual Standard Error is incomparable with previous models due to changing the response variable.

We can alternatively construct a model with the response being y, and the error term additive rather than multiplicative.
```{r}
m7 = lm(Volume~0+I(Girth^2):Height,data=trees)
summary(m7)
```

Note that the parameter estimates for the last two models are slightly different... this is due to differences in the error model.

# Section 2: Model Selection

Of the last two models, the one with the log-Normal error model would seem to have the more Normal residuals. This can be inspected by looking at diagnostic plots, by and using the `shapiro.test()`:
```{r}
fitted = fitted(m6)
resid = resid(m6)
plot(fitted, resid, xlab="Fitted values", ylab="Raw residuals")
fitted = fitted(m7)
resid = resid(m7)
plot(fitted, resid, xlab="Fitted values", ylab="Raw residuals")
shapiro.test(residuals(m6))
shapiro.test(residuals(m7))
```

The Akaike Information Criterion (AIC) can help to make decisions regarding which model is the most appropriate. Now calculate the AIC for each of the above models:
```{r}
summary(m1)
AIC(m1)
summary(m2)
AIC(m2)
summary(m3)
AIC(m3)
summary(m4)
AIC(m4)
summary(m5)
AIC(m5)
summary(m6)
AIC(m6)
summary(m7)
AIC(m7)
```

Whilst the AIC can help differentiate between similar models, it cannot help deciding between models that have different responses. Which model would you select as the most appropriate?


# Section 3: Practical Exercises

## Puromycin

The in-built R dataset `Puromycin` contains data regarding the reaction velocity versus
substrate concentration in an enzymatic reaction involving untreated cells or cells
treated with Puromycin.

- Plot `conc` (concentration) against `rate`. What is the nature of the relationship
between `conc` and `rate`?
```{r message = FALSE, warning = FALSE, echo = TRUE}
plot(conc~rate,data=Puromycin)
# There is a non-linear positive relationship between conc and rate
```

- Find a transformation that linearises the data and stabilises the variance,
making it possible to use linear regression. Create the corresponding linear
regression model. Are all terms significant?
```{r message = FALSE, warning = FALSE, echo = TRUE}
plot(log(conc)~rate,data=Puromycin)
m10 = lm(log(conc)~rate,data=Puromycin)
fitted = fitted(m10)
resid = resid(m10)
plot(fitted, resid, xlab="Fitted values", ylab="Raw residuals")
summary(m10)
# Both terms are significant
```

- Add the `state` term to the model. What type of variable is this? Is the
inclusion of this term appropriate?
```{r message = FALSE, warning = FALSE, echo = TRUE}
m11 = lm(log(conc)~rate+state,data=Puromycin)
fitted = fitted(m11)
resid = resid(m11)
plot(fitted, resid, xlab="Fitted values", ylab="Raw residuals")
summary(m11)
# `state` is a boolean factor or indicator variable
# The inclusion of `state` is appropriate, as the term is significant and the diagnostic plots look reasonable 
```

- Now add a term representing the interaction between `rate` and `state`. Are all
terms significant? What can you conclude?
```{r message = FALSE, warning = FALSE, echo = TRUE}
m12 = lm(log(conc)~rate*state,data=Puromycin)
summary(m12)
# The `state` term is not significant when the interaction between `rate` and `state` is included in the model. So it may be better to remove the `state` term from the model.
```

- Given this information, create the regression model you believe to be the most
appropriate for modelling `conc`. Regenerate the plot of `conc` against `rate`.
Draw curves corresponding to the fitted values of the final model onto this
plot (note that two separate curves should be drawn, corresponding to the
two levels of `state`).
```{r message = FALSE, warning = FALSE, echo = TRUE}
m13 = lm(log(conc)~rate+rate:state,data=Puromycin)
summary(m13)

# Solution one:
plot(conc~rate,data=Puromycin)
idx = order(Puromycin$rate)
treated = Puromycin$state[idx] == "treated"
untreated = Puromycin$state[idx] == "untreated"
lines(exp(fitted(m13))[idx][treated]~Puromycin$rate[idx][treated])
lines(exp(fitted(m13))[idx][untreated]~Puromycin$rate[idx][untreated],col="red")

# Solution two (better - more general):
plot(conc~rate,data=Puromycin)
xvals = range(Puromycin$rate)[1]:range(Puromycin$rate)[2]
lines(exp(coef(m13)[1] + coef(m13)[2]*xvals) ~ xvals)
lines(exp(coef(m13)[1] + coef(m13)[2]*xvals + coef(m13)[3]*xvals) ~ xvals, col="red")
```