library(MASS)

par(mfrow=c(1,2))
p <- seq(from=0.01,to=0.99,by=0.005)
logit <- log(p/(1-p))
plot(p,logit,xlab="p",ylab="logit(p)",type="l")
abline(h=0)
abline(v=0.5)
plot(logit,exp(logit)/(1+exp(logit)),ylab="p",type="l",xlab="logit(p)")
abline(h=0.5)
abline(v=0)
Dose <- c(1,2,4,8,16,32,1,2,4,8,16,32)
sex <- c(rep("Male",6),rep("Female",6))
sex <- as.factor(sex)
dead <- c(1,4,9,13,18,20,0,2,6,10,12,16)

Data <- data.frame(cbind(Dose,sex,dead))
Data
par(mfrow=c(1,1))
plot(c(1,32),c(0,1),type="n",xlab="Dose",ylab="Probability of death", log="x")
points(Dose[sex=="Male"],dead[sex=="Male"]/20,pch=16,col=1)
points(Dose[sex=="Female"],dead[sex=="Female"]/20,pch=16,col=2)
legend(x="bottomright",legend=rev(levels(sex)),pch=16,col=1:2)

ldose <- log2(Dose)
SF <- cbind(dead,20-dead)

contrasts(sex) <- contr.treatment(n=2,base=1)


glm.obj <- glm(SF~ldose*sex,family=binomial)
summary(glm.obj)
x <- seq(0,5,0.01)
n <- length(x)
a <- glm.obj$coef
eta <- a[1] + a[3] + (a[2] + a[4])*x
p <- exp(eta)/(1 + exp(eta))
lines(2^x,p,col=1)
eta <- a[1] + a[2]*x
p <- exp(eta)/(1 + exp(eta))
lines(2^x,p,col=2)
legend(2,1,legend=rev(levels(sex)),pch=16,col=1:2,lty=1)

ldose3 <- ldose-3
glm.obj <- glm(SF~ldose3*sex,family=binomial)
summary(glm.obj)

glm.obj <- glm(SF~sex+ldose-1,family=binomial)
summary(glm.obj)

ld<- dose.p(glm.obj,cf=c(1,3),p=1:3/4) # females
ld
# approx 95% CI for LD50 on original dose scale
2^( 3.263587-2*0.2297539)
2^3.263587
2^( 3.263587+2*0.2297539)


ld50 <- dose.p(glm.obj,cf=c(2,3),p=1:3/4) # males
ld50
2^c( 2.229262,0.2259649)


glm.obj <- glm(SF~sex*ldose,family=binomial)
summary(glm.obj)
anova(glm.obj)
1-pchisq(1.763,1)

glm.obj <- glm(SF~sex+ldose,family=binomial)
summary(glm.obj)
anova(glm.obj)
1-pchisq(10.227,1)