Kway: Fit All K-way Models in a GLM

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

View source: R/Kway.R

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

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

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

Value

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

Author(s)

Michael Friendly and Heather Turner

See Also

glmlist, Summarise (soon to be deprecated), LRstats

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
## 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)

Example output

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

vcdExtra documentation built on Jan. 25, 2021, 5:08 p.m.