context("JAGS formula")
test_that("JAGS formula works", {
# check the posterior distributions with weak priors against a maximum likelihood estimates with ML
skip_on_os(("mac", "linux", "solaris")) # multivariate sampling does not exactly match across OSes
skip_on_cran()
(1)
df_all <- (
x_cont1 = (60),
x_cont2 = (60),
x_bin = (60, 1, .5),
x_fac2o = ((("A", "B"), 30), = ("A", "B")),
x_fac2t = ((("A", "B"), 30), = ("A", "B")),
x_fac3o = ((("A", "B", "C"), 20), = ("A", "B", "C")),
x_fac3t = ((("A", "B", "C"), 20), = ("A", "B", "C")),
x_fac3i = ((("A", "B", "C"), 20), = ("A", "B", "C")),
x_fac3md= ((("A", "B", "C"), 20), = ("A", "B", "C"))
)
df_all$y <- (60, 0.1, 0.5) + 0.30 * df_all$x_cont1 - 0.15 * df_all$x_cont1 * df_all$x_cont2 + 0.2 * df_all$x_bin +
(df_all$x_fac3t == "A", 0.2, (df_all$x_fac3t == "B", -0.2, 0)) +
(df_all$x_fac3o == "A", 0.2, (df_all$x_fac3o == "B", -0.2, 0))
prior_list_all <- (
"intercept" = ("normal", (0, 5)),
"x_cont1" = ("normal", (0, 1)),
"x_cont2" = ("normal", (0, 1)),
"x_cont1:x_cont2" = ("normal", (0, 1)),
"x_fac2o" = ("mcauchy", contrast = "orthonormal", (0, 1)),
"x_fac2t" = ("normal", contrast = "treatment", (0, 1)),
"x_fac3o" = ("mnormal", contrast = "orthonormal", (0, 1)),
"x_fac3t" = ("uniform", contrast = "treatment", (-2, 2)),
"x_fac3i" = ("normal", contrast = "independent", (0, 1)),
"x_fac3md" = ("mnormal", contrast = "meandif", (0, 1)),
"x_fac2t:x_fac3o" = ("mnormal", contrast = "orthonormal", (0, 2)),
"x_fac2o:x_fac3t" = ("normal", contrast = "treatment", (0, 2)),
"x_cont1:x_fac3o" = ("mnormal", contrast = "orthonormal", (0, 2)),
"x_cont1:x_fac3t" = ("normal", contrast = "treatment", (0, 2))
)
prior_list2 <- (
"sigma" = ("cauchy", (0, 1), (0, 1))
)
model_syntax <- (
"for(i in 1:N){\n",
" y[i] ~ dnorm(mu[i], 1/pow(sigma, 2))\n",
"}\n"
)
# simple linear regression ----
formula_1 <- (~ x_cont1, parameter = "mu", = df_all,"x_cont1", = ], prior_list = prior_list_all[c("intercept", "x_cont1")])
prior_list_1 <- (formula_1$prior_list, prior_list2)
model_syntax_1 <- (("model{", formula_1$formula_syntax, model_syntax, "}"), prior_list_1)
data_1 <- (formula_1$, N = (df_all), y = (df_all$y))
model_1 <- rjags::jags.model( = (model_syntax_1), inits = (prior_list_1, chains = 2, seed = 1), = data_1, n.chains = 2, quiet = )
samples_1 <- rjags::coda.samples(model = model_1, = (prior_list_1), n.iter = 5000, quiet = , progress.bar = "none")
samples_1 <- (, samples_1)
lm_1 <- ::(y ~ x_cont1, = df_all)
vdiffr::expect_doppelganger("JAGS-formula-lm-1", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfcol = (1, 3))
(samples_1,"mu_intercept"], freq = , main = "Intercept")
((x, = (lm_1)"(Intercept)"], = (lm_1)$"(Intercept)", "Std. Error"]), add = , lwd = 2)
(samples_1,"mu_x_cont1"], freq = , main = "x_cont1")
((x, = (lm_1)"x_cont1"], = (lm_1)$"x_cont1", "Std. Error"]), add = , lwd = 2)
(samples_1,"sigma"], freq = , main = "sigma")
(v = (lm_1), lwd = 3)
})
# linear regression with two continuous predictors and their interaction ----
formula_2 <- (~ x_cont1 * x_cont2, parameter = "mu", = df_all,("x_cont1", "x_cont2")], prior_list = prior_list_all[c("intercept", "x_cont1", "x_cont2", "x_cont1:x_cont2")])
prior_list_2 <- (formula_2$prior_list, prior_list2)
model_syntax_2 <- (("model{", formula_2$formula_syntax, model_syntax, "}"), prior_list_2)
data_2 <- (formula_2$, N = (df_all), y = (df_all$y))
model_2 <- rjags::jags.model( = (model_syntax_2), inits = (prior_list_2, chains = 2, seed = 1), = data_2, n.chains = 2, quiet = )
samples_2 <- rjags::coda.samples(model = model_2, = (prior_list_2), n.iter = 5000, quiet = , progress.bar = "none")
samples_2 <- (, samples_2)
lm_2 <- ::(y ~ x_cont1 * x_cont2, = df_all)
vdiffr::expect_doppelganger("JAGS-formula-lm-2", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfcol = (1, 3))
(samples_2,"mu_x_cont1"], freq = , main = "x_cont1")
((x, = (lm_2)"x_cont1"], = (lm_2)$"x_cont1", "Std. Error"]), add = , lwd = 2)
(samples_2,"mu_x_cont2"], freq = , main = "x_cont2")
((x, = (lm_2)"x_cont2"], = (lm_2)$"x_cont2", "Std. Error"]), add = , lwd = 2)
(samples_2,"mu_x_cont1__xXx__x_cont2"], freq = , main = "x_cont1:x_cont2")
((x, = (lm_2)"x_cont1:x_cont2"], = (lm_2)$"x_cont1:x_cont2", "Std. Error"]), add = , lwd = 2)
})
# linear regression with a treatment factor (2 levels) ----
formula_3 <- (~ x_fac2t, parameter = "mu", = df_all,"x_fac2t", = ], prior_list = prior_list_all[c("intercept", "x_fac2t")])
prior_list_3 <- (formula_3$prior_list, prior_list2)
model_syntax_3 <- (("model{", formula_3$formula_syntax, model_syntax, "}"), prior_list_3)
data_3 <- (formula_3$, N = (df_all), y = (df_all$y))
model_3 <- rjags::jags.model( = (model_syntax_3), inits = (prior_list_3, chains = 2, seed = 1), = data_3, n.chains = 2, quiet = )
samples_3 <- rjags::coda.samples(model = model_3, = (prior_list_3), n.iter = 5000, quiet = , progress.bar = "none")
samples_3 <- (, samples_3)
lm_3 <- ::(y ~ x_fac2t, = df_all)
vdiffr::expect_doppelganger("JAGS-formula-lm-3", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfcol = (1, 3))
(samples_3,"mu_intercept"], freq = , main = "Intercept")
((x, = (lm_3)"(Intercept)"], = (lm_3)$"(Intercept)", "Std. Error"]), add = , lwd = 2)
(samples_3,"mu_x_fac2t"], freq = , main = "x_fac2t")
((x, = (lm_3)"x_fac2tB"], = (lm_3)$"x_fac2tB", "Std. Error"]), add = , lwd = 2)
(samples_3,"sigma"], freq = , main = "sigma")
(v = (lm_3), lwd = 3)
})
# linear regression with an orthonormal factor (2 levels) ----
formula_4 <- (~ x_fac2o, parameter = "mu", = df_all,"x_fac2o", = ], prior_list = prior_list_all[c("intercept", "x_fac2o")])
prior_list_4 <- (formula_4$prior_list, prior_list2)
model_syntax_4 <- (("model{", formula_4$formula_syntax, model_syntax, "}"), prior_list_4)
data_4 <- (formula_4$, N = (df_all), y = (df_all$y))
model_4 <- rjags::jags.model( = (model_syntax_4), inits = (prior_list_4, chains = 2, seed = 1), = data_4, n.chains = 2, quiet = )
samples_4 <- rjags::coda.samples(model = model_4, = (prior_list_4), n.iter = 5000, quiet = , progress.bar = "none")
samples_4 <- (, samples_4)
df_4 <- df_all
(df_4$x_fac2o) <- ((df_4$x_fac2o))
lm_4 <- ::(y ~ x_fac2o, = df_4)
vdiffr::expect_doppelganger("JAGS-formula-lm-4", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfcol = (1, 3))
(samples_4,"mu_intercept"], freq = , main = "Intercept")
((x, = (lm_4)"(Intercept)"], = (lm_4)$"(Intercept)", "Std. Error"]), add = , lwd = 2)
(samples_4,"mu_x_fac2o"], freq = , main = "x_fac2o")
((x, = (lm_4)"x_fac2o1"], = (lm_4)$"x_fac2o1", "Std. Error"]), add = , lwd = 2)
(samples_4,"sigma"], freq = , main = "sigma")
(v = (lm_4), lwd = 3)
})
# linear regression with a treatment factor (3 levels) ----
formula_5 <- (~ x_fac3t, parameter = "mu", = df_all,"x_fac3t", = ], prior_list = prior_list_all[c("intercept", "x_fac3t")])
prior_list_5 <- (formula_5$prior_list, prior_list2)
model_syntax_5 <- (("model{", formula_5$formula_syntax, model_syntax, "}"), prior_list_5)
data_5 <- (formula_5$, N = (df_all), y = (df_all$y))
model_5 <- rjags::jags.model( = (model_syntax_5), inits = (prior_list_5, chains = 2, seed = 1), = data_5, n.chains = 2, quiet = )
samples_5 <- rjags::coda.samples(model = model_5, = (prior_list_5), n.iter = 5000, quiet = , progress.bar = "none")
samples_5 <- (, samples_5)
lm_5 <- ::(y ~ x_fac3t, = df_all)
vdiffr::expect_doppelganger("JAGS-formula-lm-5", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfcol = (1, 3))
(samples_5,"mu_intercept"], freq = , main = "Intercept")
((x, = (lm_5)"(Intercept)"], = (lm_5)$"(Intercept)", "Std. Error"]), add = , lwd = 2)
(samples_5,"mu_x_fac3t[1]"], freq = , main = "x_fac3t[1]")
((x, = (lm_5)"x_fac3tB"], = (lm_5)$"x_fac3tB", "Std. Error"]), add = , lwd = 2)
(samples_5,"mu_x_fac3t[2]"], freq = , main = "x_fac3t[2]")
((x, = (lm_5)"x_fac3tC"], = (lm_5)$"x_fac3tC", "Std. Error"]), add = , lwd = 2)
})
# linear regression with an orthonormal factor (3 levels) ----
formula_6 <- (~ x_fac3o, parameter = "mu", = df_all,"x_fac3o", = ], prior_list = prior_list_all[c("intercept", "x_fac3o")])
prior_list_6 <- (formula_6$prior_list, prior_list2)
model_syntax_6 <- (("model{", formula_6$formula_syntax, model_syntax, "}"), prior_list_6)
data_6 <- (formula_6$, N = (df_all), y = (df_all$y))
model_6 <- rjags::jags.model( = (model_syntax_6), inits = (prior_list_6, chains = 2, seed = 1), = data_6, n.chains = 2, quiet = )
samples_6 <- rjags::coda.samples(model = model_6, = (prior_list_6), n.iter = 5000, quiet = , progress.bar = "none")
samples_6 <- (, samples_6)
df_6 <- df_all
(df_6$x_fac3o) <- ((df_6$x_fac3o))
lm_6 <- ::(y ~ x_fac3o, = df_6)
vdiffr::expect_doppelganger("JAGS-formula-lm-6", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfcol = (1, 3))
(samples_6,"mu_intercept"], freq = , main = "Intercept")
((x, = (lm_6)"(Intercept)"], = (lm_6)$"(Intercept)", "Std. Error"]), add = , lwd = 2)
(samples_6,"mu_x_fac3o[1]"], freq = , main = "x_fac3o")
((x, = (lm_6)"x_fac3o1"], = (lm_6)$"x_fac3o1", "Std. Error"]), add = , lwd = 2)
(samples_6,"mu_x_fac3o[2]"], freq = , main = "x_fac3o")
((x, = (lm_6)"x_fac3o2"], = (lm_6)$"x_fac3o2", "Std. Error"]), add = , lwd = 2)
})
# linear regression with an orthonormal interaction between factors ----
formula_7 <- (~ x_fac2t * x_fac3o, parameter = "mu", = df_all,("x_fac2t", "x_fac3o")], prior_list = prior_list_all[c("intercept", "x_fac2t", "x_fac3o", "x_fac2t:x_fac3o")])
prior_list_7 <- (formula_7$prior_list, prior_list2)
model_syntax_7 <- (("model{", formula_7$formula_syntax, model_syntax, "}"), prior_list_7)
data_7 <- (formula_7$, N = (df_all), y = (df_all$y))
model_7 <- rjags::jags.model( = (model_syntax_7), inits = (prior_list_7, chains = 2, seed = 1), = data_7, n.chains = 2, quiet = )
samples_7 <- rjags::coda.samples(model = model_7, = (prior_list_7), n.iter = 5000, quiet = , progress.bar = "none")
samples_7 <- (, samples_7)
df_7 <- df_all
(df_7$x_fac3o) <- ((df_7$x_fac3o))
lm_7 <- ::(y ~ x_fac2t * x_fac3o, = df_7)
vdiffr::expect_doppelganger("JAGS-formula-lm-7", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfrow = (2, 3))
(samples_7,"mu_intercept"], freq = , main = "Intercept")
((x, = (lm_7)"(Intercept)"], = (lm_7)$"(Intercept)", "Std. Error"]), add = , lwd = 2)
(samples_7,"mu_x_fac3o[1]"], freq = , main = "x_fac3o")
((x, = (lm_7)"x_fac3o1"], = (lm_7)$"x_fac3o1", "Std. Error"]), add = , lwd = 2)
(samples_7,"mu_x_fac3o[2]"], freq = , main = "x_fac3o")
((x, = (lm_7)"x_fac3o2"], = (lm_7)$"x_fac3o2", "Std. Error"]), add = , lwd = 2)
(samples_7,"mu_x_fac2t"], freq = , main = "x_fac2t")
((x, = (lm_7)"x_fac2tB"], = (lm_7)$"x_fac2tB", "Std. Error"]), add = , lwd = 2)
(samples_7,"mu_x_fac2t__xXx__x_fac3o[1]"], freq = , main = "x_fac2t:x_fac3o")
((x, = (lm_7)"x_fac2tB:x_fac3o1"], = (lm_7)$"x_fac2tB:x_fac3o1", "Std. Error"]), add = , lwd = 2)
(samples_7,"mu_x_fac2t__xXx__x_fac3o[2]"], freq = , main = "x_fac2t:x_fac3o")
((x, = (lm_7)"x_fac2tB:x_fac3o2"], = (lm_7)$"x_fac2tB:x_fac3o2", "Std. Error"]), add = , lwd = 2)
})
# linear regression with a treatment interaction between factors ----
formula_8 <- (~ x_fac2o * x_fac3t, parameter = "mu", = df_all,("x_fac2o", "x_fac3t")], prior_list = prior_list_all[c("intercept", "x_fac2o", "x_fac3t", "x_fac2o:x_fac3t")])
prior_list_8 <- (formula_8$prior_list, prior_list2)
model_syntax_8 <- (("model{", formula_8$formula_syntax, model_syntax, "}"), prior_list_8)
data_8 <- (formula_8$, N = (df_all), y = (df_all$y))
model_8 <- rjags::jags.model( = (model_syntax_8), inits = (prior_list_8, chains = 2, seed = 1), = data_8, n.chains = 2, quiet = )
samples_8 <- rjags::coda.samples(model = model_8, = (prior_list_8), n.iter = 5000, quiet = , progress.bar = "none")
samples_8 <- (, samples_8)
df_8 <- df_all
(df_8$x_fac2o) <- ((df_8$x_fac2o))
lm_8 <- ::(y ~ x_fac2o * x_fac3t, = df_8)
vdiffr::expect_doppelganger("JAGS-formula-lm-8", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfrow = (2, 3))
(samples_8,"mu_intercept"], freq = , main = "Intercept")
((x, = (lm_8)"(Intercept)"], = (lm_8)$"(Intercept)", "Std. Error"]), add = , lwd = 2)
(samples_8,"mu_x_fac3t[1]"], freq = , main = "x_fac3t")
((x, = (lm_8)"x_fac3tB"], = (lm_8)$"x_fac3tB", "Std. Error"]), add = , lwd = 2)
(samples_8,"mu_x_fac3t[2]"], freq = , main = "x_fac3t")
((x, = (lm_8)"x_fac3tC"], = (lm_8)$"x_fac3tC", "Std. Error"]), add = , lwd = 2)
(samples_8,"mu_x_fac2o"], freq = , main = "x_fac2o")
((x, = (lm_8)"x_fac2o1"], = (lm_8)$"x_fac2o1", "Std. Error"]), add = , lwd = 2)
(samples_8,"mu_x_fac2o__xXx__x_fac3t[1]"], freq = , main = "x_fac2o:fac3t")
((x, = (lm_8)"x_fac2o1:x_fac3tB"], = (lm_8)$"x_fac2o1:x_fac3tB", "Std. Error"]), add = , lwd = 2)
(samples_8,"mu_x_fac2o__xXx__x_fac3t[2]"], freq = , main = "x_fac2o:fac3t")
((x, = (lm_8)"x_fac2o1:x_fac3tC"], = (lm_8)$"x_fac2o1:x_fac3tC", "Std. Error"]), add = , lwd = 2)
})
# linear regression with an interaction between continuous variable and orthonormal factor ----
formula_9 <- (~ x_cont1 * x_fac3o , parameter = "mu", = df_all,("x_cont1", "x_fac3o")], prior_list = prior_list_all[c("intercept", "x_cont1", "x_fac3o", "x_cont1:x_fac3o")])
prior_list_9 <- (formula_9$prior_list, prior_list2)
model_syntax_9 <- (("model{", formula_9$formula_syntax, model_syntax, "}"), prior_list_9)
data_9 <- (formula_9$, N = (df_all), y = (df_all$y))
model_9 <- rjags::jags.model( = (model_syntax_9), inits = (prior_list_9, chains = 2, seed = 1), = data_9, n.chains = 2, quiet = )
samples_9 <- rjags::coda.samples(model = model_9, = (prior_list_9), n.iter = 5000, quiet = , progress.bar = "none")
samples_9 <- (, samples_9)
df_9 <- df_all
(df_9$x_fac3o) <- ((df_9$x_fac3o))
lm_9 <- ::(y ~ x_cont1 * x_fac3o, = df_9)
vdiffr::expect_doppelganger("JAGS-formula-lm-9", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfrow = (2, 3))
(samples_9,"mu_intercept"], freq = , main = "Intercept")
((x, = (lm_9)"(Intercept)"], = (lm_9)$"(Intercept)", "Std. Error"]), add = , lwd = 2)
(samples_9,"mu_x_fac3o[1]"], freq = , main = "x_fac3o")
((x, = (lm_9)"x_fac3o1"], = (lm_9)$"x_fac3o1", "Std. Error"]), add = , lwd = 2)
(samples_9,"mu_x_fac3o[2]"], freq = , main = "x_fac3o")
((x, = (lm_9)"x_fac3o2"], = (lm_9)$"x_fac3o2", "Std. Error"]), add = , lwd = 2)
(samples_9,"mu_x_cont1"], freq = , main = "x_cont1")
((x, = (lm_9)"x_cont1"], = (lm_9)$"x_cont1", "Std. Error"]), add = , lwd = 2)
(samples_9,"mu_x_cont1__xXx__x_fac3o[1]"], freq = , main = "x_cont1:x_fac3o")
((x, = (lm_9)"x_cont1:x_fac3o1"], = (lm_9)$"x_cont1:x_fac3o1", "Std. Error"]), add = , lwd = 2)
(samples_9,"mu_x_cont1__xXx__x_fac3o[2]"], freq = , main = "x_cont1:x_fac3o")
((x, = (lm_9)"x_cont1:x_fac3o2"], = (lm_9)$"x_cont1:x_fac3o2", "Std. Error"]), add = , lwd = 2)
})
# linear regression with an interaction between continuous variable and orthonormal factor ----
formula_10 <- (~ x_cont1 * x_fac3t , parameter = "mu", = df_all,("x_cont1", "x_fac3t")], prior_list = prior_list_all[c("intercept", "x_cont1", "x_fac3t", "x_cont1:x_fac3t")])
prior_list_10 <- (formula_10$prior_list, prior_list2)
model_syntax_10 <- (("model{", formula_10$formula_syntax, model_syntax, "}"), prior_list_10)
data_10 <- (formula_10$, N = (df_all), y = (df_all$y))
model_10 <- rjags::jags.model( = (model_syntax_10), inits = (prior_list_10, chains = 2, seed = 1), = data_10, n.chains = 2, quiet = )
samples_10 <- rjags::coda.samples(model = model_10, = (prior_list_10), n.iter = 5000, quiet = , progress.bar = "none")
samples_10 <- (, samples_10)
lm_10 <- ::(y ~ x_cont1 * x_fac3t, = df_all)
vdiffr::expect_doppelganger("JAGS-formula-lm-10", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfrow = (2, 3))
(samples_10,"mu_intercept"], freq = , main = "Intercept")
((x, = (lm_10)"(Intercept)"], = (lm_10)$"(Intercept)", "Std. Error"]), add = , lwd = 2)
(samples_10,"mu_x_fac3t[1]"], freq = , main = "x_fac3t")
((x, = (lm_10)"x_fac3tB"], = (lm_10)$"x_fac3tB", "Std. Error"]), add = , lwd = 2)
(samples_10,"mu_x_fac3t[2]"], freq = , main = "x_fac3t")
((x, = (lm_10)"x_fac3tC"], = (lm_10)$"x_fac3tC", "Std. Error"]), add = , lwd = 2)
(samples_10,"mu_x_cont1"], freq = , main = "x_cont1")
((x, = (lm_10)"x_cont1"], = (lm_10)$"x_cont1", "Std. Error"]), add = , lwd = 2)
(samples_10,"mu_x_cont1__xXx__x_fac3t[1]"], freq = , main = "x_cont1:x_fac3t")
((x, = (lm_10)"x_cont1:x_fac3tB"], = (lm_10)$"x_cont1:x_fac3tB", "Std. Error"]), add = , lwd = 2)
(samples_10,"mu_x_cont1__xXx__x_fac3t[2]"], freq = , main = "x_cont1:x_fac3t")
((x, = (lm_10)"x_cont1:x_fac3tC"], = (lm_10)$"x_cont1:x_fac3tC", "Std. Error"]), add = , lwd = 2)
})
# scaling formula parameters by another parameter works ----
prior_list_1s <- prior_list_all[c("intercept", "x_cont1")]
(prior_list_1s$x_cont1, "multiply_by") <- "sigma"
formula_1s <- (~ x_cont1, parameter = "mu", = df_all,"x_cont1", = ], prior_list = prior_list_1s)
prior_list_1s <- (formula_1s$prior_list, prior_list2)
model_syntax_1s<- (("model{", formula_1s$formula_syntax, model_syntax, "}"), prior_list_1s)
data_1 <- (formula_1$, N = (df_all), y = (df_all$y))
model_1s <- rjags::jags.model( = (model_syntax_1s), inits = (prior_list_1s, chains = 2, seed = 1), = data_1, n.chains = 2, quiet = )
samples_1s <- rjags::coda.samples(model = model_1s, = (prior_list_1s), n.iter = 5000, quiet = , progress.bar = "none")
samples_1s <- (, samples_1s)
expect_equal(formula_1s$formula_syntax, "for(i in 1:N_mu){\n mu[i] = mu_intercept + sigma * mu_x_cont1 * mu_data_x_cont1[i]\n}\n")
lm_1s <- ::(y ~ ((y) * x_cont1), = df_all)
vdiffr::expect_doppelganger("JAGS-formula-lm-1s", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfcol = (1, 3))
(samples_1s,"mu_intercept"], freq = , main = "Intercept")
((x, = (lm_1s)"(Intercept)"], = (lm_1s)$"(Intercept)", "Std. Error"]), add = , lwd = 2)
(samples_1s,"mu_x_cont1"], freq = , main = "I(sd(y) * x_cont1)")
((x, = (lm_1s)"I(sd(y) * x_cont1)"], = (lm_1s)$"I(sd(y) * x_cont1)", "Std. Error"]), add = , lwd = 2)
(samples_1s,"sigma"], freq = , main = "sigma")
(v = (lm_1s), lwd = 3)
})
# input checks work
expect_error((~ x_cont1 , parameter = "mu", = df_all,("x_cont1"), = ], prior_list = prior_list_all[c("x_cont1")]),
"The 'intercept' objects are missing in the 'prior_list' argument.")
expect_error((~ x_cont1 , parameter = "mu", = df_all,("x_cont1"), = ], prior_list = prior_list_all[c("intercept")]),
"The 'x_cont1' objects are missing in the 'prior_list' argument.")
expect_error((~ x_fac2t , parameter = "mu", = df_all,("x_cont1"), = ], prior_list = prior_list_all[c("intercept", "x_fac2t")]),
"The 'x_fac2t' predictor variable is missing in the data set.")
expect_error((~ x_fac2t , parameter = "mu", = (df_all), prior_list = prior_list_all[c("intercept", "x_fac2t")]),
"'data' must be a data.frame")
expect_error((~ x_fac2t , parameter = "mu", = df_all, prior_list = (
"intercept" = ("normal", (0, 1)),
"x_fac2t" = ("normal", (0, 1))
)), "Unsupported prior distribution defined for 'x_fac2t' factor variable")
expect_error((~ x_cont1 , parameter = "mu", = df_all, prior_list = (
"intercept" = ("normal", (0, 1)),
"x_cont1" = ("normal", (0, 1), contrast = "treatment")
)), "Unsupported prior distribution defined for 'x_cont1' continuous variable.")
# linear regression with an independent factor (3 levels) ----
formula_11 <- (~ x_fac3i - 1, parameter = "mu", = df_all,"x_fac3i", = ], prior_list = prior_list_all[c("x_fac3i")])
prior_list_11 <- (formula_11$prior_list, prior_list2)
model_syntax_11 <- (("model{", formula_11$formula_syntax, model_syntax, "}"), prior_list_11)
data_11 <- (formula_11$, N = (df_all), y = (df_all$y))
model_11 <- rjags::jags.model( = (model_syntax_11), inits = (prior_list_11, chains = 2, seed = 1), = data_11, n.chains = 2, quiet = )
samples_11 <- rjags::coda.samples(model = model_11, = (prior_list_11), n.iter = 5000, quiet = , progress.bar = "none")
samples_11 <- (, samples_11)
lm_11 <- ::(y ~ x_fac3i - 1, = df_all)
vdiffr::expect_doppelganger("JAGS-formula-lm-11", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfcol = (1, 3))
(samples_11,"mu_x_fac3i[1]"], freq = , main = "x_fac3i[1]")
((x, = (lm_11)"x_fac3iA"], = (lm_11)$"x_fac3iA", "Std. Error"]), add = , lwd = 2)
(samples_11,"mu_x_fac3i[2]"], freq = , main = "x_fac3i[2]")
((x, = (lm_11)"x_fac3iB"], = (lm_11)$"x_fac3iB", "Std. Error"]), add = , lwd = 2)
(samples_11,"mu_x_fac3i[3]"], freq = , main = "x_fac3i[3]")
((x, = (lm_11)"x_fac3iC"], = (lm_11)$"x_fac3iC", "Std. Error"]), add = , lwd = 2)
})
# linear regression with a meandif factor (3 levels) ----
formula_12 <- (~ x_fac3md, parameter = "mu", = df_all,"x_fac3md", = ], prior_list = prior_list_all[c("intercept", "x_fac3md")])
prior_list_12 <- (formula_12$prior_list, prior_list2)
model_syntax_12 <- (("model{", formula_12$formula_syntax, model_syntax, "}"), prior_list_12)
data_12 <- (formula_12$, N = (df_all), y = (df_all$y))
model_12 <- rjags::jags.model( = (model_syntax_12), inits = (prior_list_12, chains = 2, seed = 1), = data_12, n.chains = 2, quiet = )
samples_12 <- rjags::coda.samples(model = model_12, = (prior_list_12), n.iter = 5000, quiet = , progress.bar = "none")
samples_12 <- (, samples_12)
df_12 <- df_all
(df_12$x_fac3md) <- contr.meandif((df_12$x_fac3o))
lm_12 <- ::(y ~ x_fac3md, = df_12)
vdiffr::expect_doppelganger("JAGS-formula-lm-12", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfcol = (1, 3))
(samples_12,"mu_intercept"], freq = , main = "Intercept")
((x, = (lm_12)"(Intercept)"], = (lm_12)$"(Intercept)", "Std. Error"]), add = , lwd = 2)
(samples_12,"mu_x_fac3md[1]"], freq = , main = "x_fac3md")
((x, = (lm_12)"x_fac3md1"], = (lm_12)$"x_fac3md1", "Std. Error"]), add = , lwd = 2)
(samples_12,"mu_x_fac3md[2]"], freq = , main = "x_fac3md")
((x, = (lm_12)"x_fac3md2"], = (lm_12)$"x_fac3md2", "Std. Error"]), add = , lwd = 2)
})
# linear regression with a spike independent factor (3 levels) ----
prior_list_13 <- ("x_fac3i" = ("spike", contrast = "independent", (1.5)))
formula_13 <- (~ x_fac3i - 1, parameter = "mu", = df_all,"x_fac3i", = ], prior_list = prior_list_13)
prior_list_13 <- (formula_13$prior_list, prior_list2)
model_syntax_13 <- (("model{", formula_13$formula_syntax, model_syntax, "}"), prior_list_13)
data_13 <- (formula_13$, N = (df_all), y = (df_all$y))
model_13 <- rjags::jags.model( = (model_syntax_13), inits = (prior_list_13, chains = 2, seed = 1), = data_13, n.chains = 2, quiet = )
samples_13 <- rjags::coda.samples(model = model_13, = (prior_list_13), n.iter = 5000, quiet = , progress.bar = "none")
samples_13 <- (, samples_13)
expect_equal((3), (1:3))
vdiffr::expect_doppelganger("JAGS-formula-lm-13", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfcol = (1, 3))
(samples_13,"mu_x_fac3i[1]"], freq = , main = "x_fac3i[1]")
(samples_13,"mu_x_fac3i[2]"], freq = , main = "x_fac3i[2]")
(samples_13,"mu_x_fac3i[3]"], freq = , main = "x_fac3i[3]")
})
# linear regression with a meandif spike factor (3 levels) ----
prior_list_14 <- ("intercept" = prior_list_all$intercept, "x_fac3md" = ("spike", contrast = "meandif", (0)))
formula_14 <- (~ x_fac3md, parameter = "mu", = df_all,"x_fac3md", = ], prior_list = prior_list_14)
prior_list_14 <- (formula_14$prior_list, prior_list2)
model_syntax_14 <- (("model{", formula_14$formula_syntax, model_syntax, "}"), prior_list_14)
data_14 <- (formula_14$, N = (df_all), y = (df_all$y))
model_14 <- rjags::jags.model( = (model_syntax_14), inits = (prior_list_14, chains = 2, seed = 1), = data_14, n.chains = 2, quiet = )
samples_14 <- rjags::coda.samples(model = model_14, = (prior_list_14), n.iter = 5000, quiet = , progress.bar = "none")
samples_14 <- (, samples_14)
df_14 <- df_all
(df_14$x_fac3md) <- contr.meandif((df_14$x_fac3o))
lm_14 <- ::(y ~ 1, = df_14)
vdiffr::expect_doppelganger("JAGS-formula-lm-14", (){
oldpar <- ::(no.readonly = )
(::(mfcol = oldpar"mfcol"]]))
(mfcol = (1, 3))
(samples_14,"mu_intercept"], freq = , main = "Intercept")
((x, = (lm_14)"(Intercept)"], = (lm_14)$"(Intercept)", "Std. Error"]), add = , lwd = 2)
(samples_14,"mu_x_fac3md[1]"], freq = , main = "x_fac3md")
(samples_14,"mu_x_fac3md[2]"], freq = , main = "x_fac3md")
})
})
test_that("JAGS formula tools work", {
# additional tools work
expect_equal(
(("mu_x_cont", "mu_x_fac3t", "mu_x_fac3t__xXx__x_cont")),
("mu_x_cont", "mu_x_fac3t", "mu_x_fac3t:x_cont")
)
expect_equal(
(("mu_x_cont", "mu_x_fac3t", "mu_x_fac3t__xXx__x_cont"), formula_parameters = "mu"),
("(mu) x_cont", "(mu) x_fac3t", "(mu) x_fac3t:x_cont")
)
expect_equal(
(("mu_x_cont", "mu_x_fac3t", "mu_x_fac3t__xXx__x_cont"), formula_parameters = "mu", formula_prefix = ),
("x_cont", "x_fac3t", "x_fac3t:x_cont")
)
expect_equal(
(("x_cont", "x_fac3t", "x_fac3t:x_cont")),
("x_cont", "x_fac3t", "x_fac3t__xXx__x_cont")
)
expect_equal(
(("x_cont", "x_fac3t", "x_fac3t:x_cont"), formula_parameter = "mu"),
("mu_x_cont", "mu_x_fac3t", "mu_x_fac3t__xXx__x_cont")
)
})
test_that("JAGS evaluate formula works", {
# check the posterior distributions with weak priors against a maximum likelihood estimates with ML
skip_on_os(("mac", "linux", "solaris")) # multivariate sampling does not exactly match across OSes
skip_on_cran()
# complex formula including scaling
(1)
df_all <- (
x_cont1 = (60),
x_cont2 = (60),
x_fac2t = ((("A", "B"), 30), = ("A", "B")),
x_fac3o = ((("A", "B", "C"), 20), = ("A", "B", "C"))
)
df_all$y <- (60, 0.1, 0.5) + 0.30 * df_all$x_cont1 - 0.15 * df_all$x_cont1 * df_all$x_cont2 +
(df_all$x_fac3o == "A", 0.2, (df_all$x_fac3o == "B", -0.2, 0))
prior_list_all <- (
"intercept" = ("normal", (0, 5)),
"x_cont1" = ("normal", (0, 1)),
"x_cont2" = ("normal", (0, 1)),
"x_fac2t" = ("normal", contrast = "treatment", (0, 1)),
"x_fac3o" = ("mnormal", contrast = "orthonormal", (0, 1)),
"x_cont1:x_fac3o" = ("mnormal", contrast = "orthonormal", (0, 2))
)
prior_list2 <- (
"sigma" = ("cauchy", (0, 1), (0, 1))
)
model_syntax <- (
"for(i in 1:N){\n",
" y[i] ~ dnorm(mu[i], 1/pow(sigma, 2))\n",
"}\n"
)
<- (~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, parameter = "mu", = df_all, prior_list = prior_list_all)
prior_list <- ($prior_list, prior_list2)
model_syntax <- (("model{", $formula_syntax, model_syntax, "}"), prior_list)
<- ($, N = (df_all), y = (df_all$y))
model <- rjags::jags.model( = (model_syntax), inits = (prior_list, chains = 1, seed = 1), = , n.chains = 1, quiet = )
samples <- rjags::coda.samples(model = model, = (prior_list), n.iter = 10, quiet = , progress.bar = "none")
new_data <- (
x_cont1 = (0, 0, 1, 1),
x_cont2 = (0, 1, 0, 1),
x_fac2t = (("A", "B", "A", "B"), = ("A", "B")),
x_fac3o = (("A", "B", "C", "A"), = ("A", "B", "C"))
)
# test the results against the lm function (by passing the ML estimates)
(df_all$x_fac3o) <- ((df_all$x_fac3o))
fit_lm <- ::(y~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, = df_all)
samples_new <- ((fit_lm), = (fit_lm))[c("(Intercept)","x_cont1","x_cont1:x_fac3o1","x_cont1:x_fac3o2","x_cont2","x_fac2tB","x_fac3o1","x_fac3o2","sigma")]
samples_new <- (samples_new, = 1)
(samples_new) <- (samples[1]])
samples_new <- coda::as.mcmc.list(coda::as.mcmc(samples_new))
expect_equal((fit_lm, newdata = new_data), (samples_new, ~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, "mu", new_data, prior_list),1])
# for a posterior samples matrix
samples_new <- ((fit_lm), = (fit_lm))[c("(Intercept)","x_cont1","x_cont1:x_fac3o1","x_cont1:x_fac3o2","x_cont2","x_fac2tB","x_fac3o1","x_fac3o2","sigma")]
samples_new <- (samples_new, = 5, = (samples_new), byrow = )
(samples_new) <- (samples[1]])
samples_new <- coda::as.mcmc.list(coda::as.mcmc(samples_new))
expect_equal(((fit_lm, newdata = new_data), = 4, = 5), ((samples_new, ~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, "mu", new_data, prior_list)))
# check filling in missing or miss ordered factor levels
samples_new <- ((fit_lm), = (fit_lm))[c("(Intercept)","x_cont1","x_cont1:x_fac3o1","x_cont1:x_fac3o2","x_cont2","x_fac2tB","x_fac3o1","x_fac3o2","sigma")]
samples_new <- (samples_new, = 1)
(samples_new) <- (samples[1]])
samples_new <- coda::as.mcmc.list(coda::as.mcmc(samples_new))
new_data2 <- new_data
new_data2$x_fac2t <- ((new_data2$x_fac2t), = ("B", "A"))
expect_equal((fit_lm, newdata = new_data), (samples_new, ~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, "mu", new_data2, prior_list),1])
new_data3 <- new_data
new_data3$x_fac3o <- (("A", "B", "A", "B"), = ("B", "A"))
expect_equal((fit_lm, newdata = new_data3), (samples_new, ~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, "mu", new_data3, prior_list),1])
new_data4 <- new_data
new_data4$x_fac3o <- ("A", "B", "A", "B")
expect_equal((fit_lm, newdata = new_data3), (samples_new, ~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, "mu", new_data4, prior_list),1])
# check scaling works (by multiplying be zero)
prior_list2 <- prior_list
(prior_list2$mu_x_cont2, "multiply_by") <- 0
(prior_list2$mu_x_fac2t, "multiply_by") <- 0
samples_new2 <- ((fit_lm), = (fit_lm))[c("(Intercept)","x_cont1","x_cont1:x_fac3o1","x_cont1:x_fac3o2","x_cont2","x_fac2tB","x_fac3o1","x_fac3o2","sigma")]
samples_new2 <- (samples_new2, = 1)
(samples_new2) <- (samples[1]])
samples_new2,"mu_x_cont2"] <- 0
samples_new2,"mu_x_fac2t"] <- 0
samples_new2 <- coda::as.mcmc.list(coda::as.mcmc(samples_new2))
expect_equal((samples_new, ~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, "mu", new_data, prior_list2),1],
(samples_new2, ~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, "mu", new_data, prior_list),1])
# check scaling by another parameter works
prior_list2 <- prior_list
(prior_list2$mu_x_cont2, "multiply_by") <- "sigma"
expect_equal((((samples, ~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, "mu", new_data, prior_list2),1])),
(0.4436353, -0.0658681, 0.1870391, 0.8548012), tolerance = 1e-5)
### test input tests
expect_error((samples_new, ~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, "mu", new_data,1:3], prior_list),
"The 'x_fac3o' predictor variable is missing in the data.")
expect_error((samples_new, ~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, "mu", new_data, prior_list-1]),
"The prior distribution for the 'x_fac2t' term is missing in the prior_list.")
bad_data <- new_data
bad_data$x_fac2t <- (("C", "B", "C", "B"), = ("B", "C"))
expect_error((samples_new, ~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, "mu", bad_data, prior_list),
"Levels specified in the 'x_fac2t' factor variable do not match the levels used for model specification.")
bad_data2 <- new_data
bad_data2$x_fac2t <- ("C", "B", "C", "B")
expect_error((samples_new, ~ x_fac2t + x_cont2 + x_cont1 * x_fac3o, "mu", bad_data2, prior_list),
"Levels specified in the 'x_fac2t' factor variable do not match the levels used for model specification.")
# evaluate formula with spike prior distributions ----
(1)
df_all <- (
x_fac2i = ((("A", "B"), 30), = ("A", "B")),
x_fac3o = ((("A", "B", "C"), 60, = ), = ("A", "B", "C")),
x_fac3t = ((("A", "B", "C"), 60, = ), = ("A", "B", "C")),
x_fac3md = ((("A", "B", "C"), 60, = ), = ("A", "B", "C"))
)
df_all$y <- (60, 0.1, 0.5)
prior_list_all <- (
"intercept" = ("normal", (0, 5)),
"x_fac2i" = ("spike", contrast = "independent", (1)),
"x_fac3o" = ("spike", contrast = "orthonormal", (0)),
"x_fac3t" = ("spike", contrast = "treatment", (2)),
"x_fac3md" = ("spike", contrast = "meandif", (0))
)
prior_list2 <- (
"sigma" = ("cauchy", (0, 1), (0, 1))
)
model_syntax <- (
"model{",
"for(i in 1:N){\n",
" y[i] ~ dnorm(mu[i], 1/pow(sigma, 2))\n",
"}\n",
"}"
)
fit1 <- (
model_syntax = model_syntax,
formula_list = (mu = ~ x_fac2i + x_fac3o + x_fac3t + x_fac3md),
= (y = df_all$y, N = (df_all)),
prior_list = prior_list2,
formula_data_list = (mu = df_all),
formula_prior_list = (mu = prior_list_all))
new_data <- (
x_fac2i = (("A", "B", "A"), = ("A", "B")),
x_fac3o = (("A", "A", "B"), = ("A", "B", "C")),
x_fac3t = (("A", "B", "C"), = ("A", "B", "C")),
x_fac3md = (("B", "B", "C"), = ("A", "B", "C"))
)
new_samples <- (fit1, ~ x_fac2i + x_fac3o + x_fac3t + x_fac3md, "mu", new_data, (fit1, "prior_list"))
new_samples <- (new_samples, 1, )
intercept_estimate <- (fit1)"(mu) intercept", "Mean"]
expect_equivalent(intercept_estimate + 1, new_samples[1])
expect_equivalent(intercept_estimate + 1 + 2, new_samples[2])
expect_equivalent(intercept_estimate + 1 + 2, new_samples[3])
})
