[R] Generating model formulas for all k-way terms

[Michael Friendly]

For the vcdExtra package, I'm exploring methods to generate and fit 
collections of glm models,
and handling lists of such model objects, of class "glmlist".  The 
simplest example is fitting all
k-way models, from k=0 (null model) to the model with the highest-order 
interaction.  I'm
having trouble writing a function, Kway (below) to do what is done in 
the example below

 > factors <- expand.grid(A=1:3, B=1:2, C=1:3)
 > Freq <- rpois(nrow(factors), lambda=40)
 > df <- cbind(factors, Freq)
 >
 > mod.0 <- glm(Freq ~ 1, data=df, family=poisson)
 > mod.1 <- update(mod.0, . ~ A + B + C)
 > mod.2 <- update(mod.1, . ~ .^2)
 > mod.3 <- update(mod.1, . ~ .^3)
 >
 > library(vcdExtra)
 > summarize(glmlist(mod.0, mod.1, mod.2, mod.3))
Model Summary:
      LR Chisq Df Pr(>Chisq)      AIC
mod.0  21.3310 17     0.2118 -12.6690
mod.1  21.1390 14     0.0981  -6.8610
mod.2  15.8044 11     0.1485  -6.1956
mod.3  14.7653 10     0.1409  -5.2347

# Generate and fit all 1-way, 2-way, ... k-way terms in a glm

Kway <- function(formula, family, data, ..., order=nt) {
    models <- list()
    mod <- glm(formula, family, data, ...)
    terms <- terms(formula)
    tl <- attr(terms, "term.labels")
    nt <- length(tl)
    models[[1]] <- mod   
    for(i in 2:order) {
        models[[i]] <- update(mod, .~.^i)
    }   
    # null model
    mod0 <- update(mod, .~1)
    models <- c(mod0, models)   
    class(models) <- "glmlist"
    models
}

 > mods <- Kway(Freq ~ A + B + C, data=df, family=poisson)
Error in terms.formula(tmp, simplify = TRUE) : invalid power in formula

I still don't understand how to manipulate formulas in functions.  Can 
someone help?

-- 
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