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("catspec-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('catspec') > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "ctab" > > ### * ctab > > flush(stderr()); flush(stdout()) > > ### Name: ctab > ### Title: Percentage tables > ### Aliases: ctab print.ctab summary.ctab ctab0 > ### Keywords: category > > ### ** Examples > > ftable(Titanic) Survived No Yes Class Sex Age 1st Male Child 0 5 Adult 118 57 Female Child 0 1 Adult 4 140 2nd Male Child 0 11 Adult 154 14 Female Child 0 13 Adult 13 80 3rd Male Child 35 13 Adult 387 75 Female Child 17 14 Adult 89 76 Crew Male Child 0 0 Adult 670 192 Female Child 0 0 Adult 3 20 > ctab(Titanic) # same output Survived No Yes Class Sex Age 1st Male Child 0 5 Adult 118 57 Female Child 0 1 Adult 4 140 2nd Male Child 0 11 Adult 154 14 Female Child 0 13 Adult 13 80 3rd Male Child 35 13 Adult 387 75 Female Child 17 14 Adult 89 76 Crew Male Child 0 0 Adult 670 192 Female Child 0 0 Adult 3 20 > ctab(Titanic,type="r") Survived No Yes Class Sex Age 1st Male Child 0.00 100.00 Adult 67.43 32.57 Female Child 0.00 100.00 Adult 2.78 97.22 2nd Male Child 0.00 100.00 Adult 91.67 8.33 Female Child 0.00 100.00 Adult 13.98 86.02 3rd Male Child 72.92 27.08 Adult 83.77 16.23 Female Child 54.84 45.16 Adult 53.94 46.06 Crew Male Child NaN NaN Adult 77.73 22.27 Female Child NaN NaN Adult 13.04 86.96 > ctab(Titanic,type=c("n","r"),addmargins=TRUE) Survived No Yes Sum Class Sex Age 1st Male Child Count 0.00 5.00 5.00 Row % 0.00 100.00 100.00 Adult Count 118.00 57.00 175.00 Row % 67.43 32.57 100.00 Sum Count 118.00 62.00 180.00 Row % 67.43 132.57 200.00 Female Child Count 0.00 1.00 1.00 Row % 0.00 100.00 100.00 Adult Count 4.00 140.00 144.00 Row % 2.78 97.22 100.00 Sum Count 4.00 141.00 145.00 Row % 2.78 197.22 200.00 2nd Male Child Count 0.00 11.00 11.00 Row % 0.00 100.00 100.00 Adult Count 154.00 14.00 168.00 Row % 91.67 8.33 100.00 Sum Count 154.00 25.00 179.00 Row % 91.67 108.33 200.00 Female Child Count 0.00 13.00 13.00 Row % 0.00 100.00 100.00 Adult Count 13.00 80.00 93.00 Row % 13.98 86.02 100.00 Sum Count 13.00 93.00 106.00 Row % 13.98 186.02 200.00 3rd Male Child Count 35.00 13.00 48.00 Row % 72.92 27.08 100.00 Adult Count 387.00 75.00 462.00 Row % 83.77 16.23 100.00 Sum Count 422.00 88.00 510.00 Row % 156.68 43.32 200.00 Female Child Count 17.00 14.00 31.00 Row % 54.84 45.16 100.00 Adult Count 89.00 76.00 165.00 Row % 53.94 46.06 100.00 Sum Count 106.00 90.00 196.00 Row % 108.78 91.22 200.00 Crew Male Child Count 0.00 0.00 0.00 Row % 0.00 0.00 0.00 Adult Count 670.00 192.00 862.00 Row % 77.73 22.27 100.00 Sum Count 670.00 192.00 862.00 Row % 77.73 22.27 100.00 Female Child Count 0.00 0.00 0.00 Row % 0.00 0.00 0.00 Adult Count 3.00 20.00 23.00 Row % 13.04 86.96 100.00 Sum Count 3.00 20.00 23.00 Row % 13.04 86.96 100.00 > ctab(Titanic,type=c("n","c","t","r"),style="w") Count Column % Total % Row % Survived No Yes No Yes No Yes No Yes Class Sex Age 1st Male Child 0.00 5.00 0.00 8.06 0.00 2.78 0.00 100.00 Adult 118.00 57.00 100.00 91.94 65.56 31.67 67.43 32.57 Female Child 0.00 1.00 0.00 0.71 0.00 0.69 0.00 100.00 Adult 4.00 140.00 100.00 99.29 2.76 96.55 2.78 97.22 2nd Male Child 0.00 11.00 0.00 44.00 0.00 6.15 0.00 100.00 Adult 154.00 14.00 100.00 56.00 86.03 7.82 91.67 8.33 Female Child 0.00 13.00 0.00 13.98 0.00 12.26 0.00 100.00 Adult 13.00 80.00 100.00 86.02 12.26 75.47 13.98 86.02 3rd Male Child 35.00 13.00 8.29 14.77 6.86 2.55 72.92 27.08 Adult 387.00 75.00 91.71 85.23 75.88 14.71 83.77 16.23 Female Child 17.00 14.00 16.04 15.56 8.67 7.14 54.84 45.16 Adult 89.00 76.00 83.96 84.44 45.41 38.78 53.94 46.06 Crew Male Child 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Adult 670.00 192.00 100.00 100.00 77.73 22.27 77.73 22.27 Female Child 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Adult 3.00 20.00 100.00 100.00 13.04 86.96 13.04 86.96 > mytab<-ftable(Titanic,row.vars=c(1,3),type="r") > mytab Sex Male Female Survived No Yes No Yes Class Age 1st Child 0 5 0 1 Adult 118 57 4 140 2nd Child 0 11 0 13 Adult 154 14 13 80 3rd Child 35 13 17 14 Adult 387 75 89 76 Crew Child 0 0 0 0 Adult 670 192 3 20 > ctab(mytab) Sex Male Female Survived No Yes No Yes Class Age 1st Child 0 5 0 1 Adult 118 57 4 140 2nd Child 0 11 0 13 Adult 154 14 13 80 3rd Child 35 13 17 14 Adult 387 75 89 76 Crew Child 0 0 0 0 Adult 670 192 3 20 > newtab<-ctab(mytab,type="r") > newtab Sex Male Female Survived No Yes No Yes Class Age 1st Child 0.00 100.00 0.00 100.00 Adult 67.43 32.57 2.78 97.22 2nd Child 0.00 100.00 0.00 100.00 Adult 91.67 8.33 13.98 86.02 3rd Child 72.92 27.08 54.84 45.16 Adult 83.77 16.23 53.94 46.06 Crew Child NaN NaN NaN NaN Adult 77.73 22.27 13.04 86.96 > summary(newtab) Number of cases in table: 2201 Number of factors: 4 Test for independence of all factors: Chisq = 1637.4, df = 25, p-value = 0 Chi-squared approximation may be incorrect > ctab0(Titanic,type="r") For future use: ctab(Titanic, type = "r", row.vars=c(4,3,1), col.vars=2) Sex Male Female Survived Age Class No Child 1st NaN NaN 2nd NaN NaN 3rd 67.31 32.69 Crew NaN NaN Adult 1st 96.72 3.28 2nd 92.22 7.78 3rd 81.30 18.70 Crew 99.55 0.45 Yes Child 1st 83.33 16.67 2nd 45.83 54.17 3rd 48.15 51.85 Crew NaN NaN Adult 1st 28.93 71.07 2nd 14.89 85.11 3rd 49.67 50.33 Crew 90.57 9.43 > > #second example using a data frame rather than table data > data(logan) > class(logan) #"data.frame" [1] "data.frame" > ctab(occ) Count Total % occ farm 19.00 2.27 operatives 217.00 25.89 craftsmen 202.00 24.11 sales 105.00 12.53 professional 295.00 35.20 > ctab(occ,addmargins=TRUE) Count Total % occ farm 19.00 2.27 operatives 217.00 25.89 craftsmen 202.00 24.11 sales 105.00 12.53 professional 295.00 35.20 Sum 838.00 100.00 > ctab(occ,style="w",type="c") occ farm operatives craftsmen sales professional 2.27 25.89 24.11 12.53 35.20 > ctab(occ,style="l",type="n") occ farm 19 operatives 217 craftsmen 202 sales 105 professional 295 > z<-ctab(occ,addmargins=TRUE,style="l") > z occ farm Count 19.00 Total % 2.27 operatives Count 217.00 Total % 25.89 craftsmen Count 202.00 Total % 24.11 sales Count 105.00 Total % 12.53 professional Count 295.00 Total % 35.20 Sum Count 838.00 Total % 100.00 > print(z,addmargins=FALSE,dec.places=5) occ farm Count 19.00000 Total % 2.26730 operatives Count 217.00000 Total % 25.89499 craftsmen Count 202.00000 Total % 24.10501 sales Count 105.00000 Total % 12.52983 professional Count 295.00000 Total % 35.20286 > summary(z) Number of cases in table: 838 Number of factors: 1 > > t<-ctab(focc,occ,type=c("n","r","c")) > t occ farm operatives craftsmen sales professional focc farm Count 15.00 29.00 26.00 7.00 15.00 Row % 16.30 31.52 28.26 7.61 16.30 Column % 78.95 13.36 12.87 6.67 5.08 operatives Count 2.00 94.00 54.00 27.00 58.00 Row % 0.85 40.00 22.98 11.49 24.68 Column % 10.53 43.32 26.73 25.71 19.66 craftsmen Count 1.00 55.00 79.00 27.00 70.00 Row % 0.43 23.71 34.05 11.64 30.17 Column % 5.26 25.35 39.11 25.71 23.73 sales Count 0.00 15.00 16.00 8.00 43.00 Row % 0.00 18.29 19.51 9.76 52.44 Column % 0.00 6.91 7.92 7.62 14.58 professional Count 1.00 24.00 27.00 36.00 109.00 Row % 0.51 12.18 13.71 18.27 55.33 Column % 5.26 11.06 13.37 34.29 36.95 > summary(t) Number of cases in table: 838 Number of factors: 2 Test for independence of all factors: Chisq = 201.44, df = 16, p-value = 4.061e-34 Chi-squared approximation may be incorrect > > > > cleanEx(); ..nameEx <- "mcl" > > ### * mcl > > flush(stderr()); flush(stdout()) > > ### Name: mclgen > ### Title: Restructure a data-frame as a "person-choice" file > ### Aliases: mclgen > ### Keywords: multivariate > > ### ** Examples > > ## Example 1 > # data from the Data from the 1972-78 GSS used by Logan (1983) > data(logan) > > # create the "person-choice" file > pc<-mclgen(logan,occ) > summary(pc) occ.1 focc.1 educ.1 black.1 farm : 95 farm : 460 Min. : 2.00 non-black:3820 operatives :1085 operatives :1175 1st Qu.:12.00 black : 370 craftsmen :1010 craftsmen :1160 Median :13.00 sales : 525 sales : 410 Mean :13.58 professional:1475 professional: 985 3rd Qu.:16.00 Max. :20.00 occ.2 focc.2 educ.2 black.2 farm : 95 farm : 460 Min. : 2.00 non-black:3820 operatives :1085 operatives :1175 1st Qu.:12.00 black : 370 craftsmen :1010 craftsmen :1160 Median :13.00 sales : 525 sales : 410 Mean :13.58 professional:1475 professional: 985 3rd Qu.:16.00 Max. :20.00 occ.3 focc.3 educ.3 black.3 farm : 95 farm : 460 Min. : 2.00 non-black:3820 operatives :1085 operatives :1175 1st Qu.:12.00 black : 370 craftsmen :1010 craftsmen :1160 Median :13.00 sales : 525 sales : 410 Mean :13.58 professional:1475 professional: 985 3rd Qu.:16.00 Max. :20.00 occ.4 focc.4 educ.4 black.4 farm : 95 farm : 460 Min. : 2.00 non-black:3820 operatives :1085 operatives :1175 1st Qu.:12.00 black : 370 craftsmen :1010 craftsmen :1160 Median :13.00 sales : 525 sales : 410 Mean :13.58 professional:1475 professional: 985 3rd Qu.:16.00 Max. :20.00 newy occ focc educ Min. :1 farm :838 farm : 460 Min. : 2.00 1st Qu.:2 operatives :838 operatives :1175 1st Qu.:12.00 Median :3 craftsmen :838 craftsmen :1160 Median :13.00 Mean :3 sales :838 sales : 410 Mean :13.58 3rd Qu.:4 professional:838 professional: 985 3rd Qu.:16.00 Max. :5 Max. :20.00 black id depvar non-black:3820 Min. : 1.0 Min. :0.0 black : 370 1st Qu.:210.0 1st Qu.:0.0 Median :419.5 Median :0.0 Mean :419.5 Mean :0.2 3rd Qu.:629.0 3rd Qu.:0.0 Max. :838.0 Max. :1.0 > attach(pc) > > library(survival) Loading required package: splines > # The following specification will work but R won't drop > # cl.lr<-clogit(depvar~occ+occ:educ+occ:black+strata(id),data=pc) > # However, R won't drop the first category of "occ" > # in the interaction effects. The last category will be omitted > # instead due to linear dependence within strata. > # Fix for the problem, create dummies manually for "occ" > occ.X<-model.matrix(~pc$occ) > occ.X<-occ.X[,attributes(occ.X)$assign==1] > cl.lr<-clogit(depvar~occ.X+occ.X:educ+occ.X:black+strata(id),data=pc) > summary(cl.lr) Call: coxph(formula = Surv(rep(1, 4190), depvar) ~ occ.X + occ.X:educ + occ.X:black + strata(id), data = pc, method = "exact") n= 4190 coef exp(coef) se(coef) z p occ.Xpc$occoperatives 2.9135 18.42203 1.374 2.121 3.4e-02 occ.Xpc$occcraftsmen 1.8433 6.31713 1.382 1.334 1.8e-01 occ.Xpc$occsales -3.1389 0.04333 1.476 -2.127 3.3e-02 occ.Xpc$occprofessional -6.1314 0.00217 1.441 -4.254 2.1e-05 occ.Xpc$occoperatives:educ -0.0505 0.95074 0.111 -0.456 6.5e-01 occ.Xpc$occcraftsmen:educ 0.0387 1.03943 0.111 0.348 7.3e-01 occ.Xpc$occsales:educ 0.3692 1.44659 0.116 3.175 1.5e-03 occ.Xpc$occprofessional:educ 0.6440 1.90399 0.114 5.671 1.4e-08 occ.Xpc$occoperatives:blackblack 1.3052 3.68827 1.043 1.251 2.1e-01 occ.Xpc$occcraftsmen:blackblack 0.6284 1.87454 1.055 0.595 5.5e-01 occ.Xpc$occsales:blackblack 0.3326 1.39462 1.103 0.301 7.6e-01 occ.Xpc$occprofessional:blackblack -0.2259 0.79781 1.093 -0.207 8.4e-01 exp(coef) exp(-coef) lower .95 upper .95 occ.Xpc$occoperatives 18.42203 0.0543 1.247069 272.1349 occ.Xpc$occcraftsmen 6.31713 0.1583 0.421248 94.7331 occ.Xpc$occsales 0.04333 23.0783 0.002402 0.7815 occ.Xpc$occprofessional 0.00217 460.0593 0.000129 0.0366 occ.Xpc$occoperatives:educ 0.95074 1.0518 0.765108 1.1814 occ.Xpc$occcraftsmen:educ 1.03943 0.9621 0.835972 1.2924 occ.Xpc$occsales:educ 1.44659 0.6913 1.151724 1.8169 occ.Xpc$occprofessional:educ 1.90399 0.5252 1.524078 2.3786 occ.Xpc$occoperatives:blackblack 3.68827 0.2711 0.477308 28.5000 occ.Xpc$occcraftsmen:blackblack 1.87454 0.5335 0.236896 14.8331 occ.Xpc$occsales:blackblack 1.39462 0.7170 0.160386 12.1267 occ.Xpc$occprofessional:blackblack 0.79781 1.2534 0.093626 6.7983 Rsquare= 0.15 (max possible= 0.475 ) Likelihood ratio test= 683 on 12 df, p=0 Wald test = 421 on 12 df, p=0 Score (logrank) test = 703 on 12 df, p=0 > > # Estimate a "quasi-uniform association" loglinear model for "focc" and "occ" > # with "educ" and "black" as covariates at the respondent level > cl.qu<-clogit(depvar~occ.X+occ.X:educ+occ.X:black+ + mob.qi(focc,occ)+mob.unif(focc,occ)+strata(id),data=pc) > summary(cl.qu) Call: coxph(formula = Surv(rep(1, 4190), depvar) ~ occ.X + occ.X:educ + occ.X:black + mob.qi(focc, occ) + mob.unif(focc, occ) + strata(id), data = pc, method = "exact") n= 4190 coef exp(coef) se(coef) z p occ.Xpc$occoperatives 4.9676 143.6760 1.7264 2.877 4.0e-03 occ.Xpc$occcraftsmen 3.6728 39.3635 1.7509 2.098 3.6e-02 occ.Xpc$occsales -1.0365 0.3547 1.8184 -0.570 5.7e-01 occ.Xpc$occprofessional -4.3603 0.0128 1.8238 -2.391 1.7e-02 mob.qi(focc, occ)2 2.9339 18.8005 0.6043 4.855 1.2e-06 mob.qi(focc, occ)3 0.2347 1.2646 0.2036 1.153 2.5e-01 mob.qi(focc, occ)4 0.5355 1.7083 0.1891 2.832 4.6e-03 mob.qi(focc, occ)5 -0.4538 0.6352 0.3947 -1.150 2.5e-01 mob.qi(focc, occ)6 -0.3347 0.7156 0.2680 -1.249 2.1e-01 mob.unif(focc, occ) 0.1118 1.1183 0.0403 2.775 5.5e-03 occ.Xpc$occoperatives:educ -0.1379 0.8712 0.1295 -1.065 2.9e-01 occ.Xpc$occcraftsmen:educ -0.0646 0.9374 0.1294 -0.499 6.2e-01 occ.Xpc$occsales:educ 0.2356 1.2657 0.1327 1.776 7.6e-02 occ.Xpc$occprofessional:educ 0.5109 1.6668 0.1301 3.927 8.6e-05 occ.Xpc$occoperatives:blackblack 1.7547 5.7818 1.0686 1.642 1.0e-01 occ.Xpc$occcraftsmen:blackblack 1.3101 3.7064 1.0827 1.210 2.3e-01 occ.Xpc$occsales:blackblack 0.9633 2.6203 1.1321 0.851 3.9e-01 occ.Xpc$occprofessional:blackblack 0.5038 1.6550 1.1241 0.448 6.5e-01 exp(coef) exp(-coef) lower .95 upper .95 occ.Xpc$occoperatives 143.6760 0.00696 4.873848 4235.419 occ.Xpc$occcraftsmen 39.3635 0.02540 1.272614 1217.559 occ.Xpc$occsales 0.3547 2.81946 0.010047 12.520 occ.Xpc$occprofessional 0.0128 78.28074 0.000358 0.456 mob.qi(focc, occ)2 18.8005 0.05319 5.751157 61.459 mob.qi(focc, occ)3 1.2646 0.79078 0.848410 1.885 mob.qi(focc, occ)4 1.7083 0.58537 1.179281 2.475 mob.qi(focc, occ)5 0.6352 1.57435 0.293067 1.377 mob.qi(focc, occ)6 0.7156 1.39747 0.423184 1.210 mob.unif(focc, occ) 1.1183 0.89421 1.033396 1.210 occ.Xpc$occoperatives:educ 0.8712 1.14785 0.675927 1.123 occ.Xpc$occcraftsmen:educ 0.9374 1.06677 0.727379 1.208 occ.Xpc$occsales:educ 1.2657 0.79010 0.975897 1.641 occ.Xpc$occprofessional:educ 1.6668 0.59996 1.291635 2.151 occ.Xpc$occoperatives:blackblack 5.7818 0.17296 0.712004 46.950 occ.Xpc$occcraftsmen:blackblack 3.7064 0.26980 0.443948 30.944 occ.Xpc$occsales:blackblack 2.6203 0.38163 0.284919 24.098 occ.Xpc$occprofessional:blackblack 1.6550 0.60425 0.182795 14.983 Rsquare= 0.165 (max possible= 0.475 ) Likelihood ratio test= 756 on 18 df, p=0 Wald test = 398 on 18 df, p=0 Score (logrank) test = 748 on 18 df, p=0 > > data(housing,package="MASS") > housing.prsch<-mclgen(housing,Sat) > library(survival) > # clogit doesn't support the weights argument at present > # a work-around is to call coxph directly > # coxph warns that X is singular, because the main > # effects of Infl, Type, and Cont are dropped > coxph.prsch<-coxph(Surv(rep(1, NROW(housing.prsch)), depvar) ~ + Sat+Sat*Infl+Sat*Type+Sat*Cont+strata(id), + weights = housing.prsch$Freq, data = housing.prsch) Warning in coxph(Surv(rep(1, NROW(housing.prsch)), depvar) ~ Sat + Sat * : X matrix deemed to be singular; variable 3 4 5 6 7 8 > summary(coxph.prsch) Call: coxph(formula = Surv(rep(1, NROW(housing.prsch)), depvar) ~ Sat + Sat * Infl + Sat * Type + Sat * Cont + strata(id), data = housing.prsch, weights = housing.prsch$Freq) n= 216 coef exp(coef) se(coef) z p SatMedium -0.419 0.658 0.173 -2.424 1.5e-02 SatHigh -0.139 0.870 0.159 -0.871 3.8e-01 InflMedium NA NA 0.000 NA NA InflHigh NA NA 0.000 NA NA TypeApartment NA NA 0.000 NA NA TypeAtrium NA NA 0.000 NA NA TypeTerrace NA NA 0.000 NA NA ContHigh NA NA 0.000 NA NA SatMedium:InflMedium 0.446 1.563 0.142 3.153 1.6e-03 SatHigh:InflMedium 0.735 2.085 0.137 5.366 8.0e-08 SatMedium:InflHigh 0.665 1.944 0.186 3.568 3.6e-04 SatHigh:InflHigh 1.613 5.016 0.167 9.649 0.0e+00 SatMedium:TypeApartment -0.436 0.647 0.173 -2.525 1.2e-02 SatHigh:TypeApartment -0.736 0.479 0.155 -4.738 2.2e-06 SatMedium:TypeAtrium 0.131 1.140 0.223 0.589 5.6e-01 SatHigh:TypeAtrium -0.408 0.665 0.211 -1.929 5.4e-02 SatMedium:TypeTerrace -0.667 0.513 0.206 -3.232 1.2e-03 SatHigh:TypeTerrace -1.412 0.244 0.200 -7.056 1.7e-12 SatMedium:ContHigh 0.361 1.435 0.132 2.726 6.4e-03 SatHigh:ContHigh 0.482 1.619 0.124 3.881 1.0e-04 exp(coef) exp(-coef) lower .95 upper .95 SatMedium 0.658 1.521 0.469 0.923 SatHigh 0.870 1.149 0.637 1.189 InflMedium NA NA NA NA InflHigh NA NA NA NA TypeApartment NA NA NA NA TypeAtrium NA NA NA NA TypeTerrace NA NA NA NA ContHigh NA NA NA NA SatMedium:InflMedium 1.563 0.640 1.184 2.062 SatHigh:InflMedium 2.085 0.480 1.594 2.727 SatMedium:InflHigh 1.944 0.514 1.349 2.801 SatHigh:InflHigh 5.016 0.199 3.615 6.960 SatMedium:TypeApartment 0.647 1.546 0.461 0.907 SatHigh:TypeApartment 0.479 2.087 0.353 0.650 SatMedium:TypeAtrium 1.140 0.877 0.736 1.766 SatHigh:TypeAtrium 0.665 1.504 0.439 1.007 SatMedium:TypeTerrace 0.513 1.948 0.343 0.769 SatHigh:TypeTerrace 0.244 4.106 0.165 0.361 SatMedium:ContHigh 1.435 0.697 1.107 1.860 SatHigh:ContHigh 1.619 0.618 1.269 2.065 Rsquare= 0.645 (max possible= 1 ) Likelihood ratio test= 223 on 14 df, p=0 Wald test = 213 on 14 df, p=0 Score (logrank) test = 233 on 14 df, p=0 > > # the same model using multinomial logistic regression > library(nnet) > house.mult<- multinom(Sat ~ Infl + Type + Cont, weights = Freq, + data = housing) # weights: 24 (14 variable) initial value 1846.767257 iter 10 value 1747.045232 final value 1735.041933 converged > summary(house.mult,correlation=FALSE) Call: multinom(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq) Coefficients: (Intercept) InflMedium InflHigh TypeApartment TypeAtrium TypeTerrace Medium -0.4192316 0.4464003 0.6649367 -0.4356851 0.1313663 -0.6665728 High -0.1387453 0.7348626 1.6126294 -0.7356261 -0.4079808 -1.4123333 ContHigh Medium 0.3608513 High 0.4818236 Std. Errors: (Intercept) InflMedium InflHigh TypeApartment TypeAtrium TypeTerrace Medium 0.1729344 0.1415572 0.1863374 0.1725327 0.2231065 0.2062532 High 0.1592295 0.1369380 0.1671316 0.1552714 0.2114965 0.2001496 ContHigh Medium 0.1323975 High 0.1241371 Residual Deviance: 3470.084 AIC: 3498.084 > > # compare the coefficients > m1<-coef(coxph.prsch) > m1<-m1[!is.na(m1)] > dim(m1)<-c(2,7) > m2<-coef(house.mult) > m1 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] -0.4192287 0.4463959 0.6649353 -0.4356887 0.1313703 -0.6665705 0.3608519 [2,] -0.1387428 0.7348632 1.6126311 -0.7356317 -0.4079781 -1.4123277 0.4818270 > m2 (Intercept) InflMedium InflHigh TypeApartment TypeAtrium TypeTerrace Medium -0.4192316 0.4464003 0.6649367 -0.4356851 0.1313663 -0.6665728 High -0.1387453 0.7348626 1.6126294 -0.7356261 -0.4079808 -1.4123333 ContHigh Medium 0.3608513 High 0.4818236 > m1-m2 (Intercept) InflMedium InflHigh TypeApartment TypeAtrium Medium 2.818297e-06 -4.371143e-06 -1.344318e-06 -3.570156e-06 3.972080e-06 High 2.498949e-06 6.589361e-07 1.696236e-06 -5.648022e-06 2.760696e-06 TypeTerrace ContHigh Medium 2.372064e-06 5.397481e-07 High 5.656258e-06 3.363651e-06 > max(abs(m1-m2)) [1] 5.656258e-06 > mean(abs(m1-m2)) [1] 2.947897e-06 > > > > cleanEx(); ..nameEx <- "sqtab" > > ### * sqtab > > flush(stderr()); flush(stdout()) > > ### Name: sqtab > ### Title: sqtab: models for square tables > ### Aliases: sqtab check.square fitmacro mob.cp mob.eqmain mob.qi mob.rc1 > ### mob.symint mob.unif > ### Keywords: multivariate > > ### ** Examples > > # Examples of loglinear models for square tables, > # from Hout, M. (1983). "Mobility Tables". Sage Publication 07-031 > > # Table from page 11 of "Mobility Tables" > # Original source: Featherman D.L., R.M. Hauser. (1978) "Opportunity and Change." > # New York: Academic, page 49 > > data(FHtab) > FHtab<-as.data.frame(FHtab) > attach(FHtab) > > xtabs(Freq ~ .,FHtab) OccSon OccFather Upper nonmanual Lower nonmanual Upper manual Lower manual Upper nonmanual 1414 521 302 643 Lower nonmanual 724 524 254 703 Upper manual 798 648 856 1676 Lower manual 756 914 771 3325 Farm 409 357 441 1611 OccSon OccFather Farm Upper nonmanual 40 Lower nonmanual 48 Upper manual 108 Lower manual 237 Farm 1832 > > # independence model > indep<-glm(Freq~OccFather+OccSon,family=poisson()) > summary(indep) Call: glm(formula = Freq ~ OccFather + OccSon, family = poisson()) Deviance Residuals: Min 1Q Median 3Q Max -20.371 -14.730 -2.554 4.015 44.109 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 6.39925 0.02315 276.383 <2e-16 *** OccFatherLower nonmanual -0.25932 0.02804 -9.248 <2e-16 *** OccFatherUpper manual 0.33598 0.02423 13.865 <2e-16 *** OccFatherLower manual 0.72068 0.02256 31.942 <2e-16 *** OccFatherFarm 0.46528 0.02361 19.706 <2e-16 *** OccSonLower nonmanual -0.32469 0.02411 -13.468 <2e-16 *** OccSonUpper manual -0.44653 0.02500 -17.862 <2e-16 *** OccSonLower manual 0.66295 0.01922 34.488 <2e-16 *** OccSonFarm -0.59366 0.02618 -22.677 <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: 13132.1 on 24 degrees of freedom Residual deviance: 6170.1 on 16 degrees of freedom AIC: 6390.8 Number of Fisher Scoring iterations: 6 > fitmacro(indep) deviance: 6170.130 df: 16 bic: 6011.745 aic: 6138.130 Number of parameters: 9 Number of cases: 19912 > > wt <- (OccFather != OccSon) > qi0<-glm(Freq~OccFather+OccSon,weights=wt,family=poisson()) > # A quasi-independence loglinear model, using structural zeros > # (page 23 of "Mobility Tables"). > # 0 1 1 1 1 values of variable "wt" > # 1 0 1 1 1 > # 1 1 0 1 1 > # 1 1 1 0 1 > # 1 1 1 1 0 > qi0<-glm(Freq~OccFather+OccSon,weights=wt,family=poisson()) > summary(qi0) Call: glm(formula = Freq ~ OccFather + OccSon, family = poisson(), weights = wt) Deviance Residuals: Min 1Q Median 3Q Max -8.181 -3.336 0.000 1.558 13.465 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.91935 0.03433 172.444 < 2e-16 *** OccFatherLower nonmanual 0.11956 0.03588 3.332 0.000861 *** OccFatherUpper manual 0.70710 0.03188 22.179 < 2e-16 *** OccFatherLower manual 0.92052 0.03430 26.838 < 2e-16 *** OccFatherFarm 0.44093 0.03230 13.653 < 2e-16 *** OccSonLower nonmanual -0.07868 0.02847 -2.764 0.005710 ** OccSonUpper manual -0.26496 0.03137 -8.446 < 2e-16 *** OccSonLower manual 0.77953 0.02591 30.081 < 2e-16 *** OccSonFarm -1.74558 0.05210 -33.506 < 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: 6123.71 on 19 degrees of freedom Residual deviance: 683.34 on 11 degrees of freedom AIC: 858.97 Number of Fisher Scoring iterations: 4 > fitmacro(qi0) deviance: 683.342 df: 11 bic: 574.452 aic: 661.342 Number of parameters: 9 Number of cases: 19912 > > # Quasi-independence using a "dummy factor" to create the design > # vectors for the diagonal cells (page 23). > # 1 0 0 0 0 > # 0 2 0 0 0 > # 0 0 3 0 0 > # 0 0 0 4 0 > # 0 0 0 0 5 > glm.qi<-glm(Freq~OccFather+OccSon+mob.qi(OccFather,OccSon),family=poisson()) > summary(glm.qi) Call: glm(formula = Freq ~ OccFather + OccSon + mob.qi(OccFather, OccSon), family = poisson()) Deviance Residuals: Min 1Q Median 3Q Max -8.181e+00 -3.336e+00 -4.108e-07 1.558e+00 1.347e+01 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.91935 0.03433 172.444 < 2e-16 *** OccFatherLower nonmanual 0.11956 0.03588 3.332 0.000861 *** OccFatherUpper manual 0.70710 0.03188 22.179 < 2e-16 *** OccFatherLower manual 0.92052 0.03430 26.838 < 2e-16 *** OccFatherFarm 0.44093 0.03230 13.653 < 2e-16 *** OccSonLower nonmanual -0.07868 0.02847 -2.764 0.005710 ** OccSonUpper manual -0.26496 0.03137 -8.446 < 2e-16 *** OccSonLower manual 0.77953 0.02591 30.081 < 2e-16 *** OccSonFarm -1.74558 0.05210 -33.506 < 2e-16 *** mob.qi(OccFather, OccSon)2 1.33483 0.04342 30.741 < 2e-16 *** mob.qi(OccFather, OccSon)3 0.30126 0.05512 5.465 4.62e-08 *** mob.qi(OccFather, OccSon)4 0.39079 0.04695 8.323 < 2e-16 *** mob.qi(OccFather, OccSon)5 0.48983 0.03330 14.711 < 2e-16 *** mob.qi(OccFather, OccSon)6 2.89847 0.05774 50.197 < 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: 13132.14 on 24 degrees of freedom Residual deviance: 683.34 on 11 degrees of freedom AIC: 914.05 Number of Fisher Scoring iterations: 4 > fitmacro(glm.qi) deviance: 683.342 df: 11 bic: 574.452 aic: 661.342 Number of parameters: 14 Number of cases: 19912 > > # Quasi-independence constrained (QPM-C, page 31) > # Single immobility parameter > # 1 0 0 0 0 > # 0 1 0 0 0 > # 0 0 1 0 0 > # 0 0 0 1 0 > # 0 0 0 0 1 > glm.q0<-glm(Freq~OccFather+OccSon+mob.qi(OccFather,OccSon,constrained=TRUE),family=poisson()) > # slightly different results than Hout also found in Stata: L2=2567.658, q0=0.964 > summary(glm.q0) Call: glm(formula = Freq ~ OccFather + OccSon + mob.qi(OccFather, OccSon, constrained = TRUE), family = poisson()) Deviance Residuals: Min 1Q Median 3Q Max -17.784 -8.712 1.196 4.545 20.115 Coefficients: Estimate Std. Error z value (Intercept) 6.16634 0.02193 281.135 OccFatherLower nonmanual -0.18999 0.02866 -6.628 OccFatherUpper manual 0.45477 0.02500 18.194 OccFatherLower manual 0.54946 0.02351 23.369 OccFatherFarm 0.61469 0.02445 25.142 OccSonLower nonmanual -0.29151 0.02462 -11.840 OccSonUpper manual -0.54963 0.02574 -21.350 OccSonLower manual 0.53357 0.02001 26.670 OccSonFarm -0.74220 0.02705 -27.441 mob.qi(OccFather, OccSon, constrained = TRUE)2 0.96449 0.01550 62.221 Pr(>|z|) (Intercept) < 2e-16 *** OccFatherLower nonmanual 3.39e-11 *** OccFatherUpper manual < 2e-16 *** OccFatherLower manual < 2e-16 *** OccFatherFarm < 2e-16 *** OccSonLower nonmanual < 2e-16 *** OccSonUpper manual < 2e-16 *** OccSonLower manual < 2e-16 *** OccSonFarm < 2e-16 *** mob.qi(OccFather, OccSon, constrained = TRUE)2 < 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: 13132.1 on 24 degrees of freedom Residual deviance: 2567.7 on 15 degrees of freedom AIC: 2790.4 Number of Fisher Scoring iterations: 5 > fitmacro(glm.q0) deviance: 2567.658 df: 15 bic: 2419.172 aic: 2537.658 Number of parameters: 10 Number of cases: 19912 > > # Quasi-symmetry using the symmetric cross-classification (page 23) > # 0 1 2 3 4 values of variable "sym" > # 1 0 5 6 7 > # 2 5 0 8 9 > # 3 6 8 0 10 > # 4 7 9 10 0 */ > glm.qsym<- + glm(Freq~OccFather+OccSon+mob.symint(OccFather,OccSon),family=poisson()) > summary(glm.qsym) Call: glm(formula = Freq ~ OccFather + OccSon + mob.symint(OccFather, OccSon), family = poisson()) Deviance Residuals: Min 1Q Median 3Q Max -2.3247 -0.4558 0.0000 0.6836 2.0935 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 7.25418 0.02659 272.780 <2e-16 *** OccFatherLower nonmanual -0.74821 0.04311 -17.355 <2e-16 *** OccFatherUpper manual -0.55598 0.04166 -13.345 <2e-16 *** OccFatherLower manual -0.56958 0.04100 -13.891 <2e-16 *** OccFatherFarm -1.24550 0.05447 -22.868 <2e-16 *** OccSonLower nonmanual -0.89828 0.04450 -20.184 <2e-16 *** OccSonUpper manual -1.58756 0.05056 -31.397 <2e-16 *** OccSonLower manual -0.85890 0.04341 -19.786 <2e-16 *** OccSonFarm -3.51518 0.07705 -45.624 <2e-16 *** mob.symint(OccFather, OccSon)2:2 0.32690 0.03820 8.557 <2e-16 *** mob.symint(OccFather, OccSon)2:3 0.65816 0.06032 10.910 <2e-16 *** mob.symint(OccFather, OccSon)2:4 0.97607 0.05347 18.256 <2e-16 *** mob.symint(OccFather, OccSon)2:5 0.78009 0.08074 9.662 <2e-16 *** mob.symint(OccFather, OccSon)3:3 0.82082 0.03841 21.370 <2e-16 *** mob.symint(OccFather, OccSon)3:4 1.57396 0.05351 29.414 <2e-16 *** mob.symint(OccFather, OccSon)3:5 1.63239 0.07953 20.526 <2e-16 *** mob.symint(OccFather, OccSon)4:4 1.14176 0.03121 36.583 <2e-16 *** mob.symint(OccFather, OccSon)4:5 2.24279 0.06714 33.406 <2e-16 *** mob.symint(OccFather, OccSon)5:5 2.50983 0.05607 44.763 <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: 13132.143 on 24 degrees of freedom Residual deviance: 27.448 on 6 degrees of freedom AIC: 268.16 Number of Fisher Scoring iterations: 3 > fitmacro(glm.qsym) deviance: 27.448 df: 6 bic: -31.947 aic: 15.448 Number of parameters: 19 Number of cases: 19912 > > symmetry<-glm(Freq~mob.eqmain(OccFather,OccSon) + +mob.symint(OccFather,OccSon),family=poisson()) > summary(symmetry) Call: glm(formula = Freq ~ mob.eqmain(OccFather, OccSon) + mob.symint(OccFather, OccSon), family = poisson()) Deviance Residuals: Min 1Q Median 3Q Max -2.700e+01 -1.012e+01 2.859e-07 8.701e+00 2.042e+01 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 7.25418 0.02659 272.780 <2e-16 *** mob.eqmain(OccFather, OccSon)2 -0.82043 0.03886 -21.110 <2e-16 *** mob.eqmain(OccFather, OccSon)3 -0.94426 0.04020 -23.487 <2e-16 *** mob.eqmain(OccFather, OccSon)4 -0.70381 0.03771 -18.664 <2e-16 *** mob.eqmain(OccFather, OccSon)5 -1.84030 0.05417 -33.973 <2e-16 *** mob.symint(OccFather, OccSon)2:2 0.32409 0.03817 8.490 <2e-16 *** mob.symint(OccFather, OccSon)2:3 0.62198 0.05940 10.471 <2e-16 *** mob.symint(OccFather, OccSon)2:4 0.96525 0.05333 18.101 <2e-16 *** mob.symint(OccFather, OccSon)2:5 0.71730 0.07878 9.105 <2e-16 *** mob.symint(OccFather, OccSon)3:3 0.69331 0.03712 18.677 <2e-16 *** mob.symint(OccFather, OccSon)3:4 1.50336 0.05234 28.722 <2e-16 *** mob.symint(OccFather, OccSon)3:5 1.14534 0.07527 15.217 <2e-16 *** mob.symint(OccFather, OccSon)4:4 1.13134 0.03109 36.385 <2e-16 *** mob.symint(OccFather, OccSon)4:5 2.11865 0.06473 32.729 <2e-16 *** mob.symint(OccFather, OccSon)5:5 1.96979 0.05040 39.081 <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: 13132.1 on 24 degrees of freedom Residual deviance: 2804.9 on 10 degrees of freedom AIC: 3037.6 Number of Fisher Scoring iterations: 5 > fitmacro(symmetry) deviance: 2804.873 df: 10 bic: 2705.882 aic: 2784.873 Number of parameters: 15 Number of cases: 19912 > > # Crossings parameter model (page 35) > # 0 v1 v1 v1 v1 | 0 0 v2 v2 v2 | 0 0 0 v3 v3 | 0 0 0 0 v4 > # v1 0 0 0 0 | 0 0 v2 v2 v2 | 0 0 0 v3 v3 | 0 0 0 0 v4 > # v1 0 0 0 0 | v2 v2 0 0 0 | 0 0 0 v3 v3 | 0 0 0 0 v4 > # v1 0 0 0 0 | v2 v2 0 0 0 | v3 v3 v3 0 0 | 0 0 0 0 v4 > # v1 0 0 0 0 | v2 v2 0 0 0 | v3 v3 v3 0 0 | v4 v4 v4 v4 0 > glm.cp<-glm(Freq~OccFather+OccSon+mob.cp(OccFather,OccSon),family=poisson()) > summary(glm.cp) Call: glm(formula = Freq ~ OccFather + OccSon + mob.cp(OccFather, OccSon), family = poisson()) Deviance Residuals: Min 1Q Median 3Q Max -3.1888 -0.8059 0.0000 0.6484 5.2990 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 7.25418 0.02659 272.780 < 2e-16 *** OccFatherLower nonmanual -0.35986 0.02879 -12.500 < 2e-16 *** OccFatherUpper manual 0.20073 0.02568 7.816 5.45e-15 *** OccFatherLower manual 0.57421 0.02441 23.521 < 2e-16 *** OccFatherFarm 1.25914 0.03463 36.365 < 2e-16 *** OccSonLower nonmanual -0.53016 0.02700 -19.637 < 2e-16 *** OccSonUpper manual -0.79953 0.02810 -28.449 < 2e-16 *** OccSonLower manual 0.28056 0.02344 11.969 < 2e-16 *** OccSonFarm -1.00015 0.03463 -28.885 < 2e-16 *** mob.cp(OccFather, OccSon)1 -0.42562 0.02394 -17.776 < 2e-16 *** mob.cp(OccFather, OccSon)2 -0.36748 0.02115 -17.372 < 2e-16 *** mob.cp(OccFather, OccSon)3 -0.29347 0.01728 -16.980 < 2e-16 *** mob.cp(OccFather, OccSon)4 -1.40259 0.02894 -48.458 < 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: 13132.143 on 24 degrees of freedom Residual deviance: 89.914 on 12 degrees of freedom AIC: 318.63 Number of Fisher Scoring iterations: 4 > fitmacro(glm.cp) deviance: 89.914 df: 12 bic: -28.875 aic: 65.914 Number of parameters: 13 Number of cases: 19912 > > # Uniform association model: linear by linear association (page 58) > glm.unif<-glm(Freq~OccFather+OccSon+mob.unif(OccFather,OccSon),family=poisson()) > summary(glm.unif) Call: glm(formula = Freq ~ OccFather + OccSon + mob.unif(OccFather, OccSon), family = poisson()) Deviance Residuals: Min 1Q Median 3Q Max -24.103 -3.302 -1.461 4.088 25.781 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 6.980009 0.021329 327.25 <2e-16 *** OccFatherLower nonmanual -0.860564 0.029168 -29.50 <2e-16 *** OccFatherUpper manual -0.991014 0.031298 -31.66 <2e-16 *** OccFatherLower manual -1.457724 0.042152 -34.58 <2e-16 *** OccFatherFarm -2.672356 0.059456 -44.95 <2e-16 *** OccSonLower nonmanual -1.046149 0.026393 -39.64 <2e-16 *** OccSonUpper manual -2.018346 0.035704 -56.53 <2e-16 *** OccSonLower manual -1.876387 0.047874 -39.19 <2e-16 *** OccSonFarm -4.195326 0.069883 -60.03 <2e-16 *** mob.unif(OccFather, OccSon) 0.269259 0.004804 56.05 <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: 13132.1 on 24 degrees of freedom Residual deviance: 2280.7 on 15 degrees of freedom AIC: 2503.4 Number of Fisher Scoring iterations: 5 > fitmacro(glm.unif) deviance: 2280.694 df: 15 bic: 2132.207 aic: 2250.694 Number of parameters: 10 Number of cases: 19912 > > # RC model 1 (unequal row and column effects, page 58) > # Fits a uniform association parameter and row and column effect > # parameters. Row and column effect parameters have the > # restriction that the first and last categories are zero. > glm.rc1<-glm(Freq~OccFather+OccSon+mob.rc1(OccFather,OccSon),family=poisson()) > summary(glm.rc1) Call: glm(formula = Freq ~ OccFather + OccSon + mob.rc1(OccFather, OccSon), family = poisson()) Deviance Residuals: Min 1Q Median 3Q Max -15.254 -2.453 -1.352 4.013 13.785 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 6.672678 0.026313 253.585 < 2e-16 *** OccFatherLower nonmanual -1.154320 0.060824 -18.978 < 2e-16 *** OccFatherUpper manual -1.617053 0.061608 -26.247 < 2e-16 *** OccFatherLower manual -2.048505 0.067952 -30.147 < 2e-16 *** OccFatherFarm -3.636547 0.078946 -46.064 < 2e-16 *** OccSonLower nonmanual -1.125905 0.057383 -19.621 < 2e-16 *** OccSonUpper manual -1.814371 0.066842 -27.144 < 2e-16 *** OccSonLower manual -1.376182 0.060296 -22.824 < 2e-16 *** OccSonFarm -8.992250 0.214486 -41.925 < 2e-16 *** mob.rc1(OccFather, OccSon)R2 0.029654 0.020919 1.418 0.156 mob.rc1(OccFather, OccSon)R3 0.079253 0.017742 4.467 7.93e-06 *** mob.rc1(OccFather, OccSon)R4 0.002588 0.017053 0.152 0.879 mob.rc1(OccFather, OccSon)U 0.533707 0.011078 48.178 < 2e-16 *** mob.rc1(OccFather, OccSon)C2 -0.252840 0.020567 -12.293 < 2e-16 *** mob.rc1(OccFather, OccSon)C3 -0.614921 0.028553 -21.536 < 2e-16 *** mob.rc1(OccFather, OccSon)C4 -0.960196 0.034636 -27.722 < 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: 13132.1 on 24 degrees of freedom Residual deviance: 877.8 on 9 degrees of freedom AIC: 1112.5 Number of Fisher Scoring iterations: 5 > fitmacro(glm.rc1) deviance: 877.808 df: 9 bic: 788.717 aic: 859.808 Number of parameters: 16 Number of cases: 19912 > > # Homogeneous row and column effects model 1 (page 58) > # An equality restriction is placed on the row and column effects > glm.hrc1<-glm(Freq~OccFather+OccSon+mob.rc1(OccFather,OccSon,equal=TRUE),family=poisson()) > # Results differ from those in Hout, replicated by other programs > summary(glm.hrc1) Call: glm(formula = Freq ~ OccFather + OccSon + mob.rc1(OccFather, OccSon, equal = TRUE), family = poisson()) Deviance Residuals: Min 1Q Median 3Q Max -15.801 -3.195 -1.188 3.070 17.434 Coefficients: Estimate Std. Error z value (Intercept) 6.798987 0.025260 269.157 OccFatherLower nonmanual -0.982508 0.047712 -20.592 OccFatherUpper manual -1.022509 0.048903 -20.909 OccFatherLower manual -1.141834 0.050288 -22.706 OccFatherFarm -3.752391 0.083624 -44.872 OccSonLower nonmanual -1.170357 0.050894 -22.996 OccSonUpper manual -2.049682 0.057562 -35.608 OccSonLower manual -1.570274 0.055556 -28.264 OccSonFarm -5.599423 0.100750 -55.577 mob.rc1(OccFather, OccSon, equal = TRUE)RC2 -0.035255 0.012639 -2.789 mob.rc1(OccFather, OccSon, equal = TRUE)RC3 -0.141634 0.011321 -12.510 mob.rc1(OccFather, OccSon, equal = TRUE)RC4 -0.313500 0.011305 -27.730 mob.rc1(OccFather, OccSon, equal = TRUE)U 0.392262 0.007312 53.649 Pr(>|z|) (Intercept) < 2e-16 *** OccFatherLower nonmanual < 2e-16 *** OccFatherUpper manual < 2e-16 *** OccFatherLower manual < 2e-16 *** OccFatherFarm < 2e-16 *** OccSonLower nonmanual < 2e-16 *** OccSonUpper manual < 2e-16 *** OccSonLower manual < 2e-16 *** OccSonFarm < 2e-16 *** mob.rc1(OccFather, OccSon, equal = TRUE)RC2 0.00528 ** mob.rc1(OccFather, OccSon, equal = TRUE)RC3 < 2e-16 *** mob.rc1(OccFather, OccSon, equal = TRUE)RC4 < 2e-16 *** mob.rc1(OccFather, OccSon, equal = TRUE)U < 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: 13132.1 on 24 degrees of freedom Residual deviance: 1451.3 on 12 degrees of freedom AIC: 1680.1 Number of Fisher Scoring iterations: 5 > fitmacro(glm.hrc1) deviance: 1451.347 df: 12 bic: 1332.558 aic: 1427.347 Number of parameters: 13 Number of cases: 19912 > > > > ### *