# chex97.r                 J F Monahan, June 2001, October 2007, 2009
 #
 # IRWLS for Cox logistic regression example
 # Model has y as binomial( m, p ) with
 #  logit(p) = b[1] + b[2]*x
 #
 X <- matrix( c(rep(1,4), 7, 14, 27, 51),4,2)
 m <- c( 55, 157, 159, 16)    # trials
 y <- c( 0, 2, 7, 3)          # failures
 #
 # initial values
 beta <- c( log(sum(y)/sum(m-y)), 0 )
 #
 # start looping
 for (j in 1:10) {                      # ten iterations are enough
   cat("***************** iteration",j,"\n")
   cat("beta",beta,"\n")
   gam <- as.vector(X %*% beta)
   p <- exp(gam)/(1+exp(gam))           # probabilities
   omp <- 1/(1+exp(gam))                # 1 - Pr
   rll <- sum(y*gam) + sum(m*log(omp))  # unlikelihood
   cat("log-like",rll,"\n")
   del <- t(X) %*% (y-m*p)              # gradient function
   cat("gradient",del,"\n")
   w <- m*p*omp                         # weights
 #      drop negative from Hessian
   mh <-  t(X) %*% ( X * w )            #***goofy***t(X)%*%diag(w)%*%X
   print("-Hessian")
   print(mh)
   L <- t(chol(mh))
   step <- backsolve( t(L), forwardsolve(L,del) )
   cat("Newton/Scoring Step",step,"\n")
   beta <- beta + step                  # update
                             }          # end iteration loop on j
   print("Covariance Matrix")
   backsolve( t(L), forwardsolve(L,diag(1,nrow(L))) )
 #  done
  rm(list=ls())
  q()