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

#### Statistics for Laboratory Scientists ( 140.615 )

## Prediction and Calibration

#### Example: 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)
dat <- data.frame(conc=h2o2, od=pf3d7)
```

The regression of optical density on hydrogen peroxide concentration.

```{r}
lm.fit <- lm(od ~ conc, data=dat)
lm.out <- summary(lm.fit)
lm.out$coef
```

Predict the optical density at H202 concentration 30.

```{r}
b0.hat <- lm.out$coef[1,1]
b0.hat
b1.hat <- lm.out$coef[2,1]
b1.hat
y.hat.30 <- b0.hat+b1.hat*30
y.hat.30
```

Get the 95% confidence interval for the mean optical density at H202 concentration 30.

```{r}
n <- length(dat$conc)
n
m <- mean(dat$conc)
m
SXX <- sum((dat$conc-m)^2)
SXX
s.hat <- lm.out$sigma
s.hat
y.hat.30 + c(-1,1)*qt(0.975,n-2)*s.hat*sqrt(1/n+(30-m)^2/SXX)
```

Easier:

```{r}
predict(lm.fit, data.frame(conc=30), interval="confidence")
```

Get the 95% prediction interval for a new observation at H202 concentration 30.

```{r}
y.hat.30 + c(-1,1)*qt(0.975,n-2)*s.hat*sqrt(1+1/n+(30-m)^2/SXX)
```

Easier:

```{r}
predict(lm.fit, data.frame(conc=30), interval="prediction")
```

#### Confidence and prediction intervals: a simple example

Generate some random data, 20 design points equally spaced on the x-axis, and assume E[y|x] = 10 + 5x. Assume Gaussian noise with mean zero and standard deviation 10.

```{r}
set.seed(1)
x <- 1:20
x
eps <- rnorm(20, sd=10)
eps
y <- 10 + 5*x + eps
y
par(las=1)
plot(x,y)
```

Fit a linear regression model.

```{r}
xydat <- data.frame(x,y)
head(xydat)
lm.fit <- lm(y ~ x, data=xydat)
summary(lm.fit)
confint(lm.fit)
```

The fitted values.

```{r}
lm.fit$fitted
fitted(lm.fit)
predict(lm.fit)
```

Make a picture with the regression line and the fitted values.

```{r}
par(las=1)
plot(xydat$x, xydat$y, xlab="x", ylab="y")
abline(lsfit(xydat$x, xydat$y))
points(xydat$x, predict(lm.fit), pch=21, bg="grey")
```

You can use the 'predict' function to get predicted values and confidence limits for the mean response for any value of x.

```{r}
predict(lm.fit, data.frame(x=c(5,10,20)), interval="confidence")
predict(lm.fit, data.frame(x=c(5,10,20)), interval="confidence", level=0.99)
```

Plot the confidence band for the mean response.

```{r}
xx <- seq(1, 20, by=0.1)
predict.mean <- predict(lm.fit, data.frame(x=xx), interval="confidence")
par(las=1)
plot(xydat$x, xydat$y, xlab="x", ylab="y", ylim=range(predict.mean), col="lightgrey")
lines(xx, predict.mean[,1], lwd=2)
lines(xx, predict.mean[,2], lty=2)
lines(xx, predict.mean[,3], lty=2)
```

You can also use the 'predict' function to get prediction intervals for new observations.

```{r}
predict(lm.fit, data.frame(x=c(5,10,20)), interval="prediction")
```

Plot the confidence band for the mean response, and the prediction interval for new observations.

```{r}
predict.new <- predict(lm.fit, data.frame(x=xx), interval="prediction")
par(las=1)
plot(xydat$x, xydat$y, xlab="x", ylab="y", ylim=range(predict.new), col="lightgrey")
lines(xx, predict.mean[,1], lwd=2)
lines(xx, predict.mean[,2], lty=2)
lines(xx, predict.mean[,3], lty=2)
lines(xx, predict.new[,2], lty=2, col="blue")
lines(xx, predict.new[,3], lty=2, col="blue")
```