Michael Friendly
2010-Mar-01 13:19 UTC
[R] MASS::loglm - exploring a collection of models with add1, drop1
I'd like to fit and explore a collection of hierarchical loglinear models that might range from the independence model, ~ 1 + 2 + 3 + 4 to the saturated model, ~ 1 * 2 * 3 * 4 I can use add1 starting with a baseline model or drop1 starting with the saturated model, but I can't see how to get the model formulas or terms in each model as a *list* that I can work with further. Consider this 2 x 3 x 7 x 2 table: Hoyt1 <- structure(c(87, 125, 216, 136, 256, 65, 72, 233, 159, 269, 176, 125, 52, 572, 119, 635, 119, 300, 88, 351, 158, 355, 144, 147, 32, 137, 43, 131, 42, 65, 14, 155, 24, 152, 24, 67, 20, 116, 41, 106, 32, 47, 53, 96, 163, 176, 309, 144, 36, 138, 116, 308, 225, 327, 52, 598, 162, 901, 243, 711, 48, 238, 130, 414, 237, 339, 12, 116, 35, 200, 72, 136, 9, 146, 19, 209, 42, 134, 3, 95, 25, 182, 36, 126), .Dim = c(2L, 3L, 7L, 2L), .Dimnames = structure(list( Status = c("College", "Non-College"), Rank = c("Low", "Middle", "High"), Occupation = c("1", "2", "3", "4", "5", "6", "7" ), Sex = c("Male", "Female")), .Names = c("Status", "Rank", "Occupation", "Sex")), class = c("xtabs", "table"), call = xtabs(formula = as.formula(paste("freq ~", paste(tvars, collapse = "+"))), data = table)) > str(Hoyt1) xtabs [1:2, 1:3, 1:7, 1:2] 87 125 216 136 256 65 72 233 159 269 ... - attr(*, "dimnames")=List of 4 ..$ Status : chr [1:2] "College" "Non-College" ..$ Rank : chr [1:3] "Low" "Middle" "High" ..$ Occupation: chr [1:7] "1" "2" "3" "4" ... ..$ Sex : chr [1:2] "Male" "Female" - attr(*, "class")= chr [1:2] "xtabs" "table" - attr(*, "call")= language xtabs(formula = as.formula(paste("freq ~", paste(tvars, collapse = "+"))), data = table) # fit baseline log-linear model for Status as response hoyt.mod0 <- loglm(~ Status + (Sex*Rank*Occupation), data=Hoyt1) > (hoyt.add1 <- add1(hoyt.mod0, ~.^2, test="Chisq")) Single term additions Model: ~Status + (Sex * Rank * Occupation) Df AIC LRT Pr(Chi) <none> 2166.36 Status:Sex 1 2129.54 38.82 4.658e-10 *** Status:Rank 2 1430.03 740.33 < 2.2e-16 *** Status:Occupation 6 841.86 1336.50 < 2.2e-16 *** --- Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 > str(hoyt.add1) Classes ?anova? and 'data.frame': 4 obs. of 4 variables: $ Df : num NA 1 2 6 $ AIC : num 2166 2130 1430 842 $ LRT : num NA 38.8 740.3 1336.5 $ Pr(Chi): num NA 4.66e-10 1.74e-161 1.36e-285 - attr(*, "heading")= chr "Single term additions" "\nModel:" "~Status + (Sex * Rank * Occupation)" So, all I get is a data frame containing the results for added terms. I'm also not sure which add1 function to look at since I don't find an add1.loglm > methods(add1) [1] add1.default* add1.glm* add1.gnm* add1.lm* add1.mlm* add1.multinom* -Michael -- Michael Friendly Email: friendly AT yorku DOT ca Professor, Psychology Dept. York University Voice: 416 736-5115 x66249 Fax: 416 736-5814 4700 Keele Street http://www.math.yorku.ca/SCS/friendly.html Toronto, ONT M3J 1P3 CANADA
Maybe Matching Threads
- unknown dimensions for loglm
- fitting the null loglinear model with MASS::loglm??
- loglm() uses only a reference to data, and not data itsel f - is that on purpose??
- [R] possible bug in loglm (PR#122)
- Replicated LR goodness-of-fit tests, heterogeneity G, with loglm?