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?
>
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
