9

I was doing some log-linear models to test for interactions/associations in multiway contingency tables (based on the tutorial here, http://ww2.coastal.edu/kingw/statistics/R-tutorials/loglin.html). I was doing this using a Poisson glm on the observed frequencies as well as with MASS's loglm. I was just wondering though what type of hypothesis test would make most sense here, sequential type I using anova() (not good since p values there depend on the order of the factors in the model), type III test using Anova() in car (independent of the order of the factors in the model) or using drop1 starting from the most complex model?

E.g. using the Titanic passenger survival data

library(COUNT)
data(titanic)
titanic=droplevels(titanic)
head(titanic)
mytable=xtabs(~class+age+sex+survived, data=titanic)
ftable(mytable)
                       survived  no yes
class     age    sex                   
1st class child  women            0   1
                 man              0   5
          adults women            4 140
                 man            118  57
2nd class child  women            0  13
                 man              0  11
          adults women           13  80
                 man            154  14
3rd class child  women           17  14
                 man             35  13
          adults women           89  76
                 man            387  75
freqdata=data.frame(mytable)
fullmodel=glm(Freq~SITE*SEX*MORTALITY,family=poisson,data=freqdata)

Would the most sensible test for interactions between the different categorical factors then be given by type I SS as in

anova(fullmodel, test="Chisq")
Analysis of Deviance Table

Model: poisson, link: log

Response: Freq

Terms added sequentially (first to last)


                       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                                      23    2173.33              
class                   2   231.18        21    1942.15 < 2.2e-16 ***
age                     1  1072.61        20     869.54 < 2.2e-16 ***
sex                     1   137.74        19     731.80 < 2.2e-16 ***
survived                1    77.61        18     654.19 < 2.2e-16 ***
class:age               2    32.41        16     621.78 9.178e-08 ***
class:sex               2    29.61        14     592.17 3.719e-07 ***
age:sex                 1     6.09        13     586.09   0.01363 *  
class:survived          2   132.69        11     453.40 < 2.2e-16 ***
age:survived            1    25.58        10     427.81 4.237e-07 ***
sex:survived            1   312.93         9     114.89 < 2.2e-16 ***
class:age:sex           2     4.04         7     110.84   0.13250    
class:age:survived      2    35.45         5      75.39 2.002e-08 ***
class:sex:survived      2    73.71         3       1.69 < 2.2e-16 ***
age:sex:survived        1     1.69         2       0.00   0.19421    
class:age:sex:survived  2     0.00         0       0.00   1.00000

or using type III SS using car's Anova :

library(car)
library(afex)
set_sum_contrasts()
Anova(fullmodel, test="LR", type="III")
Analysis of Deviance Table (Type III tests)

Response: Freq
                       LR Chisq Df Pr(>Chisq)    
class                    37.353  2  7.744e-09 ***
age                       5.545  1  0.0185317 *  
sex                       0.000  1  0.9999999    
survived                  1.386  1  0.2390319    
class:age                 5.476  2  0.0646851 .  
class:sex                 0.000  2  1.0000000    
age:sex                   0.000  1  0.9999888    
class:survived           16.983  2  0.0002052 ***
age:survived              0.056  1  0.8126973    
sex:survived              0.000  1  0.9999953    
class:age:sex             0.000  2  1.0000000    
class:age:survived        3.461  2  0.1771673    
class:sex:survived        0.000  2  1.0000000    
age:sex:survived          0.000  1  0.9999905    
class:age:sex:survived    0.000  2  1.0000000

or using single term deletions and LRTs with drop1 :

fullmodel=glm(Freq~class+age+sex+survived+class:age+class:sex+class:survived+age:sex+age:survived+sex:survived, family=poisson, data=freqdata)
drop1(fullmodel,test="Chisq")
Single term deletions

Model:
Freq ~ class + age + sex + survived + class:age + class:sex + 
    class:survived + age:sex + age:survived + sex:survived
               Df Deviance    AIC     LRT  Pr(>Chi)    
<none>              114.89 249.01                      
class:age       2   162.76 292.89  47.877 4.016e-11 ***
class:sex       2   115.74 245.86   0.850    0.6537    
class:survived  2   230.95 361.08 116.067 < 2.2e-16 ***
age:sex         1   114.89 247.02   0.008    0.9294    
age:survived    1   134.39 266.52  19.505 1.003e-05 ***
sex:survived    1   427.81 559.94 312.927 < 2.2e-16 ***

?

[This last result appears to match that of MASS's loglm, as should be the case :

fullmodel=loglm(~class+age+sex+survived+class:age+class:sex+class:survived+age:sex+age:survived+sex:survived, mytable)
stepAIC(fullmodel) 
drop1(fullmodel,test="Chisq")
Single term deletions

Model:
~class + age + sex + survived + class:age + class:sex + class:survived + 
    age:sex + age:survived + sex:survived
               Df    AIC     LRT  Pr(>Chi)    
<none>            144.89                      
class:age       2 188.76  47.877 4.016e-11 ***
class:sex       2 141.74   0.850    0.6537    
class:survived  2 256.95 116.067 < 2.2e-16 ***
age:sex         1 142.89   0.008    0.9294    
age:survived    1 162.39  19.505 1.003e-05 ***
sex:survived    1 455.81 312.927 < 2.2e-16 ***

]

(Any other more elegant ways btw to specify a model with main effects + all first order interaction effects?)

Any thoughts what would be the best way to analyse such multiway contingency tables, and adequately test for associations for unbalanced data sets?

EDIT: based on the answer below I went for the drop1 solution :

fullmodel=glm(Freq~class+age+sex+survived+class:age+class:sex+class:survived+age:sex+age:survived+sex:survived, family=poisson, data=freqdata)
drop1(fullmodel,test="Chisq")

which is equivalent to the log-linear model in MASS :

fullmodel=loglm(~class+age+sex+survived+class:age+class:sex+class:survived+age:sex+age:survived+sex:survived, mytable)
stepAIC(fullmodel) 
drop1(fullmodel,test="Chisq")
Tom Wenseleers
  • 2,413
  • 1
  • 21
  • 39
  • 1
    Don't you need to define a study question to answer this? – DWin Nov 18 '14 at 21:28
  • What do you mean with study question? Well, basically I would just like to know which of the categorical factors significantly interact with one another... – Tom Wenseleers Nov 18 '14 at 21:29
  • So you have no hypothesis and there is no science behind this investigation? Just an exercise in curiosity? – DWin Nov 18 '14 at 21:32
  • Ha well I was using a different dataset myself, and was just using the Titanic passenger survival data as an example here, but basically I would like to know which factors are associated with each other, as in a loglinear analysis - that should be pretty standard right? As in a two-way contingency table, but then multiway... – Tom Wenseleers Nov 18 '14 at 21:38
  • I was directing my question at the data you were actually working with. I suspected it had some value and that there was some goal at improving knowledge in a particular domain where prior research had occurred. I do not believe you can do useful statistics without some effort at "design". Trying to abstract the statistics from the scientific foundation is "just wrong". – DWin Nov 18 '14 at 22:34

2 Answers2

3

If I had to choose based on how you set this up, I guess I would go with Anova(), but neither makes much sense. The order R enters the variables into the model is standardized and arbitrary. I would not use that to define the tests I would run.

Instead, use ?loglin or ?loglm in the MASS library, and then drop the specific variables / combinations that you are interested in testing. There is an R loglm tutorial using the Titanic dataset here.

As @DWin notes, the sequential vs not distinction corresponds to meaningfully different hypotheses. So that cannot be answered except by the researcher. The standard version of this point in the R world is Venables' paper, Exegesis on linear models (pdf). Given that you state you just wonder "which factors are associated with each other", that seems less like a conditional inference and more like dropping the specified association from the full model and testing that, or perhaps testing associations dropped from stripped down models where the other variables aren't included at all.

Community
  • 1
gung - Reinstate Monica
  • 132,789
  • 81
  • 357
  • 650
  • Yes I saw that tutorial many thx for that - I've added the `loglm` solution to my answer now too, but I am still puzzled by what would be the best solution. The `Anova` type III SS solution at least has the advantage that it's not dependent on the order in which the factors are added to the model, which should be much better right? In the tutorial you refer to they use type I SS with `step()` to do backward stepwise model reduction, but this is not guaranteed to find the very best model. And they use `anova()` as well together with `glm`, which is dependent on the order of the factors... – Tom Wenseleers Nov 18 '14 at 21:50
  • So if I would just consider 1st order interactions, would that correspond to `fullmodel=loglm(~class+age+sex+survived+class:age+class:sex+class:survived+age:sex+age:survived+sex:survived, mytable) drop1(fullmodel,test="Chisq")` ? – Tom Wenseleers Nov 18 '14 at 21:54
  • 1
    Which would be the "best" solution depends on the hypothesis you want to test. There is some sense of sequentiality (?) necessary in that I would not drop a 'main effect' while leaving a higher level interaction that includes it in the model. (Thus, if dropping the top interaction yields a significant reduction in fit, I would stop there, eg.) But otherwise, your tests could be conditional only, could be type 'III'-esque, or could be marginal (analogous to pairwise correlations instead of multiple regression). – gung - Reinstate Monica Nov 18 '14 at 21:58
  • Ha OK I see - this kind of sequential model reduction in which higher order interactions are removed first is done by drop1, right, and that's not dependent on the order of the variables in the model right? So I gather that the two most sensible hypothesis tests would probably be with drop1 or type III LRTs then, and not with sequential type I tests at least, right? And would type III tests be marginal, and drop1 be conditional tests? (sorry need to read up on this) – Tom Wenseleers Nov 18 '14 at 22:18
  • It's hard to say, @TomWenseleers. I would have to dig into `drop1()` to figure out what it's doing. Personally, I don't use it; I like to have clearer control over which model / variable I'm testing, so I always fit 2 models & then use `anova()` for the nested model test (nb, this can lead to a lot of models). – gung - Reinstate Monica Nov 18 '14 at 22:22
  • Yes but it's just that in my case I would like to test for all possible interactions, starting from the most complex model, and test if dropping each interaction term would or would not significantly reduce the likelihood of the model. I understand that that is what `drop1` does (basically it just automates repeated calls to `anova()`). I think I'll go with that solution then, as it makes sense to me and matches what you would do manually... Plus I would like to obtain a nice table with p values, rather than have to reconstruct them by hand with many repeated calls to `anova()`... – Tom Wenseleers Nov 18 '14 at 22:40
  • And the `drop1` result also matches the result of the `loglm` combined with `stepAIC`, which is good in terms of showing the equivalence I think... – Tom Wenseleers Nov 18 '14 at 22:48
  • @TomWenseleers, that seems reasonable to me. – gung - Reinstate Monica Nov 18 '14 at 22:58
0

The approach recommended in CAR for this type of analysis is to use Anova Type II tests, which conform to the principle of marginality.

Fit saturated model: 

library(COUNT)
data(titanic)
titanic=droplevels(titanic)
mytable=xtabs(~class+age+sex+survived, data=titanic)
freqdata=data.frame(mytable)
fullmodel=glm(Freq~class*age*sex*survived, family=poisson, data=freqdata)

Then Type II Anova: 

> Anova(fullmodel)
Analysis of Deviance Table (Type II tests)

Response: Freq
                       LR Chisq Df Pr(>Chisq)    
class                    231.18  2  < 2.2e-16 ***
age                     1072.61  1  < 2.2e-16 ***
sex                      137.74  1  < 2.2e-16 ***
survived                  77.61  1  < 2.2e-16 ***
class:age                 41.24  2  1.107e-09 ***
class:sex                  2.30  2   0.316761    
age:sex                    0.27  1   0.604214    
class:survived           114.88  2  < 2.2e-16 ***
age:survived              20.34  1  6.486e-06 ***
sex:survived             318.53  1  < 2.2e-16 ***
class:age:sex              9.78  2   0.007509 ** 
class:age:survived        37.26  2  8.101e-09 ***
class:sex:survived        64.07  2  1.220e-14 ***
age:sex:survived           1.69  1   0.194209    
class:age:sex:survived     0.00  2   1.000000    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

And you work your way around this table from bottom to top. It is clear that the terms class:age:sex:survived and age:sex:survived have large p values, so they can probably be ignored. We refit the model without those terms: 

secondmodel=glm(Freq~class*age*sex*survived - age:sex:survived - class:age:sex:survived, 
                family=poisson, data=freqdata)
Anova(secondmodel)

Yielding: 

> Anova(secondmodel)
Analysis of Deviance Table (Type II tests)

Response: Freq
                   LR Chisq Df Pr(>Chisq)    
class                231.18  2  < 2.2e-16 ***
age                 1072.61  1  < 2.2e-16 ***
sex                  137.74  1  < 2.2e-16 ***
survived              77.61  1  < 2.2e-16 ***
class:age             48.00  2  3.767e-11 ***
class:sex              0.85  2     0.6530    
age:sex                0.27  1     0.6042    
class:survived       115.42  2  < 2.2e-16 ***
age:survived          20.34  1  6.486e-06 ***
sex:survived         318.53  1  < 2.2e-16 ***
class:age:sex         20.27  2  3.971e-05 ***
class:age:survived    44.21  2  2.507e-10 ***
class:sex:survived    73.71  2  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

We can inspect the loglinear fit itself: 

> secondmodel

Call:  glm(formula = Freq ~ class * age * sex * survived - age:sex:survived - 
    class:age:sex:survived, family = poisson, data = freqdata)

Coefficients:
                         (Intercept)                        class2nd class  
                           -25.43983                               1.08137  
                      class3rd class                             ageadults  
                            28.11861                              26.82613  
                              sexman                           survivedyes  
                             5.89242                              25.43983  
            class2nd class:ageadults              class3rd class:ageadults  
                             0.09728                             -24.98930  
               class2nd class:sexman                 class3rd class:sexman  
                            -1.84450                              -4.94864  
                    ageadults:sexman            class2nd class:survivedyes  
                            -2.50803                               1.48358  
          class3rd class:survivedyes                 ageadults:survivedyes  
                           -25.31933                             -21.88449  
                  sexman:survivedyes       class2nd class:ageadults:sexman  
                            -4.28298                               0.93211  
     class3rd class:ageadults:sexman  class2nd class:ageadults:survivedyes  
                             3.00077                              -3.22185  
class3rd class:ageadults:survivedyes     class2nd class:sexman:survivedyes  
                            21.54657                               0.06801  
   class3rd class:sexman:survivedyes  
                             2.89768  

Degrees of Freedom: 23 Total (i.e. Null);  3 Residual
Null Deviance:      2173 
Residual Deviance: 1.685    AIC: 147.8
> 1 - pchisq(1.685, 3)
[1] 0.6402738

We do not reject the null for the lack-of-fit test (fitted model vs saturated model), so we conclude that the model fits the data well. You can then inspect the coefficients from summary(secondmodel).

landroni
  • 1,003
  • 15
  • 30