#
# 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(shocks)
attach(shocks)

fraction = Y / m 

postscript("../../WriteUp/Graphics/Chapter12/prob_3_scatterplot.eps", onefile=FALSE, horizontal=FALSE)
plot( Intensity, fraction, xlab="Intensity", ylab="fraction" ) 
dev.off()

logisticRegressionTarget = as.matrix( cbind( Y, m - Y ) )
shocks$lRT = logisticRegressionTarget

m0 = glm( lRT ~ Intensity, family=binomial(), data=shocks )
s0 = summary(m0)
s0

# replot with the logistic fit on top:
# 
postscript("../../WriteUp/Graphics/Chapter12/prob_3_scatterplot.eps", onefile=FALSE, horizontal=FALSE)
plot( Intensity, fraction, xlab="Intensity", ylab="fraction", lty=1 ) 
xx <- seq(min(Intensity),max(Intensity),length=100)
lines(xx,predict(m0,data.frame(Intensity=xx),type="response"),lty=2)
legend(2.5,.2,legend=c("logistic regression fit"),lty=c(2))
dev.off()