[R] MASS::loglm - exploring a collection of models with add1, drop1

[Michael Friendly]

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