---
title: "Sep222021.Rmd"
author: "Nancy"
date: "22/09/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(faraway)
data(prostate)
options(digits=4)
head(prostate)
```
Continuing with prostate example (see also Sep 15 Rmd file). Here we will look mainly at F-Tests, as well as factor variables. 

First, model 1 is the original linear model with all the covariates: 
```{r full model}
model1 <- lm(lpsa ~ ., data = prostate)
summary(model1)
summary(prostate$gleason)
prostate$gleason
``` 
We should have done some data exploration first, but speeding along I'll just note that the variable `gleason` has just 4 values. I'll make a new variable that treats this as a factor variable
```{r factors}
prostate2 <- prostate[,-7]
prostate2$gleason_factor <- as.factor(prostate$gleason)
summary(prostate2$gleason)
```
Note that this is a different type of variable, but you can't really tell this from the data frame:
```{r}
head(prostate2)
```
You can tell though with `is.factor(prostate2$gleason_factor)`, for example. 
```{r factor model}
model2 <- lm(lpsa ~ ., data = prostate2)
summary(model2)
model3 <- update(model2, . ~ . - gleason_factor)# same model as before, but minus g_f
summary(model3)
```
Note that we now have used 3 parameters to fit the factor variable; each parameter measures the difference between the baseline level (here "6") and the others. This is very suspect, there is only 1 observation with a score of 8, and very few with scores of 9. 
```{r}
anova(model2,model3) # anova(model3,model2) looks a little nicer
```
The effect of `gleason_factor` is assessed by an $F$-test, with $3$ and $88$ numerator and denominator degrees of freedom, respectively. 

This is one example of comparing two models. Another example done in class is less interesting:
```{r}
anova(model1,lm(lpsa ~ lcavol + lweight + age + svi + lbph, data = prostate ))
```
The omitted variables do not reduce the residual SS very much, but these were chosen based on the output of the first fitted model, so this test is not very useful. Neither is the test that all coefficients are equal to zero (it has a p-value of $10^{-6}$). 

Another setting where $F$-tests are needed, as opposed to just the summary output from `lm`, is when fitting "smooth" functions of covariates. As you will recall from "IJALM", these are also linear regression models, even though they're sometimes described as nonparametric. 
```{r smoothing}
require("splines")
library(splines)
ns(prostate$age,4)[1:10,] # first 10 rows of this constructed variable
model4 <- lm(lpsa ~ .-age + ns(age,4), data = prostate, )
summary(model4)
```
Note that although some of the coefficients have a "*", it doesn't make sense to look at them individually, only as a group. 
```{r}
anova(model4,lm(lpsa ~ . - age, data = prostate))
```
Have a look at `model.matrix` of model4 to see the extra 4 columns that got added with the `ns` command. 

There was a question about the difference between unordered and ordered factors. Ordered factors are treated differently in `lm` and other methods built on `lm`: the model matrix includes values that correspond to linear, quadratic, cubic, etc. polynomials in the underlying variable. The coefficients of these linear, quadratic, cubic etc. polynomials are computed on the fly. (See the documentation for `contrasts`.)
```{r ordered factors}
is.ordered(prostate2$gleason_factor)
prostate2$go <- as.ordered(prostate2$gleason_factor)
head(prostate2)
model5 <- lm(lpsa ~ . - gleason_factor, data = prostate2)
summary(model5)
X2 <- model.matrix(model5)
X2[1:10,]
contrasts(prostate2$go)
```
Note the 3 terms for `gleason_factor` are now called `.L`, `.Q`, `.C`, for linear, quadratic, cubic. These could in principle be tested one at a time, because they represent successively more complicated relationships between Gleason score and `lpsa`. But this data is not very suitable for illustrating ordered factors, so I'll come back to it when it is more relevant.