######################################################################
# Computer notes             Statistics for Laboratory Scientists (2)
# Lecture 21                                 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".)
######################################################################
##############################
# Spider mites example
##############################
# Data available at:
#  http://www.biostat.jhsph.edu/~kbroman/teaching/labstat/fourth/spiders.csv
spiders <- read.csv("spiders.csv")

# Fit the logistic regression model
glm.out <- glm(n.dead/n ~ dose, data=spiders, weights=n,
               family=binomial(link=logit))
summary(glm.out)$coef

# Plot the data with the fitted curve
par(las=1)
plot(n.dead/n ~ dose, data=spiders,
     lwd=2, ylim=c(0,1), ylab="proportion dead")
u <- par("usr")
x <- seq(u[1],u[2],len=250)
y <- predict(glm.out, data.frame(dose=x), type="response")
lines(x,y,col="blue",lwd=2)

# LD50
co <- glm.out$coef
ld50 <- -co[1]/co[2]
glm.sum <- summary(glm.out)
se.co <- glm.sum$coef[,2]
cov.co <- glm.sum$cov.scaled[1,2]
se.ld50 <- abs(ld50) * sqrt( (se.co[1]/co[1])^2 + (se.co[2]/co[2])^2 -
                             2*cov.co/(co[1]*co[2]) )
ci.ld50 <- ld50 + c(-1,1) * qnorm(0.975) * se.ld50


##############################
# Worms example
##############################
# Data available at:
#  http://www.biostat.jhsph.edu/~kbroman/teaching/labstat/fourth/worms.csv
worms <- read.csv("worms.csv")

# plot the data
par(las=1)
plot(n.dead/n ~ dose, data=worms, col=ifelse(sex=="male","blue","red"),
     lwd=2, ylim=c(0,1), ylab="proportion dead")
u <- par("usr")
legend(u[2],u[3],c("Male","Female"),pch=1,col=c("blue","red"),xjust=1,yjust=0,
       cex=1.3)

# no sex difference
glm.out <- glm(n.dead/n ~ dose, weights=n, data=worms, 
               family=binomial(link=logit))
summary(glm.out)$coef

# sexes completely different
glm.outB <- glm(n.dead/n ~ sex*dose, weights=n, data=worms,
                family=binomial(link=logit))
summary(glm.outB)$coef

# different slopes but common "intercepts"
glm.outC <- glm(n.dead/n ~ dose + sex:dose, weights=n, data=worms,
                family=binomial(link=logit))
summary(glm.outC)$coef


# plot the data with fitted curves
par(las=1)
plot(n.dead/n ~ dose, data=worms, col=ifelse(sex=="male","blue","red"),
     lwd=2, ylim=c(0,1), ylab="proportion dead")
u <- par("usr")
x <- seq(u[1],u[2],len=250)
y <- predict(glm.out, data.frame(dose=x), type="response")
lines(x,y,col="green",lwd=2)
ym <- predict(glm.outB, data.frame(dose=x,sex=factor(rep("male",length(x)))),
              type="response")
lines(x,ym,col="blue",lwd=2, lty=2)
yf <- predict(glm.outB, data.frame(dose=x,sex=factor(rep("female",length(x)))),
              type="response")
lines(x,yf,col="red",lwd=2, lty=2)
ym <- predict(glm.outC, data.frame(dose=x,sex=factor(rep("male",length(x)))),
              type="response")
lines(x,ym,col="blue",lwd=2)
yf <- predict(glm.outC, data.frame(dose=x,sex=factor(rep("female",length(x)))),
              type="response")
lines(x,yf,col="red",lwd=2)
legend(u[2], u[3], c("Common intercept","Separate intercepts","Same curve"),lwd=2,
       col=c("blue","blue","green"),lty=c(1,2,1),xjust=1,yjust=0,cex=1.3)

# convert dose to log2(dose)
worms$dose <- log2(worms$dose)

# no sex difference
glm.out <- glm(n.dead/n ~ dose, weights=n, data=worms, 
               family=binomial(link=logit))
summary(glm.out)$coef

# sexes completely different
glm.outB <- glm(n.dead/n ~ sex*dose, weights=n, data=worms,
                family=binomial(link=logit))
summary(glm.outB)$coef

# different slopes but common "intercepts"
glm.outC <- glm(n.dead/n ~ dose + sex:dose, weights=n, data=worms,
                family=binomial(link=logit))
summary(glm.outC)$coef

# plot the data with fitted curves
par(las=1)
plot(n.dead/n ~ dose, data=worms, col=ifelse(sex=="male","blue","red"),
     lwd=2, ylim=c(0,1), ylab="proportion dead",
     xlab=expression(paste(log[2]," dose")))
u <- par("usr")
x <- seq(u[1],u[2],len=250)
y <- predict(glm.out, data.frame(dose=x), type="response")
lines(x,y,col="green",lwd=2)
ym <- predict(glm.outB, data.frame(dose=x,sex=factor(rep("male",length(x)))),
              type="response")
lines(x,ym,col="blue",lwd=2, lty=2)
yf <- predict(glm.outB, data.frame(dose=x,sex=factor(rep("female",length(x)))),
              type="response")
lines(x,yf,col="red",lwd=2, lty=2)
ym <- predict(glm.outC, data.frame(dose=x,sex=factor(rep("male",length(x)))),
              type="response")
lines(x,ym,col="blue",lwd=2)
yf <- predict(glm.outC, data.frame(dose=x,sex=factor(rep("female",length(x)))),
              type="response")
lines(x,yf,col="red",lwd=2)
legend(u[2], u[3], c("Common intercept","Separate intercepts","Same curve"),lwd=2,
       col=c("blue","blue","green"),lty=c(1,2,1), xjust=1,yjust=0,cex=1.3)

##################
# End of comp21.R
##################