Zelig for NOT WTP
Thang, TB
22 March 2017
reading data
library(foreign)## Warning: package 'foreign' was built under R version 3.2.5r=read.dta("C:/Users/BINH THANG/Dropbox/Korea/STudy/Thesis/data management/DataR/dataR3.dta")
r1 <- subset(r)
attach(r1)
r1$wtp[cost1 == 0] <- 0
r1$wtp[cost1 == 1] <- 0
r1$wtp[is.na(r1$wtp)] <- 1
r1$b16a[b16a == 5] <- 0
attach(r1)## The following objects are masked from r1 (pos = 3):
##
## a1, advice, advice1, ag, age_group, anticam, anticam2, b1,
## b10, b10a1, b10a2, b10a3, b10a4, b10a5, b10a6, b10a7, b11,
## b11a, b11a2, b11a3, b12, b12a, b13, b14, b15, b16, b16a, b17,
## b18, b18a, b1a, b2, b2a1, b3, b4, b5, b5a, b6, b6a, b6a1,
## b6a2, b6a3, b6a4, b6a5, b6a6, b6a7, b6a7a, b7, b8, b9, br,
## branch, branch1, branch2, branch3, c1, c2, c23a, c3, c4, c5,
## c6, c7, c8, c9, COST, cost_inc, cost1, costincrease, ct, d1,
## d10, d11, d1a, d2, d3, d3a, d4, d5, d6, d7, d8, d9, Decision,
## e1, e2, edu, educ1, educ2, f1, ghi1, ghi2, ghiro2,
## giadinhkoUH, group_age, group_age1, h1, h10a, h10a1, h10a10,
## h10a11, h10a11a, h10a2, h10a3, h10a4, h10a5, h10a6, h10a7,
## h10a8, h10a9, h12, h12a, h12a_1, h12log, h13, h2, h3, h4, h5,
## h6, h7, h8, h9, ha000, itc1, itc2, l1, l2, l3, l4, l5, l6,
## label1, label2, moneyspent, msdt, n01, n02, n03, n05, n06,
## n07, n08, n1, n10, n100, n101, n102, n103, n11, n12, n13, n14,
## n15, n16, n1b, n2, n3, n35, n36, n37, n38, n39, n3posterb, n4,
## n40, n41, n42, n43, n44, n45, n46, n47, n48, n49, n5, n50,
## n51, n52, n53, n54, n55, n56, n57, n58, n59, n6, n60, n61,
## n61a, n62, n62a, n63, n63a, n64, n64a, n65, n65a, n66, n66a,
## n67, n67a, n68, n68a, n7, n77, n78, n7tren, n8, n88, n89,
## n8nha, n9, n97, n98, n99, n9khach, noEnvi2, occup1, policy,
## policy_a, policy2, reasons, reasons1, s1, Screening,
## selfhealth, SH1, smostt, ter_fa1, ter_in, tertile_fa,
## tertile_indi, test, unitsdiffi1, var242, w1, w2, w3, w4mytable <- table(wtp, c1)
print(mytable)## c1
## wtp 1 2
## 0 398 63
## 1 329 32prop.table(mytable)## c1
## wtp 1 2
## 0 0.48418491 0.07664234
## 1 0.40024331 0.03892944newdata2=r1
attach(newdata2 )## The following objects are masked from r1 (pos = 3):
##
## a1, advice, advice1, ag, age_group, anticam, anticam2, b1,
## b10, b10a1, b10a2, b10a3, b10a4, b10a5, b10a6, b10a7, b11,
## b11a, b11a2, b11a3, b12, b12a, b13, b14, b15, b16, b16a, b17,
## b18, b18a, b1a, b2, b2a1, b3, b4, b5, b5a, b6, b6a, b6a1,
## b6a2, b6a3, b6a4, b6a5, b6a6, b6a7, b6a7a, b7, b8, b9, br,
## branch, branch1, branch2, branch3, c1, c2, c23a, c3, c4, c5,
## c6, c7, c8, c9, COST, cost_inc, cost1, costincrease, ct, d1,
## d10, d11, d1a, d2, d3, d3a, d4, d5, d6, d7, d8, d9, Decision,
## e1, e2, edu, educ1, educ2, f1, ghi1, ghi2, ghiro2,
## giadinhkoUH, group_age, group_age1, h1, h10a, h10a1, h10a10,
## h10a11, h10a11a, h10a2, h10a3, h10a4, h10a5, h10a6, h10a7,
## h10a8, h10a9, h12, h12a, h12a_1, h12log, h13, h2, h3, h4, h5,
## h6, h7, h8, h9, ha000, itc1, itc2, l1, l2, l3, l4, l5, l6,
## label1, label2, moneyspent, msdt, n01, n02, n03, n05, n06,
## n07, n08, n1, n10, n100, n101, n102, n103, n11, n12, n13, n14,
## n15, n16, n1b, n2, n3, n35, n36, n37, n38, n39, n3posterb, n4,
## n40, n41, n42, n43, n44, n45, n46, n47, n48, n49, n5, n50,
## n51, n52, n53, n54, n55, n56, n57, n58, n59, n6, n60, n61,
## n61a, n62, n62a, n63, n63a, n64, n64a, n65, n65a, n66, n66a,
## n67, n67a, n68, n68a, n7, n77, n78, n7tren, n8, n88, n89,
## n8nha, n9, n97, n98, n99, n9khach, noEnvi2, occup1, policy,
## policy_a, policy2, reasons, reasons1, s1, Screening,
## selfhealth, SH1, smostt, ter_fa1, ter_in, tertile_fa,
## tertile_indi, test, unitsdiffi1, var242, w1, w2, w3, w4, wtp## The following objects are masked from r1 (pos = 4):
##
## a1, advice, advice1, ag, age_group, anticam, anticam2, b1,
## b10, b10a1, b10a2, b10a3, b10a4, b10a5, b10a6, b10a7, b11,
## b11a, b11a2, b11a3, b12, b12a, b13, b14, b15, b16, b16a, b17,
## b18, b18a, b1a, b2, b2a1, b3, b4, b5, b5a, b6, b6a, b6a1,
## b6a2, b6a3, b6a4, b6a5, b6a6, b6a7, b6a7a, b7, b8, b9, br,
## branch, branch1, branch2, branch3, c1, c2, c23a, c3, c4, c5,
## c6, c7, c8, c9, COST, cost_inc, cost1, costincrease, ct, d1,
## d10, d11, d1a, d2, d3, d3a, d4, d5, d6, d7, d8, d9, Decision,
## e1, e2, edu, educ1, educ2, f1, ghi1, ghi2, ghiro2,
## giadinhkoUH, group_age, group_age1, h1, h10a, h10a1, h10a10,
## h10a11, h10a11a, h10a2, h10a3, h10a4, h10a5, h10a6, h10a7,
## h10a8, h10a9, h12, h12a, h12a_1, h12log, h13, h2, h3, h4, h5,
## h6, h7, h8, h9, ha000, itc1, itc2, l1, l2, l3, l4, l5, l6,
## label1, label2, moneyspent, msdt, n01, n02, n03, n05, n06,
## n07, n08, n1, n10, n100, n101, n102, n103, n11, n12, n13, n14,
## n15, n16, n1b, n2, n3, n35, n36, n37, n38, n39, n3posterb, n4,
## n40, n41, n42, n43, n44, n45, n46, n47, n48, n49, n5, n50,
## n51, n52, n53, n54, n55, n56, n57, n58, n59, n6, n60, n61,
## n61a, n62, n62a, n63, n63a, n64, n64a, n65, n65a, n66, n66a,
## n67, n67a, n68, n68a, n7, n77, n78, n7tren, n8, n88, n89,
## n8nha, n9, n97, n98, n99, n9khach, noEnvi2, occup1, policy,
## policy_a, policy2, reasons, reasons1, s1, Screening,
## selfhealth, SH1, smostt, ter_fa1, ter_in, tertile_fa,
## tertile_indi, test, unitsdiffi1, var242, w1, w2, w3, w4define varibles
newdata2$c1=as.factor(newdata2$c1)
newdata2$h1=as.factor(newdata2$h1)
newdata2$label1=as.factor(newdata2$label1)
newdata2$educ2=as.factor(newdata2$educ2)
newdata2$age_group=as.factor(newdata2$age_group)
newdata2$d1a=as.factor(newdata2$d1a)
newdata2$selfhealth=as.factor(newdata2$selfhealth)
newdata2$b18a=as.factor(newdata2$b18a)
newdata2$b16a=as.factor(newdata2$b16a)
newdata2$b6a=as.factor(newdata2$b6a)
newdata2$smostt=as.factor(newdata2$smostt)
newdata2$reasons1=as.factor(newdata2$reasons1)
newdata2$ter_in=as.factor(newdata2$ter_in)
newdata2$group_age1=as.factor(newdata2$group_age1)
attach(newdata2)## The following objects are masked from newdata2 (pos = 3):
##
## a1, advice, advice1, ag, age_group, anticam, anticam2, b1,
## b10, b10a1, b10a2, b10a3, b10a4, b10a5, b10a6, b10a7, b11,
## b11a, b11a2, b11a3, b12, b12a, b13, b14, b15, b16, b16a, b17,
## b18, b18a, b1a, b2, b2a1, b3, b4, b5, b5a, b6, b6a, b6a1,
## b6a2, b6a3, b6a4, b6a5, b6a6, b6a7, b6a7a, b7, b8, b9, br,
## branch, branch1, branch2, branch3, c1, c2, c23a, c3, c4, c5,
## c6, c7, c8, c9, COST, cost_inc, cost1, costincrease, ct, d1,
## d10, d11, d1a, d2, d3, d3a, d4, d5, d6, d7, d8, d9, Decision,
## e1, e2, edu, educ1, educ2, f1, ghi1, ghi2, ghiro2,
## giadinhkoUH, group_age, group_age1, h1, h10a, h10a1, h10a10,
## h10a11, h10a11a, h10a2, h10a3, h10a4, h10a5, h10a6, h10a7,
## h10a8, h10a9, h12, h12a, h12a_1, h12log, h13, h2, h3, h4, h5,
## h6, h7, h8, h9, ha000, itc1, itc2, l1, l2, l3, l4, l5, l6,
## label1, label2, moneyspent, msdt, n01, n02, n03, n05, n06,
## n07, n08, n1, n10, n100, n101, n102, n103, n11, n12, n13, n14,
## n15, n16, n1b, n2, n3, n35, n36, n37, n38, n39, n3posterb, n4,
## n40, n41, n42, n43, n44, n45, n46, n47, n48, n49, n5, n50,
## n51, n52, n53, n54, n55, n56, n57, n58, n59, n6, n60, n61,
## n61a, n62, n62a, n63, n63a, n64, n64a, n65, n65a, n66, n66a,
## n67, n67a, n68, n68a, n7, n77, n78, n7tren, n8, n88, n89,
## n8nha, n9, n97, n98, n99, n9khach, noEnvi2, occup1, policy,
## policy_a, policy2, reasons, reasons1, s1, Screening,
## selfhealth, SH1, smostt, ter_fa1, ter_in, tertile_fa,
## tertile_indi, test, unitsdiffi1, var242, w1, w2, w3, w4, wtp## The following objects are masked from r1 (pos = 4):
##
## a1, advice, advice1, ag, age_group, anticam, anticam2, b1,
## b10, b10a1, b10a2, b10a3, b10a4, b10a5, b10a6, b10a7, b11,
## b11a, b11a2, b11a3, b12, b12a, b13, b14, b15, b16, b16a, b17,
## b18, b18a, b1a, b2, b2a1, b3, b4, b5, b5a, b6, b6a, b6a1,
## b6a2, b6a3, b6a4, b6a5, b6a6, b6a7, b6a7a, b7, b8, b9, br,
## branch, branch1, branch2, branch3, c1, c2, c23a, c3, c4, c5,
## c6, c7, c8, c9, COST, cost_inc, cost1, costincrease, ct, d1,
## d10, d11, d1a, d2, d3, d3a, d4, d5, d6, d7, d8, d9, Decision,
## e1, e2, edu, educ1, educ2, f1, ghi1, ghi2, ghiro2,
## giadinhkoUH, group_age, group_age1, h1, h10a, h10a1, h10a10,
## h10a11, h10a11a, h10a2, h10a3, h10a4, h10a5, h10a6, h10a7,
## h10a8, h10a9, h12, h12a, h12a_1, h12log, h13, h2, h3, h4, h5,
## h6, h7, h8, h9, ha000, itc1, itc2, l1, l2, l3, l4, l5, l6,
## label1, label2, moneyspent, msdt, n01, n02, n03, n05, n06,
## n07, n08, n1, n10, n100, n101, n102, n103, n11, n12, n13, n14,
## n15, n16, n1b, n2, n3, n35, n36, n37, n38, n39, n3posterb, n4,
## n40, n41, n42, n43, n44, n45, n46, n47, n48, n49, n5, n50,
## n51, n52, n53, n54, n55, n56, n57, n58, n59, n6, n60, n61,
## n61a, n62, n62a, n63, n63a, n64, n64a, n65, n65a, n66, n66a,
## n67, n67a, n68, n68a, n7, n77, n78, n7tren, n8, n88, n89,
## n8nha, n9, n97, n98, n99, n9khach, noEnvi2, occup1, policy,
## policy_a, policy2, reasons, reasons1, s1, Screening,
## selfhealth, SH1, smostt, ter_fa1, ter_in, tertile_fa,
## tertile_indi, test, unitsdiffi1, var242, w1, w2, w3, w4, wtp## The following objects are masked from r1 (pos = 5):
##
## a1, advice, advice1, ag, age_group, anticam, anticam2, b1,
## b10, b10a1, b10a2, b10a3, b10a4, b10a5, b10a6, b10a7, b11,
## b11a, b11a2, b11a3, b12, b12a, b13, b14, b15, b16, b16a, b17,
## b18, b18a, b1a, b2, b2a1, b3, b4, b5, b5a, b6, b6a, b6a1,
## b6a2, b6a3, b6a4, b6a5, b6a6, b6a7, b6a7a, b7, b8, b9, br,
## branch, branch1, branch2, branch3, c1, c2, c23a, c3, c4, c5,
## c6, c7, c8, c9, COST, cost_inc, cost1, costincrease, ct, d1,
## d10, d11, d1a, d2, d3, d3a, d4, d5, d6, d7, d8, d9, Decision,
## e1, e2, edu, educ1, educ2, f1, ghi1, ghi2, ghiro2,
## giadinhkoUH, group_age, group_age1, h1, h10a, h10a1, h10a10,
## h10a11, h10a11a, h10a2, h10a3, h10a4, h10a5, h10a6, h10a7,
## h10a8, h10a9, h12, h12a, h12a_1, h12log, h13, h2, h3, h4, h5,
## h6, h7, h8, h9, ha000, itc1, itc2, l1, l2, l3, l4, l5, l6,
## label1, label2, moneyspent, msdt, n01, n02, n03, n05, n06,
## n07, n08, n1, n10, n100, n101, n102, n103, n11, n12, n13, n14,
## n15, n16, n1b, n2, n3, n35, n36, n37, n38, n39, n3posterb, n4,
## n40, n41, n42, n43, n44, n45, n46, n47, n48, n49, n5, n50,
## n51, n52, n53, n54, n55, n56, n57, n58, n59, n6, n60, n61,
## n61a, n62, n62a, n63, n63a, n64, n64a, n65, n65a, n66, n66a,
## n67, n67a, n68, n68a, n7, n77, n78, n7tren, n8, n88, n89,
## n8nha, n9, n97, n98, n99, n9khach, noEnvi2, occup1, policy,
## policy_a, policy2, reasons, reasons1, s1, Screening,
## selfhealth, SH1, smostt, ter_fa1, ter_in, tertile_fa,
## tertile_indi, test, unitsdiffi1, var242, w1, w2, w3, w4Packages and library of
library(Zelig)## Warning: package 'Zelig' was built under R version 3.2.5Model1
z.out<- zelig(wtp~age_group+h1+ d1a+educ2+ter_in+selfhealth+smostt
+b16a+b18a+label1+anticam+noEnvi2+reasons,
above = 150000, thin=5, mcmc=30000,
data=newdata2, model = "logit.bayes", verbose = FALSE)## How to cite this model in Zelig:
## Ben Goodrich, and Ying Lu. 2013.
## logit-bayes: Bayesian Logistic Regression for Dichotomous Dependent Variables
## in Christine Choirat, Christopher Gandrud, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
## "Zelig: Everyone's Statistical Software," http://zeligproject.org/You can check for convergence before summarizing the estimates with three diagnostic tests. See the section Diagnostics for Zelig Models for examples of the output with interpretation:
z.out$geweke.diag()##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## (Intercept) age_groupgr3039 age_groupgr4049 age_groupgr5059
## -0.45150 -1.28943 -1.00804 -1.31650
## age_group60plus h12 d1amarried educ22
## -0.69423 -0.26735 0.95832 -1.32998
## educ23 ter_in2 ter_in3 selfhealthgood
## -0.76254 -0.44250 -1.16195 -1.13595
## smosttmedium smosttheavy b16a1 b18abad
## 1.30202 0.08466 1.45160 1.46227
## label11 anticam noEnvi2 reasons
## -0.36386 0.45640 0.86699 0.08843z.out$heidel.diag()##
## Stationarity start p-value
## test iteration
## (Intercept) passed 1 0.948
## age_groupgr3039 passed 1 0.649
## age_groupgr4049 passed 1 0.817
## age_groupgr5059 passed 1 0.277
## age_group60plus passed 1 0.963
## h12 passed 1 0.825
## d1amarried passed 1 0.783
## educ22 passed 1 0.267
## educ23 passed 1 0.405
## ter_in2 passed 1 0.449
## ter_in3 passed 1 0.805
## selfhealthgood passed 1 0.866
## smosttmedium passed 1 0.380
## smosttheavy passed 1 0.647
## b16a1 passed 1 0.580
## b18abad passed 1 0.220
## label11 passed 1 0.689
## anticam passed 1 0.379
## noEnvi2 passed 1 0.772
## reasons passed 1 0.716
##
## Halfwidth Mean Halfwidth
## test
## (Intercept) failed -0.0194 0.05159
## age_groupgr3039 passed 0.4486 0.03990
## age_groupgr4049 passed 0.5818 0.04374
## age_groupgr5059 passed 0.5858 0.04490
## age_group60plus passed 0.8566 0.05104
## h12 failed 0.0429 0.02402
## d1amarried failed -0.3117 0.03480
## educ22 failed -0.2385 0.03133
## educ23 failed 0.1850 0.04093
## ter_in2 failed 0.2920 0.02991
## ter_in3 passed 0.4843 0.02776
## selfhealthgood failed 0.1180 0.02233
## smosttmedium failed -0.1499 0.03007
## smosttheavy passed 0.3877 0.03415
## b16a1 passed -0.7619 0.02832
## b18abad passed -0.3502 0.02815
## label11 passed 0.6087 0.03591
## anticam failed -0.0316 0.03577
## noEnvi2 failed -0.1672 0.03017
## reasons failed -0.0480 0.00562z.out$raftery.diag()##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## (Intercept) 245 269370 3746 71.9
## age_groupgr3039 270 293620 3746 78.4
## age_groupgr4049 360 408455 3746 109.0
## age_groupgr5059 275 313335 3746 83.6
## age_group60plus 235 269195 3746 71.9
## h12 545 654885 3746 175.0
## d1amarried 245 269370 3746 71.9
## educ22 225 245660 3746 65.6
## educ23 460 518590 3746 138.0
## ter_in2 310 350445 3746 93.6
## ter_in3 340 384765 3746 103.0
## selfhealthgood 315 354950 3746 94.8
## smosttmedium 285 316135 3746 84.4
## smosttheavy 220 249555 3746 66.6
## b16a1 215 236780 3746 63.2
## b18abad 225 245660 3746 65.6
## label11 200 224195 3746 59.8
## anticam 270 297515 3746 79.4
## noEnvi2 270 293620 3746 78.4
## reasons 165 179805 3746 48.0summary(z.out)## Model:
##
## Iterations = 1001:30996
## Thinning interval = 5
## Number of chains = 1
## Sample size per chain = 6000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## (Intercept) -0.01936 0.34543 0.0044595 0.026322
## age_groupgr3039 0.44860 0.26577 0.0034311 0.020358
## age_groupgr4049 0.58184 0.28811 0.0037195 0.022317
## age_groupgr5059 0.58582 0.28987 0.0037423 0.022906
## age_group60plus 0.85656 0.33922 0.0043793 0.026042
## h12 0.04292 0.16392 0.0021162 0.012254
## d1amarried -0.31172 0.22938 0.0029613 0.017753
## educ22 -0.23854 0.20310 0.0026220 0.015985
## educ23 0.18495 0.25257 0.0032607 0.020880
## ter_in2 0.29195 0.19607 0.0025313 0.015262
## ter_in3 0.48428 0.19039 0.0024579 0.014162
## selfhealthgood 0.11795 0.15201 0.0019624 0.011391
## smosttmedium -0.14988 0.20889 0.0026968 0.015340
## smosttheavy 0.38766 0.22601 0.0029178 0.017423
## b16a1 -0.76190 0.19009 0.0024540 0.014450
## b18abad -0.35017 0.16847 0.0021750 0.014363
## label11 0.60873 0.23171 0.0029914 0.018320
## anticam -0.03158 0.21738 0.0028064 0.018252
## noEnvi2 -0.16721 0.20753 0.0026792 0.015394
## reasons -0.04796 0.03634 0.0004691 0.002867
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## (Intercept) -0.642399 -0.24848 -0.00481 0.19490 0.64057
## age_groupgr3039 -0.091911 0.26543 0.41992 0.61919 1.04710
## age_groupgr4049 -0.008236 0.40638 0.57387 0.77719 1.13800
## age_groupgr5059 0.017918 0.38635 0.56045 0.76872 1.16856
## age_group60plus 0.177202 0.62711 0.85999 1.09071 1.52637
## h12 -0.288542 -0.06151 0.03815 0.14725 0.37973
## d1amarried -0.769350 -0.46465 -0.30238 -0.16823 0.16865
## educ22 -0.600874 -0.39142 -0.24567 -0.10940 0.19459
## educ23 -0.364928 0.03253 0.19186 0.32090 0.72369
## ter_in2 -0.100021 0.18120 0.29182 0.43551 0.67171
## ter_in3 0.126922 0.35117 0.46892 0.60669 0.86837
## selfhealthgood -0.216456 0.03193 0.11545 0.22356 0.40254
## smosttmedium -0.635523 -0.27987 -0.12812 -0.01059 0.22307
## smosttheavy -0.076335 0.25492 0.39313 0.51244 0.88842
## b16a1 -1.140302 -0.88214 -0.74817 -0.64736 -0.36741
## b18abad -0.682144 -0.46042 -0.33508 -0.23589 -0.01542
## label11 0.146424 0.47830 0.60523 0.74747 1.08282
## anticam -0.452560 -0.15923 -0.02424 0.10755 0.39605
## noEnvi2 -0.591611 -0.30620 -0.15659 -0.05140 0.26584
## reasons -0.121594 -0.07099 -0.04686 -0.02316 0.02237
##
## Next step: Use 'setx' methodSetting values for the explanatory variables to their sample averages:
x.out <- setx(z.out)Simulating quantities of interest from the posterior distribution given: x.out
s.out1 <- sim(z.out, x = x.out)
summary(s.out1)##
## sim x :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.204286 0.05003832 0.2040698 0.1193315 0.3131794
## pv
## 0 1
## [1,] 4.798 1.202Simulating First Differences Estimating the first difference (and risk ratio) in individual’s probability of voting when education is set to be low (25th percentile) versus high (75th percentile) while all the other variables are held at their default values:
x.high <- setx(z.out, reasons = quantile(newdata2$reasons, prob = 0.8))
x.low <- setx(z.out, reasons= quantile(newdata2$reasons, prob = 0.2))s.out2 <- sim(z.out, x = x.high, x1 = x.low)
summary(s.out2)##
## sim x :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.1976573 0.04879861 0.1987198 0.1161147 0.3045151
## pv
## 0 1
## [1,] 4.8 1.2
##
## sim x1 :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.2131926 0.05244707 0.2111619 0.1238252 0.3275657
## pv
## 0 1
## [1,] 4.692 1.308
## fd
## mean sd 50% 2.5% 97.5%
## [1,] 0.01553538 0.01243249 0.01554525 -0.00753272 0.04404592Checking
z.out1 <- zelig(wtp~age_group+h1+ d1a+educ2+ter_in+selfhealth+smostt
+b16a+b18a+label1+anticam+noEnvi2+reasons, model = "logit", data = newdata2)## How to cite this model in Zelig:
## R Core Team. 2007.
## logit: Logistic Regression for Dichotomous Dependent Variables
## in Christine Choirat, Christopher Gandrud, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
## "Zelig: Everyone's Statistical Software," http://zeligproject.org/summary(z.out1)## Model:
##
## Call:
## z5$zelig(formula = wtp ~ age_group + h1 + d1a + educ2 + ter_in +
## selfhealth + smostt + b16a + b18a + label1 + anticam + noEnvi2 +
## reasons, data = newdata2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8643 -1.0122 -0.7353 1.1741 1.9335
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.00481 0.34331 -0.014 0.988822
## age_groupgr3039 0.41992 0.26203 1.603 0.109033
## age_groupgr4049 0.57387 0.27332 2.100 0.035761
## age_groupgr5059 0.56045 0.29173 1.921 0.054712
## age_group60plus 0.85999 0.36901 2.331 0.019780
## h12 0.02666 0.16221 0.164 0.869444
## d1amarried -0.30238 0.22522 -1.343 0.179407
## educ22 -0.24567 0.19830 -1.239 0.215380
## educ23 0.19186 0.23908 0.802 0.422273
## ter_in2 0.29182 0.19596 1.489 0.136440
## ter_in3 0.46892 0.19290 2.431 0.015061
## selfhealthgood 0.09294 0.15691 0.592 0.553638
## smosttmedium -0.12812 0.20824 -0.615 0.538401
## smosttheavy 0.39313 0.20801 1.890 0.058756
## b16a1 -0.74816 0.19823 -3.774 0.000161
## b18abad -0.33508 0.16185 -2.070 0.038427
## label11 0.57910 0.22818 2.538 0.011152
## anticam -0.02424 0.21069 -0.115 0.908400
## noEnvi2 -0.15659 0.20589 -0.761 0.446936
## reasons -0.04686 0.03680 -1.274 0.202831
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1081.6 on 787 degrees of freedom
## Residual deviance: 1004.7 on 768 degrees of freedom
## (35 observations deleted due to missingness)
## AIC: 1044.7
##
## Number of Fisher Scoring iterations: 4
##
## Next step: Use 'setx' methodx.out1 <- setx(z.out1, reasons = 0, age_group = "gr5059")
s.out1 <- sim(z.out1, x = x.out1)
summary(s.out1)##
## sim x :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.3302158 0.06092046 0.3285895 0.2204554 0.4599791
## pv
## 0 1
## [1,] 0.661 0.339plot(s.out2)