Surviving the Titanic with Zelig

Re-implementing computations from slides 5.2-5.18 of Lecture 6 using Zelig 5 Syntax

Previewing, recoding and cleaning up the data set

library(radiant.data)
data("Titanic")
titanic <- titanic %>%
  mutate(surv = as.integer(survived)) %>% 
  mutate(survival = sjmisc::rec(surv, rec = "2=0; 1=1")) %>% 
  select(pclass, survived, survival, everything())
head(titanic)

Interaction Model (Slide 5.2)

library(Zelig)
z5_tit <- zlogit$new()
z5_tit$zelig(survival ~ age + sex*pclass + fare, data = titanic)
summary(z5_tit)

Model: 

Call:
z5_tit$zelig(formula = survival ~ age + sex * pclass + fare, 
    data = titanic)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0634  -0.6641  -0.4943   0.4336   2.4941  

Coefficients:
                    Estimate Std. Error z value            Pr(>|z|)
(Intercept)        4.8978215  0.6131092   7.988 0.00000000000000137
age               -0.0387245  0.0068044  -5.691 0.00000001262563680
sexmale           -3.8996177  0.5015659  -7.775 0.00000000000000755
pclass2nd         -1.5923247  0.6024844  -2.643             0.00822
pclass3rd         -4.1382715  0.5601819  -7.387 0.00000000000014976
fare              -0.0009058  0.0020509  -0.442             0.65874
sexmale:pclass2nd -0.0600076  0.6372949  -0.094             0.92498
sexmale:pclass3rd  2.5019110  0.5479901   4.566 0.00000498035051247

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1409.99  on 1042  degrees of freedom
Residual deviance:  931.45  on 1035  degrees of freedom
AIC: 947.45

Number of Fisher Scoring iterations: 5

Next step: Use 'setx' method

Interaction between sex and class along with age and fare to determine the factors that have an effect on survival during the titanic.

Age Effect (Slides 5.3-5.4)

z5_tit$setrange(age=0:80)
z5_tit$sim()
z5_tit$graph()

This plot shows the effect age (independent variable) has on surviving (dependent variable) in the titanic. As we can see elderly people were less likely to survive. As age increases, the probability of surviving the titanic decreases.

Fare Effect (Slides 5.5-5.6)

z5_tit$setrange(fare=0:513)
z5_tit$sim()
z5_tit$graph()

This graph shows the effect fare has on survival. According to this graph, fare does not have much of an affect on survival because as fare increases, the probability of surviving is almost the same for everyone. In other words surviving does not depend on whether the passenger paid a low or high fare during the titanic.

Sex Difference and viewing only the simulated first difference (Slides 5.7-5.9)

z5_tit$setx(sex="male")
z5_tit$setx1(sex="female")
z5_tit$sim()
fd <- z5_tit$get_qi(xvalue="x1", qi="fd")
summary(fd)

       V1         
 Min.   :0.09044  
 1st Qu.:0.22501  
 Median :0.25545  
 Mean   :0.25504  
 3rd Qu.:0.28319  
 Max.   :0.42953  

plot(z5_tit)

This plot shows the simulated first difference between males and females in the titanic. The predicted and expected values for both are shown.

Testing Class Variations in Sex Difference Using Zelig (Slides 5.11-5.14)

First Class

z5t.1 <- zlogit$new()
z5t.1$zelig(survival ~ age +sex*pclass + fare, data=titanic)
z5t.1$setx(sex="male", pclass="1st")
z5t.1$setx1(sex="female", pclass="1st")
z5t.1$sim()
plot(z5t.1)

2nd Class

z5t.2 <- zlogit$new()
z5t.2$zelig(survival ~ age +sex*pclass + fare, data=titanic)
z5t.2$setx(sex="male", pclass="2nd")
z5t.2$setx1(sex="female", pclass="2nd")
z5t.2$sim()
plot(z5t.2)

3rd Class

z5t.3 <- zlogit$new()
z5t.3$zelig(survival ~ age +sex*pclass + fare, data=titanic)
z5t.3$setx(sex="male", pclass="3rd")
z5t.3$setx1(sex="female", pclass="3rd")
z5t.3$sim()
plot(z5t.3)

Putting them into one place (Slides 5.15-5.16)

d1 <- z5t.1$get_qi(xvalue="x1", qi="fd")
d2 <- z5t.2$get_qi(xvalue="x1", qi="fd")
d3 <- z5t.3$get_qi(xvalue="x1", qi="fd")
  
dfd <- as.data.frame(cbind(d1, d2, d3))
head(dfd)

tidd <- dfd %>% 
  gather(class, simv, 1:3)
head(tidd)

Group and plot by class (Slides 5.17-5.18)

tidd %>% 
  group_by(class) %>% 
  summarise(mean = mean(simv), sd = sd(simv))

ggplot(tidd, aes(simv)) + geom_histogram() + facet_grid(~class)