R version 3.1.1 (2014-07-10) -- "Sock it to Me" Copyright (C) 2014 The R Foundation for Statistical Computing Platform: x86_64-unknown-linux-gnu (64-bit) 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. Natural language support but running in an English locale 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 an HTML browser interface to help. Type 'q()' to quit R. > > # data from Kihlberg, Narragon, and Campbell (1964) used in Fienberg (1980, > # The Analysis of Cross-Classified Categorical Data, MIT Press), Table 5-9. > > foo <- read.table("auto.txt", header = TRUE) > names(foo) [1] "count" "type" "severity" "ejected" "weight" > foo count type severity ejected weight 1 350 collision not.severe no small 2 150 collision severe no small 3 60 rollover not.severe no small 4 112 rollover severe no small 5 26 collision not.severe yes small 6 23 collision severe yes small 7 19 rollover not.severe yes small 8 80 rollover severe yes small 9 1878 collision not.severe no standard 10 1022 collision severe no standard 11 148 rollover not.severe no standard 12 404 rollover severe no standard 13 111 collision not.severe yes standard 14 161 collision severe yes standard 15 22 rollover not.severe yes standard 16 265 rollover severe yes standard > > out1 <- glm(count ~ severity + type + ejected + weight, + family = poisson, data = foo) > summary(out1) Call: glm(formula = count ~ severity + type + ejected + weight, family = poisson, data = foo) Deviance Residuals: Min 1Q Median 3Q Max -15.577 -5.843 -3.169 3.511 19.094 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.67584 0.03863 146.927 < 2e-16 *** severitysevere -0.16473 0.02887 -5.705 1.16e-08 *** typerollover -1.20963 0.03420 -35.369 < 2e-16 *** ejectedyes -1.76355 0.04070 -43.325 < 2e-16 *** weightstandard 1.58749 0.03833 41.422 < 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: 7686.4 on 15 degrees of freedom Residual deviance: 1193.1 on 11 degrees of freedom AIC: 1309.9 Number of Fisher Scoring iterations: 5 > out2 <- glm(count ~ (severity + type + ejected + weight)^2, + family = poisson, data = foo) > summary(out2) Call: glm(formula = count ~ (severity + type + ejected + weight)^2, family = poisson, data = foo) Deviance Residuals: 1 2 3 4 5 6 7 8 -0.48389 0.81525 0.34793 -0.30773 0.46555 -0.58786 1.03833 -0.39015 9 10 11 12 13 14 15 16 0.06835 -0.11229 0.29131 -0.14306 0.37153 -0.25326 -1.93167 0.60399 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.88369 0.04970 118.380 < 2e-16 *** severitysevere -0.94036 0.08281 -11.356 < 2e-16 *** typerollover -1.83460 0.09787 -18.744 < 2e-16 *** ejectedyes -2.71830 0.12118 -22.432 < 2e-16 *** weightstandard 1.65270 0.05365 30.807 < 2e-16 *** severitysevere:typerollover 1.63871 0.08281 19.789 < 2e-16 *** severitysevere:ejectedyes 1.03060 0.09891 10.420 < 2e-16 *** severitysevere:weightstandard 0.33701 0.08607 3.915 9.03e-05 *** typerollover:ejectedyes 1.36560 0.09199 14.845 < 2e-16 *** typerollover:weightstandard -0.72862 0.09472 -7.692 1.44e-14 *** ejectedyes:weightstandard -0.14402 0.10992 -1.310 0.19 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 7686.3582 on 15 degrees of freedom Residual deviance: 7.3281 on 5 degrees of freedom AIC: 136.11 Number of Fisher Scoring iterations: 4 > out3 <- glm(count ~ (severity + type + ejected + weight)^3, + family = poisson, data = foo) > summary(out3) Call: glm(formula = count ~ (severity + type + ejected + weight)^3, family = poisson, data = foo) Deviance Residuals: 1 2 3 4 5 6 7 8 0.08607 -0.13081 -0.20574 0.15266 -0.30904 0.34340 0.37980 -0.17857 9 10 11 12 13 14 15 16 -0.03709 0.05032 0.13264 -0.07989 0.15335 -0.12629 -0.33477 0.09897 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.85333 0.05328 109.863 < 2e-16 severitysevere -0.83203 0.09551 -8.712 < 2e-16 typerollover -1.73254 0.13303 -13.024 < 2e-16 ejectedyes -2.53523 0.18223 -13.912 < 2e-16 weightstandard 1.68549 0.05795 29.087 < 2e-16 severitysevere:typerollover 1.41528 0.17351 8.157 3.44e-16 severitysevere:ejectedyes 0.57696 0.24400 2.365 0.0181 severitysevere:weightstandard 0.22116 0.10254 2.157 0.0310 typerollover:ejectedyes 1.27046 0.25261 5.029 4.92e-07 typerollover:weightstandard -0.81999 0.15417 -5.319 1.04e-07 ejectedyes:weightstandard -0.30865 0.19902 -1.551 0.1209 severitysevere:typerollover:ejectedyes 0.38569 0.22837 1.689 0.0912 severitysevere:typerollover:weightstandard 0.21469 0.19383 1.108 0.2680 severitysevere:ejectedyes:weightstandard 0.43032 0.25158 1.710 0.0872 typerollover:ejectedyes:weightstandard -0.25130 0.24125 -1.042 0.2976 (Intercept) *** severitysevere *** typerollover *** ejectedyes *** weightstandard *** severitysevere:typerollover *** severitysevere:ejectedyes * severitysevere:weightstandard * typerollover:ejectedyes *** typerollover:weightstandard *** ejectedyes:weightstandard severitysevere:typerollover:ejectedyes . severitysevere:typerollover:weightstandard severitysevere:ejectedyes:weightstandard . typerollover:ejectedyes:weightstandard --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 7686.35821 on 15 degrees of freedom Residual deviance: 0.66893 on 1 degrees of freedom AIC: 137.45 Number of Fisher Scoring iterations: 3 > anova(out1, out2, out3, test = "Chisq") Analysis of Deviance Table Model 1: count ~ severity + type + ejected + weight Model 2: count ~ (severity + type + ejected + weight)^2 Model 3: count ~ (severity + type + ejected + weight)^3 Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 11 1193.11 2 5 7.33 6 1185.78 <2e-16 *** 3 1 0.67 4 6.66 0.155 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova(out1, out2, out3, test = "LRT") Analysis of Deviance Table Model 1: count ~ severity + type + ejected + weight Model 2: count ~ (severity + type + ejected + weight)^2 Model 3: count ~ (severity + type + ejected + weight)^3 Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 11 1193.11 2 5 7.33 6 1185.78 <2e-16 *** 3 1 0.67 4 6.66 0.155 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > # so 3-way interactions are unnecessary, now how about 2-way? > # from summary(out2) it seems that all 2-way interactions are statistically > # significant except for ejected * weight, so our final model is > > out4 <- glm(count ~ (severity + type + ejected + weight)^2 - ejected : weight, + family = poisson, data = foo) > summary(out4) Call: glm(formula = count ~ (severity + type + ejected + weight)^2 - ejected:weight, family = poisson, data = foo) Deviance Residuals: Min 1Q Median 3Q Max -2.0975 -0.1681 0.1188 0.2972 1.3400 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.88939 0.04945 119.094 < 2e-16 *** severitysevere -0.92454 0.08160 -11.330 < 2e-16 *** typerollover -1.81340 0.09609 -18.872 < 2e-16 *** ejectedyes -2.83714 0.08150 -34.813 < 2e-16 *** weightstandard 1.64585 0.05338 30.835 < 2e-16 *** severitysevere:typerollover 1.63834 0.08284 19.777 < 2e-16 *** severitysevere:ejectedyes 1.02343 0.09872 10.367 < 2e-16 *** severitysevere:weightstandard 0.31941 0.08486 3.764 0.000167 *** typerollover:ejectedyes 1.38157 0.09120 15.148 < 2e-16 *** typerollover:weightstandard -0.75905 0.09167 -8.280 < 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: 7686.36 on 15 degrees of freedom Residual deviance: 9.02 on 6 degrees of freedom AIC: 135.8 Number of Fisher Scoring iterations: 4 > > # now try loglin > > bar <- xtabs(count ~ severity + type + ejected + weight, data = foo) > bar , , ejected = no, weight = small type severity collision rollover not.severe 350 60 severe 150 112 , , ejected = yes, weight = small type severity collision rollover not.severe 26 19 severe 23 80 , , ejected = no, weight = standard type severity collision rollover not.severe 1878 148 severe 1022 404 , , ejected = yes, weight = standard type severity collision rollover not.severe 111 22 severe 161 265 > lout1 <- loglin(bar, list(1, 2, 3, 4)) 2 iterations: deviation 4.547474e-13 > lout2 <- loglin(bar, + list(c(1, 2), c(1, 3), c(1, 4), c(2, 3), c(2, 4), c(3, 4))) 6 iterations: deviation 0.03456946 > lout3 <- loglin(bar, + list(c(1, 2, 3), c(1, 2, 4), c(1, 3, 4), c(2, 3, 4))) 5 iterations: deviation 0.02584502 > names(lout1) [1] "lrt" "pearson" "df" "margin" > > lout1$lrt [1] 1193.106 > lout2$lrt [1] 7.328126 > lout3$lrt [1] 0.6689364 > # appears these are differences from the saturated model > lout1$df [1] 11 > lout2$df [1] 5 > lout3$df [1] 1 > # ditto here > > # LRT > qux <- matrix(NA, 3, 5) > qux[1, 1] <- lout1$df > qux[2, 1] <- lout2$df > qux[3, 1] <- lout3$df > qux[1, 2] <- lout1$lrt > qux[2, 2] <- lout2$lrt > qux[3, 2] <- lout3$lrt > qux[2:3, 3] <- qux[1:2, 1] - qux[2:3, 1] > qux[2:3, 4] <- qux[1:2, 2] - qux[2:3, 2] > qux[2:3, 5] <- pchisq(qux[2:3, 4], df = qux[2:3, 3], lower.tail = FALSE) > colnames(qux) <- c("Resid. Df", "Resid. Dev", "Df", "Deviance", "Pr(>Chi)") > rownames(qux) <- 1:3 > class(qux) <- "anova" > print(qux) Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 11 1193.11 2 5 7.33 6 1185.78 <2e-16 *** 3 1 0.67 4 6.66 0.155 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova(out1, out2, out3, test = "LRT") Analysis of Deviance Table Model 1: count ~ severity + type + ejected + weight Model 2: count ~ (severity + type + ejected + weight)^2 Model 3: count ~ (severity + type + ejected + weight)^3 Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 11 1193.11 2 5 7.33 6 1185.78 <2e-16 *** 3 1 0.67 4 6.66 0.155 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > # Pearson Chi-square > quux <- matrix(NA, 3, 5) > quux[1, 1] <- lout1$df > quux[2, 1] <- lout2$df > quux[3, 1] <- lout3$df > quux[1, 2] <- lout1$pearson > quux[2, 2] <- lout2$pearson > quux[3, 2] <- lout3$pearson > quux[2:3, 3] <- quux[1:2, 1] - quux[2:3, 1] > quux[2:3, 4] <- quux[1:2, 2] - quux[2:3, 2] > quux[2:3, 5] <- pchisq(quux[2:3, 4], df = quux[2:3, 3], lower.tail = FALSE) > colnames(quux) <- c("Resid. Df", "Resid. Dev", "Df", "Deviance", "Pr(>Chi)") > rownames(quux) <- 1:3 > class(quux) <- "anova" > print(quux) Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 11 1601.73 2 5 7.01 6 1594.72 <2e-16 *** 3 1 0.67 4 6.34 0.175 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova(out1, out2, out3, test = "Chisq") Analysis of Deviance Table Model 1: count ~ severity + type + ejected + weight Model 2: count ~ (severity + type + ejected + weight)^2 Model 3: count ~ (severity + type + ejected + weight)^3 Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 11 1193.11 2 5 7.33 6 1185.78 <2e-16 *** 3 1 0.67 4 6.66 0.155 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > # Clearly, in the latter "Chisq" means "LRT". How quaint! > > # The help for loglin says the loglm function in the MASS package > # is more user-friendly. Let's try it. > library(MASS) > mout1 <- loglm(count ~ severity + type + ejected + weight, data = foo) > summary(mout1) Formula: count ~ severity + type + ejected + weight attr(,"variables") list(count, severity, type, ejected, weight) attr(,"factors") severity type ejected weight count 0 0 0 0 severity 1 0 0 0 type 0 1 0 0 ejected 0 0 1 0 weight 0 0 0 1 attr(,"term.labels") [1] "severity" "type" "ejected" "weight" attr(,"order") [1] 1 1 1 1 attr(,"intercept") [1] 1 attr(,"response") [1] 1 attr(,".Environment") attr(,"predvars") list(count, severity, type, ejected, weight) attr(,"dataClasses") count severity type ejected weight "numeric" "factor" "factor" "factor" "factor" Statistics: X^2 df P(> X^2) Likelihood Ratio 1193.106 11 0 Pearson 1601.729 11 0 > tout1 <- loglm( ~ severity + type + ejected + weight, data = bar) > summary(tout1) Formula: ~severity + type + ejected + weight attr(,"variables") list(severity, type, ejected, weight) attr(,"factors") severity type ejected weight severity 1 0 0 0 type 0 1 0 0 ejected 0 0 1 0 weight 0 0 0 1 attr(,"term.labels") [1] "severity" "type" "ejected" "weight" attr(,"order") [1] 1 1 1 1 attr(,"intercept") [1] 1 attr(,"response") [1] 0 attr(,".Environment") Statistics: X^2 df P(> X^2) Likelihood Ratio 1193.106 11 0 Pearson 1601.729 11 0 > ##### seems to work either way !!!!! > mout2 <- loglm(count ~ (severity + type + ejected + weight)^2, data = foo) > mout3 <- loglm(count ~ (severity + type + ejected + weight)^3, data = foo) > tout2 <- loglm( ~ (severity + type + ejected + weight)^2, data = bar) > tout3 <- loglm( ~ (severity + type + ejected + weight)^3, data = bar) > #### tests > anova(tout1, tout2, tout3, test = "LR") LR tests for hierarchical log-linear models Model 1: ~severity + type + ejected + weight Model 2: ~(severity + type + ejected + weight)^2 Model 3: ~(severity + type + ejected + weight)^3 Deviance df Delta(Dev) Delta(df) P(> Delta(Dev) Model 1 1193.1057775 11 Model 2 7.3281264 5 1185.7776511 6 0.00000 Model 3 0.6689364 1 6.6591900 4 0.15503 Saturated 0.0000000 0 0.6689364 1 0.41342 > print(qux) Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 11 1193.11 2 5 7.33 6 1185.78 <2e-16 *** 3 1 0.67 4 6.66 0.155 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova(tout1, tout2, tout3, test = "Chisq") LR tests for hierarchical log-linear models Model 1: ~severity + type + ejected + weight Model 2: ~(severity + type + ejected + weight)^2 Model 3: ~(severity + type + ejected + weight)^3 Deviance df Delta(Dev) Delta(df) P(> Delta(Dev) Model 1 1193.1057775 11 Model 2 7.3281264 5 1185.7776511 6 0.00000 Model 3 0.6689364 1 6.6591900 4 0.15503 Saturated 0.0000000 0 0.6689364 1 0.41342 > print(quux) Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 11 1601.73 2 5 7.01 6 1594.72 <2e-16 *** 3 1 0.67 4 6.34 0.175 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > ##### a bunch of not so user-friendly stuff about loglm > ##### > ##### * no help page for anova.loglm (although function exists) > ##### * summary.loglm leaves a lot to be desired > ##### * both loglm and glm heavily biased against Pearson test statistic > ##### If you want it, you have to do it yourself using loglin > ##### (Well, not quite. We saw that loglm calculates the Pearson > ##### statistic. It's anova.loglm that does not use it no matter > ##### what arguments are specified.) > > > ### the book Graphical Models with R by Hojsgaard, Edwards, and Lauritzen > ### recommends the function dmod in the R package gRim (this package is > ### a pain to install, since one must install the Bioconductor package > ### RGBL and the R function install.packages doesn't handle that > > library(gRim) Loading required package: gRbase Loading required package: gRain > > gout1 <- dmod( ~ severity + type + ejected + weight, data = bar) > > #### according for help for dmod can also do > > gout1 <- dmod(~.^., interactions = 1, data = bar) > gout1 Model: A dModel with 4 variables graphical : TRUE decomposable : TRUE -2logL : 20295.50 mdim : 4 aic : 20303.50 ideviance : -0.00 idf : 0 bic : 20329.43 deviance : 1193.11 df : 11 > > gout2 <- dmod(~.^., interactions = 2, data = bar) > gout2 Model: A dModel with 4 variables graphical : FALSE decomposable : FALSE -2logL : 19109.72 mdim : 10 aic : 19129.72 ideviance : 1185.78 idf : 6 bic : 19194.55 deviance : 7.33 df : 5 > > gout3 <- dmod(~.^., interactions = 3, data = bar) > gout3 Model: A dModel with 4 variables graphical : FALSE decomposable : FALSE -2logL : 19103.06 mdim : 14 aic : 19131.06 ideviance : 1192.44 idf : 10 bic : 19221.82 deviance : 0.67 df : 1 > > # there is no anova function in this package > # how annoying > # have to do LRT by hand > > comparemodels <- function(...) { + foo <- list(...) + bar <- sapply(foo, function(x) inherits(x, "dModel")) + if (! all(bar)) + stop("all arguments must be objects of class \"dModel\"") + if (length(foo) < 2) + stop("must have at least 2 arguments") + dev <- sapply(foo, function(x) x$fitinfo$dev) + df <- sapply(foo, function(x) x$fitinfo$dimension[1]) + if (! all(diff(df) > 0)) + stop("dimensions of models not increasing; models cannot be nested") + qux <- cbind(df, dev) + qux <- cbind(qux, c(NA, diff(df))) + qux <- cbind(qux, c(NA, - diff(dev))) + qux <- cbind(qux, pchisq(qux[ , 4], df = qux[ , 3], lower.tail = FALSE)) + colnames(qux) <- c("Resid. Df", "Resid. Dev", "Df", "Deviance", + "Pr(>Chi)") + rownames(qux) <- 1:nrow(qux) + class(qux) <- "anova" + return(qux) + } > > comparemodels(gout1, gout2, gout3) Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 4 1193.11 2 10 7.33 6 1185.78 <2e-16 *** 3 14 0.67 4 6.66 0.155 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > # as before we find no point in 3-way interactions > > # no how about graphical models ? > > bout <- backward(gout2) . BACKWARD: type=unrestricted search=all, criterion=aic(2.00), alpha=0.00 . Initial model: is graphical=FALSE is decomposable=FALSE change.AIC -2.1308 Edge deleted: ejected weight > bout Model: A dModel with 4 variables graphical : FALSE decomposable : FALSE -2logL : 19111.41 mdim : 9 aic : 19129.41 ideviance : 1184.09 idf : 5 bic : 19187.76 deviance : 9.02 df : 6 > bout$fitinfo$margin [[1]] [1] "severity" "type" [[2]] [1] "severity" "ejected" [[3]] [1] "severity" "weight" [[4]] [1] "type" "ejected" [[5]] [1] "type" "weight" > > fout <- forward(gout1) . FORWARD: type=decomposable search=all, criterion=aic(2.00), alpha=0.00 . Initial model: is graphical=TRUE is decomposable=TRUE change.AIC -599.4196 Edge added: severity type change.AIC -399.6905 Edge added: ejected type change.AIC -111.3534 Edge added: severity ejected change.AIC -50.9642 Edge added: weight type change.AIC -11.8090 Edge added: severity weight > fout Model: A dModel with 4 variables graphical : TRUE decomposable : TRUE -2logL : 19108.26 mdim : 11 aic : 19130.26 ideviance : 1187.24 idf : 7 bic : 19201.57 deviance : 5.87 df : 4 > fout$fitinfo$margin [[1]] [1] "ejected" "severity" "type" [[2]] [1] "severity" "weight" "type" > > #### clearly the forward function disagrees with our analysis up to this point > > gout4 <- dmod(~.^., data = bar) > gout4 Model: A dModel with 4 variables graphical : TRUE decomposable : TRUE -2logL : 19102.39 mdim : 15 aic : 19132.39 ideviance : 1193.11 idf : 11 bic : 19229.64 deviance : 0.00 df : 0 > > sout <- stepwise(gout4) > sout Model: A dModel with 4 variables graphical : TRUE decomposable : TRUE -2logL : 19108.26 mdim : 11 aic : 19130.26 ideviance : 1187.24 idf : 7 bic : 19201.57 deviance : 5.87 df : 4 > sout$fitinfo$margin [[1]] [1] "severity" "type" "ejected" [[2]] [1] "severity" "type" "weight" > > #### same as the model found by the forward function > #### hopefully, this is the best model by AIC > #### apparently does not have an all decomposible models function > > soutu <- stepwise(gout4, type = "unrestricted") > soutu Model: A dModel with 4 variables graphical : TRUE decomposable : TRUE -2logL : 19108.26 mdim : 11 aic : 19130.26 ideviance : 1187.24 idf : 7 bic : 19201.57 deviance : 5.87 df : 4 > soutu$fitinfo$margin [[1]] [1] "severity" "type" "ejected" [[2]] [1] "severity" "type" "weight" > > #### same model again > #### according to this analysis, the best hierarchical model is decomposible > > > proc.time() user system elapsed 2.486 0.090 3.292