unm_glm {unmconf} | R Documentation |
Fitting Multilevel Bayesian Regression Model with Unmeasured Confounders
Description
unm_glm()
fits a multilevel Bayesian regression model that accounts for
unmeasured confounders. Users can input model information into unm_glm()
in
a similar manner as they would for the standard stats::glm()
function,
providing arguments like formula
, family
, and data
. Results are stored
as MCMC iterations.
Usage
unm_glm(
form1,
form2 = NULL,
form3 = NULL,
family1 = binomial(),
family2 = NULL,
family3 = NULL,
data,
n.iter = 2000,
n.adapt = 1000,
thin = 1,
n.chains = 4,
filename = tempfile(fileext = ".jags"),
quiet = getOption("unm_quiet"),
progress.bar = getOption("unm_progress.bar"),
code_only = FALSE,
priors,
response_nuisance_priors,
response_params_to_track,
confounder1_nuisance_priors,
confounder1_params_to_track,
confounder2_nuisance_priors,
confounder2_params_to_track,
...
)
jags_code(mod)
## S3 method for class 'unm_int'
print(x, digits = 3, ..., print_call = getOption("unm_print_call"))
## S3 method for class 'unm_int'
coef(object, ...)
Arguments
form1 |
The formula specification for the response model (Level I) |
form2 |
The formula specification for the first unmeasured confounder model (Level II) |
form3 |
The formula specification for the second unmeasured confounder model (Level III) |
family1 , family2 , family3 |
The family object, communicating the types of
models to be used for response ( |
data |
The dataset containing all variables (this function currently only supports a single dataset containing internally validated data) |
n.iter |
|
n.adapt |
|
thin |
|
n.chains |
|
filename |
File name where to store jags code |
quiet |
The |
progress.bar |
The |
code_only |
Should only the code be created? |
priors |
Custom priors to use on regression coefficients, see examples. |
response_nuisance_priors , confounder1_nuisance_priors , confounder2_nuisance_priors |
JAGS code for the nuisance priors on parameters in a JAGS model (see examples) |
response_params_to_track , confounder1_params_to_track , confounder2_params_to_track |
Additional parameters to track when nuisance parameter priors are used (see examples) |
... |
Additional arguments to pass into |
mod |
The output of |
x |
Object to be printed |
digits |
Number of digits to round to; defaults to 3 |
print_call |
Should the call be printed? Defaults to |
object |
Model object for which the coefficients are desired |
Value
(Invisibly) The output of rjags::coda.samples()
, an object of class
mcmc.list
, along with attributes code
containing the jags code used and
file
containing the filename of the jags code.
See Also
Examples
# ~~ One Unmeasured Confounder Examples (II-Level Model) ~~
# normal response, normal confounder model with internally validated data
(df <- runm(20, response = "norm"))
(unm_mod <- unm_glm(
form1 = y ~ x + z1 + z2 + z3 + u1, family1 = gaussian(),
form2 = u1 ~ x + z1 + z2 + z3, family2 = gaussian(),
data = df
))
## Not run: # reduce cran check time
(unm_mod <- unm_glm(
y ~ ., family1 = gaussian(),
u1 ~ . - y, family2 = gaussian(),
data = df
))
glm(y ~ x + z1 + z2 + z3, data = df)
coef(unm_mod)
jags_code(unm_mod)
unm_glm(
y ~ .,
u1 ~ . - y,
family1 = gaussian(),
family2 = gaussian(),
data = df, code_only = TRUE
)
# a normal-normal model - external validation
(df <- runm(c(10, 10), type = "ext", response = "norm"))
unm_glm(
y ~ x + z1 + z2 + z3 + u1, family1 = gaussian(),
u1 ~ x + z1 + z2 + z3, family2 = gaussian(),
data = df
)
# setting custom priors
unm_glm(
y ~ ., family1 = gaussian(),
u1 ~ . - y, family2 = gaussian(),
data = df,
code_only = TRUE
)
unm_glm(
y ~ ., family1 = gaussian(),
u1 ~ . - y, family2 = gaussian(),
data = df,
code_only = FALSE,
priors = c("lambda[u1]" = "dnorm(1, 10)"),
response_nuisance_priors = "tau_{y} <- sigma_{y}^-2; sigma_{y} ~ dunif(0, 100)",
response_params_to_track = "sigma_{y}",
confounder1_nuisance_priors = "tau_{u1} <- sigma_{u1}^-2; sigma_{u1} ~ dunif(0, 100)",
confounder1_params_to_track = "sigma_{u1}"
)
# turn progress tracking on
options("unm_progress.bar" = "text")
# more complex functional forms _for non-confounder predictors only_
# zero-intercept model
unm_glm(
y ~ . - 1,
u1 ~ . - y,
family1 = gaussian(),
family2 = gaussian(),
data = df
)
glm(y ~ . - 1, data = df)
# polynomial model
unm_glm(
y ~ x + poly(z1, 2) + u1,
u1 ~ x + z1,
family1 = gaussian(),
family2 = gaussian(),
data = df
)
glm(y ~ x + poly(z1, 2), data = df)
# interaction model
unm_glm(
y ~ x*z1 + u1, family1 = gaussian(),
u1 ~ x*z1, family2 = gaussian(),
data = df
)
glm(y ~ x*z1, data = df)
# a binomial-binomial model
(df <- runm(
50,
missing_prop = .75,
response = "bin",
unmeasured_fam_list = list("bin"),
unmeasured_param_list = list(.5)
))
(unm_mod <- unm_glm(
y ~ .,
u1 ~ . - y,
family1 = binomial(),
family2 = binomial(),
data = df
))
glm(y ~ . - u1, family = binomial(), data = df)
# a poisson-normal model
(df <- runm(
25,
response = "pois",
response_model_coefs = c("int" = -1, "z" = .5, "u1" = .5, "x" = .5),
treatment_model_coefs = c("int" = -1, "z" = .5, "u1" = .5),
covariate_fam_list = list("norm"),
covariate_param_list = list(c(mean = 0, sd = 1)),
unmeasured_fam_list = list("norm"),
unmeasured_param_list = list(c(0, 1))
))
(unm_mod <- unm_glm(
y ~ x + z + u1 + offset(log(t)),
u1 ~ x + z,
family1 = poisson(),
family2 = gaussian(),
data = df
))
glm(y ~ x + z + offset(log(t)), family = poisson(), data = df)
# a poisson-binomial model
(df <- runm(
25,
response = "pois",
response_model_coefs = c("int" = -1, "z" = .5, "u1" = .5, "x" = .5),
treatment_model_coefs = c("int" = -1, "z" = .5, "u1" = .5),
covariate_fam_list = list("norm"),
covariate_param_list = list(c(mean = 0, sd = 1)),
unmeasured_fam_list = list("bin"),
unmeasured_param_list = list(.5)
))
(unm_mod <- unm_glm(
y ~ x + z + u1 + offset(log(t)), family1 = poisson(),
u1 ~ x + z, family2 = binomial(),
data = df
))
glm(y ~ x + z + offset(log(t)), family = poisson(), data = df)
# a gamma-normal model
(df <- runm(
25,
response = "gam",
response_model_coefs = c("int" = -1, "z" = .5, "u1" = .5, "x" = .5),
treatment_model_coefs = c("int" = -1, "z" = .5, "u1" = .5),
covariate_fam_list = list("norm"),
covariate_param_list = list(c(mean = 0, sd = 1)),
unmeasured_fam_list = list("norm"),
unmeasured_param_list = list(c(0, 1))
))
(unm_mod <- unm_glm(
y ~ x + z + u1, family1 = Gamma(),
u1 ~ x + z, family2 = gaussian(),
data = df
))
glm(y ~ x + z, family = Gamma(link = "log"), data = df)
# a gamma-binomial model
(df <- runm(
25,
response = "gam",
response_model_coefs = c("int" = -1, "z" = .5, "u1" = .5, "x" = .5),
treatment_model_coefs = c("int" = -1, "z" = .5, "u1" = .5),
covariate_fam_list = list("norm"),
covariate_param_list = list(c(mean = 0, sd = 1)),
unmeasured_fam_list = list("bin"),
unmeasured_param_list = list(.5)
))
(unm_mod <- unm_glm(
y ~ x + z + u1, family1 = Gamma(),
u1 ~ x + z, family2 = binomial(),
data = df
))
print(df, n = 25)
glm(y ~ x + z, family = Gamma(link = "log"), data = df)
# the output of unm_glm() is classed jags output
(df <- runm(20, response = "norm"))
(unm_mod <- unm_glm(
y ~ .,
u1 ~ . - y,
family1 = gaussian(),
family2 = gaussian(),
data = df)
)
class(unm_mod)
jags_code(unm_mod)
unm_glm(y ~ ., u1 ~ . - y, data = df, code_only = TRUE)
# visualizing output
library("ggplot2")
library("bayesplot"); bayesplot_theme_set(ggplot2::theme_minimal())
mcmc_hist(unm_mod, facet_args = list(labeller = label_parsed))
mcmc_hist(unm_mod)
mcmc_trace(unm_mod, facet_args = list(labeller = label_parsed))
# more extensive visualization with the tidyverse
mcmc_intervals(unm_mod, prob = .90) +
geom_point(
aes(value, name), data = tibble::enframe(attr(df, "params")),
color = "red", fill = "pink", size = 4, shape = 21
)
library("dplyr")
library("tidyr")
unm_mod %>%
as.matrix() %>%
as_tibble() %>%
pivot_longer(everything(), names_to = "var", values_to = "val") %>%
ggplot(aes("0", val)) +
geom_jitter() +
geom_point(
aes("0", value), data = tibble::enframe(attr(df, "params"), name = "var"),
color = "red", fill = "pink", size = 4, shape = 21
) +
coord_flip() +
facet_grid(var ~ ., scales = "free_y", labeller = label_parsed) +
theme_bw() +
theme(
axis.title = element_blank(),
axis.text.y = element_blank(), axis.ticks.y = element_blank(),
strip.text.y = element_text(angle = 0)
)
# getting draws out
(samps <- posterior::as_draws_df(unm_mod))
samps$`.chain`
samps$`.iteration`
samps$`.draw`
# implementation is variable-name independent
(df <- runm(100, response = "norm"))
df$ht <- df$y
df$age <- df$u1
df$biom <- df$x
(unm_mod <- unm_glm(
ht ~ x + biom + age,
age ~ x + biom,
data = df,
family1 = gaussian(),
family2 = gaussian()
))
jags_code(unm_mod)
# ~~ Two Unmeasured Confounders Examples (III-Level Model) ~~
# a normal-normal-normal model - internal validation
(df <- runm(
50,
missing_prop = .75,
response = "norm",
response_model_coefs = c("int" = -1, "z" = .5,
"u1" = .5, "u2" = .5, "x" = .5),
treatment_model_coefs = c("int" = -1, "z" = .5,
"u1" = .5, "u2" = .5),
covariate_fam_list = list("norm"),
covariate_param_list = list(c(mean = 0, sd = 1)),
unmeasured_fam_list = list("norm", "norm"),
unmeasured_param_list = list(c(0, 1), c(0, 1))
))
(unm_mod <- unm_glm(
y ~ x + z + u1 + u2,
u1 ~ x + z + u2,
u2 ~ x + z,
family1 = gaussian(),
family2 = gaussian(),
family3 = gaussian(),
data = df
))
glm(y ~ x + z, data = df)
coef(unm_mod)
unm_glm(
y ~ x + z + u1 + u2, family1 = gaussian(),
u1 ~ x + z + u2, family2 = gaussian(),
u2 ~ x + z, family3 = gaussian(),
data = df,
code_only = TRUE
)
# a normal-normal-normal model - external validation
(df <- runm(
c(20, 20),
type = "ext",
response = "norm",
response_model_coefs = c("int" = -1, "z" = .5,
"u1" = .5, "u2" = .5, "x" = .5),
treatment_model_coefs = c("int" = -1, "z" = .5,
"u1" = .5, "u2" = .5),
covariate_fam_list = list("norm"),
covariate_param_list = list(c(mean = 0, sd = 1)),
unmeasured_fam_list = list("norm", "norm"),
unmeasured_param_list = list(c(0, 1), c(0, 1))
))
(unm_mod <- unm_glm(
y ~ x + z + u1 + u2, family1 = gaussian(),
u1 ~ x + z + u2, family2 = gaussian(),
u2 ~ x + z, family3 = gaussian(),
data = df
))
# a binomial-binomial-binomial model - internal validation
(df <- runm(
25,
response = "bin",
response_model_coefs = c("int" = -1, "z" = .5,
"u1" = .5, "u2" = .5, "x" = .5),
treatment_model_coefs = c("int" = -1, "z" = .5,
"u1" = .5, "u2" = .5),
covariate_fam_list = list("norm"),
covariate_param_list = list(c(mean = 0, sd = 1)),
unmeasured_fam_list = list("bin", "bin"),
unmeasured_param_list = list(.5, .75)
))
unm_glm(y ~ x + z + u1 + u2, family1 = binomial(),
u1 ~ x + z + u2, family2 = binomial(),
u2 ~ x + z, family3 = binomial(),
data = df,
code_only = TRUE
)
## End(Not run)