library(Zelig)
## Loading required package: survival
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(sjmisc)
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
library(radiant.data)
## Loading required package: magrittr
## Loading required package: ggplot2
## 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:sjmisc':
##
## replace_na
##
## Attaching package: 'radiant.data'
## The following objects are masked from 'package:lubridate':
##
## month, wday
## The following object is masked from 'package:ggplot2':
##
## diamonds
## The following objects are masked from 'package:sjmisc':
##
## center, is_empty
## The following object is masked from 'package:dplyr':
##
## mutate_each
library(tidyr)
data(titanic)
titanic <- titanic %>%
mutate(age = as.numeric(age),
surv= as.integer(survived),
survival = ifelse(surv == 1,1,0),
sex=as.factor(sex))%>%
select(pclass, survived, age ,survival, surv, sex, fare)
head(titanic)
z1 <- zlogit$new()
z1$zelig(survival ~ age + sex*pclass + fare, data = titanic)
z1$setx(sex = "male")
z1$setx1(sex = "female")
z1$sim()
summary(z1)
##
## sim x :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.1402391 0.01975662 0.1391319 0.1037799 0.1826869
## pv
## 0 1
## [1,] 0.844 0.156
##
## sim x1 :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.397574 0.04413851 0.3983417 0.3089826 0.48479
## pv
## 0 1
## [1,] 0.592 0.408
## fd
## mean sd 50% 2.5% 97.5%
## [1,] 0.2573349 0.04459602 0.2559749 0.1700498 0.3458138
plot(z1)
z1 <- zlogit$new()
z1$zelig(survival ~ age + sex*pclass + fare, data = titanic)
summary(z1)
## Model:
##
## Call:
## z1$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
z1$setrange(age = min(titanic$age):max(titanic$age))
z1$sim()
z1$graph()
z1$setrange(fare = min(titanic$fare):max(titanic$fare))
z1$sim()
z1$graph()
z1 <- zlogit$new()
z1$zelig(survival ~ age + sex*pclass + fare, data = titanic)
z1$setx(sex = "male")
z1$setx1(sex = "female")
z1$sim()
summary(z1)
##
## sim x :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.1399452 0.01909118 0.1388476 0.1056201 0.1824067
## pv
## 0 1
## [1,] 0.874 0.126
##
## sim x1 :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.3953875 0.04379721 0.3939994 0.3163083 0.4796405
## pv
## 0 1
## [1,] 0.61 0.39
## fd
## mean sd 50% 2.5% 97.5%
## [1,] 0.2554423 0.04468766 0.2551999 0.1733354 0.3399038
fd <- z1$get_qi(xvalue="x1", qi="fd")
summary(fd)
## V1
## Min. :0.1355
## 1st Qu.:0.2228
## Median :0.2552
## Mean :0.2554
## 3rd Qu.:0.2852
## Max. :0.4243
plot(z1)
z2 <- zlogit$new()
z2$zelig(survival ~ age + sex*pclass + fare, data = titanic)
z2$setx(sex = "male", pclass = "1st")
z2$setx1(sex = "female", pclass = "1st")
z2$sim()
summary(z2)
##
## sim x :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.4553696 0.05043271 0.4525382 0.3549129 0.5509884
## pv
## 0 1
## [1,] 0.543 0.457
##
## sim x1 :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.9733808 0.01460295 0.9768674 0.9383972 0.9914779
## pv
## 0 1
## [1,] 0.021 0.979
## fd
## mean sd 50% 2.5% 97.5%
## [1,] 0.5180111 0.04970335 0.5201882 0.4183179 0.619113
plot(z2)
z3 <- zlogit$new()
z3$zelig(survival ~ age + sex*pclass + fare, data = titanic)
z3$setx(sex = "male", pclass = "2nd")
z3$setx1(sex = "female", pclass = "2nd")
z3$sim()
summary(z3)
##
## sim x :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.1391962 0.02807985 0.1367332 0.08992457 0.1998919
## pv
## 0 1
## [1,] 0.862 0.138
##
## sim x1 :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.8883352 0.0324761 0.8920388 0.8148295 0.9366309
## pv
## 0 1
## [1,] 0.117 0.883
## fd
## mean sd 50% 2.5% 97.5%
## [1,] 0.749139 0.04350102 0.7519973 0.6539815 0.8262562
plot(z3)
z4 <- zlogit$new()
z4$zelig(survival ~ age + sex*pclass + fare, data = titanic)
z4$setx(sex = "male", pclass = "3rd")
z4$setx1(sex = "female", pclass = "3rd")
z4$sim()
summary(z4)
##
## sim x :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.1404313 0.01966434 0.1386446 0.1070164 0.1833493
## pv
## 0 1
## [1,] 0.871 0.129
##
## sim x1 :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.3929052 0.0438295 0.3923728 0.3131241 0.4817764
## pv
## 0 1
## [1,] 0.622 0.378
## fd
## mean sd 50% 2.5% 97.5%
## [1,] 0.2524739 0.04393798 0.2518535 0.1681392 0.3395827
plot(z4)
d1 <- z2$get_qi(xvalue="x1", qi="fd")
d2 <- z3$get_qi(xvalue="x1", qi="fd")
d3 <- z4$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)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
