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
> > 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 >debug shows the problem is here. Examining that little section separately, it seems that update() and the for loop are not working together properly, but I am not certain how to fix it.> mod0 <- update(mod, .~1) > models <- c(mod0, models) class(models) <- "glmlist" > models > } > >-- Joshua Wiley Senior in Psychology University of California, Riverside http://www.joshuawiley.com/ [[alternative HTML version deleted]]
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?
>