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] > # chex99.r J F Monahan, June 2001, October 2007 > # > # IRWLS for Poisson regression using Frome smoking example > # > ni <- 9 # age groups > nj <- 7 # smoking rate group, nonsmoker = 1 > # read in data > All <- matrix(scan("frome.dat"),63,4,byrow=T) > I <- All[,1] > J <- All[,2] > # make factors > If <- factor(I) > Jf <- factor(J) > IJf <- factor( I+ni*(J-1) ) > Nij <- All[,3] > y <- All[,4] > y # print to check [1] 1 0 0 0 0 0 0 0 2 0 0 0 0 1 2 0 3 0 0 1 1 2 0 1 2 [26] 4 3 0 0 0 4 0 2 2 2 5 0 1 1 6 5 12 9 7 7 0 1 4 9 9 [51] 11 10 5 3 0 0 0 4 6 10 7 4 1 > # > # 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 <- bage[I] + bsmok[J] + lam <- Nij*exp(Xb) + rll <- sum( y*log(lam) - lam ) + print(k) + print(rll) + d <- y - lam + dell <- c( as.vector(tapply(d,If,sum)),as.vector(tapply(d,Jf,sum)) ) + print("gradient") # full (ni + nj) gradient + print(dell) + + H <- matrix(0,ni+nj,ni+nj) + H[1:ni,1:ni] <- diag( as.vector(tapply(lam,If,sum)),ni ) + H[ni+(1:nj),ni+(1:nj)] <- diag( as.vector(tapply(lam,Jf,sum)),nj ) + H[1:ni,ni+(1:nj)] <- matrix( as.vector(tapply(lam,IJf,sum)),ni,nj) + H[ni+(1:nj),1:ni] <- t( H[1:ni,ni+(1:nj)] ) + L <- t( chol(H[-(ni+1),-(ni+1)]) ) + step <- backsolve(t(L),forwardsolve(L,dell[-(ni+1)])) + print("step") + print(step) + beta <- beta + step + } # end of iteration loop on k [1] "betas" [1] -6.801519 -6.801519 -6.801519 -6.801519 -6.801519 -6.801519 -6.801519 [8] -6.801519 -6.801519 0.000000 0.000000 0.000000 0.000000 0.000000 [15] 0.000000 [1] 1 [1] -97.15273 [1] "gradient" [1] -33.2255309 -30.5849698 -22.0723640 0.7287428 1.6608402 24.6660997 [7] 21.4180393 19.7064750 17.7026677 -39.3415279 -11.3062682 -4.8620753 [13] -6.2430495 13.3752306 25.3555532 23.0221370 [1] "step" [1] -1.57341224 -1.65688558 -1.65312221 -0.91554304 -0.87595847 0.89396456 [7] 1.58252046 2.97286329 4.80424147 -0.06292814 0.44498659 0.63602475 [13] 1.24989813 1.69686771 3.19472699 [1] "betas" [1] -8.37493080 -8.45840414 -8.45464077 -7.71706160 -7.67747703 -5.90755401 [7] -5.21899810 -3.82865527 -1.99727710 -0.06292814 0.44498659 0.63602475 [13] 1.24989813 1.69686771 3.19472699 [1] 2 [1] -1383.044 [1] "gradient" [1] -17.40380 -17.54200 -14.61355 -16.35215 -14.23094 -104.20655 [7] -144.08898 -335.31152 -1086.90860 -156.15691 -97.04768 -124.59157 [13] -100.26508 -293.06677 -347.32037 -632.20973 [1] "step" [1] -0.99842883 -0.90162640 -0.75030264 -0.43351764 -0.44071204 -0.76947506 [7] -0.86524112 -0.97060882 -1.02379425 0.04808181 0.08058757 0.09332274 [13] 0.09559311 0.07067225 -0.01279945 [1] "betas" [1] -9.37335964 -9.36003054 -9.20494341 -8.15057924 -8.11818907 -6.67702907 [7] -6.08423922 -4.79926409 -3.02107135 -0.01484633 0.52557416 0.72934749 [13] 1.34549124 1.76753996 3.18192754 [1] 3 [1] -361.0926 [1] "gradient" [1] -6.155487 -5.754602 -4.156640 -2.876885 -2.549311 -30.415354 [7] -46.172342 -117.229349 -394.806647 -57.036661 -34.368371 -42.740379 [13] -33.707600 -97.799242 -117.194902 -227.269462 [1] "step" [1] -0.98877590 -0.77497041 -0.51309923 -0.19977859 -0.20192808 -0.53810668 [7] -0.70170564 -0.92512131 -1.05640968 0.11951961 0.19297801 0.21616393 [13] 0.21861684 0.16500238 -0.02202024 [1] "betas" [1] -10.3621355 -10.1350010 -9.7180426 -8.3503578 -8.3201171 -7.2151358 [7] -6.7859449 -5.7243854 -4.0774810 0.1046733 0.7185522 0.9455114 [13] 1.5641081 1.9325423 3.1599073 [1] 4 [1] -29.62573 [1] "gradient" [1] -2.0418396 -1.5607090 -0.7825962 -0.2828016 -0.2500570 [6] -6.1026953 -11.7776265 -37.7171175 -140.7255368 -20.8134084 [11] -11.5779305 -13.3769188 -10.1802501 -29.1769263 -35.9678623 [16] -80.1476831 [1] "step" [1] -0.9333161548 -0.5885469898 -0.3403509014 -0.2247660676 -0.2206723467 [6] -0.3473263005 -0.4943432510 -0.8221300065 -1.1000967466 0.2548152334 [11] 0.3779687639 0.3957738235 0.3910830439 0.3111303498 0.0005713531 [1] "betas" [1] -11.2954517 -10.7235479 -10.0583935 -8.5751239 -8.5407895 -7.5624621 [7] -7.2802881 -6.5465154 -5.1775778 0.3594885 1.0965209 1.3412852 [13] 1.9551911 2.2436727 3.1604787 [1] 5 [1] 70.45016 [1] "gradient" [1] -0.5746115 -0.2834900 -0.1316956 -0.3722847 -0.3140119 -0.9723892 [7] -1.9828243 -9.8761645 -47.7054099 -7.4902424 -3.5422510 -3.7137607 [13] -2.8504286 -8.1835712 -10.1259765 -26.3066513 [1] "step" [1] -0.7382572 -0.4494917 -0.3681063 -0.3524647 -0.3486711 -0.3595152 [7] -0.3998240 -0.6272443 -1.0451198 0.3873027 0.4953802 0.4803792 [13] 0.4665362 0.4094559 0.1173394 [1] "betas" [1] -12.0337089 -11.1730396 -10.4264998 -8.9275886 -8.8894606 -7.9219773 [7] -7.6801121 -7.1737597 -6.2226976 0.7467912 1.5919011 1.8216645 [13] 2.4217274 2.6531286 3.2778180 [1] 6 [1] 95.7756 [1] "gradient" [1] -0.10609207 -0.04709406 -0.07474224 -0.31141834 -0.26378365 [6] -0.48152237 -0.43247063 -1.60708798 -13.75437276 -2.38625664 [11] -0.90426859 -1.04891082 -0.89747669 -2.58969754 -2.75407153 [16] -6.49790232 [1] "step" [1] -0.4077980 -0.3261214 -0.3183460 -0.3160203 -0.3149133 -0.3156789 [7] -0.3185722 -0.3659322 -0.7023772 0.3362196 0.3682993 0.3514631 [13] 0.3443307 0.3296149 0.2095304 [1] "betas" [1] -12.441507 -11.499161 -10.744846 -9.243609 -9.204374 -8.237656 [7] -7.998684 -7.539692 -6.925075 1.083011 1.960200 2.173128 [13] 2.766058 2.982744 3.487348 [1] 7 [1] 99.49502 [1] "gradient" [1] -0.008373134 -0.008826940 -0.016291261 -0.064354329 -0.054074893 [6] -0.095171049 -0.076672769 -0.118470179 -2.474077739 -0.507757720 [11] -0.152287090 -0.230018988 -0.195342968 -0.546125461 -0.501576760 [16] -0.783203307 [1] "step" [1] -0.1257074 -0.1211118 -0.1207870 -0.1206744 -0.1205843 -0.1207180 [7] -0.1208813 -0.1226766 -0.2228144 0.1255589 0.1274447 0.1244579 [13] 0.1234820 0.1221752 0.1075272 [1] "betas" [1] -12.567214 -11.620273 -10.865633 -9.364283 -9.324958 -8.358374 [7] -8.119566 -7.662369 -7.147889 1.208570 2.087645 2.297585 [13] 2.889540 3.104919 3.594876 [1] 8 [1] 99.66071 [1] "gradient" [1] -0.0003273174 -0.0006297096 -0.0010163365 -0.0035679725 -0.0029241751 [6] -0.0049134473 -0.0040062506 -0.0048527891 -0.1313213264 -0.0377749929 [11] -0.0080983884 -0.0140263734 -0.0113962349 -0.0311603738 -0.0270363436 [16] -0.0240666175 [1] "step" [1] -0.01115782 -0.01115464 -0.01115855 -0.01116341 -0.01116280 -0.01116902 [7] -0.01117180 -0.01116746 -0.01709886 0.01139497 0.01137179 0.01127329 [13] 0.01124462 0.01121522 0.01096707 [1] "betas" [1] -12.578372 -11.631427 -10.876791 -9.375447 -9.336121 -8.369543 [7] -8.130738 -7.673536 -7.164988 1.219965 2.099017 2.308859 [13] 2.900785 3.116134 3.605843 [1] 9 [1] 99.66131 [1] "gradient" [1] -2.247634e-06 -4.567245e-06 -7.118566e-06 -2.412055e-05 -1.960370e-05 [6] -3.236381e-05 -2.653810e-05 -3.295757e-05 -4.466812e-04 -2.427635e-04 [11] -2.602936e-05 -4.685907e-05 -3.718429e-05 -1.007712e-04 -8.613333e-05 [16] -5.645754e-05 [1] "step" [1] -7.659084e-05 -7.673296e-05 -7.679047e-05 -7.683239e-05 -7.683700e-05 [6] -7.686055e-05 -7.686166e-05 -7.678560e-05 -9.603630e-05 7.757069e-05 [11] 7.736765e-05 7.714036e-05 7.707256e-05 7.700775e-05 7.684807e-05 [1] "betas" [1] -12.578449 -11.631504 -10.876868 -9.375524 -9.336198 -8.369620 [7] -8.130814 -7.673613 -7.165084 1.220042 2.099094 2.308936 [13] 2.900862 3.116211 3.605920 [1] 10 [1] 99.66131 [1] "gradient" [1] -1.048975e-10 -2.136322e-10 -3.327318e-10 -1.125863e-09 -9.148153e-10 [6] -1.508200e-09 -1.236406e-09 -1.537950e-09 -6.591768e-09 -9.914056e-09 [11] -2.717888e-10 -4.951715e-10 -3.895519e-10 -1.053192e-09 -8.975292e-10 [16] -5.449743e-10 [1] "step" [1] -3.251675e-09 -3.255407e-09 -3.256951e-09 -3.258000e-09 -3.258136e-09 [6] -3.258626e-09 -3.258578e-09 -3.256645e-09 -3.478315e-09 3.270827e-09 [11] 3.267057e-09 3.263900e-09 3.262812e-09 3.261715e-09 3.261365e-09 [1] "betas" [1] -12.578449 -11.631504 -10.876868 -9.375524 -9.336198 -8.369620 [7] -8.130814 -7.673613 -7.165084 1.220042 2.099094 2.308936 [13] 2.900862 3.116211 3.605920 [1] 11 [1] 99.66131 [1] "gradient" [1] -3.469447e-16 -4.163336e-15 -3.996803e-15 -2.220446e-15 -2.042810e-14 [6] 1.687539e-14 1.376677e-14 -6.217249e-15 -6.217249e-15 2.220446e-16 [11] 3.108624e-15 8.881784e-15 -3.996803e-15 -1.731948e-14 -1.953993e-14 [16] 1.509903e-14 [1] "step" [1] 1.148320e-16 -9.196592e-16 -2.069188e-16 3.602576e-16 -5.310008e-16 [6] 8.760025e-16 8.810182e-16 1.500027e-16 5.984777e-17 2.151501e-16 [11] 3.121329e-16 -5.892703e-16 -6.927980e-16 -7.373619e-16 9.788017e-17 [1] "betas" [1] -12.578449 -11.631504 -10.876868 -9.375524 -9.336198 -8.369620 [7] -8.130814 -7.673613 -7.165084 1.220042 2.099094 2.308936 [13] 2.900862 3.116211 3.605920 [1] 12 [1] 99.66131 [1] "gradient" [1] 2.636780e-16 4.496403e-15 1.332268e-15 1.842970e-14 -3.552714e-15 [6] 4.973799e-14 3.597123e-14 1.199041e-14 8.437695e-15 4.440892e-16 [11] 3.330669e-15 -3.996803e-15 1.110223e-14 7.105427e-14 2.930989e-14 [16] 1.598721e-14 [1] "step" [1] -3.829008e-16 8.502835e-16 -4.036084e-16 1.347966e-16 -7.587590e-16 [6] 7.337613e-16 6.198830e-16 -8.989527e-17 -1.353799e-16 4.352938e-16 [11] -4.278556e-16 5.809407e-16 1.319176e-15 3.820958e-16 3.128121e-16 > rm(list=ls()) >