# March 30, 2004   data here:  claim the colder one was 3 (not 2) failures 
# http://wps.aw.com/wps/media/objects/15/15719/projects/ch5_challenger/

# space shuttle challenger data

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)
ring = c(0,1,0,0,0,0,0,0,1,1,1,0,0,2,0,0,0,0,0,0,2,0,1)

eg.1 = function() {

   k = seq(23)[ring>0];  par(err=-1)
   x=jitter(temp[k]); y=jitter(ring[k]); xr <- range(temp[k],x); yr <- range(ring[k],y)
   plot(temp[k],ring[k],pch=16,main="Problem Flights 2 10 11 12 15 22 24",xlim=xr,ylim=yr)
          text(65,1.6,"Jan 20, 1986  Shuttle Flight 25")
   plot(x,y,pch=16,main="Jittered",xlim=xr,ylim=yr)
   plot(x,y,xlim=range(30,x),pch=16,main="Temperature is 31?"); text(31,1.5,"??",cex=2)
   plot(x,y,xlim=range(30,x),pch=16,main="Linear Prediction"); text(31,1.5,"??",cex=2)
     abline(lsfit(x,y),col=2,lwd=2)

   ans = lsfit(cbind(x,x^2),y); xx = seq(30,75,l=101)
   plot(x,y,xlim=range(30,x),pch=16,main="Quadratic Prediction"); text(31,1.5,"??",cex=2)
     lines(xx,cbind(1,xx,xx^2)%*%ans$coef,col=3,lwd=2)

   plot(x,y,xlim=range(30,x),ylim=c(-.3,10),pch=16,main="Quad again")
     text(31,1.5,"??",cex=2); lines(xx,cbind(1,xx,xx^2)%*%ans$coef,col=3,lwd=2)

   plot(temp,ring,pch=16,main="Include Shuttle Flights With No Failures")
   x=jitter(temp); y=jitter(ring);
   plot(x,y,pch=16); abline(lsfit(x,y),col=2,lwd=2)
   plot(x,y,pch=16,xlim=range(30,x)); abline(lsfit(x,y),col=2,lwd=2)

   y1 = ifelse(y>1,1,y)
   plot(x,y1,pch=16,xlim=range(30,x),main="yes/no event version")
     abline(lsfit(x,y1),col=2,lwd=2)

#       GLM
  io=order(temp); y2 = ifelse(ring>1,1,ring)[io]; x2=temp[io]
  glm.orings = glm( y2 ~ x2, family=binomial)
  lines( 31:80, predict(glm.orings,data.frame(x2=(31:80)),type="response"),lwd=3,col=3)
#  usa()
}


par(ask=T); eg.1()