#make a data frame
Temp <- c(66, 70, 69, 68, 67, 72, 73, 70, 57, 63, 70, 78, 67, 53, 67, 75, 70, 81 , 76, 79, 75, 76, 58)
Damage  <- c(0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1)
challenger <- data.frame(Temp=Temp, Damage=Damage)

attach(challenger)
plot(Temp, Damage, xlab="Temperature", ylab="Damage", main="Incidence of Booster Field Joint Damage  vs. Temperature")

#fit a SLR
c.fit1 <- lm(Damage~Temp, challenger)
summary(c.fit1)

plot(Temp, Damage, xlab="Temperature", ylab="Damage", main="SLR Fitted Line")
abline(c.fit1)

resid1 <- resid(c.fit1)
plot(Temp, resid1, xlab="Temperature", ylab="Residuals", main="SLR Residuals")

#fit a logistic regression using default link (logit if binomial)
c.glm <- glm(Damage~Temp, data=challenger, family=binomial)

predict(c.glm)   # returns predictions on the Xb (linear predictor scale)
predict(c.glm,type='response')  # returns predictions on the Y scale
plot(predict(c.glm), residuals(c.glm,type='deviance'))
plot(residuals(c.glm, type='pearson'), residuals(c.glm,type='deviance'))

# to overlay data and curve, need type='response'
plot(Temp, Damage,pch=19)
lines(30:85, predict(c.glm,newdata=data.frame(Temp=30:85),type='response'))
title("Logistic Regression Predicted Values")



# Graphs of estimated probabilities for logistic model
# set a range of values for temperature
temp.range=seq(30,80,by=0.1)

# set the X of the model
X=cbind(rep(1,length(temp.range)),temp.range)
beta.log=coefficients(c.glm)
# compute linear predictor
Xb=X%*%beta.log

# use logistic cdf to compute probabilities
prob= exp(Xb)/(1+exp(Xb))

# plot probabilities versus temperature values
plot(temp.range, prob,type='l',xlab="temperature",ylab="probabilities",col="red")