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 (form1) and unmeasured confounder (⁠form2, form3⁠) models. See stats::family() for details

data

The dataset containing all variables (this function currently only supports a single dataset containing internally validated data)

n.iter

n.iter argument of rjags::coda.samples()

n.adapt

n.adapt argument of rjags::jags.model()

thin

thin argument of rjags::coda.samples()

n.chains

n.chains argument of rjags::jags.model()

filename

File name where to store jags code

quiet

The quiet parameter of rjags::jags.model(). Defaults to TRUE, but you can change it on a per-session basis with options(unm_quiet = FALSE).

progress.bar

The progress.bar parameter of rjags::update.jags(). Defaults to "none", but you can change it on a per-session basis with options(unm_progress.bar = "text").

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 rjags::jags.model(), such as inits

mod

The output of unm_glm()

x

Object to be printed

digits

Number of digits to round to; defaults to 3

print_call

Should the call be printed? Defaults to TRUE, but can be turned off with options("unm_print_call" = FALSE)

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

runm(), rjags::dic.samples()

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)


[Package unmconf version 1.0.0 Index]