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