Week Seven Part I

Interaction Model

library(Zelig)
library(dplyr)
library(radiant.data)
library(tidyr)
library(ggplot2)
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 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

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

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)

Sex Difference

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.1408838 0.01944116 0.1396559 0.1063131 0.1824498
## pv
##          0     1
## [1,] 0.846 0.154
## 
##  sim x1 :
##  -----
## ev
##           mean         sd       50%      2.5%     97.5%
## [1,] 0.3968025 0.04359156 0.3976338 0.3136083 0.4793755
## pv
##          0     1
## [1,] 0.581 0.419
## fd
##           mean         sd      50%      2.5%     97.5%
## [1,] 0.2559187 0.04324699 0.255616 0.1750195 0.3391593

First Difference

fd <- s$get_qi(xvalue="x1", qi="fd")
summary(fd)
##        V1        
##  Min.   :0.1317  
##  1st Qu.:0.2256  
##  Median :0.2556  
##  Mean   :0.2559  
##  3rd Qu.:0.2849  
##  Max.   :0.3699

Plots

plot(s)

Test The Class Variation in Sex

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

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


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

plot(c1s)

plot(c2s)

plot(c3s)

Putting Them In One Place

d1 <- c1s$get_qi(xvalue="x1", qi="fd")
d2 <- c2s$get_qi(xvalue="x1", qi="fd")
d3 <- c3s$get_qi(xvalue="x1", qi="fd")

dfd <- as.data.frame(cbind(d1, d2, d3))
head(dfd)
##          V1        V2        V3
## 1 0.5119357 0.8152948 0.2543196
## 2 0.5020596 0.7881658 0.3006135
## 3 0.5658872 0.7961920 0.3538066
## 4 0.4655966 0.8430529 0.2067243
## 5 0.5481982 0.7521200 0.2145423
## 6 0.5265109 0.7590866 0.3116191
tidd <- dfd %>% 
  gather(class, simv, 1:3)
head(tidd)
##   class      simv
## 1    V1 0.5119357
## 2    V1 0.5020596
## 3    V1 0.5658872
## 4    V1 0.4655966
## 5    V1 0.5481982
## 6    V1 0.5265109
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.521 0.0487
## 2 V2    0.747 0.0427
## 3 V3    0.255 0.0442
ggplot(tidd, aes(simv)) + geom_histogram() + facet_grid(~class)