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

[Heather Turner]

Hi Michael,

You need to substitute the value of i for the symbol i in your formula, i.e.

update(mod, substitute(.~.^i, list(i = i)))

So, with some other tidying up:

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

mods <- Kway(Freq ~ A + B + C, data=df, family=poisson)

Best wishes,

Heather

Michael Friendly wrote:
> 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?
>