Using Zelig 5 Syntax vs Zelig 4 Syntax to analyze titanic dataset.

library(Zelig)
## Warning: package 'Zelig' was built under R version 3.5.3
## Loading required package: survival
## Warning: package 'survival' was built under R version 3.5.3
library(texreg)
## Warning: package 'texreg' was built under R version 3.5.3
## Version:  1.36.23
## Date:     2017-03-03
## Author:   Philip Leifeld (University of Glasgow)
## 
## Please cite the JSS article in your publications -- see citation("texreg").
library (radiant.data)
## Warning: package 'radiant.data' was built under R version 3.5.3
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:texreg':
## 
##     extract
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:Zelig':
## 
##     stat
## Loading required package: lubridate
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
## Loading required package: tidyr
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:magrittr':
## 
##     extract
## The following object is masked from 'package:texreg':
## 
##     extract
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## 
## Attaching package: 'radiant.data'
## The following objects are masked from 'package:lubridate':
## 
##     month, wday
## The following object is masked from 'package:ggplot2':
## 
##     diamonds
data (titanic)
head (titanic)
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)

Creating Interaction 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

Slides: Age Affect

There is a negative correlation between expected value of chances of survival and age. Younger is the age, higher is the chance of survival, however, there is a much higher uncertainly of sirvival for infants. It is much more certain that people over 80 years old had much less chances of survival.

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

Effect of The Ticket Fair: There is no significant correlation between the ticket fair and expected chances of survival. Moreover, it is slightly less for people of first class. However, the uncertainly of survival is increasing drastically with ticket fair price, and extremely varies with first class passengers.

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

Slide Effect of Gender on Survival

Females has 0.26 higher chance of survival than males on average.

z5$setx(sex="male")
z5$setx1(sex="female")
z5$sim()
fd<-z5$get_qi(xvalue="x1", qi="fd")
summary(fd)
##        V1        
##  Min.   :0.1368  
##  1st Qu.:0.2254  
##  Median :0.2545  
##  Mean   :0.2542  
##  3rd Qu.:0.2822  
##  Max.   :0.3859
z5$graph()

Class Variation in Gender Difference

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

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

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

Third class women have a .26 higher chance of survival than third class men.

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

Overall, women hav higher chances of survival. And the highest chances between gender differences of suvival is in the second class, where women has 0.769 higher chance os survival. With class, the chances of survival decreases.

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)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Now, I will use dataset from kaggle which explores Student perfomance. The variables are the following:

perform <- read.csv ("C:/Users/Marcy/Documents/soc 712/StudentsPerformance.csv")
head (perform)
library(readr)
library(dplyr)
library(magrittr)
perform=perform %>%
 mutate(test.preparation.course = sjmisc::rec(test.preparation.course, rec = "none=0; completed=1")) %>% 
  select( gender, race.ethnicity, parental.level.of.education, lunch, test.preparation.course, math.score, reading.score, everything ())
head(perform)
M1<-zlogit$new()
M1$zelig(test.preparation.course ~ gender + parental.level.of.education, data=perform)
M2<-zlogit$new()
M2$zelig(test.preparation.course ~ parental.level.of.education *lunch + gender, data=perform)
M3<-zlogit$new()
M3$zelig(test.preparation.course ~ parental.level.of.education *lunch + math.score + writing.score + gender, data=perform)
texreg::htmlreg(list(M1, M2, M3),doctype = FALSE)
## 
## <table cellspacing="0" align="center" style="border: none;">
## <caption align="bottom" style="margin-top:0.3em;">Statistical models</caption>
## <tr>
## <th style="text-align: left; border-top: 2px solid black; border-bottom: 1px solid black; padding-right: 12px;"><b></b></th>
## <th style="text-align: left; border-top: 2px solid black; border-bottom: 1px solid black; padding-right: 12px;"><b>Model 1</b></th>
## <th style="text-align: left; border-top: 2px solid black; border-bottom: 1px solid black; padding-right: 12px;"><b>Model 2</b></th>
## <th style="text-align: left; border-top: 2px solid black; border-bottom: 1px solid black; padding-right: 12px;"><b>Model 3</b></th>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">(Intercept)</td>
## <td style="padding-right: 12px; border: none;">-0.55<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">-0.52<sup style="vertical-align: 0px;">*</sup></td>
## <td style="padding-right: 12px; border: none;">-6.31<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.15)</td>
## <td style="padding-right: 12px; border: none;">(0.25)</td>
## <td style="padding-right: 12px; border: none;">(0.57)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">gendermale</td>
## <td style="padding-right: 12px; border: none;">0.03</td>
## <td style="padding-right: 12px; border: none;">0.03</td>
## <td style="padding-right: 12px; border: none;">2.26<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.13)</td>
## <td style="padding-right: 12px; border: none;">(0.13)</td>
## <td style="padding-right: 12px; border: none;">(0.25)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">parental.level.of.educationbachelor's degree</td>
## <td style="padding-right: 12px; border: none;">0.09</td>
## <td style="padding-right: 12px; border: none;">0.04</td>
## <td style="padding-right: 12px; border: none;">-0.47</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.23)</td>
## <td style="padding-right: 12px; border: none;">(0.39)</td>
## <td style="padding-right: 12px; border: none;">(0.44)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">parental.level.of.educationhigh school</td>
## <td style="padding-right: 12px; border: none;">-0.38</td>
## <td style="padding-right: 12px; border: none;">-0.15</td>
## <td style="padding-right: 12px; border: none;">0.47</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.21)</td>
## <td style="padding-right: 12px; border: none;">(0.34)</td>
## <td style="padding-right: 12px; border: none;">(0.39)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">parental.level.of.educationmaster's degree</td>
## <td style="padding-right: 12px; border: none;">-0.13</td>
## <td style="padding-right: 12px; border: none;">0.51</td>
## <td style="padding-right: 12px; border: none;">0.02</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.31)</td>
## <td style="padding-right: 12px; border: none;">(0.47)</td>
## <td style="padding-right: 12px; border: none;">(0.54)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">parental.level.of.educationsome college</td>
## <td style="padding-right: 12px; border: none;">-0.13</td>
## <td style="padding-right: 12px; border: none;">-0.21</td>
## <td style="padding-right: 12px; border: none;">-0.02</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.20)</td>
## <td style="padding-right: 12px; border: none;">(0.34)</td>
## <td style="padding-right: 12px; border: none;">(0.38)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">parental.level.of.educationsome high school</td>
## <td style="padding-right: 12px; border: none;">0.25</td>
## <td style="padding-right: 12px; border: none;">0.01</td>
## <td style="padding-right: 12px; border: none;">0.63</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.21)</td>
## <td style="padding-right: 12px; border: none;">(0.35)</td>
## <td style="padding-right: 12px; border: none;">(0.41)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">lunchstandard</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">-0.05</td>
## <td style="padding-right: 12px; border: none;">-0.02</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.29)</td>
## <td style="padding-right: 12px; border: none;">(0.33)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">parental.level.of.educationbachelor's degree:lunchstandard</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">0.07</td>
## <td style="padding-right: 12px; border: none;">0.17</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.49)</td>
## <td style="padding-right: 12px; border: none;">(0.55)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">parental.level.of.educationhigh school:lunchstandard</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">-0.38</td>
## <td style="padding-right: 12px; border: none;">-0.63</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.44)</td>
## <td style="padding-right: 12px; border: none;">(0.49)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">parental.level.of.educationmaster's degree:lunchstandard</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">-1.17</td>
## <td style="padding-right: 12px; border: none;">-1.41</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.64)</td>
## <td style="padding-right: 12px; border: none;">(0.72)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">parental.level.of.educationsome college:lunchstandard</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">0.12</td>
## <td style="padding-right: 12px; border: none;">-0.13</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.42)</td>
## <td style="padding-right: 12px; border: none;">(0.47)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">parental.level.of.educationsome high school:lunchstandard</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">0.37</td>
## <td style="padding-right: 12px; border: none;">-0.03</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.44)</td>
## <td style="padding-right: 12px; border: none;">(0.50)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">math.score</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">-0.12<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.01)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">writing.score</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">0.18<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.02)</td>
## </tr>
## <tr>
## <td style="border-top: 1px solid black;">AIC</td>
## <td style="border-top: 1px solid black;">1308.86</td>
## <td style="border-top: 1px solid black;">1313.33</td>
## <td style="border-top: 1px solid black;">1091.37</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">BIC</td>
## <td style="padding-right: 12px; border: none;">1343.22</td>
## <td style="padding-right: 12px; border: none;">1377.13</td>
## <td style="padding-right: 12px; border: none;">1164.99</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">Log Likelihood</td>
## <td style="padding-right: 12px; border: none;">-647.43</td>
## <td style="padding-right: 12px; border: none;">-643.66</td>
## <td style="padding-right: 12px; border: none;">-530.69</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">Deviance</td>
## <td style="padding-right: 12px; border: none;">1294.86</td>
## <td style="padding-right: 12px; border: none;">1287.33</td>
## <td style="padding-right: 12px; border: none;">1061.37</td>
## </tr>
## <tr>
## <td style="border-bottom: 2px solid black;">Num. obs.</td>
## <td style="border-bottom: 2px solid black;">1000</td>
## <td style="border-bottom: 2px solid black;">1000</td>
## <td style="border-bottom: 2px solid black;">1000</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;" colspan="5"><span style="font-size:0.8em"><sup style="vertical-align: 0px;">***</sup>p &lt; 0.001, <sup style="vertical-align: 0px;">**</sup>p &lt; 0.01, <sup style="vertical-align: 0px;">*</sup>p &lt; 0.05</span></td>
## </tr>
## </table>
Statistical models
Model 1 Model 2 Model 3
(Intercept) -0.55*** -0.52* -6.31***
(0.15) (0.25) (0.57)
gendermale 0.03 0.03 2.26***
(0.13) (0.13) (0.25)
parental.level.of.educationbachelor’s degree 0.09 0.04 -0.47
(0.23) (0.39) (0.44)
parental.level.of.educationhigh school -0.38 -0.15 0.47
(0.21) (0.34) (0.39)
parental.level.of.educationmaster’s degree -0.13 0.51 0.02
(0.31) (0.47) (0.54)
parental.level.of.educationsome college -0.13 -0.21 -0.02
(0.20) (0.34) (0.38)
parental.level.of.educationsome high school 0.25 0.01 0.63
(0.21) (0.35) (0.41)
lunchstandard -0.05 -0.02
(0.29) (0.33)
parental.level.of.educationbachelor’s degree:lunchstandard 0.07 0.17
(0.49) (0.55)
parental.level.of.educationhigh school:lunchstandard -0.38 -0.63
(0.44) (0.49)
parental.level.of.educationmaster’s degree:lunchstandard -1.17 -1.41
(0.64) (0.72)
parental.level.of.educationsome college:lunchstandard 0.12 -0.13
(0.42) (0.47)
parental.level.of.educationsome high school:lunchstandard 0.37 -0.03
(0.44) (0.50)
math.score -0.12***
(0.01)
writing.score 0.18***
(0.02)
AIC 1308.86 1313.33 1091.37
BIC 1343.22 1377.13 1164.99
Log Likelihood -647.43 -643.66 -530.69
Deviance 1294.86 1287.33 1061.37
Num. obs. 1000 1000 1000
p < 0.001, p < 0.01, p < 0.05 
M3$setx(gender="male")
M3$setx1(gender="female")
M3$sim()
fd <- M3$get_qi(xvalue="x1", qi="fd")
summary(fd)
##        V1         
##  Min.   :-0.6253  
##  1st Qu.:-0.4762  
##  Median :-0.4425  
##  Mean   :-0.4425  
##  3rd Qu.:-0.4070  
##  Max.   :-0.3053
M3$graph()

m3s = min(perform$math.score):max(perform$math.score)
x <- setx(M3, math.score = m3s)
m <- sim(M3, x = x)
m$graph()

M3$setx(lunch="standard")
M3$setx1(lunch="free/reduced")
M3$sim()
fd <- M3$get_qi(xvalue="x1", qi="fd")
summary(fd)
##        V1           
##  Min.   :-0.116473  
##  1st Qu.:-0.008423  
##  Median : 0.016796  
##  Mean   : 0.019078  
##  3rd Qu.: 0.044229  
##  Max.   : 0.182516
M3$graph()