stan_glm: Bayesian generalized linear models via Stan

Description

Generalized linear modeling with optional prior distributions for the coefficients, intercept, and nuisance parameter.

Usage

stan_glm(formula, family = gaussian(), data, weights, subset, na.action = NULL, offset = NULL, model = TRUE, x = FALSE, y = TRUE, contrasts = NULL, ..., prior = normal(), prior_intercept = normal(), prior_ops = prior_options(), prior_PD = FALSE, algorithm = c("sampling", "optimizing", "meanfield", "fullrank"), adapt_delta = NULL, QR = FALSE, sparse = FALSE)
stan_glm.nb(formula, data, weights, subset, na.action = NULL, offset = NULL, model = TRUE, x = FALSE, y = TRUE, contrasts = NULL, link = "log", ..., prior = normal(), prior_intercept = normal(), prior_ops = prior_options(), prior_PD = FALSE, algorithm = c("sampling", "optimizing", "meanfield", "fullrank"), adapt_delta = NULL, QR = FALSE)
stan_glm.fit(x, y, weights = rep(1, NROW(x)), offset = rep(0, NROW(x)), family = gaussian(), ..., prior = normal(), prior_intercept = normal(), prior_ops = prior_options(), group = list(), prior_PD = FALSE, algorithm = c("sampling", "optimizing", "meanfield", "fullrank"), adapt_delta = NULL, QR = FALSE, sparse = FALSE)

Arguments

formula, data, subset

Same as glm.

family

Same as glm, except negative binomial GLMs are also possible using the neg_binomial_2 family object.

na.action, contrasts

Same as glm, but rarely specified.

model, offset, weights

Same as glm.

x, y

In stan_glm, stan_glm.nb, logical scalars indicating whether to return the design matrix and response vector. In stan_glm.fit, a design matrix and response vector.

...

Further arguments passed to the function in the rstan package (sampling, vb, or optimizing), corresponding to the estimation method named by algorithm. For example, if algorithm is "sampling" it is possibly to specify iter, chains, cores, refresh, etc.

prior

The prior distribution for the regression coefficients. prior can be a call to normal, student_t, cauchy, hs or hs_plus. See priors for details. To omit a prior ---i.e., to use a flat (improper) uniform prior--- prior can be set to NULL, although this is rarely a good idea. (Note: unless QR=TRUE, if the scaled argument to prior_options is left at its default and recommended value of TRUE, then the scale(s) of prior may be modified internally based on the scales of the predictors, as in the arm package. See priors for details on the rescaling and prior_summary for a summary of the priors used for a particular model.)

prior_intercept

The prior distribution for the intercept. prior_intercept can be a call to normal, student_t or cauchy. See priors for details. To to omit a prior ---i.e., to use a flat (improper) uniform prior--- set prior_intercept to NULL. (Note: if a dense representation of the design matrix is utilized ---i.e., if the sparse argument is left at its default value of FALSE--- then the prior distribution for the intercept is set so it applies to the value when all predictors are centered.)

prior_ops

Additional options related to prior distributions. Can be NULL to omit a prior on the dispersion and see prior_options otherwise.

prior_PD

A logical scalar (defaulting to FALSE) indicating whether to draw from the prior predictive distribution instead of conditioning on the outcome.

algorithm

A string (possibly abbreviated) indicating the estimation approach to use. Can be "sampling" for MCMC (the default), "optimizing" for optimization, "meanfield" for variational inference with independent normal distributions, or "fullrank" for variational inference with a multivariate normal distribution. See rstanarm-package for more details on the estimation algorithms. NOTE: not all fitting functions support all four algorithms.

adapt_delta

Only relevant if algorithm="sampling". See adapt_delta for details.

A logical scalar (defaulting to FALSE) but if TRUE applies a scaled qr decomposition to the design matrix, $X = Q* R*$, where $Q* = Q (n-1)^0.5$ and $R* = (n-1)^(-0.5) R$. The coefficients relative to $Q*$ are obtained and then premultiplied by the inverse of $R*$ to obtain coefficients relative to the original predictors, $X$. These transformations do not change the likelihood of the data but are recommended for computational reasons when there are multiple predictors. However, because the coefficients relative to $Q*$ are not very interpretable it is hard to specify an informative prior. Setting QR=TRUE is therefore only recommended if you do not have an informative prior for the regression coefficients.

sparse

A logical scalar (defaulting to FALSE) indicating whether to use a sparse representation of the design (X) matrix. Setting this to TRUE will likely be twice as slow, even if the design matrix has a considerable number of zeros, but it may allow the model to be estimated when the computer has too little RAM to utilize a dense design matrix. If TRUE, the the design matrix is not centered (since that would destroy the sparsity) and it is not possible to specify both QR = TRUE and sparse = TRUE.

link

For stan_glm.nb only, the link function to use. See neg_binomial_2.

group

A list, possibly of length zero (the default), but otherwise having the structure of that produced by mkReTrms to indicate the group-specific part of the model. In addition, this list must have elements for the regularization, concentration shape, and scale components of a decov prior for the covariance matrices among the group-specific coefficients.

Value

A stanreg object is returned for stan_glm, stan_glm.nb.A stanfit object (or a slightly modified stanfit object) is returned if stan_glm.fit is called directly.

Details

The stan_glm function is similar in syntax to glm but rather than performing maximum likelihood estimation of generalized linear models, full Bayesian estimation is performed (if algorithm is "sampling") via MCMC. The Bayesian model adds independent priors on the coefficients of the GLM. The stan_glm function calls the workhorse stan_glm.fit function, but it is also possible to call the latter directly. The stan_glm.nb function, which takes the extra argument link, is a simple wrapper for stan_glm with

family =
  neg_binomial_2(link)

References

Gelman, A. and Hill, J. (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press, Cambridge, UK. (Ch. 3-6)

Examples

Run this code

if (!grepl("^sparc",  R.version$platform)) {
### Linear regression
fit <- stan_glm(mpg / 10 ~ ., data = mtcars, QR = TRUE,
                algorithm = "fullrank") # only to make example fast enoug
plot(fit, prob = 0.5)
plot(fit, prob = 0.5, pars = "beta")
}

### Logistic regression
head(wells)
wells$dist100 <- wells$dist / 100
fit2 <- stan_glm(
  switch ~ dist100 + arsenic, 
  data = wells, 
  family = binomial(link = "logit"), 
  prior = student_t(df = 7, location = 0, scale = 2.5), 
  prior_intercept = normal(0, 10)
)
print(fit2)
prior_summary(fit2)

plot(fit2, plotfun = "areas", prob = 0.9, # ?bayesplot::mcmc_areas
     pars = c("(Intercept)", "arsenic"))
pp_check(fit2, plotfun = "error_binned")  # ?bayesplot::ppc_error_binned


### Poisson regression (example from help("glm")) 
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
fit3 <- stan_glm(counts ~ outcome + treatment, family = poisson(link="log"),
                 prior = normal(0, 1), prior_intercept = normal(0, 5))
print(fit3)

bayesplot::color_scheme_set("green")
plot(fit3)
plot(fit3, regex_pars = c("outcome", "treatment"))
plot(fit3, plotfun = "combo", regex_pars = "treatment") # ?bayesplot::mcmc_combo

### Gamma regression (example from help("glm"))
clotting <- data.frame(log_u = log(c(5,10,15,20,30,40,60,80,100)),
                       lot1 = c(118,58,42,35,27,25,21,19,18),
                       lot2 = c(69,35,26,21,18,16,13,12,12))
fit4 <- stan_glm(lot1 ~ log_u, data = clotting, family = Gamma) 
print(fit4, digits = 2)                 
fit5 <- update(fit4, formula = lot2 ~ log_u)

Run the code above in your browser using DataCamp Workspace