---
title: "oct62021.Rmd"
author: "Nancy"
date: "06/10/2021"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## October 6 model building with `nuclear`



```{r data}
require(SMPracticals); data(nuclear); head(nuclear)
help(nuclear) 
```

You can also read about the data in Table 8.13 of Davison:

![Table 8.13 of SM](SMTable813.png)

As explained in class, and in SM, it makes sense to take logs of the continuous variables: `cost`, `cap`, `T1`, `T2`, `cum.n` (called $N$ in SM).

The full model is:
```{r modelling}
nuclear.lm <- lm(log(cost) ~ date + log(t1) + log(t2) + log(cap) + pr + ne + ct + bw + log(cum.n) + pt, data = nuclear)
summary(nuclear.lm)
```
I'll let you fill in the details, but backwards elimination, taking out the least signficant variable until the remaning variables are significant, leads to
```{r modelling part 2}
nuclear.lm2 <- lm(log(cost) ~ date + log(cap) + ne + ct + log(cum.n) + pt, data = nuclear)
summary(nuclear.lm2)
```
You can compare this to SM Table 8.14:
![Table 8.14 SM](SMTable814.png) 


Cox and Snell note that adding back in the omitted variables (in a different order) doesn't lead to an improved model. Using `step` will pick a model based on the smallest value of $AIC$:
```{r modelling part 3}
nuclear.lm2 <- step(nuclear.lm)
```
This model includes `log(T2)` and `pr` as well, and we can compare the estimated coefficients and their standard errors:
```{r}
summary(nuclear.lm2)
```
This is typical, AIC tends to pick out slightly bigger models than stepwise regression.

Another model-building trick is to fit all possible subsets of the covariates; of course we can only do this with relatively few covariates. There is a pretty fast algorithm for this in the `leaps` package called `regsubsets`: see LM-1 p.149 and LM-2 p.154. 

Finally, Cox \& Snell suggest that it would be good to see if the 6 observations with `pt`$=1$ have a different relation with cost, which can be assessed with interaction terms `pt*`$z$, where $z$ is one of the other variables retained. We can't include all the interaction terms in a single model, but we can add them one at a time. For example: 
```{r}
nuclear.lm3 <- update(nuclear.lm2, . ~ . + pt*log(cap))
summary(nuclear.lm3)
```
There's not much change in the other coefficients, so a final model could be `nuclear.lm2` or even `nuclear.lm3`. Residual and diagnostic plots don't show any surprises.

```{r plots}
par(mfrow = c(2,2))
plot(nuclear.lm2)
```

