Kway: Fit All K-way Models in a GLM

Description

Generate and fit all 0-way, 1-way, 2-way, ... k-way terms in a glm.

This function is designed mainly for hierarchical loglinear models (or glms in the poission family), where it is desired to find the highest-order terms necessary to achieve a satisfactory fit.

Using anova on the resulting glmlist object will then give sequential tests of the pooled contributions of all terms of degree \(k+1\) over and above those of degree \(k\).

This function is also intended as an example of a generating function for glmlist objects, to facilitate model comparison, extraction, summary and plotting of model components, etc., perhaps using lapply or similar.

Usage

Kway(formula, family=poisson, data, ..., order = nt, prefix = "kway")

Arguments

formula

a two-sided formula for the 1-way effects in the model. The LHS should be the response, and the RHS should be the first-order terms connected by + signs.

family

a description of the error distribution and link function to be used in the model. This can be a character string naming a family function, a family function or the result of a call to a family function. (See family for details of family functions.)

data

an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which glm is called.

…

Other arguments passed to glm

order

Highest order interaction of the models generated. Defaults to the number of terms in the model formula.

prefix

Prefix used to label the models fit in the glmlist object.

Value

An object of class glmlist, of length order+1 containing the 0-way, 1-way, ... models up to degree order.

Details

With y as the response in the formula, the 0-way (null) model is y ~ 1. The 1-way ("main effects") model is that specified in the formula argument. The k-way model is generated using the formula . ~ .^k. With the default order = nt, the final model is the saturated model.

As presently written, the function requires a two-sided formula with an explicit response on the LHS. For frequency data in table form (e.g., produced by xtabs) you the data argument is coerced to a data.frame, so you should supply the formula in the form Freq ~ ….

Examples

Run this code

# NOT RUN {
## artificial data
factors <- expand.grid(A=factor(1:3), B=factor(1:2), C=factor(1:3), D=factor(1:2))
Freq <- rpois(nrow(factors), lambda=40)
df <- cbind(factors, Freq)

mods3 <- Kway(Freq ~ A + B + C, data=df, family=poisson)
LRstats(mods3)
mods4 <- Kway(Freq ~ A + B + C + D, data=df, family=poisson)
LRstats(mods4)

# JobSatisfaction data
data(JobSatisfaction, package="vcd")
modSat <- Kway(Freq ~ management+supervisor+own, data=JobSatisfaction, 
               family=poisson, prefix="JobSat")
LRstats(modSat)
anova(modSat, test="Chisq")

# Rochdale data: very sparse, in table form
data(Rochdale, package="vcd")
# }
# NOT RUN {
modRoch <- Kway(Freq~EconActive + Age + HusbandEmployed + Child + 
                     Education + HusbandEducation + Asian + HouseholdWorking,
                data=Rochdale, family=poisson)
LRstats(modRoch)
# }

Run the code above in your browser using DataCamp Workspace