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

[Michael Friendly]

Beautiful! And thanks for the clean-up
-Michael

Heather Turner wrote:
> 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?
>>
>>     


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