Untitled
Thang, TB
21 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")
r$d4a1=log10(r$d4+1)
attach(r)
qqnorm(d4a1)
qqline(d4a1)
r1=subset(r)
attach(r1)## The following objects are masked from r:
##
## 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, d4a1, 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, w4r1$c_gr[cost_inc>60000 ] <- 1
r1$c_gr[cost_inc<=60000 ] <- 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, d4a1, 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
##
## The following objects are masked from r:
##
## 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, d4a1, 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, w4newdata2=r1define 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 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, c_gr, c1, c2, c23a, c3, c4,
## c5, c6, c7, c8, c9, COST, cost_inc, cost1, costincrease, ct,
## d1, d10, d11, d1a, d2, d3, d3a, d4, d4a1, 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## 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, d4a1, 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## The following objects are masked from r:
##
## 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, d4a1, 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(c_gr~age_group+c1+h1+ d1a+educ2+selfhealth+smostt+b16a
+b6a+b18a+label1+anticam+noEnvi2+ ter_in+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
## -2.7712 -1.4960 -1.0215 -2.7903
## age_group60plus c12 h12 d1amarried
## -0.6490 1.0502 2.8413 1.4627
## educ22 educ23 selfhealthgood smosttmedium
## 1.4546 1.0740 0.4692 0.3326
## smosttheavy b16a1 b6ayes b18abad
## 0.3911 3.7682 -0.9156 -0.3535
## label11 anticam noEnvi2 ter_in2
## 0.2934 -1.8971 -0.5085 1.1664
## ter_in3 reasons
## 0.1550 0.7990z.out$heidel.diag()##
## Stationarity start p-value
## test iteration
## (Intercept) failed NA 0.0388
## age_groupgr3039 passed 1201 0.5212
## age_groupgr4049 passed 1201 0.1144
## age_groupgr5059 passed 1201 0.0865
## age_group60plus passed 1 0.4353
## c12 passed 1 0.7717
## h12 passed 601 0.8572
## d1amarried passed 1 0.2780
## educ22 passed 1 0.8846
## educ23 passed 1 0.8112
## selfhealthgood passed 1201 0.1161
## smosttmedium passed 1 0.0841
## smosttheavy passed 1 0.1580
## b16a1 passed 601 0.1160
## b6ayes passed 1 0.1993
## b18abad passed 1 0.4297
## label11 passed 1 0.1673
## anticam passed 1 0.2021
## noEnvi2 passed 1 0.8431
## ter_in2 passed 2401 0.1572
## ter_in3 passed 1 0.8336
## reasons passed 1 0.3192
##
## Halfwidth Mean Halfwidth
## test
## (Intercept) <NA> NA NA
## age_groupgr3039 failed 0.3830 0.07504
## age_groupgr4049 failed -0.5310 0.07349
## age_groupgr5059 failed -0.3447 0.06780
## age_group60plus failed -0.2392 0.09369
## c12 failed -0.0111 0.07086
## h12 passed -0.5308 0.03795
## d1amarried failed 0.1834 0.05208
## educ22 failed 0.0418 0.04474
## educ23 passed 0.9752 0.06141
## selfhealthgood passed 0.6022 0.03954
## smosttmedium failed 0.6352 0.06874
## smosttheavy passed 0.7188 0.06232
## b16a1 failed -0.1685 0.05239
## b6ayes failed 0.5811 0.06088
## b18abad failed 0.3646 0.04209
## label11 failed -0.3452 0.05516
## anticam failed 0.1326 0.05321
## noEnvi2 failed -0.4771 0.05124
## ter_in2 failed 0.0265 0.05719
## ter_in3 failed 0.1163 0.04814
## reasons failed -0.0542 0.00856z.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) 170 184350 3746 49.2
## age_groupgr3039 470 643850 3746 172.0
## age_groupgr4049 375 408445 3746 109.0
## age_groupgr5059 270 301430 3746 80.5
## age_group60plus 495 543700 3746 145.0
## c12 615 747275 3746 199.0
## h12 335 401500 3746 107.0
## d1amarried 340 389635 3746 104.0
## educ22 775 922380 3746 246.0
## educ23 495 543700 3746 145.0
## selfhealthgood 315 359485 3746 96.0
## smosttmedium 390 441385 3746 118.0
## smosttheavy 345 399470 3746 107.0
## b16a1 345 399470 3746 107.0
## b6ayes 380 424760 3746 113.0
## b18abad 485 578130 3746 154.0
## label11 455 512015 3746 137.0
## anticam 630 702260 3746 187.0
## noEnvi2 335 375110 3746 100.0
## ter_in2 500 557960 3746 149.0
## ter_in3 555 678285 3746 181.0
## reasons 380 456350 3746 122.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.77180 0.55243 0.0071318 0.052459
## age_groupgr3039 0.31587 0.43145 0.0055700 0.042840
## age_groupgr4049 -0.56515 0.37632 0.0048582 0.033686
## age_groupgr5059 -0.39135 0.37747 0.0048731 0.031570
## age_group60plus -0.23917 0.52388 0.0067633 0.047803
## c12 -0.01106 0.42347 0.0054670 0.036154
## h12 -0.50764 0.22823 0.0029464 0.020753
## d1amarried 0.18344 0.30924 0.0039923 0.026573
## educ22 0.04179 0.28543 0.0036849 0.022827
## educ23 0.97523 0.35458 0.0045775 0.031331
## selfhealthgood 0.59069 0.19639 0.0025354 0.018107
## smosttmedium 0.63520 0.36991 0.0047755 0.035070
## smosttheavy 0.71881 0.35327 0.0045607 0.031794
## b16a1 -0.14253 0.30805 0.0039769 0.027382
## b6ayes 0.58107 0.34418 0.0044433 0.031062
## b18abad 0.36457 0.23138 0.0029871 0.021476
## label11 -0.34521 0.31899 0.0041181 0.028144
## anticam 0.13259 0.29746 0.0038402 0.027148
## noEnvi2 -0.47706 0.29787 0.0038455 0.026141
## ter_in2 0.06956 0.26459 0.0034158 0.022855
## ter_in3 0.11626 0.26755 0.0034541 0.024563
## reasons -0.05417 0.05052 0.0006523 0.004369
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## (Intercept) -1.851784 -1.109394 -0.755518 -0.38342 0.22781
## age_groupgr3039 -0.425387 -0.004539 0.331925 0.63200 1.15742
## age_groupgr4049 -1.271575 -0.874574 -0.585271 -0.30914 0.27105
## age_groupgr5059 -0.985508 -0.642926 -0.456427 -0.13070 0.47027
## age_group60plus -1.405871 -0.574256 -0.167953 0.10427 0.77353
## c12 -0.829662 -0.303540 -0.001314 0.26785 0.85868
## h12 -0.941471 -0.649461 -0.541852 -0.34544 -0.03786
## d1amarried -0.387691 -0.012676 0.123960 0.40654 0.79969
## educ22 -0.543824 -0.108967 0.021060 0.20773 0.66439
## educ23 0.299498 0.758330 0.950242 1.21621 1.67741
## selfhealthgood 0.216627 0.428077 0.606331 0.72367 0.96727
## smosttmedium -0.006531 0.369628 0.646600 0.89421 1.43210
## smosttheavy 0.075317 0.472089 0.708891 0.95365 1.46318
## b16a1 -0.694005 -0.325892 -0.157790 0.04932 0.45317
## b6ayes -0.015064 0.352534 0.570911 0.79996 1.28146
## b18abad -0.100355 0.167409 0.360860 0.53919 0.78550
## label11 -1.037852 -0.562614 -0.313627 -0.09921 0.19362
## anticam -0.428885 -0.048077 0.088132 0.31226 0.77283
## noEnvi2 -1.006960 -0.689340 -0.476710 -0.29912 0.12284
## ter_in2 -0.407713 -0.066612 0.050054 0.23359 0.56989
## ter_in3 -0.460615 -0.055691 0.139947 0.29332 0.64603
## reasons -0.148409 -0.091700 -0.050796 -0.02220 0.05292
##
## 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.542197 0.0992403 0.5453919 0.3527054 0.7340393
## pv
## 0 1
## [1,] 2.775 3.225Simulating 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.5337018 0.09986704 0.5340032 0.3500591 0.7279448
## pv
## 0 1
## [1,] 2.784 3.216
##
## sim x1 :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.5595075 0.09957209 0.5674993 0.3722139 0.7538735
## pv
## 0 1
## [1,] 2.602 3.398
## fd
## mean sd 50% 2.5% 97.5%
## [1,] 0.02580566 0.02396715 0.02433592 -0.02316867 0.07299986Checking
z.out1 <- zelig(c_gr~age_group+c1+h1+ d1a+educ2+selfhealth+smostt+b16a
+b6a+b18a+label1+anticam+noEnvi2+ + ter_in+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 = c_gr ~ age_group + c1 + h1 + d1a + educ2 +
## selfhealth + smostt + b16a + b6a + b18a + label1 + anticam +
## noEnvi2 + +ter_in + reasons, data = newdata2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0147 -1.0985 0.5976 1.0370 1.7217
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.72417 0.54075 -1.339 0.18051
## age_groupgr3039 0.31486 0.36522 0.862 0.38862
## age_groupgr4049 -0.54932 0.37066 -1.482 0.13833
## age_groupgr5059 -0.38765 0.39480 -0.982 0.32616
## age_group60plus -0.28634 0.52578 -0.545 0.58603
## c12 0.04896 0.41406 0.118 0.90588
## h12 -0.47558 0.22638 -2.101 0.03565
## d1amarried 0.14970 0.31366 0.477 0.63317
## educ22 0.05405 0.28409 0.190 0.84911
## educ23 0.96089 0.35743 2.688 0.00718
## selfhealthgood 0.55955 0.21575 2.593 0.00950
## smosttmedium 0.61669 0.34441 1.791 0.07336
## smosttheavy 0.71151 0.35518 2.003 0.04515
## b16a1 -0.18684 0.31219 -0.598 0.54952
## b6ayes 0.57894 0.34230 1.691 0.09077
## b18abad 0.35937 0.22490 1.598 0.11006
## label11 -0.33895 0.32733 -1.035 0.30044
## anticam 0.17757 0.28524 0.623 0.53361
## noEnvi2 -0.47263 0.27506 -1.718 0.08574
## ter_in2 0.09375 0.26191 0.358 0.72038
## ter_in3 0.11272 0.27159 0.415 0.67812
## reasons -0.05673 0.05177 -1.096 0.27318
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 599.58 on 433 degrees of freedom
## Residual deviance: 547.87 on 412 degrees of freedom
## (389 observations deleted due to missingness)
## AIC: 591.87
##
## 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.456811 0.08935509 0.4565996 0.2894232 0.6408521
## pv
## 0 1
## [1,] 0.529 0.471plot(s.out2)
BMA
t2<-subset(newdata2, select =c(c_gr,age_group,c1,h1, d1a,educ2,selfhealth,smostt,b16a,b6a,b18a,label1,anticam,noEnvi2, ter_in))yvar = t2[,1]
xvars = t2[,-1]
library(BMA)## Warning: package 'BMA' was built under R version 3.2.5## Loading required package: survival## Warning: package 'survival' was built under R version 3.2.5## Loading required package: leaps## Warning: package 'leaps' was built under R version 3.2.5## Loading required package: robustbase## Warning: package 'robustbase' was built under R version 3.2.5##
## Attaching package: 'robustbase'## The following object is masked from 'package:survival':
##
## heart## Loading required package: inline## Warning: package 'inline' was built under R version 3.2.5## Loading required package: rrcov## Warning: package 'rrcov' was built under R version 3.2.5## Scalable Robust Estimators with High Breakdown Point (version 1.4-3)bma.search=bic.glm(xvars, yvar, strict = F,OR=20, glm.family = "binomial")## Warning in bic.glm.data.frame(xvars, yvar, strict = F, OR = 20, glm.family
## = "binomial"): There were 43 records deleted due to NA'ssummary(bma.search)##
## Call:
## bic.glm.data.frame(x = xvars, y = yvar, glm.family = "binomial", strict = F, OR = 20)
##
##
## 17 models were selected
## Best 5 models (cumulative posterior probability = 0.6529 ):
##
## p!=0 EV SD model 1 model 2
## Intercept 100 -0.297516 0.29909 -3.001e-01 -4.370e-01
## age_group.x 0.0
## .gr3039 0.000000 0.00000 . .
## .gr4049 0.000000 0.00000 . .
## .gr5059 0.000000 0.00000 . .
## .60plus 0.000000 0.00000 . .
## c1.x 3.5
## .2 -0.014469 0.09555 . .
## h1.x 20.1
## .2 -0.086927 0.19703 . .
## d1a.x 0.0
## .married 0.000000 0.00000 . .
## educ2.x 78.3
## .2 0.172912 0.24778 2.616e-01 1.679e-01
## .3 0.791167 0.49265 1.062e+00 9.277e-01
## selfhealth.x 54.0
## .good 0.296580 0.31355 . 5.041e-01
## smostt.x 0.0
## .medium 0.000000 0.00000 . .
## .heavy 0.000000 0.00000 . .
## b16a.x 0.0
## .1 0.000000 0.00000 . .
## .5 0.000000 0.00000 . .
## b6a.x 1.2
## .yes 0.002811 0.03909 . .
## b18a.x 12.2
## .bad 0.053114 0.16033 . .
## label1.x 0.0
## .1 0.000000 0.00000 . .
## anticam.x 0.0 0.000000 0.00000 . .
## noEnvi2.x 11.9 -0.058632 0.18219 . .
## ter_in.x 0.0
## .2 0.000000 0.00000 . .
## .3 0.000000 0.00000 . .
##
## nVar 1 2
## BIC -2.289e+03 -2.288e+03
## post prob 0.224 0.178
## model 3 model 4 model 5
## Intercept -1.335e-01 -6.205e-01 4.017e-02
## age_group.x
## .gr3039 . . .
## .gr4049 . . .
## .gr5059 . . .
## .60plus . . .
## c1.x
## .2 . . .
## h1.x
## .2 . . -4.746e-01
## d1a.x
## .married . . .
## educ2.x
## .2 . 3.112e-01 .
## .3 . 1.092e+00 .
## selfhealth.x
## .good 6.170e-01 . 6.277e-01
## smostt.x
## .medium . . .
## .heavy . . .
## b16a.x
## .1 . . .
## .5 . . .
## b6a.x
## .yes . . .
## b18a.x
## .bad . 4.321e-01 .
## label1.x
## .1 . . .
## anticam.x . . .
## noEnvi2.x . . .
## ter_in.x
## .2 . . .
## .3 . . .
##
## nVar 1 2 2
## BIC -2.287e+03 -2.286e+03 -2.286e+03
## post prob 0.117 0.067 0.067Binh Thang: binhthang1001@gmail.com