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

#### 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*pf3d7) - n*xb*yb
sxy
# alternatively
sum((h2o2-xb) * (pf3d7-yb))
sxx <- sum(h2o2^2) - n*xb^2
sxx
# alternatively
sum( (h2o2-xb)^2 )
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
```

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

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

##### Getting estimated standard errors for the regression parameters

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))
```

Compare to the output from the 'lm' function.

```{r}
lm.sum$coefficients
lm.sum$coefficients[,2]
```

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

The residual degrees of freedom.

```{r}
df <- lm.out$df.residual 
df
```

The t-statistic for a 95% confidence interval.

```{r}
tmult <- qt(0.975, df)  
tmult
```

The parameter estimates.

```{r}
beta <- lm.sum$coef[,1] 
beta
```

The estimated standard errors.

```{r}
se <- lm.sum$coef[,2] 
se
```

The 95% confidence interval for the intercept.

```{r}
beta[1] - tmult*se[1]
beta[1] + tmult*se[1]
beta[1] + c(-1,1)*tmult*se[1]
```

The 95% confidence interval for the slope.

```{r}
beta[2] - tmult*se[2]
beta[2] + tmult*se[2]
beta[2] + c(-1,1)*tmult*se[2]
```

Compare to the output from the 'confint' function.

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