Cruz_MA712_HW6_Zelig5_Titanic

library(sjmisc)
library(radiant.data)
library(tidyverse)
library(Zelig)
library(texreg)
library(mvtnorm)
library(radiant.data)
library(sjmisc)
library(lattice)
library(dplyr)
library(readr)
data(titanic)
head(titanic)

titanic <- titanic %>% 
  mutate(surv = as.integer(survived))

titanic <- titanic %>% 
  select(pclass, survived, surv, everything())
head(titanic)

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

library(Zelig)
z.tit1 <- zlogit$new()
z.tit1$zelig(survival ~ age, data = titanic)
summary(z.tit1)

Model: 

Call:
z.tit1$zelig(formula = survival ~ age, data = titanic)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1211  -1.0351  -0.9736   1.3195   1.5247  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.131200   0.145095  -0.904   0.3659
age         -0.008202   0.004431  -1.851   0.0642

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1410.0  on 1042  degrees of freedom
Residual deviance: 1406.5  on 1041  degrees of freedom
AIC: 1410.5

Number of Fisher Scoring iterations: 4

Next step: Use 'setx' method

library(Zelig)
z.tit2 <- zlogit$new()
z.tit2$zelig(survival ~ sex, data = titanic)
summary(z.tit2)

Model: 

Call:
z.tit2$zelig(formula = survival ~ sex, data = titanic)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6682  -0.6783  -0.6783   0.7562   1.7790  

Coefficients:
            Estimate Std. Error z value            Pr(>|z|)
(Intercept)   1.1055     0.1177   9.389 <0.0000000000000002
sexmale      -2.4579     0.1523 -16.141 <0.0000000000000002

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1410.0  on 1042  degrees of freedom
Residual deviance: 1100.4  on 1041  degrees of freedom
AIC: 1104.4

Number of Fisher Scoring iterations: 4

Next step: Use 'setx' method

library(Zelig)
z.tit3 <- zlogit$new()
z.tit3$zelig(survival ~ sex*pclass, data = titanic)
summary(z.tit3)

Model: 

Call:
z.tit3$zelig(formula = survival ~ sex * pclass, data = titanic)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.5557  -0.6096  -0.5609   0.4753   1.9632  

Coefficients:
                  Estimate Std. Error z value            Pr(>|z|)
(Intercept)        3.22684    0.45584   7.079 0.00000000000145275
sexmale           -3.84152    0.48668  -7.893 0.00000000000000294
pclass2nd         -1.10295    0.55639  -1.982              0.0474
pclass3rd         -3.33220    0.48392  -6.886 0.00000000000574319
sexmale:pclass2nd -0.05215    0.62412  -0.084              0.9334
sexmale:pclass3rd  2.35799    0.53260   4.427 0.00000953965064616

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1409.99  on 1042  degrees of freedom
Residual deviance:  966.38  on 1037  degrees of freedom
AIC: 978.38

Number of Fisher Scoring iterations: 5

Next step: Use 'setx' method

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

Model: 

Call:
z.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

The Effect of Age

z.tit <- z.tit
z.tit$setrange(age = 0:100)
z.tit$sim()
z.tit$graph()

The Effect of Fare

z.tit$setrange(fare = 0:700)
z.tit$sim()
z.tit$graph()

Sex Effect

z.tit2$setx(sex = "male")
z.tit2$setx1(sex = "female")
z.tit2$sim()
summary(z.tit2)

 sim x :
 -----
ev
          mean         sd       50%      2.5%     97.5%
[1,] 0.2065737 0.01527193 0.2062044 0.1772126 0.2371819
pv
         0     1
[1,] 0.793 0.207

 sim x1 :
 -----
ev
          mean         sd      50%      2.5%     97.5%
[1,] 0.7511008 0.02250548 0.752001 0.7050649 0.7947439
pv
         0     1
[1,] 0.245 0.755
fd
          mean         sd       50%      2.5%     97.5%
[1,] 0.5445271 0.02719262 0.5455427 0.4909122 0.5963622

fd <- z.tit$get_qi(xvalue="x1", qi="fd")
summary(fd)

       V1        
 Min.   :0.3245  
 1st Qu.:0.4831  
 Median :0.5178  
 Mean   :0.5173  
 3rd Qu.:0.5506  
 Max.   :0.6634  

z.tit$setx(sex = "male", pclass = "1st") 
z.tit$setx1(sex = "female", pclass = "1st")
z.tit$sim()
graphics.off()
 par("mar")

[1] 5.1 4.1 2.1 2.1

 par(mar=c(1,1,1,1))
plot(z.tit)

sd1 <- z.tit$get_qi(xvalue="x1", qi="fd")
summary(sd1)

       V1        
 Min.   :0.3781  
 1st Qu.:0.4845  
 Median :0.5191  
 Mean   :0.5187  
 3rd Qu.:0.5520  
 Max.   :0.6738  

z.tit3$setx(sex = "male", pclass = "2nd") 
z.tit3$setx1(sex = "female", pclass = "2nd")
z.tit3$sim()
graphics.off()
 par("mar")

[1] 5.1 4.1 2.1 2.1

 par(mar=c(1,1,1,1))
plot(z.tit3)

sd2 <- z.tit$get_qi(xvalue="x1", qi="fd")
summary(sd2)

       V1        
 Min.   :0.3781  
 1st Qu.:0.4845  
 Median :0.5191  
 Mean   :0.5187  
 3rd Qu.:0.5520  
 Max.   :0.6738  

z.tit3$setx(sex = "male", pclass = "3rd") 
z.tit3$setx1(sex = "female", pclass = "3rd")
z.tit3$sim()
graphics.off()
 par("mar")

[1] 5.1 4.1 2.1 2.1

 par(mar=c(1,1,1,1))
plot(z.tit3)

sd3 <- z.tit$get_qi(xvalue="x1", qi="fd")
summary(sd3)

       V1        
 Min.   :0.3781  
 1st Qu.:0.4845  
 Median :0.5191  
 Mean   :0.5187  
 3rd Qu.:0.5520  
 Max.   :0.6738  

dfd <- as.data.frame(cbind(sd1, sd2, sd3))
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)