describle data
outcome: WTP ()
It results in a likelihood maximized when a parameter is extremely large, and causes trouble with ordinary maximum likelihood approached.) Another option is to use Bayesian methods. Here we focus on Markov chain Monte Carlo (MCMC) approaches to Bayesian analysis.
steps to solve
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/dataR5.dta")
r$smostt=as.factor(r$smostt)
r1 <- subset(r)
attach(r1)
r1$wtp[cost_inc <=180000] <- 0
r1$wtp[cost_inc >180000] <- 0
r1$wtp[is.na(r1$wtp)] <- 1
r1$b16a[b16a == 5] <- 0
attach(r1)## The following objects are masked from r1 (pos = 3):
##
## a, a1, advice, ag, age_group, anticam, b1, b10, b10a1, b10a2,
## b10a3, b10a4, b10a5, b10a6, b10a7, b11, b11a, b11a2, b11a3,
## b12, b12a, b13, b14, b15, b16, b16a, b17, b18, b18a, b1a,
## b1a1, b1a2, b1a3, b2, b2a1, b3, b4, b5, b5a, b6, b6_123, b6a,
## b6a1, b6a2, b6a3, b6a4, b6a5, b6a6, b6a7, b6a7a, b7, b8, b9,
## br, branch, branch1, branch2, branch3, c1, c2, c23a, c3, c4,
## c5, c6, c6a111, c7, c8, c9, clog, clog1, 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, ict1a, ict2a,
## itc1, itc2, l1, l2, l3, l4, l5, l6, label_1, label1, label1a,
## logb7, logC, logitCost, 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, noE, noEnvi2,
## occup1, poli4, poli4a, policy, policy_a, policyeffect,
## reasons, reasons1, s1, Screening, selfhealth, SH1, smostt,
## taxIn, 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.03892944r1$label2[l5==1] <- 1
r1$label2[l5>1] <- 0
r1$freeEn[c9==1] <- 1
r1$freeEn[c9>1] <- 0
newdata2=r1
attach(newdata2 )## The following objects are masked from r1 (pos = 3):
##
## a, a1, advice, ag, age_group, anticam, b1, b10, b10a1, b10a2,
## b10a3, b10a4, b10a5, b10a6, b10a7, b11, b11a, b11a2, b11a3,
## b12, b12a, b13, b14, b15, b16, b16a, b17, b18, b18a, b1a,
## b1a1, b1a2, b1a3, b2, b2a1, b3, b4, b5, b5a, b6, b6_123, b6a,
## b6a1, b6a2, b6a3, b6a4, b6a5, b6a6, b6a7, b6a7a, b7, b8, b9,
## br, branch, branch1, branch2, branch3, c1, c2, c23a, c3, c4,
## c5, c6, c6a111, c7, c8, c9, clog, clog1, 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, ict1a, ict2a,
## itc1, itc2, l1, l2, l3, l4, l5, l6, label_1, label1, label1a,
## logb7, logC, logitCost, 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, noE, noEnvi2,
## occup1, poli4, poli4a, policy, policy_a, policyeffect,
## reasons, reasons1, s1, Screening, selfhealth, SH1, smostt,
## taxIn, 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):
##
## a, a1, advice, ag, age_group, anticam, b1, b10, b10a1, b10a2,
## b10a3, b10a4, b10a5, b10a6, b10a7, b11, b11a, b11a2, b11a3,
## b12, b12a, b13, b14, b15, b16, b16a, b17, b18, b18a, b1a,
## b1a1, b1a2, b1a3, b2, b2a1, b3, b4, b5, b5a, b6, b6_123, b6a,
## b6a1, b6a2, b6a3, b6a4, b6a5, b6a6, b6a7, b6a7a, b7, b8, b9,
## br, branch, branch1, branch2, branch3, c1, c2, c23a, c3, c4,
## c5, c6, c6a111, c7, c8, c9, clog, clog1, 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, ict1a, ict2a,
## itc1, itc2, l1, l2, l3, l4, l5, l6, label_1, label1, label1a,
## logb7, logC, logitCost, 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, noE, noEnvi2,
## occup1, poli4, poli4a, policy, policy_a, policyeffect,
## reasons, reasons1, s1, Screening, selfhealth, SH1, smostt,
## taxIn, ter_fa1, ter_in, tertile_fa, tertile_indi, test,
## unitsdiffi1, var242, w1, w2, w3, w4newdata2$c1=as.factor(newdata2$c1)
newdata2$h1=as.factor(newdata2$h1)
newdata2$noEnvi2=as.factor(newdata2$noEnvi2)
newdata2$policy=as.factor(newdata2$policy)
newdata2$educ2=as.factor(newdata2$educ2)
newdata2$ter_in=as.factor(newdata2$ter_in)
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)
newdata2$policy_a=as.factor(newdata2$policy_a)
attach(newdata2 )## The following objects are masked from newdata2 (pos = 3):
##
## a, a1, advice, ag, age_group, anticam, b1, b10, b10a1, b10a2,
## b10a3, b10a4, b10a5, b10a6, b10a7, b11, b11a, b11a2, b11a3,
## b12, b12a, b13, b14, b15, b16, b16a, b17, b18, b18a, b1a,
## b1a1, b1a2, b1a3, b2, b2a1, b3, b4, b5, b5a, b6, b6_123, b6a,
## b6a1, b6a2, b6a3, b6a4, b6a5, b6a6, b6a7, b6a7a, b7, b8, b9,
## br, branch, branch1, branch2, branch3, c1, c2, c23a, c3, c4,
## c5, c6, c6a111, c7, c8, c9, clog, clog1, 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,
## freeEn, 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,
## ict1a, ict2a, itc1, itc2, l1, l2, l3, l4, l5, l6, label_1,
## label1, label1a, label2, logb7, logC, logitCost, 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, noE, noEnvi2, occup1, poli4, poli4a, policy,
## policy_a, policyeffect, reasons, reasons1, s1, Screening,
## selfhealth, SH1, smostt, taxIn, 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):
##
## a, a1, advice, ag, age_group, anticam, b1, b10, b10a1, b10a2,
## b10a3, b10a4, b10a5, b10a6, b10a7, b11, b11a, b11a2, b11a3,
## b12, b12a, b13, b14, b15, b16, b16a, b17, b18, b18a, b1a,
## b1a1, b1a2, b1a3, b2, b2a1, b3, b4, b5, b5a, b6, b6_123, b6a,
## b6a1, b6a2, b6a3, b6a4, b6a5, b6a6, b6a7, b6a7a, b7, b8, b9,
## br, branch, branch1, branch2, branch3, c1, c2, c23a, c3, c4,
## c5, c6, c6a111, c7, c8, c9, clog, clog1, 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, ict1a, ict2a,
## itc1, itc2, l1, l2, l3, l4, l5, l6, label_1, label1, label1a,
## logb7, logC, logitCost, 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, noE, noEnvi2,
## occup1, poli4, poli4a, policy, policy_a, policyeffect,
## reasons, reasons1, s1, Screening, selfhealth, SH1, smostt,
## taxIn, 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):
##
## a, a1, advice, ag, age_group, anticam, b1, b10, b10a1, b10a2,
## b10a3, b10a4, b10a5, b10a6, b10a7, b11, b11a, b11a2, b11a3,
## b12, b12a, b13, b14, b15, b16, b16a, b17, b18, b18a, b1a,
## b1a1, b1a2, b1a3, b2, b2a1, b3, b4, b5, b5a, b6, b6_123, b6a,
## b6a1, b6a2, b6a3, b6a4, b6a5, b6a6, b6a7, b6a7a, b7, b8, b9,
## br, branch, branch1, branch2, branch3, c1, c2, c23a, c3, c4,
## c5, c6, c6a111, c7, c8, c9, clog, clog1, 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, ict1a, ict2a,
## itc1, itc2, l1, l2, l3, l4, l5, l6, label_1, label1, label1a,
## logb7, logC, logitCost, 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, noE, noEnvi2,
## occup1, poli4, poli4a, policy, policy_a, policyeffect,
## reasons, reasons1, s1, Screening, selfhealth, SH1, smostt,
## taxIn, ter_fa1, ter_in, tertile_fa, tertile_indi, test,
## unitsdiffi1, var242, w1, w2, w3, w4newdata2=newdata2[sample(1:nrow(newdata2), 822, replace=F),]Train data
train=newdata2[1:658, ]
View(train)Test data
test=newdata2[659:822, ]
View(test)library("brms")## Warning: package 'brms' was built under R version 3.2.5## Loading required package: Rcpp## Warning: package 'Rcpp' was built under R version 3.2.5## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 3.2.5## Loading 'brms' package (version 1.5.1). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').library("caret")## Warning: package 'caret' was built under R version 3.2.5## Loading required package: latticelibrary("coda")## Warning: package 'coda' was built under R version 3.2.5library("rstan", lib.loc="~/R/win-library/3.2")## Warning: package 'rstan' was built under R version 3.2.5## Loading required package: StanHeaders## Warning: package 'StanHeaders' was built under R version 3.2.5## rstan (Version 2.14.1, packaged: 2016-12-28 14:55:41 UTC, GitRev: 5fa1e80eb817)## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())##
## Attaching package: 'rstan'## The following object is masked from 'package:coda':
##
## traceplotoptions(mc.cores = parallel::detectCores())prior=get_prior(formula=wtp~age_group*educ2+h1+ d1a+ter_in+selfhealth+smostt
+b18a+b6a+label2+freeEn, family="bernoulli", data=train)## Warning: Rows containing NAs were excluded from the modelset.seed(1234)
baybin=brm(formula=wtp~age_group*educ2+h1+ d1a+ter_in+selfhealth+smostt
+b18a+b6a+label2+freeEn, family="bernoulli", data=train, chains=5, iter=2000, warmup=1000, prior=prior)## Warning: Rows containing NAs were excluded from the model## Compiling the C++ model## Start samplingtest interaction baybin=brm(formula=wtp~age_group*educ2, family=“bernoulli”, data=newdata2)
summary(baybin)## Warning: The model has not converged (some Rhats are > 1.1). Do not analyse the results!
## We recommend running more iterations and/or setting stronger priors.## Family: bernoulli (logit)
## Formula: wtp ~ age_group * educ2 + h1 + d1a + ter_in + selfhealth + smostt + b18a + b6a + label2 + freeEn
## Data: train (Number of observations: 342)
## Samples: 5 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 5000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## Intercept 0.44 1.02 -1.26 2.53
## age_groupgr3039 -0.90 1.23 -3.68 0.99
## age_groupgr4049 0.40 1.10 -2.03 2.35
## age_groupgr5059 0.05 1.09 -2.41 1.86
## age_group60plus 0.02 1.34 -2.88 2.23
## educ22 -0.70 0.91 -3.00 0.79
## educ23 0.10 0.97 -1.99 1.84
## h12 0.07 0.29 -0.50 0.62
## d1amarried -0.25 0.34 -1.04 0.40
## ter_in2 0.85 0.32 0.20 1.43
## ter_in3 0.91 0.30 0.35 1.55
## selfhealthgood -0.31 0.25 -0.86 0.07
## smosttmedium -0.07 0.30 -0.69 0.53
## smosttheavy 0.06 0.31 -0.48 0.73
## b18abad -0.42 0.28 -0.95 0.21
## b6ayes -0.29 0.40 -1.05 0.47
## label2 -0.54 0.24 -1.07 -0.11
## freeEn -0.22 0.25 -0.79 0.25
## age_groupgr3039:educ22 0.87 1.20 -0.96 3.63
## age_groupgr4049:educ22 0.15 1.04 -1.66 2.59
## age_groupgr5059:educ22 -0.04 1.12 -1.90 2.61
## age_group60plus:educ22 0.43 1.63 -2.34 3.93
## age_groupgr3039:educ23 0.25 1.24 -1.75 3.12
## age_groupgr4049:educ23 -0.68 1.28 -3.07 1.87
## age_groupgr5059:educ23 2.08 1.82 -0.66 6.25
## age_group60plus:educ23 -608887.20 465490.73 -1703324.64 -98536.94
## Eff.Sample Rhat
## Intercept 4 2.11
## age_groupgr3039 4 1.88
## age_groupgr4049 4 2.10
## age_groupgr5059 4 2.09
## age_group60plus 4 1.86
## educ22 4 2.03
## educ23 4 2.12
## h12 18 1.35
## d1amarried 6 1.34
## ter_in2 21 1.29
## ter_in3 26 1.20
## selfhealthgood 11 1.30
## smosttmedium 8 1.30
## smosttheavy 22 1.30
## b18abad 15 1.31
## b6ayes 23 1.19
## label2 10 1.19
## freeEn 66 1.09
## age_groupgr3039:educ22 5 1.66
## age_groupgr4049:educ22 6 1.74
## age_groupgr5059:educ22 5 1.63
## age_group60plus:educ22 5 1.69
## age_groupgr3039:educ23 5 1.65
## age_groupgr4049:educ23 4 1.80
## age_groupgr5059:educ23 12 1.30
## age_group60plus:educ23 3 4.26
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).WAIC(baybin)## Warning: 2 (0.6%) p_waic estimates greater than 0.4.
## We recommend trying loo() instead.## WAIC SE
## 493.66 14.37Notes: If you wanna report OR and 95%CI, you might need the exp funtion (using in R or in excel)
baybin$fit## Inference for Stan model: bernoulli(logit) brms-model.
## 5 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=5000.
##
## mean se_mean sd 2.5%
## b_Intercept 0.44 0.50 1.02 -1.26
## b_age_groupgr3039 -0.90 0.59 1.23 -3.68
## b_age_groupgr4049 0.40 0.55 1.10 -2.03
## b_age_groupgr5059 0.05 0.57 1.09 -2.41
## b_age_group60plus 0.02 0.64 1.34 -2.88
## b_educ22 -0.70 0.44 0.91 -3.00
## b_educ23 0.10 0.48 0.97 -1.99
## b_h12 0.07 0.07 0.29 -0.50
## b_d1amarried -0.25 0.14 0.34 -1.04
## b_ter_in2 0.85 0.07 0.32 0.20
## b_ter_in3 0.91 0.06 0.30 0.35
## b_selfhealthgood -0.31 0.08 0.25 -0.86
## b_smosttmedium -0.07 0.11 0.30 -0.69
## b_smosttheavy 0.06 0.07 0.31 -0.48
## b_b18abad -0.42 0.07 0.28 -0.95
## b_b6ayes -0.29 0.08 0.40 -1.05
## b_label2 -0.54 0.08 0.24 -1.07
## b_freeEn -0.22 0.03 0.25 -0.79
## b_age_groupgr3039:educ22 0.87 0.55 1.20 -0.96
## b_age_groupgr4049:educ22 0.15 0.44 1.04 -1.66
## b_age_groupgr5059:educ22 -0.04 0.51 1.12 -1.90
## b_age_group60plus:educ22 0.43 0.71 1.63 -2.34
## b_age_groupgr3039:educ23 0.25 0.57 1.24 -1.75
## b_age_groupgr4049:educ23 -0.68 0.62 1.28 -3.07
## b_age_groupgr5059:educ23 2.08 0.52 1.82 -0.66
## b_age_group60plus:educ23 -608887.20 270142.76 465490.73 -1703324.64
## lp__ -230.96 1.44 4.21 -240.07
## 25% 50% 75% 97.5% n_eff
## b_Intercept -0.29 0.29 1.17 2.53 4
## b_age_groupgr3039 -1.63 -0.65 -0.02 0.99 4
## b_age_groupgr4049 -0.21 0.52 1.01 2.35 4
## b_age_groupgr5059 -0.62 0.20 0.74 1.86 4
## b_age_group60plus -0.83 0.25 0.91 2.23 4
## b_educ22 -1.21 -0.47 -0.13 0.79 4
## b_educ23 -0.53 0.19 0.71 1.84 4
## b_h12 -0.13 0.07 0.27 0.62 18
## b_d1amarried -0.47 -0.22 -0.01 0.40 6
## b_ter_in2 0.62 0.89 1.06 1.43 21
## b_ter_in3 0.71 0.88 1.11 1.55 26
## b_selfhealthgood -0.47 -0.28 -0.11 0.07 11
## b_smosttmedium -0.23 -0.08 0.10 0.53 8
## b_smosttheavy -0.16 0.00 0.26 0.73 22
## b_b18abad -0.59 -0.45 -0.28 0.21 15
## b_b6ayes -0.58 -0.31 -0.01 0.47 23
## b_label2 -0.69 -0.52 -0.37 -0.11 10
## b_freeEn -0.36 -0.21 -0.06 0.25 66
## b_age_groupgr3039:educ22 0.00 0.58 1.55 3.63 5
## b_age_groupgr4049:educ22 -0.51 0.06 0.63 2.59 6
## b_age_groupgr5059:educ22 -0.81 -0.25 0.61 2.61 5
## b_age_group60plus:educ22 -0.82 0.30 1.44 3.93 5
## b_age_groupgr3039:educ23 -0.62 0.02 1.04 3.12 5
## b_age_groupgr4049:educ23 -1.49 -0.73 0.15 1.87 4
## b_age_groupgr5059:educ23 0.82 1.83 3.12 6.25 12
## b_age_group60plus:educ23 -858985.71 -429884.84 -267820.54 -98536.94 3
## lp__ -233.54 -230.67 -227.77 -224.30 8
## Rhat
## b_Intercept 2.11
## b_age_groupgr3039 1.88
## b_age_groupgr4049 2.10
## b_age_groupgr5059 2.09
## b_age_group60plus 1.86
## b_educ22 2.03
## b_educ23 2.12
## b_h12 1.35
## b_d1amarried 1.34
## b_ter_in2 1.29
## b_ter_in3 1.20
## b_selfhealthgood 1.30
## b_smosttmedium 1.30
## b_smosttheavy 1.30
## b_b18abad 1.31
## b_b6ayes 1.19
## b_label2 1.19
## b_freeEn 1.09
## b_age_groupgr3039:educ22 1.66
## b_age_groupgr4049:educ22 1.74
## b_age_groupgr5059:educ22 1.63
## b_age_group60plus:educ22 1.69
## b_age_groupgr3039:educ23 1.65
## b_age_groupgr4049:educ23 1.80
## b_age_groupgr5059:educ23 1.30
## b_age_group60plus:educ23 4.26
## lp__ 1.22
##
## Samples were drawn using NUTS(diag_e) at Mon Mar 27 16:41:42 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).next steps, ww will check the accuracy of model using ROC
prob=as.data.frame(predict(baybin,test,type="response"))## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs produced
## Warning in rbinom(length(eta), size = 1, prob = ilink(eta, draws$f$link)):
## NAs producedlibrary(e1071)## Warning: package 'e1071' was built under R version 3.2.5pred=ifelse(prob$Estimate> 0.5, 1, 0)confusionMatrix(data=pred,reference=test$wtp)## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 46 23
## 1 14 15
##
## Accuracy : 0.6224
## 95% CI : (0.5188, 0.7184)
## No Information Rate : 0.6122
## P-Value [Acc > NIR] : 0.4618
##
## Kappa : 0.1687
## Mcnemar's Test P-Value : 0.1884
##
## Sensitivity : 0.7667
## Specificity : 0.3947
## Pos Pred Value : 0.6667
## Neg Pred Value : 0.5172
## Prevalence : 0.6122
## Detection Rate : 0.4694
## Detection Prevalence : 0.7041
## Balanced Accuracy : 0.5807
##
## 'Positive' Class : 0
##