data(nigeria)
## Define design parameters
p <- 2/3 # probability of answering honestly in Forced Response Design
p1 <- 1/6 # probability of forced 'yes'
p0 <- 1/6 # probability of forced 'no'
## run three chains with overdispersed starting values
set.seed(1)
## starting values constructed from MLE model
mle.estimates <- rrreg(rr.q1 ~ cov.asset.index + cov.married +
I(cov.age/10) + I((cov.age/10)^2) + cov.education + cov.female,
data = nigeria,
p = p, p1 = p1, p0 = p0,
design = "forced-known")
# \donttest{
library(MASS)
draws <- mvrnorm(n = 3, mu = coef(mle.estimates),
Sigma = vcov(mle.estimates) * 9)
## run three chains
bayes.1 <- rrreg.bayes(rr.q1 ~ cov.asset.index + cov.married +
I(cov.age/10) + I((cov.age/10)^2) + cov.education + cov.female,
data = nigeria, p = p, p1 = p1, p0 = p0,
beta.tune = .0001, beta.start = draws[1,],
design = "forced-known")
bayes.2 <- rrreg.bayes(rr.q1 ~ cov.asset.index + cov.married +
I(cov.age/10) + I((cov.age/10)^2) + cov.education + cov.female,
data = nigeria, p = p, p1 = p1, p0 = p0,
beta.tune = .0001, beta.start = draws[2,],
design = "forced-known")
bayes.3 <- rrreg.bayes(rr.q1 ~ cov.asset.index + cov.married +
I(cov.age/10) + I((cov.age/10)^2) + cov.education + cov.female,
data = nigeria, p = p, p1 = p1, p0 = p0,
beta.tune = .0001, beta.start = draws[3,],
design = "forced-known")
bayes <- as.list(bayes.1, bayes.2, bayes.3)
summary(bayes)
# }
Run the code above in your browser using DataLab