######################################################################
# Computer notes             Statistics for Laboratory Scientists (2)
# Lecture 15                                 Johns Hopkins University 
######################################################################
# These notes provide R code related to the course lectures, to
# further illustrate the topics in the lectures and to assist the
# students to learn R.
#
# Lines beginning with the symbol '#' are comments in R.  All other
# lines contain code.
#
# In R for Windows, you may wish to open this file from the menu bar
# (File:Display file); you can then easily copy commands into the
# command window.  (Use the mouse to highlight one or more lines; then
# right-click and select "Paste to console".)
######################################################################

###############################################
# David Sullivan's pf3d7 plasmodium data
###############################################

# read the data
h2o2 <- c(0,0,0,10,10,10,25,25,25,50,50,50)
od <- 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)

# make it a data frame
pf3d7 <- data.frame(h2o2=h2o2, od=od)

# the regression of od (optical density) on h2o2 (hydrogen
#     peroxide concentration)
lm.out <- lm(od ~ h2o2, data=pf3d7)

# summary table: the key bit is the table with
#                the est'd coefficients, SEs, and P-values
lm.sum <- summary(lm.out)

# just the table of coefficients
lm.sum$coef

# the residual SD
lm.sum$sigma

# 95% confidence intervals for the true intercept and slope
df <- lm.out$df
tmult <- qt(0.975, df)  # multiplier from the t distribution
beta <- lm.sum$coef[,1] # estimated coefficients
se <- lm.sum$coef[,2]   # estimated SEs

# 95% CI for intercept
beta[1] + c(-1,1)*tmult*se[1]

# 95% CI for slope
beta[2] + c(-1,1)*tmult*se[2]



######################################################################
# estimated SE of ratio
######################################################################
# read the data
h2o2 <- c(0,0,0,10,10,10,25,25,25,50,50,50)
od <- 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)

# function to do the work
estratio <-
function(x,y)
{  
  lm.sum <- summary(lm(y~x))

  # coefficients
  betahat <- lm.sum$coef[,1]

  # estimated ratio
  ratio <- betahat[2]/betahat[1]

  # SEs
  ses <- lm.sum$coef[,2]

  # covariance between estimated coefficients
  covar <- lm.sum$cov.unscaled[1,2]*lm.sum$sigma^2

  # SE of ratio
  se <- abs(ratio) * sqrt( sum( (ses/betahat)^2 + 2*covar/prod(betahat)) )

  print(betahat)
  print(ses)
  print(covar)

  c(ratio=ratio, se=se)
}

estratio(h2o2, od)*100


##################
# End of comp15.R
##################