Homework 6.1 - Professor Song’s Zelig 4 Model in Zelig 5

Titanic Survival: Zelig 4 vs. Zelig 5

The model simulates and extracts quantities of interest by estimating effects of age, fare, sex and class on survival probablity using the Titanic data from the radiant.data package. The class slides used a Zelig 4 syntax, the following will show the the syntax of Zelig 4 followed by Zelig 5.

Loading Packages

library(Zelig)
library(texreg)
library(car)
library(mvtnorm)

Loading Data

library(radiant.data)
data(titanic)

Creating Variable - Survival

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

Model - Surviving the titanic with Zelig 4 syntax

z_tit <- zelig(survival ~ age + sex*pclass + fare, model = "logit", data = titanic, cite = F)
summary(z_tit)

Model: 

Call:
z5$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

Model - Zelig5

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

Model: 

Call:
z_tit5$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

Age effect

a.range = min(titanic$age):max(titanic$age)
x <- setx(z_tit5, age = a.range)
s <- sim(z_tit5, x = x)
ci.plot(s)

Age effect in Zelig5

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

Model: 

Call:
z5.age$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

z5.age$setrange(age = min(titanic$age):max(titanic$age))
z5.age$sim()
ci.plot(z5.age)

Fare effect

f.range = min(titanic$fare):max(titanic$fare)
x <- setx(z_tit, fare = f.range)
s <- sim(z_tit, x = x)
ci.plot(s)

Fare effect in Zelig5

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

Model: 

Call:
z5.fare$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

z5.fare$setrange(fare = min(titanic$fare):max(titanic$fare))
z5.fare$sim()
ci.plot(z5.fare)

Sex effect

x <- setx(z_tit, sex = "male")
x1 <- setx(z_tit, sex = "female")
s <- sim(z_tit, x = x, x1 = x1)
summary(s)

 sim x :
 -----
ev
          mean         sd       50%      2.5%     97.5%
[1,] 0.1390279 0.01975385 0.1379809 0.1054139 0.1787193
pv
         0     1
[1,] 0.865 0.135

 sim x1 :
 -----
ev
         mean         sd       50%      2.5%     97.5%
[1,] 0.396558 0.04308903 0.3949084 0.3112113 0.4778312
pv
         0     1
[1,] 0.578 0.422
fd
          mean         sd       50%      2.5%     97.5%
[1,] 0.2575301 0.04390234 0.2597158 0.1735156 0.3418706

Sex effect in Zelig5

z_tit5$setx(sex = "male")
z_tit5$setx1(sex = "female")
z_tit5$sim()
summary(z_tit5)

 sim x :
 -----
ev
         mean         sd       50%      2.5%     97.5%
[1,] 0.141008 0.01892692 0.1401055 0.1065529 0.1825971
pv
         0     1
[1,] 0.847 0.153

 sim x1 :
 -----
ev
          mean         sd       50%      2.5%     97.5%
[1,] 0.3948343 0.04429635 0.3940561 0.3077843 0.4820824
pv
         0     1
[1,] 0.605 0.395
fd
          mean        sd       50%      2.5%     97.5%
[1,] 0.2538263 0.0448389 0.2531833 0.1695305 0.3436792

First difference

fd <- s$get_qi(xvalue="x1", qi="fd")
summary(fd)

       V1        
 Min.   :0.1329  
 1st Qu.:0.2235  
 Median :0.2597  
 Mean   :0.2575  
 3rd Qu.:0.2885  
 Max.   :0.3983  

First difference in Zelig5

fd <- z_tit5$get_qi(xvalue="x1", qi="fd")
summary(fd)

       V1        
 Min.   :0.1131  
 1st Qu.:0.2246  
 Median :0.2532  
 Mean   :0.2538  
 3rd Qu.:0.2811  
 Max.   :0.4029  

Plot- Zelig 5

plot(z_tit5)

How to test the class variations in sex difference using Zelig?

First Class

c1x <- setx(z_tit, sex = "male", pclass = "1st")
c1x1 <- setx(z_tit, sex = "female", pclass = "1st")
c1s <- sim(z_tit, x = c1x, x1 = c1x1)

First Class in Zelig 5

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

 sim x :
 -----
ev
          mean         sd       50%      2.5%     97.5%
[1,] 0.4533515 0.05137299 0.4537543 0.3558283 0.5509001
pv
         0     1
[1,] 0.561 0.439

 sim x1 :
 -----
ev
          mean         sd       50%      2.5%     97.5%
[1,] 0.9730933 0.01344537 0.9760756 0.9402655 0.9909094
pv
        0    1
[1,] 0.03 0.97
fd
          mean        sd       50%      2.5%     97.5%
[1,] 0.5197418 0.0509812 0.5214277 0.4188457 0.6179676

Second Class

c2x <- setx(z_tit, sex = "male", pclass = "2nd")
c2x1 <- setx(z_tit, sex = "female", pclass = "2nd")
c2s <- sim(z_tit, x = c2x, x1 = c2x1)

Second Class in Zelig 5

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

 sim x :
 -----
ev
          mean        sd       50%       2.5%    97.5%
[1,] 0.1383051 0.0279243 0.1360758 0.08944321 0.197538
pv
        0    1
[1,] 0.87 0.13

 sim x1 :
 -----
ev
          mean         sd       50%      2.5%     97.5%
[1,] 0.8902234 0.03228658 0.8948994 0.8166275 0.9403089
pv
         0     1
[1,] 0.087 0.913
fd
          mean         sd      50%      2.5%     97.5%
[1,] 0.7519182 0.04177891 0.754372 0.6662137 0.8210498

Third Class

c3x <- setx(z_tit, sex = "male", pclass = "3rd")
c3x1 <- setx(z_tit, sex = "female", pclass = "3rd")
c3s <- sim(z_tit, x = c3x, x1 = c3x1)

Third Class in Zelig

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

 sim x :
 -----
ev
          mean         sd       50%      2.5%     97.5%
[1,] 0.1399685 0.01972803 0.1392533 0.1029011 0.1858248
pv
         0     1
[1,] 0.859 0.141

 sim x1 :
 -----
ev
          mean         sd       50%      2.5%    97.5%
[1,] 0.3956632 0.04307623 0.3945303 0.3145296 0.483128
pv
         0     1
[1,] 0.592 0.408
fd
          mean         sd       50%      2.5%     97.5%
[1,] 0.2556948 0.04355542 0.2564176 0.1668021 0.3404756

plot(z_firstclass)

plot(z_secondclass)

plot(z_thirdclass)

Putting them together - in Zelig 5 syntax

d1 <- z_firstclass$get_qi(xvalue="x1", qi="fd")
d2 <- z_secondclass$get_qi(xvalue="x1", qi="fd")
d3 <- z_thirdclass$get_qi(xvalue="x1", qi="fd")
dfd <- as.data.frame(cbind(d1, d2, d3))
head(dfd)

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

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

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