glm_nab: Fit GLM with the 'naive adaptive bayes' prior

Description Usage Arguments Value Examples

View source: R/glm_nab.R

Description

Program for fitting a GLM equipped with the 'naive adaptive bayes' prior evaluated in the manuscript.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
glm_nab(
  stan_fit = stanmodels$NAB_Stable,
  y,
  x_standardized,
  alpha_prior_mean,
  alpha_prior_cov,
  phi_mean,
  phi_sd,
  beta_orig_scale,
  beta_aug_scale,
  beta_aug_scale_tilde,
  local_dof = 1,
  global_dof = 1,
  slab_precision = (1/15)^2,
  only_prior = F,
  mc_warmup = 50,
  mc_iter_after_warmup = 50,
  mc_chains = 1,
  mc_thin = 1,
  mc_stepsize = 0.1,
  mc_adapt_delta = 0.9,
  mc_max_treedepth = 15,
  ntries = 1,
  return_as_stanfit = FALSE,
  eigendecomp_hist_var = NULL,
  scale_to_variance225 = NULL
)

Arguments

stan_fit

an R object of class stanfit, which allows the function to run without recompiling the stan code.

y

(vector) outcomes corresponding to the type of glm desired. This should match whatever datatype is expected by the stan program.

x_standardized

(matrix) matrix of numeric values with number of rows equal to the length of y and number of columns equal to p+q. It is assumed without verification that each column is standardized to whatever scale the prior expects - in Boonstra and Barbaro, all predictors are marginally generated to have mean zero and unit variance, so no standardization is conducted. In practice, all data should be standardized to have a common scale before model fitting. If regression coefficients on the natural scale are desired, they can be easily obtained through unstandardizing.

alpha_prior_mean

(vector) p-length vector giving the mean of alpha from the historical analysis, corresponds to m_alpha in Boonstra and Barbaro

alpha_prior_cov

(matrix) pxp positive definite matrix giving the variance of alpha from the historical analysis, corresponds to S_alpha in Boonstra and Barbaro

phi_mean

(real) mean of phi corresponding to a truncated normal distribution. Since the support of the distribution is truncated to [0,1], it would make sense, but is not required, that 'phi_mean' itself also be in [0,1]

phi_sd

(pos. real) sd of phi corresponding to a truncated normal distribution.

beta_aug_scale

(pos. real) constants indicating the prior scale of the horseshoe. Both values correspond to 'c' in the notation of Boonstra and Barbaro, because that paper never considers beta_orig_scale!=beta_aug_scale

beta_aug_scale_tilde

(pos. real) constant indicating the prior scale of the horseshoe for the augmented covariates when phi = 1, i.e. when the historical analysis is fully used. This corresponds to tilde_c in Boonstra and Barbaro

local_dof

(pos. integer) numbers indicating the degrees of freedom for lambda_j and tau, respectively. Boonstra, et al. never considered local_dof != 1 or global_dof != 1.

global_dof

(pos. integer) numbers indicating the degrees of freedom for lambda_j and tau, respectively. Boonstra, et al. never considered local_dof != 1 or global_dof != 1.

slab_precision

(pos. real) the slab-part of the regularized horseshoe, this is equivalent to (1/d)^2 in the notation of Boonstra and Barbaro

only_prior

(logical) should all data be ignored, sampling only from the prior?

mc_warmup

number of MCMC warm-up iterations

mc_iter_after_warmup

number of MCMC iterations after warm-up

mc_chains

number of MCMC chains

mc_thin

every nth draw to keep

mc_stepsize

positive stepsize

mc_adapt_delta

between 0 and 1

mc_max_treedepth

max tree depth

ntries

(pos. integer) the stan function will run up to this many times, stopping either when the number of divergent transitions* is zero or when ntries has been reached. The reported fit will be that with the fewest number of divergent iterations.

return_as_stanfit

(logical) should the function return the stanfit object asis or should a summary of stanfit be returned as a regular list

eigendecomp_hist_var:

R object of class 'eigen' containing a pxp matrix of eigenvectors in each row (equivalent to v_0 in Boonstra and Barbaro) and a p-length vector of eigenvalues. This is by default equal to eigen(alpha_prior_cov)

scale_to_variance225:

a vector assumed to be such that, when multiplied by the diagonal elements of alpha_prior_cov, the result is a vector of elements each equal to 225. This is explicitly calculated if it is not provided

Value

list object containing the draws and other information.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
data(current)

alpha_prior_cov = matrix(data = c(0.02936, -0.02078, 0.00216, -0.00637,
                                  -0.02078, 0.03192, -0.01029, 0.00500,
                                  0.00216, -0.01029, 0.01991, -0.00428,
                                  -0.00637, 0.00500, -0.00428, 0.01650),
                         byrow = FALSE, nrow = 4);
scale_to_variance225 = diag(alpha_prior_cov) / 225;
eigendecomp_hist_var = eigen(alpha_prior_cov);

foo = glm_nab(y = current$y_curr,
              x_standardized = current[,2:11],
              alpha_prior_mean = c(1.462, -1.660, 0.769, -0.756),
              alpha_prior_cov = alpha_prior_cov,
              phi_mean = 1,
              phi_sd = 0.25,
              beta_orig_scale = 0.0223,
              beta_aug_scale = 0.0223,
              beta_aug_scale_tilde = 0.05,
              local_dof = 1,
              global_dof = 1,
              slab_precision = 0.00444,
              only_prior = 0,
              mc_warmup = 1000,
              mc_iter_after_warmup = 1000,
              mc_chains = 2,
              mc_thin = 1,
              mc_stepsize = 0.1,
              mc_adapt_delta = 0.999,
              mc_max_treedepth = 15,
              ntries = 2,
              eigendecomp_hist_var = eigendecomp_hist_var,
              scale_to_variance225 = scale_to_variance225);

umich-biostatistics/AdaptiveBayesianUpdates documentation built on July 29, 2021, 3:06 a.m.