sink("c:/friday/logistic/1_1splusout.txt")
# All this could be done using the GUI.

# Bring in data; this can be done in various ways, and the
# variable names can be assigned in various ways.
# The following command is one of the logically simplest.

icu <- read.table("c:/friday/logistic/icu.dat",
       col.names=c("ID", "STA", "AGE", "SEX", "RACE", "SER",
       "CAN", "CRN", "INF", "CPR", "SYS", "HRA", "PRE", "TYP", 
       "FRA", "PO2", "PH", "PCO", "BIC", "CRE", "LOC"))

glmout <- glm(formula = STA ~ AGE, family = binomial(link = logit), 
          data = icu, na.action = na.exclude, 
          control = list(epsilon = 0.0001, maxit = 50, trace = F))

summary(glmout)

# get stuff for AGE=60

logest60 <- glmout$coefficients[1] + 60*glmout$coefficients[2]
pred60 <- exp(logest60) / (1 + exp(logest60))
# ??????????????????????????????????????????
# logestse60 <- sqrt(Var(bhat0) +
#                    x^2*Var(bhat1) + 
#                    2*x*Cov(bhat0,bhat1))
# ??????????????????????????????????????????
logestse60 <- sqrt(0.69514564^2 +
                   60^2*0.01055256^2 -
                   2*60*0.9658631*(0.69514564*0.01055256))

lpred60 <- exp(logest60 - qnorm(.975)*logestse60) / 
           (1 + exp(logest60 - qnorm(.975)*logestse60))
upred60 <- exp(logest60 + qnorm(.975)*logestse60) / 
           (1 + exp(logest60 + qnorm(.975)*logestse60))
pred60
lpred60
upred60
logest60
logestse60
out60<-c(pred60,lpred60,upred60,logest60,logestse60)
names(out60)<-c("pred60","lpred60","upred60","logest60","logestse60")
print (out60)

attach(icu)

logest <- glmout$coefficients[1] + AGE*glmout$coefficients[2]
# logestse <- sqrt(Var(bhat0)+x^2*Var(bhat1) + 2*x*Cov(bhat0,bhat1))
logestse <- sqrt(0.69514564^2 +
                   AGE^2*0.01055256^2 -
                   2*AGE*0.9658631*(0.69514564*0.01055256))

pred <- glmout$fitted.values
lpred <- exp(logest - qnorm(.975)*logestse) / 
           (1 + exp(logest - qnorm(.975)*logestse))
upred <- exp(logest + qnorm(.975)*logestse) / 
           (1 + exp(logest + qnorm(.975)*logestse))

ord<-order(AGE)
plot(AGE[ord],pred[ord],xlab='AGE',ylab='Pr (STA | AGE)',type='l',
     ylim=c(min(lpred),max(upred)))
lines(AGE[ord], lpred[ord])
lines(AGE[ord], upred[ord])