Kway: Fit All K-way Models in a GLM
In vcdExtra: 'vcd' Extensions and Additions

Description Usage Arguments Details Value Author(s) See Also Examples

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.

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

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

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

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

Michael Friendly and Heather Turner

glmlist, Summarise (soon to be deprecated), LRstats

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

## End(Not run)

Loading required package: vcd
Loading required package: grid
Loading required package: gnm
Likelihood summary table:
          AIC    BIC LR Chisq Df Pr(>Chisq)
kway.0 236.43 238.02   35.431 35     0.4479
kway.1 240.10 249.60   29.103 30     0.5122
kway.2 249.35 271.52   22.353 22     0.4390
kway.3 253.47 281.98   18.470 18     0.4251
Likelihood summary table:
          AIC    BIC LR Chisq Df Pr(>Chisq)
kway.0 236.43 238.02   35.431 35     0.4479
kway.1 239.20 250.28   26.194 29     0.6151
kway.2 251.80 283.47   12.798 16     0.6875
kway.3 264.38 315.05    1.379  4     0.8479
kway.4 271.00 328.01    0.000  0     1.0000
Likelihood summary table:
             AIC     BIC LR Chisq Df Pr(>Chisq)    
JobSat.0 260.251 260.330  208.775  7     <2e-16 ***
JobSat.1 175.472 175.790  117.997  4     <2e-16 ***
JobSat.2  63.541  64.097    0.065  1     0.7989    
JobSat.3  65.476  66.111    0.000  0     1.0000    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Deviance Table

Model 1: Freq ~ 1
Model 2: Freq ~ management + supervisor + own
Model 3: Freq ~ management + supervisor + own + management:supervisor + 
    management:own + supervisor:own
Model 4: Freq ~ management + supervisor + own + management:supervisor + 
    management:own + supervisor:own + management:supervisor:own
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)    
1         7    208.775                         
2         4    117.997  3   90.778   <2e-16 ***
3         1      0.065  3  117.932   <2e-16 ***
4         0      0.000  1    0.065   0.7989    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Warning messages:
1: glm.fit: fitted rates numerically 0 occurred 
2: glm.fit: algorithm did not converge 
3: glm.fit: fitted rates numerically 0 occurred 
4: glm.fit: algorithm did not converge 
5: glm.fit: fitted rates numerically 0 occurred 
6: glm.fit: fitted rates numerically 0 occurred 
Likelihood summary table:
           AIC     BIC LR Chisq  Df Pr(>Chisq)    
kway.0 2620.23 2623.78  2332.66 255     <2e-16 ***
kway.1 1155.95 1187.86   852.37 247     <2e-16 ***
kway.2  504.14  635.31   144.56 219          1    
kway.3  534.29  863.99    62.71 163          1    
kway.4  611.58 1189.44     0.00  93          1    
kway.5  723.58 1499.97     0.00  37          1    
kway.6  779.58 1655.24     0.00   9          1    
kway.7  795.58 1699.60     0.00   1          1    
kway.8  797.58 1705.15     0.00   0     <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1