Loglinear models in Splus

Mike Meyer

December 19, 1997

Abstract The Splus lm command and related glm commands are very powerful. Most students and instructors know how to use lm and the associated modelling and plotting commands, but many do not know that the glm command is equally powerful and versatile. This handout gives a few examples of fitting loglinear and logistic models in Splus.

1 Overview

Full documentation of the Splus glm command is available in **Chambers and Hastie**, and good examples are in **Venables and Ripley**. In order to avoid being creative, I have taken data examples from the Fienberg Xmas Tree book. The rest of this handout is little more than an annotated session using Splus to do mechanical analyses of some data. I will not even pretend that these analyses are definitive, or even good.

2 Examples

Consider the following data (from Problem 5.3 in Fienberg). In a study of the relationship between car size and accident injuries, accidents were classified according to the type of accident, severity of accident, and whether or not the driver was ejected. The data are presented below.

Accident Type

Collision  Rollover
Severity  Severity
Car Driver Not Not
weight ejected severe Severe severe Severe
Small  No  350  150  60  112
Yes  26  23  19  80
Standard No  1878 1022  148  404
Yes  111  161  22  265

Our mission is to find an appropriate model for these data, but in order to do that we need to;

• Enter the data into Splus and organize it like a table,

• Find a way of looking at the data,

• Coerce Splus into fitting some appropriate loglinear models,

• Look at some diagnostic measures.

2.1 Entering the Data

Assume that data is in a nice simple format, with only the numbers and not the cell labels. This is most common. Here is the data, as stored in the file called glm.dat.

350    150   60  112
26      23   19   80
1878  1022  148  404
111    161   22  265

Just read this as a list of numbers into Splus,

> mydat <- scan("glm.dat")

Now you need set figure out the “design” of the data. Splus has a useful command called fac.design which will generate factorial designs. Here is an example of using it. For documentation use Splus’s help facility.

> fac.design(c(2,2,2,2),factor.names=c("Type","Sev","Eject","Weight"))
    Type  Sev Eject Weight
 1 Type1 Sev1 Ejct1  Wght1
 2 Type2 Sev1 Ejct1  Wght1
 3 Type1 Sev2 Ejct1  Wght1
 4 Type2 Sev2 Ejct1  Wght1
 5 Type1 Sev1 Ejct2  Wght1
 6 Type2 Sev1 Ejct2  Wght1
 7 Type1 Sev2 Ejct2  Wght1
 8 Type2 Sev2 Ejct2  Wght1
 9 Type1 Sev1 Ejct1  Wght2
10 Type2 Sev1 Ejct1  Wght2
11 Type1 Sev2 Ejct1  Wght2
12 Type2 Sev2 Ejct1  Wght2
13 Type1 Sev1 Ejct2  Wght2
14 Type2 Sev1 Ejct2  Wght2
15 Type1 Sev2 Ejct2  Wght2
16 Type2 Sev2 Ejct2  Wght2

As you can see, this generated the correct design. Using the fac.design command makes sure that the classifying variables really are treated as Splus factors.

Now all we need to do is marry the data with the design in a data frame,

> data.frame(Resp=mydat,fac.design(c(2,2,2,2),
+  factor.names=c("Type","Sev","Eject","Weight")))
> mydat
   Resp  Type  Sev Eject Weight
 1  350 Type1 Sev1 Ejct1  Wght1
 2  150 Type2 Sev1 Ejct1  Wght1
 3   60 Type1 Sev2 Ejct1  Wght1
 4  112 Type2 Sev2 Ejct1  Wght1

 5   26 Type1 Sev1 Ejct2  Wght1
 6   23 Type2 Sev1 Ejct2  Wght1
 7   19 Type1 Sev2 Ejct2  Wght1
 8   80 Type2 Sev2 Ejct2  Wght1
 9 1878 Type1 Sev1 Ejct1  Wght2
10 1022 Type2 Sev1 Ejct1  Wght2
11  148 Type1 Sev2 Ejct1  Wght2
12  404 Type2 Sev2 Ejct1  Wght2
13  111 Type1 Sev1 Ejct2  Wght2
14  161 Type2 Sev1 Ejct2  Wght2
15   22 Type1 Sev2 Ejct2  Wght2
16  265 Type2 Sev2 Ejct2  Wght2

2.2 Looking at the data

I wish I knew good ways of doing this, but I don’t.

2.3 Fitting Models

This is the fun part. The Splus glm command allows you to fit all sorts of hierarchical loglinear modes. There is also the older loglin command which is sometimes useful for very large data sets, but otherwise is slower and less versatile than glm. I am assuming that you already know something about the lm command and the modelling language. If you don’t, please see the handout entitled Linear Models in Splus.

We might as well start out by fitting just the null model. This is not very interesting, but provides a useful beginning position.

> mymod <- glm( Resp ~ 1, family=poisson, data=mydat)

This should be read as “The model mymod gets the output from the result of fitting a generalized linear model with the formula Resp ~ 1 (i.e., the null model), using a poisson regression model and the default (log) link, and drawing the data from the data frame called mymod.”

There are now many questions one can ask of the resulting mymod, including questions about the deviance, residuals, fitted and estimated values, etc. For example the deviance (or likelihood ratio goodness of fit statistic, or G^{2}) is

> deviance(mymod)
[1] 7686.358

and a summary of the fitted model is,

> summary(mymod)
Call: glm(formula = Resp ~ 1, family = poisson, data = mydat)
Deviance Residuals:
       Min        1Q   Median         3Q      Max
 -21.46568 -17.92141 -11.1995 -0.9543187 60.93338
Coefficients:
              Value Std. Error  t value

(Intercept) 5.71022   0.404268 14.12484
(Dispersion Parameter for Gaussian family taken to be 789.5429 )
    Null Deviance: 7686.358 on 15 degrees of freedom
Residual Deviance: 7686.358 on 15 degrees of freedom
Number of Fisher Scoring Iterations: 0

Similarly,

> family(mymod)
Poisson family
    link function:   Log: log(mu)
variance function:   Identity: mu

and so on.

We can manually fit alternate models, either by entering a completely new model via glm, or more handily using the update function, like,

> mymod <- update(mymod, . ~ Type + Sev + Eject + Weight)

which should be read as “The new value of the model structure mymod should be updated from the old value by changing the model formula, but leaving the data and class of model (poisson) the same. The new model formula has the same left hand side as the old one, but the right hand side is now Type + Sev + Eject + Weight, which represents the complete independence model.” If the mymod structure represents an interesting model, then on might consider assigning the output from update to a new structure, e.g.,

> mymod2 <- update(mymod, . ~ Type + Sev + Eject + Weight)

The same summaries as for linear models can now be produced, such as, the analysis of deviance table,

> anova(mymod)
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
       Df Deviance Resid. Df Resid. Dev
  NULL                    15   7686.358
  Type  1   32.661        14   7653.697
   Sev  1 1489.432        13   6164.265
 Eject  1 2674.713        12   3489.551
Weight  1 2296.445        11   1193.106

and the standard summary, including parameter estimates,

> summary(mymod)
Call: glm(formula = Resp ~ Type + Sev + Eject + Weight, family = poisson, data =
        mydat)
Deviance Residuals:
       Min        1Q    Median       3Q      Max
 -15.57728 -5.843355 -3.169368 3.511405 19.09357
Coefficients:
                  Value Std. Error     t value
(Intercept)  4.90063958 0.02567420  190.877988
       Type -0.08236347 0.01442982   -5.707864
        Sev -0.60481489 0.01708364  -35.403166
      Eject -0.88176920 0.02031366  -43.407704
     Weight  0.79374577 0.01915400   41.440201
(Dispersion Parameter for Poisson family taken to be 1 )
    Null Deviance: 7686.358 on 15 degrees of freedom
Residual Deviance: 1193.106 on 11 degrees of freedom
Number of Fisher Scoring Iterations: 4
Correlation of Coefficients:
       (Intercept)       Type        Sev      Eject
  Type  0.0462656
   Sev  0.3591861   0.0000000
 Eject  0.5588092   0.0000000  0.0000000
Weight -0.4927811   0.0000000  0.0000000  0.0000000

2.4 Things still to do

There is a lot more to write here, but this is at least a start. All the standard lm things, like stepwise selection, can be done. Most of the lm plots can be done, but not all make complete sense.