context("Tests Bayesian regression code (rrreg.bayes)")
## 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'
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")
library(MASS)
draws <- mvrnorm(n = 3, mu = coef(mle.estimates),
Sigma = vcov(mle.estimates) * 9)
test_that("basic code runs with all defaults", {
bayes <- 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,
design = "forced-known")
expect_is(bayes, "rrreg.bayes")
expect_equal(ncol(bayes$beta), 7)
})
test_that("basic code runs varying the burnin / thin / n.draws options", {
expect_is(
bayes <- 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,
n.draws = 250, burnin = 100,
design = "forced-known"),
"rrreg.bayes"
)
expect_is(
bayes <- 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,
n.draws = 250, burnin = 100, thin = 10,
design = "forced-known"),
"rrreg.bayes"
)
})
test_that("basic code runs setting customized Metropolis options", {
expect_is(
bayes <- 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,],
n.draws = 500, burnin = 250, thin = 1,
design = "forced-known"),
"rrreg.bayes"
)
expect_is(
bayes <- 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,], beta.mu0 = rep(0, 7),
beta.A0 = diag(7) * 1,
n.draws = 500, burnin = 250, thin = 1,
design = "forced-known"),
"rrreg.bayes"
)
})
test_that("standard functions coef, sd.rrreg.bayes, summary work", {
bayes <- 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,
n.draws = 500, burnin = 250, thin = 1,
design = "forced-known")
expect_equal(
length(coef(bayes)$beta), 7
)
expect_is(
coef(bayes)$beta, "numeric"
)
expect_is(
vcov(bayes), "matrix"
)
expect_equal(
length(sd.rrreg.bayes(bayes)$beta), 7
)
expect_output(
print(bayes), "Randomized Response Technique Bayesian Regression"
)
expect_output(
print(summary(bayes)), "Individual-level predictors"
)
expect_is(
as.list(bayes), "rrreg.bayes.list"
)
})
test_that("doing multiple chains works as well as coef, sd.rrreg.bayes, summary work", {
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,],
n.draws = 500, burnin = 250, thin = 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,],
n.draws = 500, burnin = 250, thin = 1,
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,],
n.draws = 500, burnin = 250, thin = 1,
design = "forced-known")
bayes <- as.list(bayes.1, bayes.2, bayes.3)
expect_is(bayes, "rrreg.bayes.list")
expect_equal(
length(coef(bayes)$beta), 7
)
expect_is(
coef(bayes)$beta, "numeric"
)
expect_is(
vcov(bayes), "matrix"
)
expect_equal(
length(sd.rrreg.bayes.list(bayes)$beta), 7
)
expect_output(
print(bayes), "Randomized Response Technique Bayesian Regression"
)
summ <- summary(bayes)
expect_output(
print(summ), "Individual-level predictors"
)
expect_output(
print(summ), "Gelman-Rubin statistics"
)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.