R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. [Previously saved workspace restored] > # chex98.r J F Monahan, June 2001, October 2007 > # > # IRWLS for Finney logistic regression example > # Model has y as binomial( m, p ) with > nobs <- 39 > p <- 3 > # read data in goofy order > X <- scan("finney.dat") > X <- matrix( X, p, nobs) > y <- X[3,] > y [1] 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 0 1 0 0 0 0 1 0 1 0 1 0 1 0 0 1 1 1 0 0 [39] 1 > X <- matrix( c(rep(1,nobs),log(X[1,]),log(X[2,])),nobs, 3 ) > X [,1] [,2] [,3] [1,] 1 1.30833282 -0.19237189 [2,] 1 1.25276297 0.08617770 [3,] 1 0.22314355 0.91629073 [4,] 1 -0.28768207 0.40546511 [5,] 1 -0.22314355 1.16315081 [6,] 1 -0.35667494 1.25276297 [7,] 1 -0.51082562 -0.28768207 [8,] 1 0.09531018 0.53062825 [9,] 1 -0.10536052 -0.28768207 [10,] 1 -0.10536052 -0.79850770 [11,] 1 -0.22314355 -0.56211892 [12,] 1 -0.59783700 1.01160091 [13,] 1 -0.51082562 1.09861229 [14,] 1 0.33647224 0.84586827 [15,] 1 -0.28768207 1.32175584 [16,] 1 0.83290912 0.49469624 [17,] 1 1.16315081 0.47000363 [18,] 1 -0.16251893 0.34712953 [19,] 1 0.53062825 0.05826891 [20,] 1 0.58778666 0.58778666 [21,] 1 -0.91629073 0.69314718 [22,] 1 -0.05129329 0.30748470 [23,] 1 0.30010459 0.30010459 [24,] 1 0.40546511 0.30748470 [25,] 1 0.47000363 0.57661336 [26,] 1 -0.51082562 0.40546511 [27,] 1 0.58778666 0.40546511 [28,] 1 -0.05129329 0.64185389 [29,] 1 0.64185389 -0.05129329 [30,] 1 0.47000363 -0.91629073 [31,] 1 0.99325177 -0.28768207 [32,] 1 0.85441533 -3.50655790 [33,] 1 0.09531018 0.60431597 [34,] 1 0.09531018 0.78845736 [35,] 1 0.18232156 0.69314718 [36,] 1 -0.22314355 1.20297230 [37,] 1 -0.05129329 0.64185389 [38,] 1 -0.28768207 0.64185389 [39,] 1 0.26236426 0.48550782 > # initial values > beta <- c( log(sum(y)/(nobs-sum(y))), 0 , 0) > m <- rep(1,nobs) # all m = 1 > # > # > for (j in 1:10) { # ten iterations + cat("***************** iteration",j,"\n") + print('beta') + print(beta) + gam <- X %*% beta + 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 + # take negative from Hessian + mh <- t(X) %*% diag(w[,1],nrow=length(m)) %*% X + print("Hessian") + print(mh) + L <- t(chol(mh)) # cholesky factor + step <- backsolve( t(L), forwardsolve(L,del) ) + print("Newton/Scoring Step") + print(step) + beta <- beta + step # update + } # end iteration loop on j ***************** iteration 1 [1] "beta" [1] 0.05129329 0.00000000 0.00000000 [1] "log-like" [1] -27.01992 [1] "gradient" [,1] [1,] -1.554312e-15 [2,] 4.203881e+00 [3,] 5.155115e+00 [1] "Hessian" [,1] [,2] [,3] [1,] 9.743590 1.555429 3.096897 [2,] 1.555429 2.992714 -1.079171 [3,] 3.096897 -1.079171 7.501865 [1] "Newton/Scoring Step" [,1] [1,] -0.7960684 [2,] 2.3042844 [3,] 1.3472879 ***************** iteration 2 [1] "beta" [,1] [1,] -0.7447751 [2,] 2.3042844 [3,] 1.3472879 [1] "log-like" [1] -17.47327 [1] "gradient" [,1] [1,] -0.2710879 [2,] 0.7050777 [3,] 1.7856513 [1] "Hessian" [,1] [,2] [,3] [1,] 7.364347 0.7208890 3.5511332 [2,] 0.720889 1.5913588 -0.5273007 [3,] 3.551133 -0.5273007 4.1252068 [1] "Newton/Scoring Step" [,1] [1,] -0.7430715 [2,] 1.1852657 [3,] 1.2240328 ***************** iteration 3 [1] "beta" [,1] [1,] -1.487847 [2,] 3.489550 [3,] 2.571321 [1] "log-like" [1] -15.42521 [1] "gradient" [,1] [1,] -0.3655699 [2,] 0.1474687 [3,] 0.5986119 [1] "Hessian" [,1] [,2] [,3] [1,] 5.9503027 0.5320558 3.3471170 [2,] 0.5320558 1.1271756 -0.4333254 [3,] 3.3471170 -0.4333254 3.1454677 [1] "Newton/Scoring Step" [,1] [1,] -0.8305949 [2,] 0.9881674 [3,] 1.2102836 ***************** iteration 4 [1] "beta" [,1] [1,] -2.318441 [2,] 4.477717 [3,] 3.781604 [1] "log-like" [1] -14.71815 [1] "gradient" [,1] [1,] -0.14245209 [2,] 0.05679688 [3,] 0.13121700 [1] "Hessian" [,1] [,2] [,3] [1,] 5.1634950 0.5456021 2.9853562 [2,] 0.5456021 0.9642083 -0.3781120 [3,] 2.9853562 -0.3781120 2.6641084 [1] "Newton/Scoring Step" [,1] [1,] -0.4702846 [2,] 0.5834657 [3,] 0.6590571 ***************** iteration 5 [1] "beta" [,1] [1,] -2.788726 [2,] 5.061183 [3,] 4.440661 [1] "log-like" [1] -14.61608 [1] "gradient" [,1] [1,] -0.01613479 [2,] 0.01446453 [3,] 0.01412735 [1] "Hessian" [,1] [,2] [,3] [1,] 4.8174428 0.5557891 2.7726685 [2,] 0.5557891 0.8940811 -0.3507086 [3,] 2.7726685 -0.3507086 2.4479173 [1] "Newton/Scoring Step" [,1] [1,] -0.08451294 [2,] 0.11498850 [3,] 0.11797015 ***************** iteration 6 [1] "beta" [,1] [1,] -2.873239 [2,] 5.176172 [3,] 4.558632 [1] "log-like" [1] -14.61369 [1] "gradient" [,1] [1,] -0.0003080518 [2,] 0.0005088565 [3,] 0.0002798141 [1] "Hessian" [,1] [,2] [,3] [1,] 4.7558844 0.5552568 2.7341913 [2,] 0.5552568 0.8800155 -0.3459285 [3,] 2.7341913 -0.3459285 2.4110002 [1] "Newton/Scoring Step" [,1] [1,] -0.002181250 [2,] 0.003150191 [3,] 0.003041688 ***************** iteration 7 [1] "beta" [,1] [1,] -2.875420 [2,] 5.179322 [3,] 4.561673 [1] "log-like" [1] -14.61369 [1] "gradient" [,1] [1,] -1.833721e-07 [2,] 3.969222e-07 [3,] 1.529669e-07 [1] "Hessian" [,1] [,2] [,3] [1,] 4.7542293 0.5551989 2.7331848 [2,] 0.5551989 0.8796031 -0.3458038 [3,] 2.7331848 -0.3458038 2.4100488 [1] "Newton/Scoring Step" [,1] [1,] -1.456580e-06 [2,] 2.167250e-06 [3,] 2.026312e-06 ***************** iteration 8 [1] "beta" [,1] [1,] -2.875422 [2,] 5.179324 [3,] 4.561675 [1] "log-like" [1] -14.61369 [1] "gradient" [,1] [1,] -8.215650e-14 [2,] 1.950523e-13 [3,] 5.959122e-14 [1] "Hessian" [,1] [,2] [,3] [1,] 4.7542281 0.5551989 2.7331841 [2,] 0.5551989 0.8796028 -0.3458038 [3,] 2.7331841 -0.3458038 2.4100482 [1] "Newton/Scoring Step" [,1] [1,] -6.666329e-13 [2,] 1.006221e-12 [3,] 9.251171e-13 ***************** iteration 9 [1] "beta" [,1] [1,] -2.875422 [2,] 5.179324 [3,] 4.561675 [1] "log-like" [1] -14.61369 [1] "gradient" [,1] [1,] -2.220446e-16 [2,] 2.220446e-16 [3,] -7.216450e-16 [1] "Hessian" [,1] [,2] [,3] [1,] 4.7542281 0.5551989 2.7331841 [2,] 0.5551989 0.8796028 -0.3458038 [3,] 2.7331841 -0.3458038 2.4100482 [1] "Newton/Scoring Step" [,1] [1,] 8.043840e-16 [2,] -7.753718e-16 [3,] -1.322920e-15 ***************** iteration 10 [1] "beta" [,1] [1,] -2.875422 [2,] 5.179324 [3,] 4.561675 [1] "log-like" [1] -14.61369 [1] "gradient" [,1] [1,] -2.831069e-15 [2,] 1.110223e-16 [3,] -1.887379e-15 [1] "Hessian" [,1] [,2] [,3] [1,] 4.7542281 0.5551989 2.7331841 [2,] 0.5551989 0.8796028 -0.3458038 [3,] 2.7331841 -0.3458038 2.4100482 [1] "Newton/Scoring Step" [,1] [1,] -8.866352e-16 [2,] 8.195116e-16 [3,] 3.399718e-16 > print("Covariance Matrix") [1] "Covariance Matrix" > # invert -Hessian > backsolve( t(L), forwardsolve(L,diag(1,nrow(L))) ) [,1] [,2] [,3] [1,] 1.744495 -1.991213 -2.264102 [2,] -1.991213 3.477664 2.757182 [3,] -2.264102 2.757182 3.378211 > # done > rm(list=ls()) > q()