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)
newdata2$freeEn=as.factor(newdata2$freeEn)
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, w4`` newdata2=newdata2[sample(1:nrow(newdata2), 822, replace=F),]
###Create new data frame (90%-10%) (N=822)
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=newdata2)## 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=newdata2, chains=5, iter=2000, warmup=1000, prior=prior)## Warning: Rows containing NAs were excluded from the model## Compiling the C++ model## Start samplingsummary(baybin)## Family: bernoulli (logit)
## Formula: wtp ~ age_group + educ2 + h1 + d1a + ter_in + selfhealth + smostt + b18a + b6a * label2 + freeEn
## Data: newdata2 (Number of observations: 439)
## 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 Eff.Sample Rhat
## Intercept 0.51 0.43 -0.32 1.39 5000 1
## age_groupgr3039 -0.05 0.37 -0.77 0.68 5000 1
## age_groupgr4049 0.44 0.37 -0.27 1.19 3692 1
## age_groupgr5059 0.17 0.39 -0.59 0.95 3807 1
## age_group60plus 0.26 0.51 -0.74 1.26 5000 1
## educ22 -0.60 0.27 -1.16 -0.08 5000 1
## educ23 -0.10 0.33 -0.74 0.54 5000 1
## h12 0.03 0.23 -0.42 0.50 5000 1
## d1amarried -0.41 0.31 -1.00 0.20 3955 1
## ter_in2 0.55 0.27 0.02 1.10 5000 1
## ter_in3 0.70 0.27 0.18 1.24 5000 1
## selfhealthgood -0.32 0.22 -0.73 0.12 5000 1
## smosttmedium -0.18 0.29 -0.75 0.40 5000 1
## smosttheavy 0.16 0.29 -0.41 0.74 5000 1
## b18abad -0.30 0.22 -0.75 0.14 5000 1
## b6ayes -1.12 0.78 -2.76 0.29 4576 1
## label2 -0.46 0.23 -0.92 -0.02 5000 1
## freeEn1 -0.16 0.27 -0.67 0.37 5000 1
## b6ayes:label2 0.75 0.85 -0.87 2.50 4408 1
##
## 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).test interaction baybin=brm(formula=wtp~age_group*educ2, family=“bernoulli”, data=newdata2)
`` summary(baybin)
WAIC(baybin) ```
Notes: 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% 25% 50% 75%
## b_Intercept 0.51 0.01 0.43 -0.32 0.23 0.51 0.80
## b_age_groupgr3039 -0.05 0.01 0.37 -0.77 -0.31 -0.06 0.20
## b_age_groupgr4049 0.44 0.01 0.37 -0.27 0.19 0.44 0.69
## b_age_groupgr5059 0.17 0.01 0.39 -0.59 -0.09 0.17 0.43
## b_age_group60plus 0.26 0.01 0.51 -0.74 -0.08 0.27 0.59
## b_educ22 -0.60 0.00 0.27 -1.16 -0.79 -0.60 -0.42
## b_educ23 -0.10 0.00 0.33 -0.74 -0.32 -0.10 0.12
## b_h12 0.03 0.00 0.23 -0.42 -0.12 0.03 0.18
## b_d1amarried -0.41 0.00 0.31 -1.00 -0.62 -0.41 -0.20
## b_ter_in2 0.55 0.00 0.27 0.02 0.36 0.54 0.72
## b_ter_in3 0.70 0.00 0.27 0.18 0.51 0.70 0.88
## b_selfhealthgood -0.32 0.00 0.22 -0.73 -0.46 -0.32 -0.17
## b_smosttmedium -0.18 0.00 0.29 -0.75 -0.37 -0.18 0.01
## b_smosttheavy 0.16 0.00 0.29 -0.41 -0.04 0.16 0.35
## b_b18abad -0.30 0.00 0.22 -0.75 -0.45 -0.30 -0.16
## b_b6ayes -1.12 0.01 0.78 -2.76 -1.62 -1.08 -0.58
## b_label2 -0.46 0.00 0.23 -0.92 -0.62 -0.46 -0.30
## b_freeEn1 -0.16 0.00 0.27 -0.67 -0.34 -0.16 0.03
## b_b6ayes:label2 0.75 0.01 0.85 -0.87 0.16 0.71 1.32
## lp__ -292.27 0.07 3.11 -299.21 -294.14 -291.95 -290.06
## 97.5% n_eff Rhat
## b_Intercept 1.39 5000 1
## b_age_groupgr3039 0.68 5000 1
## b_age_groupgr4049 1.19 3692 1
## b_age_groupgr5059 0.95 3807 1
## b_age_group60plus 1.26 5000 1
## b_educ22 -0.08 5000 1
## b_educ23 0.54 5000 1
## b_h12 0.50 5000 1
## b_d1amarried 0.20 3955 1
## b_ter_in2 1.10 5000 1
## b_ter_in3 1.24 5000 1
## b_selfhealthgood 0.12 5000 1
## b_smosttmedium 0.40 5000 1
## b_smosttheavy 0.74 5000 1
## b_b18abad 0.14 5000 1
## b_b6ayes 0.29 4576 1
## b_label2 -0.02 5000 1
## b_freeEn1 0.37 5000 1
## b_b6ayes:label2 2.50 4408 1
## lp__ -287.18 1943 1
##
## Samples were drawn using NUTS(diag_e) at Mon Mar 27 23:18:05 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"))
library(e1071)
pred=ifelse(prob$Estimate> 0.5, 1, 0)
`` confusionMatrix(data=pred,reference=test$wtp) ```
Calculate Bayes factors (BF)
BAYES FACTOR BF_{10} LABEL > 100 Extreme evidence for H1 30 - 100 Very strong evidence for H1 10 - 30 Strong evidence for H1 3 - 10 Moderate evidence for H1 1 - 3 Anecdotal evidence for H1 1 No evidence 1/3 - 1 Anecdotal evidence for H0 1/3 - 1/10 Moderate evidence for H0 1/10 - 1/30 Strong evidence for H0 1/30 - 1/100 Very strong evidence for H0 < 1/100 Extreme evidence for H0
hypothesis(baybin,"exp(age_groupgr3039 )>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(age_groupgr3... > 0 0.02 0.4 -0.48 Inf 0.8
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(age_groupgr4049 )>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(age_groupgr4... > 0 0.67 0.65 -0.14 Inf 7.47
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(age_groupgr5059 )>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(age_groupgr5... > 0 0.28 0.52 -0.38 Inf 1.99
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(age_group60plus)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(age_group60p... > 0 0.48 0.8 -0.44 Inf 2.32
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(educ23)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(educ23))-(1) > 0 -0.05 0.32 -0.47 Inf 0.6
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(educ22)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(educ22))-(1) > 0 -0.43 0.16 -0.66 Inf 0.01
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(d1amarried)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(d1amarried))... > 0 -0.3 0.22 -0.6 Inf 0.1
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(ter_in2)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(ter_in2))-(1) > 0 0.79 0.5 0.12 Inf 46.62 *
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(ter_in3)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(ter_in3))-(1) > 0 1.08 0.57 0.29 Inf 249 *
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(selfhealthgood)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(selfhealthgo... > 0 -0.25 0.17 -0.49 Inf 0.08
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(smosttmedium)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(smosttmedium... > 0 -0.13 0.26 -0.48 Inf 0.36
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(smosttheavy)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(smosttheavy)... > 0 0.22 0.37 -0.28 Inf 2.39
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(b18abad)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(b18abad))-(1) > 0 -0.24 0.17 -0.5 Inf 0.1
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(b6ayes)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(b6ayes))-(1) > 0 -0.57 0.35 -0.92 Inf 0.07
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(h12)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(h12))-(1) > 0 0.06 0.25 -0.3 Inf 1.24
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(label2)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(label2))-(1) > 0 -0.35 0.15 -0.57 Inf 0.02
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(freeEn1)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(freeEn1))-(1) > 0 -0.12 0.24 -0.45 Inf 0.39
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.hypothesis(baybin,"exp(b6ayes:label2)>1",alpha=0.05)## Hypothesis Tests for class b:
## Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(b6ayes:label... > 0 2.07 3.46 -0.46 Inf 4.28
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.