Bayesian regularized linear but big models via Stan
Source:R/stan_biglm.R, R/stan_biglm.fit.R
stan_biglm.Rd
This is the same model as with stan_lm but it utilizes the
output from biglm in the biglm package in order to
proceed when the data is too large to fit in memory.
stan_biglm( biglm, xbar, ybar, s_y, ..., prior = R2(stop("'location' must be specified")), prior_intercept = NULL, prior_PD = FALSE, algorithm = c("sampling", "meanfield", "fullrank"), adapt_delta = NULL ) stan_biglm.fit( b, R, SSR, N, xbar, ybar, s_y, has_intercept = TRUE, ..., prior = R2(stop("'location' must be specified")), prior_intercept = NULL, prior_PD = FALSE, algorithm = c("sampling", "meanfield", "fullrank", "optimizing"), adapt_delta = NULL, importance_resampling = TRUE, keep_every = 1 )
Arguments
| biglm | The list output by |
|---|---|
| xbar | A numeric vector of column means in the implicit design matrix excluding the intercept for the observations included in the model. |
| ybar | A numeric scalar indicating the mean of the outcome for the observations included in the model. |
| s_y | A numeric scalar indicating the unbiased sample standard deviation of the outcome for the observations included in the model. |
| ... | Further arguments passed to the function in the rstan
package ( |
| prior | Must be a call to |
| prior_intercept | Either Note: If using a dense representation of the design matrix
---i.e., if the |
| prior_PD | A logical scalar (defaulting to |
| algorithm | A string (possibly abbreviated) indicating the
estimation approach to use. Can be |
| adapt_delta | Only relevant if |
| b | A numeric vector of OLS coefficients, excluding the intercept |
| R | A square upper-triangular matrix from the QR decomposition of the design matrix, excluding the intercept |
| SSR | A numeric scalar indicating the sum-of-squared residuals for OLS |
| N | A integer scalar indicating the number of included observations |
| has_intercept | A logical scalar indicating whether to add an intercept to the model when estimating it. |
| importance_resampling | Logical scalar indicating whether to use
importance resampling when approximating the posterior distribution with
a multivariate normal around the posterior mode, which only applies
when |
| keep_every | Positive integer, which defaults to 1, but can be higher
in order to thin the importance sampling realizations and also only
apples when |
Value
The output of both stan_biglm and stan_biglm.fit is an
object of stanfit-class rather than
stanreg-objects, which is more limited and less convenient
but necessitated by the fact that stan_biglm does not bring the full
design matrix into memory. Without the full design matrix,some of the
elements of a stanreg-objects object cannot be calculated,
such as residuals. Thus, the functions in the rstanarm package that
input stanreg-objects, such as
posterior_predict cannot be used.
Details
The stan_biglm function is intended to be used in the same
circumstances as the biglm function in the biglm
package but with an informative prior on the R2 of the regression.
Like biglm, the memory required to estimate the model
depends largely on the number of predictors rather than the number of
observations. However, stan_biglm and stan_biglm.fit have
additional required arguments that are not necessary in
biglm, namely xbar, ybar, and s_y.
If any observations have any missing values on any of the predictors or the
outcome, such observations do not contribute to these statistics.
Examples
# create inputs ols <- lm(mpg ~ wt + qsec + am, data = mtcars, # all row are complete so ... na.action = na.exclude) # not necessary in this case b <- coef(ols)[-1] R <- qr.R(ols$qr)[-1,-1] SSR <- crossprod(ols$residuals)[1] not_NA <- !is.na(fitted(ols)) N <- sum(not_NA) xbar <- colMeans(mtcars[not_NA,c("wt", "qsec", "am")]) y <- mtcars$mpg[not_NA] ybar <- mean(y) s_y <- sd(y) post <- stan_biglm.fit(b, R, SSR, N, xbar, ybar, s_y, prior = R2(.75), # the next line is only to make the example go fast chains = 1, iter = 500, seed = 12345)#> #> SAMPLING FOR MODEL 'lm' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 1.9e-05 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 500 [ 0%] (Warmup) #> Chain 1: Iteration: 50 / 500 [ 10%] (Warmup) #> Chain 1: Iteration: 100 / 500 [ 20%] (Warmup) #> Chain 1: Iteration: 150 / 500 [ 30%] (Warmup) #> Chain 1: Iteration: 200 / 500 [ 40%] (Warmup) #> Chain 1: Iteration: 250 / 500 [ 50%] (Warmup) #> Chain 1: Iteration: 251 / 500 [ 50%] (Sampling) #> Chain 1: Iteration: 300 / 500 [ 60%] (Sampling) #> Chain 1: Iteration: 350 / 500 [ 70%] (Sampling) #> Chain 1: Iteration: 400 / 500 [ 80%] (Sampling) #> Chain 1: Iteration: 450 / 500 [ 90%] (Sampling) #> Chain 1: Iteration: 500 / 500 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 0.636559 seconds (Warm-up) #> Chain 1: 0.289041 seconds (Sampling) #> Chain 1: 0.9256 seconds (Total) #> Chain 1:#> Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See #> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup#> Warning: Examine the pairs() plot to diagnose sampling problems#> Warning: The largest R-hat is 1.06, indicating chains have not mixed. #> Running the chains for more iterations may help. See #> http://mc-stan.org/misc/warnings.html#r-hat#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mc-stan.org/misc/warnings.html#bulk-ess#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. #> Running the chains for more iterations may help. See #> http://mc-stan.org/misc/warnings.html#tail-ess#> lm stan_lm #> wt -3.916504 -3.776792 #> qsec 1.225886 1.218030 #> am 2.935837 3.002232