Importing Necessary Packages

library(Zelig)
library(texreg)
library(mvtnorm)
library(car)
library(sjmisc)
library(radiant.data)
data(titanic)
head(titanic)
## # A tibble: 6 x 10
##   pclass survived sex       age sibsp parch  fare name       cabin embarked
##   <fct>  <fct>    <fct>   <dbl> <int> <int> <dbl> <chr>      <chr> <fct>   
## 1 1st    Yes      female 29         0     0 211.  Allen, Mi~ B5    Southam~
## 2 1st    Yes      male    0.917     1     2 152.  Allison, ~ C22 ~ Southam~
## 3 1st    No       female  2         1     2 152.  Allison, ~ C22 ~ Southam~
## 4 1st    No       male   30         1     2 152.  Allison, ~ C22 ~ Southam~
## 5 1st    No       female 25         1     2 152.  Allison, ~ C22 ~ Southam~
## 6 1st    Yes      male   48         0     0  26.5 Anderson,~ E12   Southam~

Creating the Binomial Variable Survival from Survived

titanic <- titanic %>%
  mutate(survint = as.integer(survived)) %>% 
  mutate(survival = sjmisc::rec(survint, rec = "2=0; 1=1")) %>% 
  select(pclass, survived, survival, everything())
head(titanic)
## # A tibble: 6 x 12
##   pclass survived survival sex      age sibsp parch  fare name  cabin
##   <fct>  <fct>       <dbl> <fct>  <dbl> <int> <int> <dbl> <chr> <chr>
## 1 1st    Yes             1 fema~ 29         0     0 211.  Alle~ B5   
## 2 1st    Yes             1 male   0.917     1     2 152.  Alli~ C22 ~
## 3 1st    No              0 fema~  2         1     2 152.  Alli~ C22 ~
## 4 1st    No              0 male  30         1     2 152.  Alli~ C22 ~
## 5 1st    No              0 fema~ 25         1     2 152.  Alli~ C22 ~
## 6 1st    Yes             1 male  48         0     0  26.5 Ande~ E12  
## # ... with 2 more variables: embarked <fct>, survint <int>

Logit Model

z5 <- zlogit$new()
z5$zelig(survival ~ age + sex*pclass + fare, data = titanic)
summary(z5)
## 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 1.37e-15
## age               -0.0387245  0.0068044  -5.691 1.26e-08
## sexmale           -3.8996177  0.5015659  -7.775 7.55e-15
## pclass2nd         -1.5923247  0.6024844  -2.643  0.00822
## pclass3rd         -4.1382715  0.5601819  -7.387 1.50e-13
## 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 4.98e-06
## 
## (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

z5$setrange(age=min(titanic$age):max(titanic$age))
z5$sim()
z5$graph()

Fare Effect

z5$setrange(fare=min(titanic$fare):max(titanic$fare))
z5$sim()
z5$graph()

Sex Difference

z.outsex <- zlogit$new()
z.outsex$zelig(survival ~ age + sex*pclass + fare,  data = titanic)
z.outsex$setx(sex = "male")
z.outsex$setx1(sex = "female")
z.outsex$sim()
summary(z.outsex)
## 
##  sim x :
##  -----
## ev
##           mean         sd      50%      2.5%   97.5%
## [1,] 0.1396722 0.01971422 0.138052 0.1047187 0.18554
## pv
##          0     1
## [1,] 0.869 0.131
## 
##  sim x1 :
##  -----
## ev
##           mean         sd       50%      2.5%     97.5%
## [1,] 0.3958124 0.04205127 0.3921881 0.3192209 0.4751355
## pv
##          0     1
## [1,] 0.611 0.389
## fd
##           mean         sd       50%      2.5%     97.5%
## [1,] 0.2561401 0.04298339 0.2543503 0.1775146 0.3436489

First Difference: Sex Difference for Survival

fd <- z.outsex$get_qi(xvalue="x1", qi="fd")
summary(fd)
##        V1        
##  Min.   :0.1302  
##  1st Qu.:0.2263  
##  Median :0.2544  
##  Mean   :0.2561  
##  3rd Qu.:0.2847  
##  Max.   :0.4050
z.outsex$graph()

Sex DIfference in Class Difference for Survival

First Class

z5one <- zlogit$new()
z5one$zelig(survival ~ age + sex*pclass + fare,  data = titanic)
z5one$setx(sex = "male", pclass = "1st")
z5one$setx1(sex = "female", pclass = "1st")
z5one$sim()
summary(z5one)
## 
##  sim x :
##  -----
## ev
##           mean         sd       50%      2.5%    97.5%
## [1,] 0.4536915 0.04994006 0.4519512 0.3599637 0.548139
## pv
##          0     1
## [1,] 0.546 0.454
## 
##  sim x1 :
##  -----
## ev
##           mean         sd       50%      2.5%     97.5%
## [1,] 0.9732642 0.01344564 0.9766458 0.9412318 0.9906823
## pv
##          0     1
## [1,] 0.025 0.975
## fd
##           mean        sd       50%      2.5%     97.5%
## [1,] 0.5195727 0.0491683 0.5218508 0.4201201 0.6103604
z5one$graph()

Second Class

z5two <- zlogit$new()
z5two$zelig(survival ~ age + sex*pclass + fare,  data = titanic)
z5two$setx(sex = "male", pclass = "2nd")
z5two$setx1(sex = "female", pclass = "2nd")
z5two$sim()
summary(z5two)
## 
##  sim x :
##  -----
## ev
##           mean         sd       50%       2.5%     97.5%
## [1,] 0.1389105 0.02733794 0.1361191 0.09079693 0.1984747
## pv
##          0     1
## [1,] 0.868 0.132
## 
##  sim x1 :
##  -----
## ev
##           mean         sd       50%      2.5%     97.5%
## [1,] 0.8895154 0.03232589 0.8924987 0.8131954 0.9389934
## pv
##          0     1
## [1,] 0.118 0.882
## fd
##          mean         sd       50%     2.5%     97.5%
## [1,] 0.750605 0.04170076 0.7531786 0.661502 0.8264599
z5two$graph()

Third Class

z5three <- zlogit$new()
z5three$zelig(survival ~ age + sex*pclass + fare,  data = titanic)
z5three$setx(sex = "male", pclass = "3rd")
z5three$setx1(sex = "female", pclass = "3rd")
z5three$sim()
summary(z5three)
## 
##  sim x :
##  -----
## ev
##           mean         sd       50%      2.5%     97.5%
## [1,] 0.1397389 0.01887632 0.1394167 0.1045682 0.1808309
## pv
##          0     1
## [1,] 0.859 0.141
## 
##  sim x1 :
##  -----
## ev
##           mean        sd      50%      2.5%     97.5%
## [1,] 0.3972308 0.0432744 0.397811 0.3151352 0.4840281
## pv
##          0     1
## [1,] 0.631 0.369
## fd
##           mean         sd       50%      2.5%     97.5%
## [1,] 0.2574919 0.04498704 0.2572497 0.1725971 0.3450113
z5three$graph()

Combining the Results in a Single Table

d1<-z5one$get_qi(xvalue="x1", qi="fd")
d2<-z5two$get_qi(xvalue="x1", qi="fd")
d3<-z5three$get_qi(xvalue="x1", qi="fd")
dfd <- as.data.frame(cbind(d1, d2, d3))
head(dfd)
##          V1        V2        V3
## 1 0.5125521 0.7116227 0.2969740
## 2 0.4025412 0.7986892 0.3332708
## 3 0.5846946 0.7759991 0.2724982
## 4 0.4615134 0.7615697 0.3029418
## 5 0.5639483 0.7541183 0.2945763
## 6 0.5506270 0.8470110 0.2760863

Tables and GGPLOTS

library(tidyr)

tidd <- dfd %>% 
  gather(class, simv, 1:3)
head(tidd)
##   class      simv
## 1    V1 0.5125521
## 2    V1 0.4025412
## 3    V1 0.5846946
## 4    V1 0.4615134
## 5    V1 0.5639483
## 6    V1 0.5506270
library(dplyr)
tidd %>% 
  group_by(class) %>% 
  summarise(mean = mean(simv), sd = sd(simv))
## # A tibble: 3 x 3
##   class  mean     sd
##   <chr> <dbl>  <dbl>
## 1 V1    0.520 0.0492
## 2 V2    0.751 0.0417
## 3 V3    0.257 0.0450
library(ggplot2)
ggplot(tidd, aes(simv)) + geom_histogram() + facet_grid(~class)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.