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

#### Statistics for Laboratory Scientists ( 140.615 )

## Error Propagation

David Sullivan's pf3d7 plasmodium data.

```{r}
h2o2 <- c(0,0,0, 10,10,10, 25,25,25, 50,50,50)
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)
```

Plot the data.

```{r}
par(las=1)
plot(h2o2, pf3d7, xaxt="n", xlab="H2O2 concentration", ylab="OD")
axis(1,c(0,10,25,50))
abline(lsfit(h2o2,pf3d7), lwd=2, col="blue", lty=2)
```

What is the percent change per unit H202 concentration?

```{r}
lm.fit <- lm(pf3d7 ~ h2o2)
betahat <- lm.fit$coefficients
betahat
ratio <- betahat[2]/betahat[1]
ratio
round(100*ratio,1)
```

Plot the data as percentages of optical density at zero concentration.

```{r}
par(las=1)
plot(h2o2, pf3d7, xaxt="n", yaxt="n", xlab="H2O2 concentration", ylab="OD")
axis(1,c(0,10,25,50))
abline(lsfit(h2o2,pf3d7), lwd=2, col="blue", lty=2)
axis(2, betahat[1]*seq(0.5,1,by=0.1), paste0(seq(50,100,10),"%"))
abline(h=betahat[1]*seq(0.5,1,by=0.1), lty="dotted", col="lightgrey")
```

Side note: what H2O2 concentration leads to an estimated 50% reduction in optical density?

```{r}
-0.5/ratio
```

What is the standard error of the percent change per unit H202 concentration? Here is a function to do all the work.

```{r}
estratio <- function(x,y){  
  lm.sum <- summary(lm(y~x))
  # the coefficients
  betahat <- lm.sum$coef[,1]
  # the estimated ratio
  ratio <- betahat[2]/betahat[1]
  # the estimated standard errors
  ses <- lm.sum$coef[,2]
  # the covariance between estimated coefficients
  covar <- lm.sum$cov.unscaled[1,2]*lm.sum$sigma^2
  # the estimated standard error 
  se <- abs(ratio) * sqrt(sum((ses/betahat)^2+2*covar/prod(betahat)))
  # return the estimated ratio and its estimated standard error
  return(list(ratio=ratio,se=se))
}
```

Return the percent change per unit H202 concentration, and its standard error.

```{r}
est <- estratio(h2o2, pf3d7)
est
```

The estimated percent change per unit H202 concentration.

```{r}
round(100*est$ratio,2)
```

The estimated standard error.

```{r}
round(100*est$se,3)
```

The 95% confidence interval for the estimated percent change per unit H202 concentration.

```{r}
round(100 * (est$ratio + c(-1,1)*1.96*est$se) ,2)
```