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