>  # 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
***************** iteration 1 
beta -3.442019 0 
log-like -53.49422 
gradient -1.776357e-15 131.4884 
[1] "-Hessian"
          [,1]      [,2]
[1,]  11.62791  231.1159
[2,] 231.11592 5738.6575
Newton/Scoring Step -2.282487 0.1148365 
***************** iteration 2 
beta -5.724507 0.1148365 
log-like -52.45622 
gradient -10.19443 -393.2244 
[1] "-Hessian"
          [,1]       [,2]
[1,]  16.87900   511.2036
[2,] 511.20364 18170.9802
Newton/Scoring Step 0.3476281 -0.03142005 
***************** iteration 3 
beta -5.376879 0.08341648 
log-like -47.77531 
gradient -1.373675 -46.4419 
[1] "-Hessian"
          [,1]       [,2]
[1,]  12.08995   359.1453
[2,] 359.14530 12849.6330
Newton/Scoring Step -0.0368593 -0.002584047 
***************** iteration 4 
beta -5.413738 0.08083243 
log-like -47.68746 
gradient -0.06045993 -2.064531 
[1] "-Hessian"
          [,1]       [,2]
[1,]  11.05433   327.4891
[2,] 327.48911 11711.0813
Newton/Scoring Step -0.001438136 -0.0001360726 
***************** iteration 5 
beta -5.415176 0.08069636 
log-like -47.68728 
gradient -0.0001367550 -0.004815942 
[1] "-Hessian"
          [,1]       [,2]
[1,]  11.00619   325.9887
[2,] 325.98874 11656.1738
Newton/Scoring Step -1.094257e-06 -3.825635e-07 
***************** iteration 6 
beta -5.415177 0.08069598 
log-like -47.68728 
gradient -7.212755e-10 -2.670454e-08 
[1] "-Hessian"
          [,1]       [,2]
[1,]  11.00608   325.9853
[2,] 325.98528 11656.0441
Newton/Scoring Step 1.353568e-11 -2.6696e-12 
***************** iteration 7 
beta -5.415177 0.08069598 
log-like -47.68728 
gradient -5.551115e-16 -1.065814e-14 
[1] "-Hessian"
          [,1]       [,2]
[1,]  11.00608   325.9853
[2,] 325.98528 11656.0441
Newton/Scoring Step -1.360527e-16 2.890606e-18 
***************** iteration 8 
beta -5.415177 0.08069598 
log-like -47.68728 
gradient -5.551115e-16 -1.065814e-14 
[1] "-Hessian"
          [,1]       [,2]
[1,]  11.00608   325.9853
[2,] 325.98528 11656.0441
Newton/Scoring Step -1.360527e-16 2.890606e-18 
***************** iteration 9 
beta -5.415177 0.08069598 
log-like -47.68728 
gradient -5.551115e-16 -1.065814e-14 
[1] "-Hessian"
          [,1]       [,2]
[1,]  11.00608   325.9853
[2,] 325.98528 11656.0441
Newton/Scoring Step -1.360527e-16 2.890606e-18 
***************** iteration 10 
beta -5.415177 0.08069598 
log-like -47.68728 
gradient -5.551115e-16 -1.065814e-14 
[1] "-Hessian"
          [,1]       [,2]
[1,]  11.00608   325.9853
[2,] 325.98528 11656.0441
Newton/Scoring Step -1.360527e-16 2.890606e-18 
>    print("Covariance Matrix")
[1] "Covariance Matrix"
>    backsolve( t(L), forwardsolve(L,diag(1,nrow(L))) )
            [,1]          [,2]
[1,]  0.52931658 -0.0148034282
[2,] -0.01480343  0.0004998008
>  #  done
>   rm(list=ls())
>   q()