-
-
Notifications
You must be signed in to change notification settings - Fork 132
/
stan_lm.R
178 lines (174 loc) · 8.48 KB
/
stan_lm.R
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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
# Part of the rstanarm package for estimating model parameters
# Copyright (C) 2013, 2014, 2015, 2016, 2017 Trustees of Columbia University
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#' Bayesian regularized linear models via Stan
#'
#' \if{html}{\figure{stanlogo.png}{options: width="25" alt="https://mc-stan.org/about/logo/"}}
#' Bayesian inference for linear modeling with regularizing priors on the model
#' parameters that are driven by prior beliefs about \eqn{R^2}, the proportion
#' of variance in the outcome attributable to the predictors. See
#' \code{\link{priors}} for an explanation of this critical point.
#' \code{\link{stan_glm}} with \code{family="gaussian"} also estimates a linear
#' model with normally-distributed errors and allows for various other priors on
#' the coefficients.
#'
#' @export
#' @templateVar fun stan_lm, stan_aov
#' @templateVar fitfun stan_lm.fit or stan_lm.wfit
#' @templateVar pkg stats
#' @templateVar pkgfun lm
#' @templateVar rareargs model,offset,weights
#' @templateVar rareargs2 na.action,singular.ok,contrasts
#' @template return-stanreg-object
#' @template return-stanfit-object
#' @template args-formula-data-subset
#' @template args-same-as-rarely
#' @template args-same-as-rarely-2
#' @template args-x-y
#' @template args-dots
#' @template args-prior_PD
#' @template args-algorithm
#' @template args-adapt_delta
#'
#' @param w Same as in \code{lm.wfit} but rarely specified.
#' @param prior Must be a call to \code{\link{R2}} with its
#' \code{location} argument specified or \code{NULL}, which would
#' indicate a standard uniform prior for the \eqn{R^2}.
#' @param prior_intercept Either \code{NULL} (the default) or a call to
#' \code{\link{normal}}. If a \code{\link{normal}} prior is specified
#' without a \code{scale}, then the standard deviation is taken to be
#' the marginal standard deviation of the outcome divided by the square
#' root of the sample size, which is legitimate because the marginal
#' standard deviation of the outcome is a primitive parameter being
#' estimated.
#'
#' \strong{Note:} If using a dense representation of the design matrix
#' ---i.e., if the \code{sparse} argument is left at its default value of
#' \code{FALSE}--- then the prior distribution for the intercept is set so it
#' applies to the value \emph{when all predictors are centered}. If you prefer
#' to specify a prior on the intercept without the predictors being
#' auto-centered, then you have to omit the intercept from the
#' \code{\link[stats]{formula}} and include a column of ones as a predictor,
#' in which case some element of \code{prior} specifies the prior on it,
#' rather than \code{prior_intercept}. Regardless of how
#' \code{prior_intercept} is specified, the reported \emph{estimates} of the
#' intercept always correspond to a parameterization without centered
#' predictors (i.e., same as in \code{glm}).
#'
#'
#' @details The \code{stan_lm} function is similar in syntax to the
#' \code{\link[stats]{lm}} function but rather than choosing the parameters to
#' minimize the sum of squared residuals, samples from the posterior
#' distribution are drawn using MCMC (if \code{algorithm} is
#' \code{"sampling"}). The \code{stan_lm} function has a formula-based
#' interface and would usually be called by users but the \code{stan_lm.fit}
#' and \code{stan_lm.wfit} functions might be called by other functions that
#' parse the data themselves and are analogous to \code{lm.fit}
#' and \code{lm.wfit} respectively.
#'
#' In addition to estimating \code{sigma} --- the standard deviation of the
#' normally-distributed errors --- this model estimates a positive parameter
#' called \code{log-fit_ratio}. If it is positive, the marginal posterior
#' variance of the outcome will exceed the sample variance of the outcome
#' by a multiplicative factor equal to the square of \code{fit_ratio}.
#' Conversely if \code{log-fit_ratio} is negative, then the model underfits.
#' Given the regularizing nature of the priors, a slight underfit is good.
#'
#' Finally, the posterior predictive distribution is generated with the
#' predictors fixed at their sample means. This quantity is useful for
#' checking convergence because it is reasonably normally distributed
#' and a function of all the parameters in the model.
#'
#' The \code{stan_aov} function is similar to \code{\link[stats]{aov}}, but
#' does a Bayesian analysis of variance that is basically equivalent to
#' \code{stan_lm} with dummy variables. \code{stan_aov} has a somewhat
#' customized \code{\link{print}} method that prints an ANOVA-like table in
#' addition to the output printed for \code{stan_lm} models.
#'
#'
#' @references
#' Lewandowski, D., Kurowicka D., and Joe, H. (2009). Generating random
#' correlation matrices based on vines and extended onion method.
#' \emph{Journal of Multivariate Analysis}. \strong{100}(9), 1989--2001.
#'
#' @seealso
#' The vignettes for \code{stan_lm} and \code{stan_aov}, which have more
#' thorough descriptions and examples.
#' \url{https://mc-stan.org/rstanarm/articles/}
#'
#' Also see \code{\link{stan_glm}}, which --- if \code{family =
#' gaussian(link="identity")} --- also estimates a linear model with
#' normally-distributed errors but specifies different priors.
#'
#'
#' @examples
#' if (.Platform$OS.type != "windows" || .Platform$r_arch !="i386") {
#' (fit <- stan_lm(mpg ~ wt + qsec + am, data = mtcars, prior = R2(0.75),
#' # the next line is only to make the example go fast enough
#' chains = 1, iter = 300, seed = 12345, refresh = 0))
#' plot(fit, "hist", pars = c("wt", "am", "qsec", "sigma"),
#' transformations = list(sigma = "log"))
#' }
stan_lm <- function(formula, data, subset, weights, na.action,
model = TRUE, x = FALSE, y = FALSE,
singular.ok = TRUE, contrasts = NULL, offset, ...,
prior = R2(stop("'location' must be specified")),
prior_intercept = NULL,
prior_PD = FALSE,
algorithm = c("sampling", "meanfield", "fullrank"),
adapt_delta = NULL) {
algorithm <- match.arg(algorithm)
validate_glm_formula(formula)
data <- validate_data(data, if_missing = environment(formula))
call <- match.call(expand.dots = TRUE)
mf <- match.call(expand.dots = FALSE)
mf[[1L]] <- as.name("lm")
mf$data <- data
mf$x <- mf$y <- mf$singular.ok <- TRUE
mf$qr <- FALSE
mf$prior <- mf$prior_intercept <- mf$prior_PD <- mf$algorithm <-
mf$adapt_delta <- NULL
mf$method <- "model.frame"
modelframe <- suppressWarnings(eval(mf, parent.frame()))
mt <- attr(modelframe, "terms")
Y <- model.response(modelframe, "numeric")
X <- model.matrix(mt, modelframe, contrasts)
w <- as.vector(model.weights(modelframe))
offset <- as.vector(model.offset(modelframe))
stanfit <- stan_lm.wfit(y = Y, x = X, w, offset, singular.ok = singular.ok,
prior = prior, prior_intercept = prior_intercept,
prior_PD = prior_PD,
algorithm = algorithm, adapt_delta = adapt_delta,
...)
if (algorithm != "optimizing" && !is(stanfit, "stanfit")) return(stanfit)
fit <- nlist(stanfit, family = gaussian(), formula, offset, weights = w,
x = X[,intersect(colnames(X), dimnames(stanfit)[[3]]), drop = FALSE],
y = Y,
data = data,
prior.info = prior,
algorithm, call, terms = mt,
model = if (model) modelframe else NULL,
na.action = attr(modelframe, "na.action"),
contrasts = attr(X, "contrasts"),
stan_function = "stan_lm")
out <- stanreg(fit)
out$xlevels <- .getXlevels(mt, modelframe)
if (!x)
out$x <- NULL
if (!y)
out$y <- NULL
return(out)
}