```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=TRUE, fig.align="center")
options(width=110)
```

#### Statistics for Laboratory Scientists ( 140.615 )

## Linear Regression - Advanced

#### Example from class: David Sullivan's heme data

```{r}
h2o2 <- rep(c(0,10,25,50),each=3)
pf3d7 <- c(0.3399,0.3563,0.3538,
           0.3168,0.3054,0.3174,
           0.2460,0.2618,0.2848,
           0.1535,0.1613,0.1525)
```

##### The regression of optical density on hydrogen peroxide concentration, from scratch

```{r}
n <- length(h2o2)
n
yb <- mean(pf3d7)
yb
xb <- mean(h2o2)
xb
sxy <- sum((h2o2-xb)*(pf3d7-yb))
sxy
sxx <- sum((h2o2-xb)^2)
sxx
b1hat <- sxy/sxx
b1hat
b0hat <- yb-b1hat*xb
b0hat
yhat <- b0hat + b1hat*h2o2
yhat
rss <- sum((pf3d7-yhat)^2)
rss
sigmahat <- sqrt(rss/(n-2))
sigmahat
```

##### Testing parameters

Testing whether the intercept is equal to zero.

```{r}
se.b0hat <- sigmahat*sqrt(1/n+xb^2/sxx)
se.b0hat
t.stat <- b0hat/se.b0hat
t.stat
p <- 2*pt(-abs(t.stat),n-2)
p
```

Testing whether the slope is equal to zero.

```{r}
se.b1hat <- sigmahat/sqrt(sxx)
se.b1hat
t.stat <- b1hat/se.b1hat
t.stat
p <- 2*pt(-abs(t.stat),n-2)
p
```

Compare to the results from the built-in 'lm' function.

```{r}
lm.out <- lm(pf3d7 ~ h2o2)
lm.sum <- summary(lm.out)
lm.sum$coef
lm.sum$coef[,1]
lm.sum$coef[,2]
lm.sum$coef[,3]
lm.sum$coef[,4]
```

##### The variance-covariance matrix

The parameter estimate variance-covariance matrix (not including sigma).

```{r}
lm.sum$cov.unscaled   
```

The square root of the diagonal elements.

```{r}
diag(lm.sum$cov)
sqrt(diag(lm.sum$cov))
```

The estimated parameter standard errors.

```{r}
lm.sum$sigma * sqrt(diag(lm.sum$cov))
```

##### Building confidence intervals for the regression parameters

```{r}
b0hat + c(-1,1)*qt(0.975,n-2)*se.b0hat
b1hat + c(-1,1)*qt(0.975,n-2)*se.b1hat
```

Compare to the output from the 'confint' function.

```{r}
confint(lm.out)
```

##### Checking model assumptions

```{r}
lm.out$fitted
lm.out$residuals
```

Alternatively, there are also generic functions.

```{r}
fitted(lm.out)
residuals(lm.out)
```

The residual qq plot.

```{r}
qqnorm(lm.out$residuals, main="")
qqline(lm.out$residuals, col="blue", lty=2)
```

Fitted values versus residuals.

```{r}
plot(lm.out$fitted, lm.out$residuals, pch=1, xlab="fitted values", ylab="residuals")
abline(h=0, col="blue", lty=2)
```