File: loo_model_weights.Rd

package info (click to toggle)
r-cran-loo 2.9.0-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 4,836 kB
  • sloc: sh: 15; makefile: 2
file content (237 lines) | stat: -rw-r--r-- 9,878 bytes parent folder | download
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
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/loo_model_weights.R
\name{loo_model_weights}
\alias{loo_model_weights}
\alias{loo_model_weights.default}
\alias{stacking_weights}
\alias{pseudobma_weights}
\title{Model averaging/weighting via stacking or pseudo-BMA weighting}
\usage{
loo_model_weights(x, ...)

\method{loo_model_weights}{default}(
  x,
  ...,
  method = c("stacking", "pseudobma"),
  optim_method = "BFGS",
  optim_control = list(),
  BB = TRUE,
  BB_n = 1000,
  alpha = 1,
  r_eff_list = NULL,
  cores = getOption("mc.cores", 1)
)

stacking_weights(lpd_point, optim_method = "BFGS", optim_control = list())

pseudobma_weights(lpd_point, BB = TRUE, BB_n = 1000, alpha = 1)
}
\arguments{
\item{x}{A list of \code{"psis_loo"} objects (objects returned by \code{\link[=loo]{loo()}}) or
pointwise log-likelihood matrices or , one for each model. If the list
elements are named the names will be used to label the models in the
results. Each matrix/object should have dimensions \eqn{S} by \eqn{N},
where \eqn{S} is the size of the posterior sample (with all chains merged)
and \eqn{N} is the number of data points. If \code{x} is a list of
log-likelihood matrices then \code{\link[=loo]{loo()}} is called internally on each matrix.
Currently the \code{loo_model_weights()} function is not implemented to be used
with results from K-fold CV, but you can still obtain weights using K-fold
CV results by calling the \code{stacking_weights()} or \code{pseudobma_weights()}
function directly.}

\item{...}{Unused, except for the generic to pass arguments to individual
methods.}

\item{method}{Either \code{"stacking"} (the default) or \code{"pseudobma"}, indicating which method
to use for obtaining the weights. \code{"stacking"} refers to stacking of
predictive distributions and  \code{"pseudobma"} refers to pseudo-BMA+ weighting
(or plain pseudo-BMA weighting if argument \code{BB} is \code{FALSE}).}

\item{optim_method}{If \code{method="stacking"}, a string passed to the \code{method}
argument of \code{\link[stats:constrOptim]{stats::constrOptim()}} to specify the optimization algorithm.
The default is \code{optim_method="BFGS"}, but other options are available (see
\code{\link[stats:optim]{stats::optim()}}).}

\item{optim_control}{If \code{method="stacking"}, a list of control parameters for
optimization passed to the \code{control} argument of \code{\link[stats:constrOptim]{stats::constrOptim()}}.}

\item{BB}{Logical used when \code{"method"}=\code{"pseudobma"}. If
\code{TRUE} (the default), the Bayesian bootstrap will be used to adjust
the pseudo-BMA weighting, which is called pseudo-BMA+ weighting. It helps
regularize the weight away from 0 and 1, so as to reduce the variance.}

\item{BB_n}{For pseudo-BMA+ weighting only, the number of samples to use for
the Bayesian bootstrap. The default is \code{BB_n=1000}.}

\item{alpha}{Positive scalar shape parameter in the Dirichlet distribution
used for the Bayesian bootstrap. The default is \code{alpha=1}, which
corresponds to a uniform distribution on the simplex space.}

\item{r_eff_list}{Optionally, a list of relative effective sample size
estimates for the likelihood \code{(exp(log_lik))} of each observation in
each model. See \code{\link[=psis]{psis()}} and  \code{\link[=relative_eff]{relative_eff()}} helper
function for computing \code{r_eff}. If \code{x} is a list of \code{"psis_loo"}
objects then \code{r_eff_list} is ignored.}

\item{cores}{The number of cores to use for parallelization. This defaults to
the option \code{mc.cores} which can be set for an entire R session by
\code{options(mc.cores = NUMBER)}. The old option \code{loo.cores} is now
deprecated but will be given precedence over \code{mc.cores} until
\code{loo.cores} is removed in a future release. \strong{As of version
2.0.0 the default is now 1 core if \code{mc.cores} is not set}, but we
recommend using as many (or close to as many) cores as possible.
\itemize{
\item Note for Windows 10 users: it is \strong{strongly}
\href{https://github.com/stan-dev/loo/issues/94}{recommended} to avoid using
the \code{.Rprofile} file to set \code{mc.cores} (using the \code{cores} argument or
setting \code{mc.cores} interactively or in a script is fine).
}}

\item{lpd_point}{If calling \code{stacking_weights()} or \code{pseudobma_weights()}
directly, a matrix of pointwise leave-one-out (or K-fold) log likelihoods
evaluated for different models. It should be a \eqn{N} by \eqn{K}  matrix
where \eqn{N} is sample size and \eqn{K} is the number of models. Each
column corresponds to one model. These values can be calculated
approximately using \code{\link[=loo]{loo()}} or by running exact leave-one-out or K-fold
cross-validation.}
}
\value{
A numeric vector containing one weight for each model.
}
\description{
Model averaging via stacking of predictive distributions, pseudo-BMA
weighting or pseudo-BMA+ weighting with the Bayesian bootstrap. See Yao et
al. (2018), Vehtari, Gelman, and Gabry (2017), and Vehtari, Simpson,
Gelman, Yao, and Gabry (2024) for background.
}
\details{
\code{loo_model_weights()} is a wrapper around the \code{stacking_weights()} and
\code{pseudobma_weights()} functions that implements stacking, pseudo-BMA, and
pseudo-BMA+ weighting for combining multiple predictive distributions. We can
use approximate or exact leave-one-out cross-validation (LOO-CV) or K-fold CV
to estimate the expected log predictive density (ELPD).

The stacking method (\code{method="stacking"}), which is the default for
\code{loo_model_weights()}, combines all models by maximizing the leave-one-out
predictive density of the combination distribution. That is, it finds the
optimal linear combining weights for maximizing the leave-one-out log score.

The pseudo-BMA method (\code{method="pseudobma"}) finds the relative weights
proportional to the ELPD of each model. However, when
\code{method="pseudobma"}, the default is to also use the Bayesian bootstrap
(\code{BB=TRUE}), which corresponds to the pseudo-BMA+ method. The Bayesian
bootstrap  takes into account the uncertainty of finite data points and
regularizes the weights away from the extremes of 0 and 1.

In general, we recommend stacking for averaging predictive distributions,
while pseudo-BMA+ can serve as a computationally easier alternative.
}
\examples{
\dontrun{
### Demonstrating usage after fitting models with RStan
library(rstan)

# generate fake data from N(0,1).
N <- 100
y <- rnorm(N, 0, 1)

# Suppose we have three models: N(-1, sigma), N(0.5, sigma) and N(0.6,sigma).
stan_code <- "
  data {
    int N;
    vector[N] y;
    real mu_fixed;
  }
  parameters {
    real<lower=0> sigma;
  }
  model {
    sigma ~ exponential(1);
    y ~ normal(mu_fixed, sigma);
  }
  generated quantities {
    vector[N] log_lik;
    for (n in 1:N) log_lik[n] = normal_lpdf(y[n]| mu_fixed, sigma);
  }"

mod <- stan_model(model_code = stan_code)
fit1 <- sampling(mod, data=list(N=N, y=y, mu_fixed=-1))
fit2 <- sampling(mod, data=list(N=N, y=y, mu_fixed=0.5))
fit3 <- sampling(mod, data=list(N=N, y=y, mu_fixed=0.6))
model_list <- list(fit1, fit2, fit3)
log_lik_list <- lapply(model_list, extract_log_lik)

# optional but recommended
r_eff_list <- lapply(model_list, function(x) {
  ll_array <- extract_log_lik(x, merge_chains = FALSE)
  relative_eff(exp(ll_array))
})

# stacking method:
wts1 <- loo_model_weights(
  log_lik_list,
  method = "stacking",
  r_eff_list = r_eff_list,
  optim_control = list(reltol=1e-10)
)
print(wts1)

# can also pass a list of psis_loo objects to avoid recomputing loo
loo_list <- lapply(1:length(log_lik_list), function(j) {
  loo(log_lik_list[[j]], r_eff = r_eff_list[[j]])
})

wts2 <- loo_model_weights(
  loo_list,
  method = "stacking",
  optim_control = list(reltol=1e-10)
)
all.equal(wts1, wts2)

# can provide names to be used in the results
loo_model_weights(setNames(loo_list, c("A", "B", "C")))


# pseudo-BMA+ method:
set.seed(1414)
loo_model_weights(loo_list, method = "pseudobma")

# pseudo-BMA method (set BB = FALSE):
loo_model_weights(loo_list, method = "pseudobma", BB = FALSE)

# calling stacking_weights or pseudobma_weights directly
lpd1 <- loo(log_lik_list[[1]], r_eff = r_eff_list[[1]])$pointwise[,1]
lpd2 <- loo(log_lik_list[[2]], r_eff = r_eff_list[[2]])$pointwise[,1]
lpd3 <- loo(log_lik_list[[3]], r_eff = r_eff_list[[3]])$pointwise[,1]
stacking_weights(cbind(lpd1, lpd2, lpd3))
pseudobma_weights(cbind(lpd1, lpd2, lpd3))
pseudobma_weights(cbind(lpd1, lpd2, lpd3), BB = FALSE)
}

}
\references{
Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model
evaluation using leave-one-out cross-validation and WAIC.
\emph{Statistics and Computing}. 27(5), 1413--1432. doi:10.1007/s11222-016-9696-4
(\href{https://link.springer.com/article/10.1007/s11222-016-9696-4}{journal version},
\href{https://arxiv.org/abs/1507.04544}{preprint arXiv:1507.04544}).

Vehtari, A., Simpson, D., Gelman, A., Yao, Y., and Gabry, J. (2024).
Pareto smoothed importance sampling. \emph{Journal of Machine Learning Research},
25(72):1-58.
\href{https://jmlr.org/papers/v25/19-556.html}{PDF}

Yao, Y., Vehtari, A., Simpson, D., and Gelman, A. (2018) Using
stacking to average Bayesian predictive distributions.
\emph{Bayesian Analysis}, advance publication,  doi:10.1214/17-BA1091.
(\href{https://projecteuclid.org/euclid.ba/1516093227}{online}).
}
\seealso{
\itemize{
\item The \strong{loo} package \href{https://mc-stan.org/loo/articles/}{vignettes}, particularly
\href{https://mc-stan.org/loo/articles/loo2-weights.html}{Bayesian Stacking and Pseudo-BMA weights using the \strong{loo} package}.
\item \code{\link[=loo]{loo()}} for details on leave-one-out ELPD estimation.
\item \code{\link[=constrOptim]{constrOptim()}} for the choice of optimization methods and control-parameters.
\item \code{\link[=relative_eff]{relative_eff()}} for computing \code{r_eff}.
}
}