#
# Epage 266
#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

if(is.null(version$language) == FALSE){ 
  require(alr3)
}else{
  library(alr3)
}

if(is.null(version$language) == FALSE) data(donner)
donner = na.omit(donner)
attach(donner)

indM = donner$Sex == "Male" 
nM   = sum(indM)
indF = donner$Sex == "Female" 
nF   = sum(indF)

fMSurvived = sum(donner$Outcome[indM])/nM
fFSurvived = sum(donner$Outcome[indF])/nF

postscript("../../WriteUp/Graphics/Chapter12/prob_4_scatterplot.eps", onefile=FALSE, horizontal=FALSE)
plot( Age, Outcome, xlab="Age", ylab="Outcome" ) 
dev.off()

m0 = glm( Outcome ~ Age, family=binomial(), data=donner )
s0 = summary(m0)
s0


# replot the scatter plot with the logistic fit and a loess smooth on top:
# 
postscript("../../WriteUp/Graphics/Chapter12/prob_4_scatterplot.eps", onefile=FALSE, horizontal=FALSE)

plot( Age, Outcome, xlab="Age", ylab="Outcome", lty=1 ) 

xx <- seq(min(Age),max(Age),length=100)
lines(xx,predict(m0,data.frame(Age=xx),type="response"),lty=2)

lo.fit <- loess( Outcome~Age,data=donner,degree=1)
lines(xx,predict(lo.fit,data.frame(Age=xx)),lty=3)

legend(40,0.8,legend=c("logistic regression fit","loess smooth"),lty=c(2,3))

dev.off()


# create a new model with a quadratic value of Age
m1 = glm( Outcome ~ Age + I(Age^2), family=binomial(), data=donner ) 
s1 = summary(m1)
s1


# replot the scatter plot with the logistic fit (both models) and a loess smooth on top:
# 
postscript("../../WriteUp/Graphics/Chapter12/prob_4_scatterplot.eps", onefile=FALSE, horizontal=FALSE)

plot( Age, Outcome, xlab="Age", ylab="Outcome", lty=1 ) 

xx <- seq(min(Age),max(Age),length=100)
lines(xx,predict(m0,data.frame(Age=xx),type="response"),lty=2)

lo.fit <- loess( Outcome~Age,data=donner,degree=1)
lines(xx,predict(lo.fit,data.frame(Age=xx)),lty=3)

lines(xx,predict(m1,data.frame(Age=xx),lty=4))

legend(40,0.8,legend=c("logistic linear fit","loess smooth","logistic quadratic fit"),lty=c(2,3,4))

dev.off()

# fit the more complicated model suggested in the text:
#
SexF = factor(Sex,ordered=FALSE)
donner$SexF = SexF 
StatusF = factor(Status,ordered=FALSE)
donner$StatusF = StatusF

m2 = update(m1, ~.+SexF+StatusF )
s2 = summary(m2)
s2

anova(m0,m1,m2,test="Chisq")