## Exercises: survival on the titanic, using loglm()

library(MASS)     # for loglm()
library(vcd)      # for mosaic, aka plot.loglm()

data(Titanic, package="datasets")  # effects::Titanic gives a case-form version

## Baseline model

Titanic <- Titanic + 0.5   # adjust for 0 cells
titanic.mod1 <- loglm(~ (Class * Age * Sex) + Survived, data=Titanic)
titanic.mod1
## Call:
## loglm(formula = ~(Class * Age * Sex) + Survived, data = Titanic)
##
## Statistics:
##                       X^2 df P(> X^2)
## Likelihood Ratio 659.3162 15        0
## Pearson          643.1600 15        0
plot(titanic.mod1, main="Model [AGC][S]")

## Associations with survival

titanic.mod2 <- loglm(~ (Class * Age * Sex) + Survived*(Class + Age + Sex), data=Titanic)
titanic.mod2
## Call:
## loglm(formula = ~(Class * Age * Sex) + Survived * (Class + Age +
##     Sex), data = Titanic)
##
## Statistics:
##                       X^2 df P(> X^2)
## Likelihood Ratio 105.5531 10        0
## Pearson          101.4396 10        0
plot(titanic.mod2,  main="Model [AGC][AS][GS][CS]")

## update() is a simpler way:

update(titanic.mod1, . ~ .+ Survived*(Class+Age+Sex))
## Call:
## loglm(formula = . ~ Class + Age + Sex + Survived + Class:Age +
##     Class:Sex + Age:Sex + Class:Survived + Age:Survived + Sex:Survived +
##     Class:Age:Sex, data = Titanic)
##
## Statistics:
##                       X^2 df P(> X^2)
## Likelihood Ratio 105.5531 10        0
## Pearson          101.4396 10        0
titanic.mod3 <- loglm(~ (Class * Age * Sex) + Survived*(Class + Age * Sex), data=Titanic)
titanic.mod3
## Call:
## loglm(formula = ~(Class * Age * Sex) + Survived * (Class + Age *
##     Sex), data = Titanic)
##
## Statistics:
##                       X^2 df     P(> X^2)
## Likelihood Ratio 85.51508  9 1.287859e-14
## Pearson          78.83273  9 2.755574e-13
plot(titanic.mod3, main="Model [AGC][AS][GS][CS][AGS]")

## compare models

anova(titanic.mod1, titanic.mod2, titanic.mod3, test="chisq")
## LR tests for hierarchical log-linear models
##
## Model 1:
##  ~(Class * Age * Sex) + Survived
## Model 2:
##  ~(Class * Age * Sex) + Survived * (Class + Age + Sex)
## Model 3:
##  ~(Class * Age * Sex) + Survived * (Class + Age * Sex)
##
##            Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   659.31623 15
## Model 2   105.55308 10  553.76314         5          0e+00
## Model 3    85.51508  9   20.03800         1          1e-05
## Saturated   0.00000  0   85.51508         9          0e+00
IycgLS0tDQojJyB0aXRsZTogIlN1cnZpdmFsIG9uIHRoZSB0aXRhbmljLCB1c2luZyBsb2dsbSgpIg0KIycgYXV0aG9yOiAiTWljaGFlbCBGcmllbmRseSINCiMnIGRhdGU6ICJgciBmb3JtYXQoU3lzLkRhdGUoKSlgIg0KIycgb3V0cHV0Og0KIycgICBodG1sX2RvY3VtZW50Og0KIycgICAgIHRoZW1lOiByZWFkYWJsZQ0KIycgICAgIGNvZGVfZG93bmxvYWQ6IHRydWUNCiMnIC0tLQ0KDQoNCiMnICMjIEV4ZXJjaXNlczogc3Vydml2YWwgb24gdGhlIHRpdGFuaWMsIHVzaW5nIGxvZ2xtKCkNCg0KbGlicmFyeShNQVNTKSAgICAgIyBmb3IgbG9nbG0oKQ0KbGlicmFyeSh2Y2QpICAgICAgIyBmb3IgbW9zYWljLCBha2EgcGxvdC5sb2dsbSgpDQoNCmRhdGEoVGl0YW5pYywgcGFja2FnZT0iZGF0YXNldHMiKSAgIyBlZmZlY3RzOjpUaXRhbmljIGdpdmVzIGEgY2FzZS1mb3JtIHZlcnNpb24NCg0KIycgIyMgQmFzZWxpbmUgbW9kZWwgDQpUaXRhbmljIDwtIFRpdGFuaWMgKyAwLjUgICAjIGFkanVzdCBmb3IgMCBjZWxscw0KdGl0YW5pYy5tb2QxIDwtIGxvZ2xtKH4gKENsYXNzICogQWdlICogU2V4KSArIFN1cnZpdmVkLCBkYXRhPVRpdGFuaWMpDQp0aXRhbmljLm1vZDENCnBsb3QodGl0YW5pYy5tb2QxLCBtYWluPSJNb2RlbCBbQUdDXVtTXSIpDQoNCiMnICMjIEFzc29jaWF0aW9ucyB3aXRoIHN1cnZpdmFsIA0KdGl0YW5pYy5tb2QyIDwtIGxvZ2xtKH4gKENsYXNzICogQWdlICogU2V4KSArIFN1cnZpdmVkKihDbGFzcyArIEFnZSArIFNleCksIGRhdGE9VGl0YW5pYykNCnRpdGFuaWMubW9kMg0KcGxvdCh0aXRhbmljLm1vZDIsICBtYWluPSJNb2RlbCBbQUdDXVtBU11bR1NdW0NTXSIpDQoNCiMnICMjIHVwZGF0ZSgpIGlzIGEgc2ltcGxlciB3YXk6DQp1cGRhdGUodGl0YW5pYy5tb2QxLCAuIH4gLisgU3Vydml2ZWQqKENsYXNzK0FnZStTZXgpKQ0KDQp0aXRhbmljLm1vZDMgPC0gbG9nbG0ofiAoQ2xhc3MgKiBBZ2UgKiBTZXgpICsgU3Vydml2ZWQqKENsYXNzICsgQWdlICogU2V4KSwgZGF0YT1UaXRhbmljKQ0KdGl0YW5pYy5tb2QzDQpwbG90KHRpdGFuaWMubW9kMywgbWFpbj0iTW9kZWwgW0FHQ11bQVNdW0dTXVtDU11bQUdTXSIpDQoNCiMnICMjIGNvbXBhcmUgbW9kZWxzDQphbm92YSh0aXRhbmljLm1vZDEsIHRpdGFuaWMubW9kMiwgdGl0YW5pYy5tb2QzLCB0ZXN0PSJjaGlzcSIpDQoNCg0K