# chex99.s                         J F Monahan, June 2001
 # IRWLS for Poisson regression using Frome smoking example
 ni <- 9              # age groups
 nj <- 7              # smoking rate group, nonsmoker = 1
 All <- scan("frome.dat")
 All <- t(matrix(All,4,63))
 I <- matrix(All[,1],ni,nj)
 J <- matrix(All[,2],ni,nj)
 Nij <- matrix(All[,3],ni,nj)
 y <- matrix(All[,4],ni,nj)
 # beta = ( age effects (1:ni), smoke effects (2:nj) )
 #                 so nonsmoke effect is zero
 # initial beta 
 bage <-  rep(log(sum(y)/sum(Nij)),ni)
 bsmok <- rep( 0, nj)
 beta <- c( bage, bsmok[2:nj] )
 for( k in 1:12 ) {
 print("betas")
 print(beta)
 bage <- beta[1:ni]
 bsmok <- c( 0, beta[(ni-1)+(2:nj)] )
 Xb <- matrix( rep(bage,nj),ni,nj) + 
      + t(matrix( rep(bsmok,ni),nj,ni) )
 lam <- exp(Xb)
 rll <- sum( y*log(Nij*lam) - Nij*lam )
 print(k); print(rll)
 d <- y - Nij*lam
 dell <- c( (d%*%rep(1,nj)), (rep(1,ni) %*% d)[2:nj] )
 print("gradient")
 print(dell)
 lam <- Nij*lam
 H <- matrix(0,ni+nj-1,ni+nj-1)
 H[1:ni,1:ni] <- diag( (lam%*%rep(1,nj))[,1],ni )
 H[(ni-1)+(2:nj),(ni-1)+(2:nj)] <- diag( (rep(1,ni)%*%lam[,2:nj])[1,], nj-1 )
 H[1:ni,(ni-1)+(2:nj)] <- lam[,2:nj]
 H[(ni-1)+(2:nj),1:ni] <- t(lam[,2:nj])
 step <- solve(H,dell)
 print("step")
 print(step)
 beta <- beta + step 
                  }           # end of iteration loop on k
 rm(list=ls())