# chex97.s                         J F Monahan, June 2001
 # IRWLS for Cox logistic regression example
 # Model has y as binomial( m, p ) with
 #  logit(p) = b[1] + b[2]*x
 #
 x <- c( 7, 14, 27, 51)
 X <- matrix( c(rep(1,4),x),4,2)
 m <- c( 55, 157, 159, 16)
 y <- c( 0, 2, 7, 3)
 # initial values
 b <- c( log(12/365), 0 )
 for (j in 1:10) {          # ten iterations
 cat("***************** iteration",j,"\n")
 print('beta')
 print(b)
 gam <- X %*% b
 p <- exp(gam)/(1+exp(gam))
 omp <- 1/(1+exp(gam))
 rll <- sum(y*gam) + sum(m*log(omp))
 print("log-like")
 print(rll)
 del <- t(X) %*% (y-m*p)
 print("gradient")
 print(del)
 w <- m*p*omp
 h <- - t(X) %*% diag(w[,1],nrow=length(m)) %*% X
 print("Hessian")
 print(h)
 s <- solve(h,del)
 print("Newton/Scoring Step")
 print(s) 
 b <- b - s                 # update
                 }          # end iteration loop on j
 print("Covariance Matrix")
 solve(-h)
 rm(list=ls())