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. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("gllm-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('gllm') This is gllm 0.27 > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "boot.gllm" > > ### * boot.gllm > > flush(stderr()); flush(stdout()) > > ### Name: boot.gllm > ### Title: Bootstrap for generalized log-linear modelling > ### Aliases: boot.gllm > ### Keywords: category models > > ### ** Examples > > # > # Fit Hochberg 1977 double sampling data > # 2x2 table of imprecise measures and 2x2x2x2 reliability data > # > # 2x2 table of imprecise measures > # > y1 <-c(1196, 13562, + 7151, 58175) > a2<- 2-as.integer(gl(2,1,4)) > b2<- 2-as.integer(gl(2,2,4)) > set1<-data.frame(y1,a2,b2) > # > # 2x2x2x2 reliability data > # > y2<-c(17, 3, 10, 258, + 3, 4, 4, 25, + 16, 3, 25, 197, + 100, 13, 107, 1014) > > a <- 2-as.integer(gl(2,1,16)) > a2<- 2-as.integer(gl(2,2,16)) > b <- 2-as.integer(gl(2,4,16)) > b2<- 2-as.integer(gl(2,8,16)) > > set2<-data.frame(y2,a,a2,b,b2) > # > # Combined analysis > # > y<-c(y1,y2) > # > # Map observed table onto underlying 2x2x2x2x2 table > # > s <-c(1, 1, 2, 2, 1, 1, 2, 2, 3, 3, 4, 4, 3, 3, 4, 4, + 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20) > # > # Model combining the tables is A*A2*B*B2 + L (dummy study variable) > # > a <- 2-as.integer(gl(2,1,32)) > a2<- 2-as.integer(gl(2,2,32)) > b <- 2-as.integer(gl(2,4,32)) > b2<- 2-as.integer(gl(2,8,32)) > l <- 2-as.integer(gl(2,16,32)) > > X <- model.matrix( ~ a*a2*b*b2+l) > > # > # Table 1 using unreliable measure > # > res1<-glm(y1 ~ a2*b2, family=poisson(),data=set1) > print(summary(res1)) Call: glm(formula = y1 ~ a2 * b2, family = poisson(), data = set1) Deviance Residuals: [1] 0 0 0 0 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 10.971211 0.004146 2646.20 <2e-16 *** a2 -2.096204 0.012531 -167.28 <2e-16 *** b2 -1.456184 0.009535 -152.71 <2e-16 *** a2:b2 -0.332086 0.032663 -10.17 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 9.2078e+04 on 3 degrees of freedom Residual deviance: 6.8967e-13 on 0 degrees of freedom AIC: 51.8 Number of Fisher Scoring iterations: 2 > # > # Table 2 using reliable measure > # > res2a<-glm(y2 ~ a*b, family=poisson(),data=set2) > print(summary(res2a)) Call: glm(formula = y2 ~ a * b, family = poisson(), data = set2) Deviance Residuals: Min 1Q Median 3Q Max -22.056 -14.234 -1.043 5.856 35.058 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.57595 0.03077 181.198 <2e-16 *** a -1.59627 0.07497 -21.293 <2e-16 *** b -0.82885 0.05582 -14.848 <2e-16 *** a:b -0.31762 0.14998 -2.118 0.0342 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 4575.5 on 15 degrees of freedom Residual deviance: 3324.8 on 12 degrees of freedom AIC: 3412.6 Number of Fisher Scoring iterations: 6 > # > # Table 2 demonstrating complex relationship between gold standard and > # unreliable measure > # > res2b<-glm(y2 ~ a*a2*b*b2, family=poisson(),data=set2) > print(summary(res2b)) Call: glm(formula = y2 ~ a * a2 * b * b2, family = poisson(), data = set2) Deviance Residuals: [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 6.92166 0.03140 220.409 < 2e-16 *** a -2.24883 0.10165 -22.124 < 2e-16 *** a2 -4.35671 0.27912 -15.609 < 2e-16 *** b -1.63845 0.07786 -21.043 < 2e-16 *** b2 -3.70278 0.20245 -18.290 < 2e-16 *** a:a2 4.28905 0.31186 13.753 < 2e-16 *** a:b 0.18450 0.23539 0.784 0.43315 a2:b 0.17212 0.64523 0.267 0.78966 a:b2 0.41625 0.54803 0.760 0.44753 a2:b2 2.52413 0.60656 4.161 3.16e-05 *** b:b2 3.97254 0.22347 17.777 < 2e-16 *** a:a2:b -0.55075 0.73360 -0.751 0.45280 a:a2:b2 -2.74415 0.98518 -2.785 0.00535 ** a:b:b2 -1.60229 0.67029 -2.390 0.01683 * a2:b:b2 -2.79388 1.02153 -2.735 0.00624 ** a:a2:b:b2 3.99082 1.38113 2.890 0.00386 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 4.5755e+03 on 15 degrees of freedom Residual deviance: 2.1982e-14 on 0 degrees of freedom AIC: 111.74 Number of Fisher Scoring iterations: 3 > # > # Combined analysis > # > require(gllm) [1] TRUE > res12<-gllm(y,s,X) > print(summary.gllm(res12)) Call: scoregllm(y = y, s = s, X = X, m = as.array(emgllm(y, s, X, maxit = em.maxit, tol = tol)$full.table)) No. cells in observed table: 20 No. cells in complete table: 32 Mean observed cell size: 4094.15 Model Deviance (df): 6.49 (3) Estimate S.E. exp(Estimate) Lower 95% CL Upper 95% CL (Intercept) 6.8949706 0.02843553 987.29670802 9.337764e+02 1.043885e+03 a -2.2488293 0.10301182 0.10552268 8.623040e-02 1.291312e-01 a2 -4.1375855 0.24002110 0.01596134 9.971474e-03 2.554933e-02 b -1.6384545 0.07890692 0.19428008 1.664414e-01 2.267750e-01 b2 -3.6511972 0.18991911 0.02596003 1.789134e-02 3.766757e-02 l 3.7958451 0.02384012 44.51584213 4.248362e+01 4.664528e+01 a:a2 4.2890502 0.28691184 72.89719693 4.154192e+01 1.279190e+02 a:b 0.1845014 0.23855135 1.20261872 7.534752e-01 1.919495e+00 a2:b 0.1721182 0.58708171 1.18781827 3.758512e-01 3.753912e+00 a:b2 0.4162479 0.54173821 1.51626168 5.243693e-01 4.384409e+00 a2:b2 2.2752713 0.55518222 9.73055880 3.277609e+00 2.888806e+01 b:b2 3.9725382 0.22141218 53.11918782 3.441768e+01 8.198252e+01 a:a2:b -0.5507467 0.67465107 0.57651916 1.536523e-01 2.163158e+00 a:a2:b2 -2.7441508 0.97537067 0.06430288 9.505570e-03 4.349935e-01 a:b:b2 -1.6022945 0.66415152 0.20143380 5.480184e-02 7.404053e-01 a2:b:b2 -2.7938841 0.98672320 0.06118311 8.845366e-03 4.232016e-01 a:a2:b:b2 3.9908229 1.34891173 54.09938762 3.845701e+00 7.610430e+02 > # > # Bootstrap the collapsed table to get estimated OR for reliable measures > # > # a and b are binary vectors the length of the *full* table > # and define the variables for which the odds ratio is to be > # estimated, here the reliable measure of injury and seatbelt usage. > # > boot.hochberg <- function (y,s,X,nrep,a,b) { + z<-boot.gllm(y,s,X,R=nrep) + boot.tab<-cbind(apply(z[,a & b],1,sum), + apply(z[,!a & b],1,sum), + apply(z[,a & !b],1,sum), + apply(z[,!a & !b],1,sum)) + oddsr<-boot.tab[,1]*boot.tab[,4]/boot.tab[,2]/boot.tab[,3] + hochberg.tab<-data.frame( c("yes","yes","no","no"), + c("yes","no","yes","no"), + boot.tab[1,], + apply(boot.tab[2:(1+nrep),],2,sd)) + colnames(hochberg.tab)<-c("Precise Injury","Precise seatbelt usage", + "Estimated Count","Bootstrap S.E.") + print(hochberg.tab) + cat("\nEstimated OR=",oddsr[1],"\n") + cat(" Bias=",oddsr[1]-mean(oddsr[2:(1+nrep)]),"\n") + cat("Bootstrap SE=",sd(oddsr[2:(1+nrep)]),"\n\nQuantiles\n\n") + quantile(oddsr[2:(1+nrep)],c(0.025,0.50,0.975)) + } > boot.hochberg(y,s,X,nrep=20,a,b) It 0 : 1196 13562 7151 58175 17 3 10 258 3 4 4 25 16 3 25 197 100 13 107 1014 It 1 : 1185 13711 7171 57952 25 4 11 258 3 0.5 9 33 8 1 33 172 130 14 118 1045 It 2 : 1190 13680 7133 58128 10 2 6 256 1 2 3 21 15 5 21 171 116 20 84 1019 It 3 : 1150 13384 7248 58306 16 2 13 246 2 2 5 30 14 1 18 189 92 12 114 1039 It 4 : 1196 13490 7270 58043 18 7 12 227 4 4 1 26 13 5 21 214 111 16 124 1081 It 5 : 1204 13547 7076 58267 13 6 15 255 3 3 3 20 9 4 24 181 92 8 118 1035 It 6 : 1167 13685 7270 58023 15 2 5 245 1 4 3 24 13 0.5 24 197 93 12 93 1007 It 7 : 1219 13644 7076 58164 11 4 9 277 4 6 3 26 11 2 26 234 108 12 99 948 It 8 : 1163 13514 7245 58174 13 2 6 247 4 3 2 16 11 3 29 202 106 12 109 1022 It 9 : 1186 13574 7051 58319 21 2 15 254 3 2 2 24 13 4 24 213 86 12 106 972 It 10 : 1161 13622 7254 58058 14 0.5 10 270 5 2 5 31 12 0.5 24 210 90 10 115 990 It 11 : 1170 13583 7105 58203 17 3 8 263 7 5 7 39 12 3 32 212 96 9 96 1013 It 12 : 1189 13705 7183 58059 15 2 9 244 6 5 2 17 20 3 23 202 110 22 122 945 It 13 : 1159 13414 7186 58311 19 1 12 275 4 5 7 25 17 4 27 201 100 11 99 1006 It 14 : 1194 13588 7162 58094 20 2 6 305 2 4 6 24 13 5 21 192 88 9 108 1040 It 15 : 1206 13636 7166 58032 16 2 6 250 8 2 6 16 14 3 24 188 107 9 98 1094 It 16 : 1187 13719 6977 58240 23 5 11 246 2 5 1 27 16 3 24 196 109 14 98 980 It 17 : 1232 13547 7103 58275 17 1 8 244 1 4 3 30 14 3 20 175 104 19 110 973 It 18 : 1296 13445 7226 58182 16 3 18 236 1 3 6 18 14 1 27 189 105 11 88 998 It 19 : 1192 13623 7184 58101 18 4 12 267 3 4 2 29 17 1 29 190 103 11 107 986 It 20 : 1230 13457 7150 58215 21 2 7 256 1 3 3 30 13 3 21 233 102 17 104 1015 Precise Injury Precise seatbelt usage Estimated Count Bootstrap S.E. 1 yes yes 3227.388 298.0281 2 yes no 21071.032 873.8611 3 no yes 10581.906 542.3849 4 no no 47002.674 860.0072 Estimated OR= 0.6803368 Bias= 0.03399542 Bootstrap SE= 0.08177061 Quantiles 2.5% 50% 97.5% 0.5270613 0.6499786 0.8022396 > > > > cleanEx(); ..nameEx <- "boot.table" > > ### * boot.table > > flush(stderr()); flush(stdout()) > > ### Name: boot.table > ### Title: Produce one bootstrap replicate of a vector of counts > ### Aliases: boot.table > ### Keywords: category manip > > ### ** Examples > > boot.table(c(1,3,4,2)) [1] 1 3 3 3 > ## 0.5 2.0 5.0 3.0 > boot.table(c(1,3,4,2),c(1,2,1,2)) [1] 1 2 4 3 > ## 2 1 3 4 > > > > cleanEx(); ..nameEx <- "emgllm" > > ### * emgllm > > flush(stderr()); flush(stdout()) > > ### Name: emgllm > ### Title: Generalized log-linear modelling by EM and iterative > ### proportional fitting > ### Aliases: emgllm > ### Keywords: category models > > ### ** Examples > > # > # latent class analysis: two latent classes > # > # Data matrix 2x2x2x2x2 table of responses to five binary items > # > y<-c( 3, 6, 2, 11, 1, 1, 3, 4, + 1, 8, 0, 16, 0, 3, 2, 15, + 10, 29, 14, 81, 3, 28, 15, 80, + 16, 56, 21, 173, 11, 61, 28, 298) > # > # Scatter matrix: full table is 2x2x2x2x2x2 > # > s<- c(1:32,1:32) > # > # Design matrix: x is the latent variable (2 levels), > # a-e are the observed variables > # > i<-rep(1,64) > x<-as.integer(gl(2,32,64))-1 > a<-as.integer(gl(2,16,64))-1 > b<-as.integer(gl(2,8 ,64))-1 > c<-as.integer(gl(2,4 ,64))-1 > d<-as.integer(gl(2,2 ,64))-1 > e<-as.integer(gl(2,1 ,64))-1 > X<-cbind(i,x,a,b,c,d,e,x*cbind(a,b,c,d,e)) > colnames(X)<-c("Int","X","A","B","C","D","E","AX","BX","CX","DX","EX") > res<-emgllm(y,s,X, tol=0.01) > res $deviance [1] 27.23202 $observed.values [1] 3 6 2 11 1 1 3 4 1 8 0 16 0 3 2 15 10 29 14 [20] 81 3 28 15 80 16 56 21 173 11 61 28 298 $fitted.values [1] 1.0015893 4.3465254 1.7638292 8.2845186 0.4672803 2.3582133 [7] 1.0383673 6.3909342 1.2274303 5.8545843 2.5058860 14.1891177 [13] 0.7531835 4.7651603 2.3024399 17.6717690 7.7492596 35.5975116 [19] 14.9305569 79.1472895 4.2884664 25.2370732 11.8739035 86.1337653 [25] 10.5723133 57.4147498 26.1409004 175.6330729 8.8765666 65.8591376 [31] 33.4611571 282.4254320 $full.table [1] 0.96260048 3.98704955 1.57107362 6.50731896 0.38402739 [6] 1.59062484 0.62677644 2.59608090 1.07567370 4.45539394 [11] 1.75562199 7.27171035 0.42913771 1.77746983 0.70040162 [16] 2.90103320 7.05993152 29.24193091 11.52261231 47.72616168 [21] 2.81654449 11.66600541 4.59692139 19.04024950 7.88923632 [26] 32.67687552 12.87613220 53.33238250 3.14739385 13.03636918 [31] 5.13690523 21.27683924 0.03898885 0.35947582 0.19275558 [36] 1.77719961 0.08325287 0.76758850 0.41159088 3.79485328 [41] 0.15175660 1.39919040 0.75026400 6.91740731 0.32404577 [46] 2.98769045 1.60203827 14.77073579 0.68932812 6.35558070 [51] 3.40794460 31.42112778 1.47192193 13.57106777 7.27698212 [56] 67.09351581 2.68307699 24.73787425 13.26476823 122.30069039 [61] 5.72917273 52.82276844 28.32425185 261.14859276 > # > # Obtain standard errors for parameter estimates > # > summary(scoregllm(y,s,X,as.array(res$full.table))) Call: scoregllm(y = y, s = s, X = X, m = as.array(res$full.table)) No. cells in observed table: 32 No. cells in complete table: 64 Mean observed cell size: 31.25 Model Deviance (df): 22.73 (20) Estimate S.E. exp(Estimate) Lower 95% CL Upper 95% CL Int 0.47512255 0.3553893 1.60821127 0.801363559 3.22742836 X -4.50412533 0.9717774 0.01106326 0.001646985 0.07431506 A 1.71057269 0.2849818 5.53212876 3.164541901 9.67105181 B 0.07796292 0.2868302 1.08108257 0.616175049 1.89676541 C -0.88062908 0.4193757 0.41452206 0.182208181 0.94303417 D 0.41663179 0.2652291 1.51684390 0.901931075 2.55098807 E 1.21264502 0.2560096 3.36236644 2.035754035 5.55347448 AX 1.56636861 0.5116575 4.78922503 1.756845496 13.05560248 BX 1.34898782 0.3146584 3.85352310 2.079769979 7.14003973 CX 1.66505560 0.3978561 5.28596713 2.423608874 11.52886045 DX 1.28246674 0.3144935 3.60552266 1.946551610 6.67837092 EX 1.24354205 0.3712358 3.46787513 1.675178910 7.17902895 > > > > cleanEx(); ..nameEx <- "gllm" > > ### * gllm > > flush(stderr()); flush(stdout()) > > ### Name: gllm > ### Title: Generalized log-linear modelling > ### Aliases: gllm > ### Keywords: category models > > ### ** Examples > > # > # latent class analysis: two latent classes > # > # Data matrix 2x2x2x2x2 table of responses to five binary items > # (items 11-15 of sections 6-7 of the Law School Admission Test) > # > y<-c( 3, 6, 2, 11, 1, 1, 3, 4, + 1, 8, 0, 16, 0, 3, 2, 15, + 10, 29, 14, 81, 3, 28, 15, 80, + 16, 56, 21, 173, 11, 61, 28, 298) > # > # Scatter matrix: full table is 2x2x2x2x2x2 > # > s<- c(1:32,1:32) > # > # Design matrix: x is the latent variable (2 levels), > # a-e are the observed variables > # > x<-as.integer(gl(2,32,64))-1 > a<-as.integer(gl(2,16,64))-1 > b<-as.integer(gl(2,8 ,64))-1 > c<-as.integer(gl(2,4 ,64))-1 > d<-as.integer(gl(2,2 ,64))-1 > e<-as.integer(gl(2,1 ,64))-1 > > res1<-gllm(y,s,~x*(a+b+c+d+e),method="em",tol=0.01) > res1 $call scoregllm(y = y, s = s, X = X, m = as.array(emgllm(y, s, X, maxit = 10000, tol = tol)$full.table)) $formula ~x * (a + b + c + d + e) $iter [1] 10 $deviance [1] 22.73392 $df [1] 20 $coefficients (Intercept) x a b c d 0.47512255 -4.50412533 1.71057269 0.07796292 -0.88062908 0.41663179 e x:a x:b x:c x:d x:e 1.21264502 1.56636861 1.34898782 1.66505560 1.28246674 1.24354205 $se [1] 0.3553893 0.9717774 0.2849818 0.2868302 0.4193757 0.2652291 0.2560096 [8] 0.5116575 0.3146584 0.3978561 0.3144935 0.3712358 $V [,1] [,2] [,3] [,4] [,5] (Intercept) 0.12630158 0.084105878 -0.081750172 -0.065361761 -0.080277929 x 0.08410588 0.944351233 -0.078503760 -0.155657916 -0.278851486 a -0.08175017 -0.078503760 0.081214643 0.040123513 0.056338303 b -0.06536176 -0.155657916 0.040123513 0.082271537 0.061023950 c -0.08027793 -0.278851486 0.056338303 0.061023950 0.175875953 d -0.06264655 -0.131710089 0.036176458 0.040590500 0.056087816 e -0.06682881 -0.091799480 0.031756560 0.035804015 0.050321760 x:a 0.04205725 -0.277408034 -0.071484329 0.022759930 0.047555088 x:b 0.02652651 0.004759244 -0.006330577 -0.063376746 0.005295488 x:c 0.03085246 0.139384321 -0.017774786 -0.014778111 -0.139980101 x:d 0.02623051 -0.026001915 -0.002777393 -0.001193245 0.008783732 x:e 0.03588475 -0.106897225 0.003620931 0.005812664 0.017725593 [,6] [,7] [,8] [,9] [,10] (Intercept) -0.062646546 -0.066828811 0.042057253 0.026526511 0.03085246 x -0.131710089 -0.091799480 -0.277408034 0.004759244 0.13938432 a 0.036176458 0.031756560 -0.071484329 -0.006330577 -0.01777479 b 0.040590500 0.035804015 0.022759930 -0.063376746 -0.01477811 c 0.056087816 0.050321760 0.047555088 0.005295488 -0.13998010 d 0.070346457 0.032296335 0.019666229 -0.005160153 -0.01538891 e 0.032296335 0.065540896 0.016381655 -0.005331130 -0.01561746 x:a 0.019666229 0.016381655 0.261793418 -0.011462403 -0.04123333 x:b -0.005160153 -0.005331130 -0.011462403 0.099009927 -0.02430249 x:c -0.015388908 -0.015617464 -0.041233332 -0.024302493 0.15828945 x:d -0.055506004 -0.002185207 -0.007869360 -0.006617990 -0.02188062 x:e 0.004678263 -0.059252865 -0.003120034 -0.006513638 -0.02244201 [,11] [,12] (Intercept) 0.026230514 0.035884749 x -0.026001915 -0.106897225 a -0.002777393 0.003620931 b -0.001193245 0.005812664 c 0.008783732 0.017725593 d -0.055506004 0.004678263 e -0.002185207 -0.059252865 x:a -0.007869360 -0.003120034 x:b -0.006617990 -0.006513638 x:c -0.021880620 -0.022442007 x:d 0.098906187 -0.004964923 x:e -0.004964923 0.137815986 $observed.values [1] 3 6 2 11 1 1 3 4 1 8 0 16 0 3 2 15 10 29 14 [20] 81 3 28 15 80 16 56 21 173 11 61 28 298 $fitted.values [1] 1.6259847 5.6148086 2.5366968 9.3367469 0.7056316 2.6960863 [7] 1.2244110 5.8861018 1.8127196 6.7100896 3.0425644 13.5939519 [13] 0.8831119 4.3170074 1.9814144 14.0326731 9.3682345 35.4110259 [19] 16.0732141 75.4363788 4.7208888 24.4440853 11.2429592 84.6767459 [25] 11.5820526 55.2385542 25.3294544 174.2875230 8.2899863 63.5797328 [31] 29.5806758 294.7384889 $residuals [1] 1.07753990 0.16255826 -0.33697272 0.54432783 0.35043077 -1.03295419 [7] 1.60464458 -0.77741200 -0.60363642 0.49796106 -1.74429481 0.65257652 [13] -0.93974032 -0.63386506 0.01320349 0.25822783 0.20640823 -1.07735363 [19] -0.51712172 0.64057079 -0.79202872 0.71922445 1.12048382 -0.50823163 [25] 1.29815865 0.10245138 -0.86024124 -0.09752632 0.94122743 -0.32353061 [31] -0.29062887 0.18997674 $full.table [1] 1.60819354 5.40735622 2.43939581 8.20217329 0.66664834 [6] 2.24152438 1.01120861 3.40006663 1.73860132 5.84583660 [11] 2.63720545 8.86728429 0.72070659 2.42328871 1.09320712 [16] 3.67577669 8.89686341 29.91462689 13.49524842 45.37614018 [21] 3.68803820 12.40058227 5.59421780 18.80988058 9.61830656 [26] 32.34039219 14.58957280 49.05567349 3.98709976 13.40614060 [31] 6.04785072 20.33516637 0.01779112 0.20745243 0.09730102 [36] 1.13457365 0.03898323 0.45456196 0.21320234 2.48603514 [41] 0.07411832 0.86425299 0.40535895 4.72666760 0.16240529 [46] 1.89371868 0.88820729 10.35689642 0.47137110 5.49639903 [51] 2.57796564 30.06023866 1.03285065 12.04350300 5.64874144 [56] 65.86686532 1.96374605 22.89816205 10.73988165 125.23184952 [61] 4.30288656 50.17359219 23.53282507 274.40332254 attr(,"class") [1] "gllm" > # > # An example of model fitting: gametic association between two diallelic loci > # > # Data matrix > # > y<-c( 187,386,156, + 352,310,20, + 136,0 ,0) > # > # Scatter matrix > # > s<- c( 1, 2, 2, 3, + 4, 5, 5, 6, + 4, 5, 5, 6, + 7, 8, 8, 9) > # > # Design matrix > # > X<- matrix(c( 1,0,0,0,0,0,1, + 1,0,1,0,0,0,0, + 1,0,1,0,0,0,0, + 1,0,2,0,1,0,0, + 1,1,0,0,0,0,0, + 1,1,1,0,0,1,0, + 1,1,1,0,0,0,1, + 1,1,2,0,1,1,1, + 1,1,0,0,0,0,0, + 1,1,1,0,0,0,1, + 1,1,1,0,0,1,0, + 1,1,2,0,1,1,1, + 1,2,0,1,0,0,0, + 1,2,1,1,0,1,1, + 1,2,1,1,0,1,1, + 1,2,2,1,1,2,2), byrow=TRUE, ncol=7) > > colnames(X)<-c("Intercept", "A", "B", "P1", "P2", "Delta", "Epsilon") > res2<-gllm(y,s,X[,c(1:6)],method="hybrid",em.maxit=1,tol=0.00001) > res2 $call scoregllm(y = y, s = s, X = X, m = as.array(emgllm(y, s, X, maxit = em.maxit, tol = tol)$full.table)) $formula NULL $iter [1] 7 $deviance [1] 28.2759 $df [1] 3 $coefficients Intercept A B P1 P2 Delta 5.303144151 -0.147717735 -0.102441387 -0.161218170 0.009094938 -3.289880831 $se [1] 0.05686545 0.05852344 0.05855679 0.12289564 0.11562174 0.24066719 $V [,1] [,2] [,3] [,4] [,5] Intercept 0.003233679 -0.0022075544 -0.0022482743 0.0011369438 0.0012182714 A -0.002207554 0.0034249932 0.0006762788 -0.0044734274 0.0008501124 B -0.002248274 0.0006762788 0.0034288977 0.0008713895 -0.0044672913 P1 0.001136944 -0.0044734274 0.0008713895 0.0151033396 -0.0026201559 P2 0.001218271 0.0008501124 -0.0044672913 -0.0026201559 0.0133683870 Delta 0.002911724 -0.0032866559 -0.0029818410 0.0000859077 -0.0003355251 [,6] Intercept 0.0029117238 A -0.0032866559 B -0.0029818410 P1 0.0000859077 P2 -0.0003355251 Delta 0.0579206969 $observed.values [1] 187 386 156 352 310 20 136 0 0 $fitted.values [1] 200.9676871 362.7993566 165.2329563 346.7394381 324.6388029 10.6217590 [7] 127.2928748 8.5618404 0.1452848 $residuals [1] -0.9852839 1.2180546 -0.7182778 0.2825078 -0.8124663 2.8775519 0.7717431 [8] -2.9260623 -0.3811624 $full.table [1] 200.9676871 181.3996783 181.3996783 165.2329563 173.3697191 5.8305065 [7] 156.4888949 5.3108795 173.3697191 156.4888949 5.8305065 5.3108795 [13] 127.2928748 4.2809202 4.2809202 0.1452848 attr(,"class") [1] "gllm" > # > > > > > cleanEx(); ..nameEx <- "hildesheim" > > ### * hildesheim > > flush(stderr()); flush(stdout()) > > ### Name: hildesheim > ### Title: Invasive Cervical Cancer v exposure to Herpes Simplex Virus > ### Aliases: hildesheim > ### Keywords: datasets > > ### ** Examples > > data(hildesheim) > ftable(xtabs(Freq ~ case+HSV.inac+HSV.gold, hildesheim)) HSV.gold ? Neg Pos case HSV.inac case Neg 318 13 5 Pos 375 3 18 control Neg 701 33 16 Pos 535 11 16 > fisher.test(xtabs(Freq ~ case+HSV.inac, hildesheim)) Fisher's Exact Test for Count Data data: xtabs(Freq ~ case + HSV.inac, hildesheim) p-value = 1.178e-06 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.5278139 0.7658777 sample estimates: odds ratio 0.6359578 > fisher.test(xtabs(Freq ~ case+HSV.gold, hildesheim, subset=HSV.gold!="?")) Fisher's Exact Test for Count Data data: p-value = 0.1148 alternative hypothesis: two.sided > > # > # Combined analysis (ordered as incomplete then complete data) > # > y<-hildesheim$Freq[c(3,9,6,12,1,2,7,8,4,5,10,11)] > # > # Map observed table onto underlying 2x2x2x2 table > # > s <-c(1, 1, 2, 2, 3, 3, 4, 4, + 5, 6, 7, 8, 9, 10, 11, 12) > # > substudy <- 2-as.integer(gl(2,8,16)) > hsv.inac <- 2-as.integer(gl(2,4,16)) > hsv.gold <- 2-as.integer(gl(2,2,16)) > cancer <- 2-as.integer(gl(2,1,16)) > > require(gllm) [1] TRUE > res<-gllm(y,s, ~substudy+hsv.inac*hsv.gold*cancer) > print(summary.gllm(res)) Call: scoregllm(y = y, s = s, X = X, m = as.array(emgllm(y, s, X, maxit = em.maxit, tol = tol)$full.table)) No. cells in observed table: 12 No. cells in complete table: 16 Mean observed cell size: 170.3333 Model Deviance (df): 1.98 (3) Estimate S.E. exp(Estimate) Lower 95% CL (Intercept) 3.3470286 0.14501705 28.4181665 21.3872549 substudy 2.8198249 0.09598987 16.7739130 13.8971632 hsv.inac -0.7912001 0.24617997 0.4533005 0.2797909 hsv.gold -0.7330717 0.19011922 0.4804310 0.3309773 cancer -0.7239187 0.32827671 0.4848485 0.2547835 hsv.inac:hsv.gold -0.6649846 0.59645056 0.5142815 0.1597688 hsv.inac:cancer 1.0986121 0.48863346 2.9999994 1.1512986 hsv.gold:cancer -0.2315927 0.60946249 0.7932691 0.2402346 hsv.inac:hsv.gold:cancer 1.6486588 0.93219691 5.2000010 0.8365683 Upper 95% CL (Intercept) 37.7604415 substudy 20.2461577 hsv.inac 0.7344103 hsv.gold 0.6973708 cancer 0.9226583 hsv.inac:hsv.gold 1.6554263 hsv.inac:cancer 7.8172565 hsv.gold:cancer 2.6194222 hsv.inac:hsv.gold:cancer 32.3225386 > # > # Bootstrap the collapsed table to get estimated OR for reliable measures > # > # a and b are binary vectors the length of the *full* table > # and define the variables for which the odds ratio is to be > # estimated, here the reliable measure of HSV exposure and Ca Cx > # > boot.hildesheim <- function (y,s,X,nrep,a,b) { + z<-boot.gllm(y,s,X,R=nrep) + boot.tab<-cbind(apply(z[,a & b],1,sum), + apply(z[,!a & b],1,sum), + apply(z[,a & !b],1,sum), + apply(z[,!a & !b],1,sum)) + oddsr<-boot.tab[,1]*boot.tab[,4]/boot.tab[,2]/boot.tab[,3] + hildesheim.tab<-data.frame( c("yes","yes","no","no"), + c("yes","no","yes","no"), + boot.tab[1,], + apply(boot.tab[2:(1+nrep),],2,sd)) + colnames(hildesheim.tab)<-c("Precise HSV","Cervical Cancer", + "Estimated Count","Bootstrap S.E.") + print(hildesheim.tab) + cat("\nEstimated OR=",oddsr[1],"\n") + cat(" Bias=",oddsr[1]-mean(oddsr[2:(1+nrep)]),"\n") + cat("Bootstrap SE=",sd(oddsr[2:(1+nrep)]),"\n\nQuantiles\n\n") + print(quantile(oddsr[2:(1+nrep)],c(0.025,0.50,0.975))) + + b<-mean(log(oddsr[2:(1+nrep)])) + se<-sd(log(oddsr[2:(1+nrep)])) + ztest<-b/se + cat("\n Estimated log(OR)=",log(oddsr[1]),"\n", + "Bootstrap mean log(OR)=",b,"\n", + " Bootstrap SE=",se,"\n", + " Wald Z=",ztest," (P=",2*pnorm(ztest,lower=FALSE),")\n") + } > boot.hildesheim(y,s,~substudy+hsv.inac*hsv.gold*cancer,nrep=50,cancer,hsv.gold) It 0 : 375 535 318 701 18 3 16 11 5 13 16 33 It 1 : 417 537 297 673 15 5 19 15 4 13 13 36 It 2 : 400 547 307 664 21 2 26 4 3 11 23 36 It 3 : 338 534 331 714 26 3 13 18 6 13 15 33 It 4 : 364 539 328 694 15 4 19 8 3 16 27 27 It 5 : 391 538 300 700 18 4 19 7 5 14 17 31 It 6 : 358 556 311 706 20 1 16 8 11 15 13 29 It 7 : 387 517 323 712 17 2 10 6 3 14 13 40 It 8 : 354 548 334 690 22 2 18 7 4 14 20 31 It 9 : 391 542 301 698 15 7 12 17 2 10 16 33 It 10 : 370 514 328 727 20 2 8 10 7 9 15 34 It 11 : 389 528 307 719 16 0.5 18 8 6 11 13 29 It 12 : 396 543 312 669 16 1 18 13 8 12 22 34 It 13 : 349 580 326 683 15 1 12 11 7 17 10 33 It 14 : 357 558 318 687 17 5 17 12 8 11 17 37 It 15 : 387 561 293 696 19 2 17 11 0.5 8 17 33 It 16 : 372 510 312 721 23 0.5 19 12 10 12 12 41 It 17 : 343 532 327 729 14 2 19 11 2 10 15 40 It 18 : 390 521 309 722 12 1 9 8 7 13 16 36 It 19 : 379 526 353 687 14 0.5 15 6 4 10 16 34 It 20 : 367 517 332 708 17 2 17 10 3 15 18 38 It 21 : 376 487 326 728 20 2 25 9 5 12 16 38 It 22 : 375 544 336 673 13 1 21 16 8 11 11 35 It 23 : 395 517 327 694 11 1 14 7 5 18 13 42 It 24 : 367 499 322 736 24 3 10 14 7 14 19 29 It 25 : 369 505 335 728 19 4 12 9 3 13 18 29 It 26 : 402 543 335 646 24 5 11 9 4 16 13 36 It 27 : 356 544 297 722 16 4 15 17 3 13 12 45 It 28 : 378 524 299 704 17 3 12 10 4 17 26 50 It 29 : 409 510 299 689 18 3 20 11 5 14 23 43 It 30 : 366 521 326 719 15 2 21 9 5 12 15 33 It 31 : 349 538 338 714 24 4 13 11 5 10 11 27 It 32 : 379 518 307 711 22 1 26 8 7 10 19 36 It 33 : 396 521 290 723 16 4 21 10 3 18 11 31 It 34 : 384 550 347 647 23 3 16 9 5 16 12 32 It 35 : 369 527 312 710 22 1 10 18 6 19 19 31 It 36 : 398 495 309 738 15 3 14 10 3 8 19 32 It 37 : 380 533 313 714 14 4 11 6 7 11 21 30 It 38 : 404 530 316 672 17 2 25 9 7 12 17 33 It 39 : 388 533 298 714 15 3 21 11 3 15 13 30 It 40 : 405 541 307 688 13 3 6 11 2 17 18 33 It 41 : 401 557 305 672 21 5 10 8 3 14 13 35 It 42 : 366 528 304 720 24 3 17 11 4 8 20 39 It 43 : 366 548 305 703 24 3 13 10 8 10 16 38 It 44 : 390 538 300 690 19 3 18 10 5 18 15 38 It 45 : 391 540 311 699 23 2 12 13 6 5 15 27 It 46 : 347 536 335 706 23 1 13 14 3 10 22 34 It 47 : 381 529 310 712 18 3 13 14 9 12 9 34 It 48 : 393 524 313 689 22 1 16 5 10 13 19 39 It 49 : 384 543 336 676 12 5 9 11 5 14 14 35 It 50 : 397 524 306 713 17 2 14 10 5 15 17 24 Precise HSV Cervical Cancer Estimated Count Bootstrap S.E. 1 yes yes 432.7619 52.73142 2 yes no 299.2381 52.11745 3 no yes 577.9348 80.65161 4 no no 734.0652 83.07459 Estimated OR= 1.836910 Bias= -0.3118162 Bootstrap SE= 0.8480334 Quantiles 2.5% 50% 97.5% 1.080421 1.859414 4.385092 Estimated log(OR)= 0.6080849 Bootstrap mean log(OR)= 0.6949742 Bootstrap SE= 0.3760445 Wald Z= 1.848117 (P= 0.06458543 ) > > > > > cleanEx(); ..nameEx <- "ld2" > > ### * ld2 > > flush(stderr()); flush(stdout()) > > ### Name: ld2 > ### Title: Estimate linkage disequilibrium between two codominant autosomal > ### loci > ### Aliases: ld2 print.ld2 > ### Keywords: category models > > ### ** Examples > > MNS<-matrix(c(91,32,5,147,78,17,85,75,7), nr=3) > colnames(MNS)<-c("S/S","S/s","s/s") > rownames(MNS)<-c("M/M","M/N","N/N") > class(MNS)<-"table" > print(MNS) S/S S/s s/s M/M 91 147 85 M/N 32 78 75 N/N 5 17 7 > res<-ld2(MNS) > print(res) Assuming HWE Call: scoregllm(y = y, s = s, X = X, m = as.array(emgllm(y, s, X, maxit = em.maxit, tol = tol)$full.table)) No. cells in observed table: 9 No. cells in complete table: 16 Mean observed cell size: 59.66667 Model Deviance (df): 12.81 (5) Estimate S.E. exp(Estimate) Lower 95% CL Upper 95% CL Int 4.3622906 0.09434106 78.4365988 65.1949596 94.3677252 a1.2 -1.5452570 0.14863599 0.2132571 0.1593609 0.2853809 a2.2 0.0242353 0.07620604 1.0245314 0.8823835 1.1895786 d.22 0.5501531 0.21141168 1.7335184 1.1454370 2.6235280 Model Resid. df Resid. Dev Pr(GOFChi) Test Df LR stat. Pr(Chi) 1 x$m0 6 20.82376 0.001973257 NA NA NA 2 x$m1 5 12.81211 0.025204571 1 vs 2 1 8.011653 0.004647728 Modelling HWD Call: scoregllm(y = y, s = s, X = X, m = as.array(emgllm(y, s, X, maxit = em.maxit, tol = tol)$full.table)) No. cells in observed table: 9 No. cells in complete table: 16 Mean observed cell size: 59.66667 Model Deviance (df): 7.95 (3) Estimate S.E. exp(Estimate) Lower 95% CL Upper 95% CL Int 4.4692797 0.1022111 87.2938223 71.4462939 106.6564967 a1.2 -1.5638403 0.1529943 0.2093306 0.1550963 0.2825298 a2.2 -0.1768361 0.1174477 0.8379171 0.6656218 1.0548108 p1.22 0.0838282 0.2433765 1.0874421 0.6749006 1.7521546 p2.22 0.3776944 0.1741681 1.4589170 1.0369928 2.0525109 d.22 0.5488856 0.2010023 1.7313226 1.1675658 2.5672883 Model Resid. df Resid. Dev Pr(GOFChi) Test Df LR stat. Pr(Chi) 1 x$m0 6 20.823765 0.0019732568 NA NA NA 2 x$m2 5 20.686490 0.0009282968 1 vs 2 1 0.1372742 0.711005960 3 x$m3 4 15.945790 0.0030927748 2 vs 3 1 4.7407005 0.029457075 4 x$m4 3 7.953696 0.0469784279 3 vs 4 1 7.9920940 0.004698205 > > > > cleanEx(); ..nameEx <- "ld2.model" > > ### * ld2.model > > flush(stderr()); flush(stdout()) > > ### Name: ld2.model > ### Title: Write design and filter matrices for log-linear model of linkage > ### disequilibrium between two codominant autosomal loci > ### Aliases: ld2.model > ### Keywords: category models > > ### ** Examples > > > > > cleanEx(); ..nameEx <- "scatter" > > ### * scatter > > flush(stderr()); flush(stdout()) > > ### Name: scatter > ### Title: Create a filter matrix from a summary array of indices > ### Aliases: scatter > ### Keywords: manip > > ### ** Examples > > y<-double(3) > z<-1:5 > z %*% scatter(y,c(1,1,2,3,3)) [,1] [,2] [,3] [1,] 3 3 9 > ## 1+2, 3, 4+5 > > > > cleanEx(); ..nameEx <- "scoregllm" > > ### * scoregllm > > flush(stderr()); flush(stdout()) > > ### Name: scoregllm > ### Title: Generalized log-linear modelling via Fisher scoring > ### Aliases: scoregllm > ### Keywords: models category > > ### ** Examples > > # > # An example of model fitting: gametic association between two diallelic loci > # Data matrix > # > y<-c( 187,386,156, + 352,310,20, + 136,0 ,0) > # > # Scatter matrix > # > s<- c( 1, 2, 2, 3, + 4, 5, 5, 6, + 4, 5, 5, 6, + 7, 8, 8, 9) > # > # Design matrix > # > X<- matrix(c( 1,0,0,0,0,0,1, + 1,0,1,0,0,0,0, + 1,0,1,0,0,0,0, + 1,0,2,0,1,0,0, + 1,1,0,0,0,0,0, + 1,1,1,0,0,1,0, + 1,1,1,0,0,0,1, + 1,1,2,0,1,1,1, + 1,1,0,0,0,0,0, + 1,1,1,0,0,0,1, + 1,1,1,0,0,1,0, + 1,1,2,0,1,1,1, + 1,2,0,1,0,0,0, + 1,2,1,1,0,1,1, + 1,2,1,1,0,1,1, + 1,2,2,1,1,2,2), byrow=TRUE, ncol=7) > > colnames(X)<-c("Intercept", "A", "B", "P1", "P2", "Delta", "Epsilon") > res<-scoregllm(y,s,X[,c(1:6)], + c(255,176,176,121,164,37,113,25,164,113,37,25,90,20,20,5)) > summary(res) Call: scoregllm(y = y, s = s, X = X[, c(1:6)], m = c(255, 176, 176, 121, 164, 37, 113, 25, 164, 113, 37, 25, 90, 20, 20, 5)) No. cells in observed table: 9 No. cells in complete table: 16 Mean observed cell size: 171.8889 Model Deviance (df): 28.28 (3) Estimate S.E. exp(Estimate) Lower 95% CL Upper 95% CL Intercept 5.303144151 0.05686545 200.96769050 179.77172982 224.66275796 A -0.147717735 0.05852344 0.86267458 0.76918507 0.96752715 B -0.102441387 0.05855679 0.90263105 0.80475878 1.01240625 P1 -0.161218170 0.12289564 0.85110636 0.66891811 1.08291588 P2 0.009094938 0.11562174 1.00913642 0.80450856 1.26581167 Delta -3.289880832 0.24066719 0.03725829 0.02324678 0.05971493 > # > > > > > ### *