File: loo2-weights.html

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 (713 lines) | stat: -rw-r--r-- 50,470 bytes
<!DOCTYPE html>

<html>

<head>

<meta charset="utf-8" />
<meta name="generator" content="pandoc" />
<meta http-equiv="X-UA-Compatible" content="IE=EDGE" />

<meta name="viewport" content="width=device-width, initial-scale=1" />

<meta name="author" content="Aki Vehtari and Jonah Gabry" />

<meta name="date" content="2025-12-22" />

<title>Bayesian Stacking and Pseudo-BMA weights using the loo package</title>

<script>// Pandoc 2.9 adds attributes on both header and div. We remove the former (to
// be compatible with the behavior of Pandoc < 2.8).
document.addEventListener('DOMContentLoaded', function(e) {
  var hs = document.querySelectorAll("div.section[class*='level'] > :first-child");
  var i, h, a;
  for (i = 0; i < hs.length; i++) {
    h = hs[i];
    if (!/^h[1-6]$/i.test(h.tagName)) continue;  // it should be a header h1-h6
    a = h.attributes;
    while (a.length > 0) h.removeAttribute(a[0].name);
  }
});
</script>

<style type="text/css">
  code{white-space: pre-wrap;}
  span.smallcaps{font-variant: small-caps;}
  span.underline{text-decoration: underline;}
  div.column{display: inline-block; vertical-align: top; width: 50%;}
  div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
  ul.task-list{list-style: none;}
    </style>



<style type="text/css">
  code {
    white-space: pre;
  }
  .sourceCode {
    overflow: visible;
  }
</style>
<style type="text/css" data-origin="pandoc">
pre > code.sourceCode { white-space: pre; position: relative; }
pre > code.sourceCode > span { display: inline-block; line-height: 1.25; }
pre > code.sourceCode > span:empty { height: 1.2em; }
code.sourceCode > span { color: inherit; text-decoration: inherit; }
div.sourceCode { margin: 1em 0; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
pre > code.sourceCode { white-space: pre-wrap; }
pre > code.sourceCode > span { text-indent: -5em; padding-left: 5em; }
}
pre.numberSource code
  { counter-reset: source-line 0; }
pre.numberSource code > span
  { position: relative; left: -4em; counter-increment: source-line; }
pre.numberSource code > span > a:first-child::before
  { content: counter(source-line);
    position: relative; left: -1em; text-align: right; vertical-align: baseline;
    border: none; display: inline-block;
    -webkit-touch-callout: none; -webkit-user-select: none;
    -khtml-user-select: none; -moz-user-select: none;
    -ms-user-select: none; user-select: none;
    padding: 0 4px; width: 4em;
    color: #aaaaaa;
  }
pre.numberSource { margin-left: 3em; border-left: 1px solid #aaaaaa;  padding-left: 4px; }
div.sourceCode
  {   }
@media screen {
pre > code.sourceCode > span > a:first-child::before { text-decoration: underline; }
}
code span.al { color: #ff0000; font-weight: bold; } /* Alert */
code span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code span.at { color: #7d9029; } /* Attribute */
code span.bn { color: #40a070; } /* BaseN */
code span.bu { } /* BuiltIn */
code span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code span.ch { color: #4070a0; } /* Char */
code span.cn { color: #880000; } /* Constant */
code span.co { color: #60a0b0; font-style: italic; } /* Comment */
code span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code span.do { color: #ba2121; font-style: italic; } /* Documentation */
code span.dt { color: #902000; } /* DataType */
code span.dv { color: #40a070; } /* DecVal */
code span.er { color: #ff0000; font-weight: bold; } /* Error */
code span.ex { } /* Extension */
code span.fl { color: #40a070; } /* Float */
code span.fu { color: #06287e; } /* Function */
code span.im { } /* Import */
code span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
code span.kw { color: #007020; font-weight: bold; } /* Keyword */
code span.op { color: #666666; } /* Operator */
code span.ot { color: #007020; } /* Other */
code span.pp { color: #bc7a00; } /* Preprocessor */
code span.sc { color: #4070a0; } /* SpecialChar */
code span.ss { color: #bb6688; } /* SpecialString */
code span.st { color: #4070a0; } /* String */
code span.va { color: #19177c; } /* Variable */
code span.vs { color: #4070a0; } /* VerbatimString */
code span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */

</style>
<script>
// apply pandoc div.sourceCode style to pre.sourceCode instead
(function() {
  var sheets = document.styleSheets;
  for (var i = 0; i < sheets.length; i++) {
    if (sheets[i].ownerNode.dataset["origin"] !== "pandoc") continue;
    try { var rules = sheets[i].cssRules; } catch (e) { continue; }
    var j = 0;
    while (j < rules.length) {
      var rule = rules[j];
      // check if there is a div.sourceCode rule
      if (rule.type !== rule.STYLE_RULE || rule.selectorText !== "div.sourceCode") {
        j++;
        continue;
      }
      var style = rule.style.cssText;
      // check if color or background-color is set
      if (rule.style.color === '' && rule.style.backgroundColor === '') {
        j++;
        continue;
      }
      // replace div.sourceCode by a pre.sourceCode rule
      sheets[i].deleteRule(j);
      sheets[i].insertRule('pre.sourceCode{' + style + '}', j);
    }
  }
})();
</script>




<style type="text/css">body {
background-color: #fff;
margin: 1em auto;
max-width: 700px;
overflow: visible;
padding-left: 2em;
padding-right: 2em;
font-family: "Open Sans", "Helvetica Neue", Helvetica, Arial, sans-serif;
font-size: 14px;
line-height: 1.35;
}
#TOC {
clear: both;
margin: 0 0 10px 10px;
padding: 4px;
width: 400px;
border: 1px solid #CCCCCC;
border-radius: 5px;
background-color: #f6f6f6;
font-size: 13px;
line-height: 1.3;
}
#TOC .toctitle {
font-weight: bold;
font-size: 15px;
margin-left: 5px;
}
#TOC ul {
padding-left: 40px;
margin-left: -1.5em;
margin-top: 5px;
margin-bottom: 5px;
}
#TOC ul ul {
margin-left: -2em;
}
#TOC li {
line-height: 16px;
}
table {
margin: 1em auto;
border-width: 1px;
border-color: #DDDDDD;
border-style: outset;
border-collapse: collapse;
}
table th {
border-width: 2px;
padding: 5px;
border-style: inset;
}
table td {
border-width: 1px;
border-style: inset;
line-height: 18px;
padding: 5px 5px;
}
table, table th, table td {
border-left-style: none;
border-right-style: none;
}
table thead, table tr.even {
background-color: #f7f7f7;
}
p {
margin: 0.5em 0;
}
blockquote {
background-color: #f6f6f6;
padding: 0.25em 0.75em;
}
hr {
border-style: solid;
border: none;
border-top: 1px solid #777;
margin: 28px 0;
}
dl {
margin-left: 0;
}
dl dd {
margin-bottom: 13px;
margin-left: 13px;
}
dl dt {
font-weight: bold;
}
ul {
margin-top: 0;
}
ul li {
list-style: circle outside;
}
ul ul {
margin-bottom: 0;
}
pre, code {
background-color: #f7f7f7;
border-radius: 3px;
color: #333;
white-space: pre-wrap; 
}
pre {
border-radius: 3px;
margin: 5px 0px 10px 0px;
padding: 10px;
}
pre:not([class]) {
background-color: #f7f7f7;
}
code {
font-family: Consolas, Monaco, 'Courier New', monospace;
font-size: 85%;
}
p > code, li > code {
padding: 2px 0px;
}
div.figure {
text-align: center;
}
img {
background-color: #FFFFFF;
padding: 2px;
border: 1px solid #DDDDDD;
border-radius: 3px;
border: 1px solid #CCCCCC;
margin: 0 5px;
}
h1 {
margin-top: 0;
font-size: 35px;
line-height: 40px;
}
h2 {
border-bottom: 4px solid #f7f7f7;
padding-top: 10px;
padding-bottom: 2px;
font-size: 145%;
}
h3 {
border-bottom: 2px solid #f7f7f7;
padding-top: 10px;
font-size: 120%;
}
h4 {
border-bottom: 1px solid #f7f7f7;
margin-left: 8px;
font-size: 105%;
}
h5, h6 {
border-bottom: 1px solid #ccc;
font-size: 105%;
}
a {
color: #0033dd;
text-decoration: none;
}
a:hover {
color: #6666ff; }
a:visited {
color: #800080; }
a:visited:hover {
color: #BB00BB; }
a[href^="http:"] {
text-decoration: underline; }
a[href^="https:"] {
text-decoration: underline; }

code > span.kw { color: #555; font-weight: bold; } 
code > span.dt { color: #902000; } 
code > span.dv { color: #40a070; } 
code > span.bn { color: #d14; } 
code > span.fl { color: #d14; } 
code > span.ch { color: #d14; } 
code > span.st { color: #d14; } 
code > span.co { color: #888888; font-style: italic; } 
code > span.ot { color: #007020; } 
code > span.al { color: #ff0000; font-weight: bold; } 
code > span.fu { color: #900; font-weight: bold; } 
code > span.er { color: #a61717; background-color: #e3d2d2; } 
</style>




</head>

<body>




<h1 class="title toc-ignore">Bayesian Stacking and Pseudo-BMA weights using the loo package</h1>
<h4 class="author">Aki Vehtari and Jonah Gabry</h4>
<h4 class="date">2025-12-22</h4>


<div id="TOC">
<ul>
<li><a href="#introduction">Introduction</a></li>
<li><a href="#setup">Setup</a></li>
<li><a href="#example-primate-milk">Example: Primate milk</a></li>
<li><a href="#example-oceanic-tool-complexity">Example: Oceanic tool complexity</a></li>
<li><a href="#simpler-coding-using-loo_model_weights-function">Simpler coding using <code>loo_model_weights</code> function</a></li>
<li><a href="#references">References</a></li>
</ul>
</div>

<!--
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{Bayesian Stacking and Pseudo-BMA weights}
-->
<div id="introduction" class="section level1">
<h1>Introduction</h1>
<p>This vignette demonstrates the new functionality in <strong>loo</strong> v2.0.0 for Bayesian stacking and Pseudo-BMA weighting. In this vignette we can’t provide all of the necessary background on this topic, so we encourage readers to refer to the paper</p>
<ul>
<li>Yao, Y., Vehtari, A., Simpson, D., and Gelman, A. (2018). Using stacking to average Bayesian predictive distributions. In Bayesian Analysis, :10.1214/17-BA1091. <a href="https://projecteuclid.org/euclid.ba/1516093227">Online</a></li>
</ul>
<p>which provides important details on the methods demonstrated in this vignette. Here we just quote from the abstract of the paper:</p>
<blockquote>
<p><strong>Abstract</strong>: Bayesian model averaging is flawed in the <span class="math inline">\(\mathcal{M}\)</span>-open setting in which the true data-generating process is not one of the candidate models being fit. We take the idea of stacking from the point estimation literature and generalize to the combination of predictive distributions. We extend the utility function to any proper scoring rule and use Pareto smoothed importance sampling to efficiently compute the required leave-one-out posterior distributions. We compare stacking of predictive distributions to several alternatives: stacking of means, Bayesian model averaging (BMA), Pseudo-BMA, and a variant of Pseudo-BMA that is stabilized using the Bayesian bootstrap. Based on simulations and real-data applications, we recommend stacking of predictive distributions, with bootstrapped-Pseudo-BMA as an approximate alternative when computation cost is an issue.</p>
</blockquote>
<p>Ideally, we would avoid the Bayesian model combination problem by extending the model to include the separate models as special cases, and preferably as a continuous expansion of the model space. For example, instead of model averaging over different covariate combinations, all potentially relevant covariates should be included in a predictive model (for causal analysis more care is needed) and a prior assumption that only some of the covariates are relevant can be presented with regularized horseshoe prior (Piironen and Vehtari, 2017a). For variable selection we recommend projective predictive variable selection (Piironen and Vehtari, 2017a; <a href="https://cran.r-project.org/package=projpred"><strong>projpred</strong> package</a>).</p>
<p>To demonstrate how to use <strong>loo</strong> package to compute Bayesian stacking and Pseudo-BMA weights, we repeat two simple model averaging examples from Chapters 6 and 10 of <em>Statistical Rethinking</em> by Richard McElreath. In <em>Statistical Rethinking</em> WAIC is used to form weights which are similar to classical “Akaike weights”. Pseudo-BMA weighting using PSIS-LOO for computation is close to these WAIC weights, but named after the Pseudo Bayes Factor by Geisser and Eddy (1979). As discussed below, in general we prefer using stacking rather than WAIC weights or the similar pseudo-BMA weights.</p>
</div>
<div id="setup" class="section level1">
<h1>Setup</h1>
<p>In addition to the <strong>loo</strong> package we will also load the <strong>rstanarm</strong> package for fitting the models.</p>
<div class="sourceCode" id="cb1"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(rstanarm)</span>
<span id="cb1-2"><a href="#cb1-2" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(loo)</span></code></pre></div>
</div>
<div id="example-primate-milk" class="section level1">
<h1>Example: Primate milk</h1>
<p>In <em>Statistical Rethinking</em>, McElreath describes the data for the primate milk example as follows:</p>
<blockquote>
<p>A popular hypothesis has it that primates with larger brains produce more energetic milk, so that brains can grow quickly. … The question here is to what extent energy content of milk, measured here by kilocalories, is related to the percent of the brain mass that is neocortex. … We’ll end up needing female body mass as well, to see the masking that hides the relationships among the variables.</p>
</blockquote>
<div class="sourceCode" id="cb2"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true" tabindex="-1"></a><span class="fu">data</span>(milk)</span>
<span id="cb2-2"><a href="#cb2-2" aria-hidden="true" tabindex="-1"></a>d <span class="ot">&lt;-</span> milk[<span class="fu">complete.cases</span>(milk),]</span>
<span id="cb2-3"><a href="#cb2-3" aria-hidden="true" tabindex="-1"></a>d<span class="sc">$</span>neocortex <span class="ot">&lt;-</span> d<span class="sc">$</span>neocortex.perc <span class="sc">/</span><span class="dv">100</span></span>
<span id="cb2-4"><a href="#cb2-4" aria-hidden="true" tabindex="-1"></a><span class="fu">str</span>(d)</span></code></pre></div>
<pre><code>&#39;data.frame&#39;:   17 obs. of  9 variables:
 $ clade         : Factor w/ 4 levels &quot;Ape&quot;,&quot;New World Monkey&quot;,..: 4 2 2 2 2 2 2 2 3 3 ...
 $ species       : Factor w/ 29 levels &quot;A palliata&quot;,&quot;Alouatta seniculus&quot;,..: 11 2 1 6 27 5 3 4 21 19 ...
 $ kcal.per.g    : num  0.49 0.47 0.56 0.89 0.92 0.8 0.46 0.71 0.68 0.97 ...
 $ perc.fat      : num  16.6 21.2 29.7 53.4 50.6 ...
 $ perc.protein  : num  15.4 23.6 23.5 15.8 22.3 ...
 $ perc.lactose  : num  68 55.2 46.9 30.8 27.1 ...
 $ mass          : num  1.95 5.25 5.37 2.51 0.68 0.12 0.47 0.32 1.55 3.24 ...
 $ neocortex.perc: num  55.2 64.5 64.5 67.6 68.8 ...
 $ neocortex     : num  0.552 0.645 0.645 0.676 0.688 ...</code></pre>
<p>We repeat the analysis in Chapter 6 of <em>Statistical Rethinking</em> using the following four models (here we use the default weakly informative priors in <strong>rstanarm</strong>, while flat priors were used in <em>Statistical Rethinking</em>).</p>
<div class="sourceCode" id="cb4"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb4-1"><a href="#cb4-1" aria-hidden="true" tabindex="-1"></a>fit1 <span class="ot">&lt;-</span> <span class="fu">stan_glm</span>(kcal.per.g <span class="sc">~</span> <span class="dv">1</span>, <span class="at">data =</span> d, <span class="at">seed =</span> <span class="dv">2030</span>)</span>
<span id="cb4-2"><a href="#cb4-2" aria-hidden="true" tabindex="-1"></a>fit2 <span class="ot">&lt;-</span> <span class="fu">update</span>(fit1, <span class="at">formula =</span> kcal.per.g <span class="sc">~</span> neocortex)</span>
<span id="cb4-3"><a href="#cb4-3" aria-hidden="true" tabindex="-1"></a>fit3 <span class="ot">&lt;-</span> <span class="fu">update</span>(fit1, <span class="at">formula =</span> kcal.per.g <span class="sc">~</span> <span class="fu">log</span>(mass))</span>
<span id="cb4-4"><a href="#cb4-4" aria-hidden="true" tabindex="-1"></a>fit4 <span class="ot">&lt;-</span> <span class="fu">update</span>(fit1, <span class="at">formula =</span> kcal.per.g <span class="sc">~</span> neocortex <span class="sc">+</span> <span class="fu">log</span>(mass))</span></code></pre></div>
<p>McElreath uses WAIC for model comparison and averaging, so we’ll start by also computing WAIC for these models so we can compare the results to the other options presented later in the vignette. The <strong>loo</strong> package provides <code>waic</code> methods for log-likelihood arrays, matrices and functions. Since we fit our model with rstanarm we can use the <code>waic</code> method provided by the <strong>rstanarm</strong> package (a wrapper around <code>waic</code> from the <strong>loo</strong> package), which allows us to just pass in our fitted model objects instead of first extracting the log-likelihood values.</p>
<div class="sourceCode" id="cb5"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb5-1"><a href="#cb5-1" aria-hidden="true" tabindex="-1"></a>waic1 <span class="ot">&lt;-</span> <span class="fu">waic</span>(fit1)</span>
<span id="cb5-2"><a href="#cb5-2" aria-hidden="true" tabindex="-1"></a>waic2 <span class="ot">&lt;-</span> <span class="fu">waic</span>(fit2)</span>
<span id="cb5-3"><a href="#cb5-3" aria-hidden="true" tabindex="-1"></a>waic3 <span class="ot">&lt;-</span> <span class="fu">waic</span>(fit3)</span></code></pre></div>
<pre><code>Warning: 
1 (5.9%) p_waic estimates greater than 0.4. We recommend trying loo instead.</code></pre>
<div class="sourceCode" id="cb7"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb7-1"><a href="#cb7-1" aria-hidden="true" tabindex="-1"></a>waic4 <span class="ot">&lt;-</span> <span class="fu">waic</span>(fit4)</span></code></pre></div>
<pre><code>Warning: 
2 (11.8%) p_waic estimates greater than 0.4. We recommend trying loo instead.</code></pre>
<div class="sourceCode" id="cb9"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb9-1"><a href="#cb9-1" aria-hidden="true" tabindex="-1"></a>waics <span class="ot">&lt;-</span> <span class="fu">c</span>(</span>
<span id="cb9-2"><a href="#cb9-2" aria-hidden="true" tabindex="-1"></a>  waic1<span class="sc">$</span>estimates[<span class="st">&quot;elpd_waic&quot;</span>, <span class="dv">1</span>],</span>
<span id="cb9-3"><a href="#cb9-3" aria-hidden="true" tabindex="-1"></a>  waic2<span class="sc">$</span>estimates[<span class="st">&quot;elpd_waic&quot;</span>, <span class="dv">1</span>],</span>
<span id="cb9-4"><a href="#cb9-4" aria-hidden="true" tabindex="-1"></a>  waic3<span class="sc">$</span>estimates[<span class="st">&quot;elpd_waic&quot;</span>, <span class="dv">1</span>],</span>
<span id="cb9-5"><a href="#cb9-5" aria-hidden="true" tabindex="-1"></a>  waic4<span class="sc">$</span>estimates[<span class="st">&quot;elpd_waic&quot;</span>, <span class="dv">1</span>]</span>
<span id="cb9-6"><a href="#cb9-6" aria-hidden="true" tabindex="-1"></a>)</span></code></pre></div>
<p>We get some warnings when computing WAIC for models 3 and 4, indicating that we shouldn’t trust the WAIC weights we will compute later. Following the recommendation in the warning, we next use the <code>loo</code> methods to compute PSIS-LOO instead. The <strong>loo</strong> package provides <code>loo</code> methods for log-likelihood arrays, matrices, and functions, but since we fit our model with <strong>rstanarm</strong> we can just pass the fitted model objects directly and <strong>rstanarm</strong> will extract the needed values to pass to the <strong>loo</strong> package. (Like <strong>rstanarm</strong>, some other R packages for fitting Stan models, e.g. <strong>brms</strong>, also provide similar methods for interfacing with the <strong>loo</strong> package.)</p>
<div class="sourceCode" id="cb10"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb10-1"><a href="#cb10-1" aria-hidden="true" tabindex="-1"></a><span class="co"># note: the loo function accepts a &#39;cores&#39; argument that we recommend specifying</span></span>
<span id="cb10-2"><a href="#cb10-2" aria-hidden="true" tabindex="-1"></a><span class="co"># when working with bigger datasets</span></span>
<span id="cb10-3"><a href="#cb10-3" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb10-4"><a href="#cb10-4" aria-hidden="true" tabindex="-1"></a>loo1 <span class="ot">&lt;-</span> <span class="fu">loo</span>(fit1)</span>
<span id="cb10-5"><a href="#cb10-5" aria-hidden="true" tabindex="-1"></a>loo2 <span class="ot">&lt;-</span> <span class="fu">loo</span>(fit2)</span>
<span id="cb10-6"><a href="#cb10-6" aria-hidden="true" tabindex="-1"></a>loo3 <span class="ot">&lt;-</span> <span class="fu">loo</span>(fit3)</span>
<span id="cb10-7"><a href="#cb10-7" aria-hidden="true" tabindex="-1"></a>loo4 <span class="ot">&lt;-</span> <span class="fu">loo</span>(fit4)</span>
<span id="cb10-8"><a href="#cb10-8" aria-hidden="true" tabindex="-1"></a>lpd_point <span class="ot">&lt;-</span> <span class="fu">cbind</span>(</span>
<span id="cb10-9"><a href="#cb10-9" aria-hidden="true" tabindex="-1"></a>  loo1<span class="sc">$</span>pointwise[,<span class="st">&quot;elpd_loo&quot;</span>], </span>
<span id="cb10-10"><a href="#cb10-10" aria-hidden="true" tabindex="-1"></a>  loo2<span class="sc">$</span>pointwise[,<span class="st">&quot;elpd_loo&quot;</span>],</span>
<span id="cb10-11"><a href="#cb10-11" aria-hidden="true" tabindex="-1"></a>  loo3<span class="sc">$</span>pointwise[,<span class="st">&quot;elpd_loo&quot;</span>], </span>
<span id="cb10-12"><a href="#cb10-12" aria-hidden="true" tabindex="-1"></a>  loo4<span class="sc">$</span>pointwise[,<span class="st">&quot;elpd_loo&quot;</span>]</span>
<span id="cb10-13"><a href="#cb10-13" aria-hidden="true" tabindex="-1"></a>)</span></code></pre></div>
<p>With <code>loo</code> we don’t get any warnings for models 3 and 4, but for illustration of good results, we display the diagnostic details for these models anyway.</p>
<div class="sourceCode" id="cb11"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb11-1"><a href="#cb11-1" aria-hidden="true" tabindex="-1"></a><span class="fu">print</span>(loo3)</span></code></pre></div>
<pre><code>
Computed from 4000 by 17 log-likelihood matrix.

         Estimate  SE
elpd_loo      4.5 2.3
p_loo         2.1 0.5
looic        -9.1 4.6
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.0]).

All Pareto k estimates are good (k &lt; 0.7).
See help(&#39;pareto-k-diagnostic&#39;) for details.</code></pre>
<div class="sourceCode" id="cb13"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb13-1"><a href="#cb13-1" aria-hidden="true" tabindex="-1"></a><span class="fu">print</span>(loo4)</span></code></pre></div>
<pre><code>
Computed from 4000 by 17 log-likelihood matrix.

         Estimate  SE
elpd_loo      8.4 2.8
p_loo         3.3 0.9
looic       -16.8 5.5
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.0]).

All Pareto k estimates are good (k &lt; 0.7).
See help(&#39;pareto-k-diagnostic&#39;) for details.</code></pre>
<p>One benefit of PSIS-LOO over WAIC is better diagnostics. Here for both models 3 and 4 all <span class="math inline">\(k&lt;0.7\)</span> and the Monte Carlo SE of <code>elpd_loo</code> is 0.1 or less, and we can expect the model comparison to be reliable.</p>
<p>Next we compute and compare 1) WAIC weights, 2) Pseudo-BMA weights without Bayesian bootstrap, 3) Pseudo-BMA+ weights with Bayesian bootstrap, and 4) Bayesian stacking weights.</p>
<div class="sourceCode" id="cb15"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb15-1"><a href="#cb15-1" aria-hidden="true" tabindex="-1"></a>waic_wts <span class="ot">&lt;-</span> <span class="fu">exp</span>(waics) <span class="sc">/</span> <span class="fu">sum</span>(<span class="fu">exp</span>(waics))</span>
<span id="cb15-2"><a href="#cb15-2" aria-hidden="true" tabindex="-1"></a>pbma_wts <span class="ot">&lt;-</span> <span class="fu">pseudobma_weights</span>(lpd_point, <span class="at">BB=</span><span class="cn">FALSE</span>)</span>
<span id="cb15-3"><a href="#cb15-3" aria-hidden="true" tabindex="-1"></a>pbma_BB_wts <span class="ot">&lt;-</span> <span class="fu">pseudobma_weights</span>(lpd_point) <span class="co"># default is BB=TRUE</span></span>
<span id="cb15-4"><a href="#cb15-4" aria-hidden="true" tabindex="-1"></a>stacking_wts <span class="ot">&lt;-</span> <span class="fu">stacking_weights</span>(lpd_point)</span>
<span id="cb15-5"><a href="#cb15-5" aria-hidden="true" tabindex="-1"></a><span class="fu">round</span>(<span class="fu">cbind</span>(waic_wts, pbma_wts, pbma_BB_wts, stacking_wts), <span class="dv">2</span>)</span></code></pre></div>
<pre><code>       waic_wts pbma_wts pbma_BB_wts stacking_wts
model1     0.01     0.02        0.07         0.01
model2     0.01     0.01        0.04         0.00
model3     0.02     0.02        0.04         0.00
model4     0.96     0.95        0.85         0.99</code></pre>
<p>With all approaches Model 4 with <code>neocortex</code> and <code>log(mass)</code> gets most of the weight. Based on theory, Pseudo-BMA weights without Bayesian bootstrap should be close to WAIC weights, and we can also see that here. Pseudo-BMA+ weights with Bayesian bootstrap provide more cautious weights further away from 0 and 1 (see Yao et al. (2018) for a discussion of why this can be beneficial and results from related experiments). In this particular example, the Bayesian stacking weights are not much different from the other weights.</p>
<p>One of the benefits of stacking is that it manages well if there are many similar models. Consider for example that there could be many irrelevant covariates that when included would produce a similar model to one of the existing models. To emulate this situation here we simply copy the first model a bunch of times, but you can imagine that instead we would have ten alternative models with about the same predictive performance. WAIC weights for such a scenario would be close to the following:</p>
<div class="sourceCode" id="cb17"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb17-1"><a href="#cb17-1" aria-hidden="true" tabindex="-1"></a>waic_wts_demo <span class="ot">&lt;-</span> </span>
<span id="cb17-2"><a href="#cb17-2" aria-hidden="true" tabindex="-1"></a>  <span class="fu">exp</span>(waics[<span class="fu">c</span>(<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">2</span>,<span class="dv">3</span>,<span class="dv">4</span>)]) <span class="sc">/</span></span>
<span id="cb17-3"><a href="#cb17-3" aria-hidden="true" tabindex="-1"></a>  <span class="fu">sum</span>(<span class="fu">exp</span>(waics[<span class="fu">c</span>(<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">2</span>,<span class="dv">3</span>,<span class="dv">4</span>)]))</span>
<span id="cb17-4"><a href="#cb17-4" aria-hidden="true" tabindex="-1"></a><span class="fu">round</span>(waic_wts_demo, <span class="dv">3</span>)</span></code></pre></div>
<pre><code> [1] 0.013 0.013 0.013 0.013 0.013 0.013 0.013 0.013 0.013 0.013 0.006 0.016
[13] 0.847</code></pre>
<p>Notice how much the weight for model 4 is lowered now that more models similar to model 1 (or in this case identical) have been added. Both WAIC weights and Pseudo-BMA approaches first estimate the predictive performance separately for each model and then compute weights based on estimated relative predictive performances. Similar models share similar weights so the weights of other models must be reduced for the total sum of the weights to remain the same.</p>
<p>On the other hand, stacking optimizes the weights <em>jointly</em>, allowing for the very similar models (in this toy example repeated models) to share their weight while more unique models keep their original weights. In our example we can see this difference clearly:</p>
<div class="sourceCode" id="cb19"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb19-1"><a href="#cb19-1" aria-hidden="true" tabindex="-1"></a><span class="fu">stacking_weights</span>(lpd_point[,<span class="fu">c</span>(<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">1</span>,<span class="dv">2</span>,<span class="dv">3</span>,<span class="dv">4</span>)])</span></code></pre></div>
<pre><code>Method: stacking
------
        weight
model1  0.001 
model2  0.001 
model3  0.001 
model4  0.001 
model5  0.001 
model6  0.001 
model7  0.001 
model8  0.001 
model9  0.001 
model10 0.001 
model11 0.000 
model12 0.000 
model13 0.987 </code></pre>
<p>Using stacking, the weight for the best model stays essentially unchanged.</p>
</div>
<div id="example-oceanic-tool-complexity" class="section level1">
<h1>Example: Oceanic tool complexity</h1>
<p>Another example we consider is the Kline oceanic tool complexity data, which McElreath describes as follows:</p>
<blockquote>
<p>Different historical island populations possessed tool kits of different size. These kits include fish hooks, axes, boats, hand plows, and many other types of tools. A number of theories predict that larger populations will both develop and sustain more complex tool kits. … It’s also suggested that contact rates among populations effectively increases population [sic, probably should be tool kit] size, as it’s relevant to technological evolution.</p>
</blockquote>
<p>We build models predicting the total number of tools given the log population size and the contact rate (high vs. low).</p>
<div class="sourceCode" id="cb21"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb21-1"><a href="#cb21-1" aria-hidden="true" tabindex="-1"></a><span class="fu">data</span>(Kline)</span>
<span id="cb21-2"><a href="#cb21-2" aria-hidden="true" tabindex="-1"></a>d <span class="ot">&lt;-</span> Kline</span>
<span id="cb21-3"><a href="#cb21-3" aria-hidden="true" tabindex="-1"></a>d<span class="sc">$</span>log_pop <span class="ot">&lt;-</span> <span class="fu">log</span>(d<span class="sc">$</span>population)</span>
<span id="cb21-4"><a href="#cb21-4" aria-hidden="true" tabindex="-1"></a>d<span class="sc">$</span>contact_high <span class="ot">&lt;-</span> <span class="fu">ifelse</span>(d<span class="sc">$</span>contact<span class="sc">==</span><span class="st">&quot;high&quot;</span>, <span class="dv">1</span>, <span class="dv">0</span>)</span>
<span id="cb21-5"><a href="#cb21-5" aria-hidden="true" tabindex="-1"></a><span class="fu">str</span>(d)</span></code></pre></div>
<pre><code>&#39;data.frame&#39;:   10 obs. of  7 variables:
 $ culture     : Factor w/ 10 levels &quot;Chuuk&quot;,&quot;Hawaii&quot;,..: 4 7 6 10 3 9 1 5 8 2
 $ population  : int  1100 1500 3600 4791 7400 8000 9200 13000 17500 275000
 $ contact     : Factor w/ 2 levels &quot;high&quot;,&quot;low&quot;: 2 2 2 1 1 1 1 2 1 2
 $ total_tools : int  13 22 24 43 33 19 40 28 55 71
 $ mean_TU     : num  3.2 4.7 4 5 5 4 3.8 6.6 5.4 6.6
 $ log_pop     : num  7 7.31 8.19 8.47 8.91 ...
 $ contact_high: num  0 0 0 1 1 1 1 0 1 0</code></pre>
<p>We start with a Poisson regression model with the log population size, the contact rate, and an interaction term between them (priors are informative priors as in <em>Statistical Rethinking</em>).</p>
<div class="sourceCode" id="cb23"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb23-1"><a href="#cb23-1" aria-hidden="true" tabindex="-1"></a>fit10 <span class="ot">&lt;-</span></span>
<span id="cb23-2"><a href="#cb23-2" aria-hidden="true" tabindex="-1"></a>  <span class="fu">stan_glm</span>(</span>
<span id="cb23-3"><a href="#cb23-3" aria-hidden="true" tabindex="-1"></a>    total_tools <span class="sc">~</span> log_pop <span class="sc">+</span> contact_high <span class="sc">+</span> log_pop <span class="sc">*</span> contact_high,</span>
<span id="cb23-4"><a href="#cb23-4" aria-hidden="true" tabindex="-1"></a>    <span class="at">family =</span> <span class="fu">poisson</span>(<span class="at">link =</span> <span class="st">&quot;log&quot;</span>),</span>
<span id="cb23-5"><a href="#cb23-5" aria-hidden="true" tabindex="-1"></a>    <span class="at">data =</span> d,</span>
<span id="cb23-6"><a href="#cb23-6" aria-hidden="true" tabindex="-1"></a>    <span class="at">prior =</span> <span class="fu">normal</span>(<span class="dv">0</span>, <span class="dv">1</span>, <span class="at">autoscale =</span> <span class="cn">FALSE</span>),</span>
<span id="cb23-7"><a href="#cb23-7" aria-hidden="true" tabindex="-1"></a>    <span class="at">prior_intercept =</span> <span class="fu">normal</span>(<span class="dv">0</span>, <span class="dv">100</span>, <span class="at">autoscale =</span> <span class="cn">FALSE</span>),</span>
<span id="cb23-8"><a href="#cb23-8" aria-hidden="true" tabindex="-1"></a>    <span class="at">seed =</span> <span class="dv">2030</span></span>
<span id="cb23-9"><a href="#cb23-9" aria-hidden="true" tabindex="-1"></a>  )</span></code></pre></div>
<p>Before running other models, we check whether Poisson is good choice as the conditional observation model.</p>
<div class="sourceCode" id="cb24"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb24-1"><a href="#cb24-1" aria-hidden="true" tabindex="-1"></a>loo10 <span class="ot">&lt;-</span> <span class="fu">loo</span>(fit10)</span>
<span id="cb24-2"><a href="#cb24-2" aria-hidden="true" tabindex="-1"></a><span class="fu">print</span>(loo10)</span></code></pre></div>
<pre><code>
Computed from 4000 by 10 log-likelihood matrix.

         Estimate   SE
elpd_loo    -40.2  5.9
p_loo         5.0  1.7
looic        80.5 11.9
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 0.7]).

All Pareto k estimates are good (k &lt; 0.7).
See help(&#39;pareto-k-diagnostic&#39;) for details.</code></pre>
<p>We get at least one observation with <span class="math inline">\(k&gt;0.7\)</span> and the estimated effective number of parameters <code>p_loo</code> is larger than the total number of parameters in the model. This indicates that Poisson might be too narrow. A negative binomial model might be better, but with so few observations it is not so clear.</p>
<p>We can compute LOO more accurately by running Stan again for the leave-one-out folds with high <span class="math inline">\(k\)</span> estimates. When using <strong>rstanarm</strong> this can be done by specifying the <code>k_threshold</code> argument:</p>
<div class="sourceCode" id="cb26"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb26-1"><a href="#cb26-1" aria-hidden="true" tabindex="-1"></a>loo10 <span class="ot">&lt;-</span> <span class="fu">loo</span>(fit10, <span class="at">k_threshold=</span><span class="fl">0.7</span>)</span></code></pre></div>
<pre><code>All pareto_k estimates below user-specified threshold of 0.7. 
Returning loo object.</code></pre>
<div class="sourceCode" id="cb28"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb28-1"><a href="#cb28-1" aria-hidden="true" tabindex="-1"></a><span class="fu">print</span>(loo10)</span></code></pre></div>
<pre><code>
Computed from 4000 by 10 log-likelihood matrix.

         Estimate   SE
elpd_loo    -40.2  5.9
p_loo         5.0  1.7
looic        80.5 11.9
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 0.7]).

All Pareto k estimates are good (k &lt; 0.7).
See help(&#39;pareto-k-diagnostic&#39;) for details.</code></pre>
<p>In this case we see that there is not much difference, and thus it is relatively safe to continue.</p>
<p>As a comparison we also compute WAIC:</p>
<div class="sourceCode" id="cb30"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb30-1"><a href="#cb30-1" aria-hidden="true" tabindex="-1"></a>waic10 <span class="ot">&lt;-</span> <span class="fu">waic</span>(fit10)</span></code></pre></div>
<pre><code>Warning: 
3 (30.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.</code></pre>
<div class="sourceCode" id="cb32"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb32-1"><a href="#cb32-1" aria-hidden="true" tabindex="-1"></a><span class="fu">print</span>(waic10)</span></code></pre></div>
<pre><code>
Computed from 4000 by 10 log-likelihood matrix.

          Estimate   SE
elpd_waic    -39.9  5.9
p_waic         4.7  1.7
waic          79.8 11.8

3 (30.0%) p_waic estimates greater than 0.4. We recommend trying loo instead. </code></pre>
<p>The WAIC computation is giving warnings and the estimated ELPD is slightly more optimistic. We recommend using the PSIS-LOO results instead.</p>
<p>To assess whether the contact rate and interaction term are useful, we can make a comparison to models without these terms.</p>
<div class="sourceCode" id="cb34"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb34-1"><a href="#cb34-1" aria-hidden="true" tabindex="-1"></a>fit11 <span class="ot">&lt;-</span> <span class="fu">update</span>(fit10, <span class="at">formula =</span> total_tools <span class="sc">~</span> log_pop <span class="sc">+</span> contact_high)</span>
<span id="cb34-2"><a href="#cb34-2" aria-hidden="true" tabindex="-1"></a>fit12 <span class="ot">&lt;-</span> <span class="fu">update</span>(fit10, <span class="at">formula =</span> total_tools <span class="sc">~</span> log_pop)</span></code></pre></div>
<div class="sourceCode" id="cb35"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb35-1"><a href="#cb35-1" aria-hidden="true" tabindex="-1"></a>(loo11 <span class="ot">&lt;-</span> <span class="fu">loo</span>(fit11))</span></code></pre></div>
<pre><code>
Computed from 4000 by 10 log-likelihood matrix.

         Estimate   SE
elpd_loo    -39.7  5.8
p_loo         4.4  1.6
looic        79.4 11.6
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.0]).

All Pareto k estimates are good (k &lt; 0.7).
See help(&#39;pareto-k-diagnostic&#39;) for details.</code></pre>
<div class="sourceCode" id="cb37"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb37-1"><a href="#cb37-1" aria-hidden="true" tabindex="-1"></a>(loo12 <span class="ot">&lt;-</span> <span class="fu">loo</span>(fit12))</span></code></pre></div>
<pre><code>Warning: Found 1 observation(s) with a pareto_k &gt; 0.7. We recommend calling &#39;loo&#39; again with argument &#39;k_threshold = 0.7&#39; in order to calculate the ELPD without the assumption that these observations are negligible. This will refit the model 1 times to compute the ELPDs for the problematic observations directly.</code></pre>
<pre><code>
Computed from 4000 by 10 log-likelihood matrix.

         Estimate  SE
elpd_loo    -42.5 4.7
p_loo         4.1 1.1
looic        85.0 9.4
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 0.6]).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     9     90.0%   653     
   (0.7, 1]   (bad)      1     10.0%   &lt;NA&gt;    
   (1, Inf)   (very bad) 0      0.0%   &lt;NA&gt;    
See help(&#39;pareto-k-diagnostic&#39;) for details.</code></pre>
<div class="sourceCode" id="cb40"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb40-1"><a href="#cb40-1" aria-hidden="true" tabindex="-1"></a>loo11 <span class="ot">&lt;-</span> <span class="fu">loo</span>(fit11, <span class="at">k_threshold=</span><span class="fl">0.7</span>)</span></code></pre></div>
<pre><code>All pareto_k estimates below user-specified threshold of 0.7. 
Returning loo object.</code></pre>
<div class="sourceCode" id="cb42"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb42-1"><a href="#cb42-1" aria-hidden="true" tabindex="-1"></a>loo12 <span class="ot">&lt;-</span> <span class="fu">loo</span>(fit12, <span class="at">k_threshold=</span><span class="fl">0.7</span>)</span></code></pre></div>
<pre><code>1 problematic observation(s) found.
Model will be refit 1 times.</code></pre>
<pre><code>
Fitting model 1 out of 1 (leaving out observation 10)</code></pre>
<div class="sourceCode" id="cb45"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb45-1"><a href="#cb45-1" aria-hidden="true" tabindex="-1"></a>lpd_point <span class="ot">&lt;-</span> <span class="fu">cbind</span>(</span>
<span id="cb45-2"><a href="#cb45-2" aria-hidden="true" tabindex="-1"></a>  loo10<span class="sc">$</span>pointwise[, <span class="st">&quot;elpd_loo&quot;</span>], </span>
<span id="cb45-3"><a href="#cb45-3" aria-hidden="true" tabindex="-1"></a>  loo11<span class="sc">$</span>pointwise[, <span class="st">&quot;elpd_loo&quot;</span>], </span>
<span id="cb45-4"><a href="#cb45-4" aria-hidden="true" tabindex="-1"></a>  loo12<span class="sc">$</span>pointwise[, <span class="st">&quot;elpd_loo&quot;</span>]</span>
<span id="cb45-5"><a href="#cb45-5" aria-hidden="true" tabindex="-1"></a>)</span></code></pre></div>
<p>For comparison we’ll also compute WAIC values for these additional models:</p>
<div class="sourceCode" id="cb46"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb46-1"><a href="#cb46-1" aria-hidden="true" tabindex="-1"></a>waic11 <span class="ot">&lt;-</span> <span class="fu">waic</span>(fit11)</span></code></pre></div>
<pre><code>Warning: 
3 (30.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.</code></pre>
<div class="sourceCode" id="cb48"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb48-1"><a href="#cb48-1" aria-hidden="true" tabindex="-1"></a>waic12 <span class="ot">&lt;-</span> <span class="fu">waic</span>(fit12)</span></code></pre></div>
<pre><code>Warning: 
5 (50.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.</code></pre>
<div class="sourceCode" id="cb50"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb50-1"><a href="#cb50-1" aria-hidden="true" tabindex="-1"></a>waics <span class="ot">&lt;-</span> <span class="fu">c</span>(</span>
<span id="cb50-2"><a href="#cb50-2" aria-hidden="true" tabindex="-1"></a>  waic10<span class="sc">$</span>estimates[<span class="st">&quot;elpd_waic&quot;</span>, <span class="dv">1</span>], </span>
<span id="cb50-3"><a href="#cb50-3" aria-hidden="true" tabindex="-1"></a>  waic11<span class="sc">$</span>estimates[<span class="st">&quot;elpd_waic&quot;</span>, <span class="dv">1</span>], </span>
<span id="cb50-4"><a href="#cb50-4" aria-hidden="true" tabindex="-1"></a>  waic12<span class="sc">$</span>estimates[<span class="st">&quot;elpd_waic&quot;</span>, <span class="dv">1</span>]</span>
<span id="cb50-5"><a href="#cb50-5" aria-hidden="true" tabindex="-1"></a>)</span></code></pre></div>
<p>The WAIC computation again gives warnings, and we recommend using PSIS-LOO instead.</p>
<p>Finally, we compute 1) WAIC weights, 2) Pseudo-BMA weights without Bayesian bootstrap, 3) Pseudo-BMA+ weights with Bayesian bootstrap, and 4) Bayesian stacking weights.</p>
<div class="sourceCode" id="cb51"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb51-1"><a href="#cb51-1" aria-hidden="true" tabindex="-1"></a>waic_wts <span class="ot">&lt;-</span> <span class="fu">exp</span>(waics) <span class="sc">/</span> <span class="fu">sum</span>(<span class="fu">exp</span>(waics))</span>
<span id="cb51-2"><a href="#cb51-2" aria-hidden="true" tabindex="-1"></a>pbma_wts <span class="ot">&lt;-</span> <span class="fu">pseudobma_weights</span>(lpd_point, <span class="at">BB=</span><span class="cn">FALSE</span>)</span>
<span id="cb51-3"><a href="#cb51-3" aria-hidden="true" tabindex="-1"></a>pbma_BB_wts <span class="ot">&lt;-</span> <span class="fu">pseudobma_weights</span>(lpd_point) <span class="co"># default is BB=TRUE</span></span>
<span id="cb51-4"><a href="#cb51-4" aria-hidden="true" tabindex="-1"></a>stacking_wts <span class="ot">&lt;-</span> <span class="fu">stacking_weights</span>(lpd_point)</span>
<span id="cb51-5"><a href="#cb51-5" aria-hidden="true" tabindex="-1"></a><span class="fu">round</span>(<span class="fu">cbind</span>(waic_wts, pbma_wts, pbma_BB_wts, stacking_wts), <span class="dv">2</span>)</span></code></pre></div>
<pre><code>       waic_wts pbma_wts pbma_BB_wts stacking_wts
model1     0.38     0.36        0.31          0.0
model2     0.58     0.63        0.53          0.8
model3     0.04     0.02        0.16          0.2</code></pre>
<p>All weights favor the second model with the log population and the contact rate. WAIC weights and Pseudo-BMA weights (without Bayesian bootstrap) are similar, while Pseudo-BMA+ is more cautious and closer to stacking weights.</p>
<p>It may seem surprising that Bayesian stacking is giving zero weight to the first model, but this is likely due to the fact that the estimated effect for the interaction term is close to zero and thus models 1 and 2 give very similar predictions. In other words, incorporating the model with the interaction (model 1) into the model average doesn’t improve the predictions at all and so model 1 is given a weight of 0. On the other hand, models 2 and 3 are giving slightly different predictions and thus their combination may be slightly better than either alone. This behavior is related to the repeated similar model illustration in the milk example above.</p>
</div>
<div id="simpler-coding-using-loo_model_weights-function" class="section level1">
<h1>Simpler coding using <code>loo_model_weights</code> function</h1>
<p>Although in the examples above we called the <code>stacking_weights</code> and <code>pseudobma_weights</code> functions directly, we can also use the <code>loo_model_weights</code> wrapper, which takes as its input either a list of pointwise log-likelihood matrices or a list of precomputed loo objects. There are also <code>loo_model_weights</code> methods for stanreg objects (fitted model objects from <strong>rstanarm</strong>) as well as fitted model objects from other packages (e.g. <strong>brms</strong>) that do the preparation work for the user (see, e.g., the examples at <code>help(&quot;loo_model_weights&quot;, package = &quot;rstanarm&quot;)</code>).</p>
<div class="sourceCode" id="cb53"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb53-1"><a href="#cb53-1" aria-hidden="true" tabindex="-1"></a><span class="co"># using list of loo objects</span></span>
<span id="cb53-2"><a href="#cb53-2" aria-hidden="true" tabindex="-1"></a>loo_list <span class="ot">&lt;-</span> <span class="fu">list</span>(loo10, loo11, loo12)</span>
<span id="cb53-3"><a href="#cb53-3" aria-hidden="true" tabindex="-1"></a><span class="fu">loo_model_weights</span>(loo_list)</span></code></pre></div>
<pre><code>Method: stacking
------
      weight
fit10 0.000 
fit11 0.802 
fit12 0.198 </code></pre>
<div class="sourceCode" id="cb55"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb55-1"><a href="#cb55-1" aria-hidden="true" tabindex="-1"></a><span class="fu">loo_model_weights</span>(loo_list, <span class="at">method =</span> <span class="st">&quot;pseudobma&quot;</span>)</span></code></pre></div>
<pre><code>Method: pseudo-BMA+ with Bayesian bootstrap
------
      weight
fit10 0.310 
fit11 0.539 
fit12 0.151 </code></pre>
<div class="sourceCode" id="cb57"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb57-1"><a href="#cb57-1" aria-hidden="true" tabindex="-1"></a><span class="fu">loo_model_weights</span>(loo_list, <span class="at">method =</span> <span class="st">&quot;pseudobma&quot;</span>, <span class="at">BB =</span> <span class="cn">FALSE</span>)</span></code></pre></div>
<pre><code>Method: pseudo-BMA
------
      weight
fit10 0.356 
fit11 0.629 
fit12 0.015 </code></pre>
</div>
<div id="references" class="section level1">
<h1>References</h1>
<p>McElreath, R. (2016). <em>Statistical rethinking: A Bayesian course with examples in R and Stan</em>. Chapman &amp; Hall/CRC. <a href="http://xcelab.net/rm/statistical-rethinking/" class="uri">http://xcelab.net/rm/statistical-rethinking/</a></p>
<p>Piironen, J. and Vehtari, A. (2017a). Sparsity information and regularization in the horseshoe and other shrinkage priors. In Electronic Journal of Statistics, 11(2):5018-5051. <a href="https://projecteuclid.org/euclid.ejs/1513306866">Online</a>.</p>
<p>Piironen, J. and Vehtari, A. (2017b). Comparison of Bayesian predictive methods for model selection. Statistics and Computing, 27(3):711-735. :10.1007/s11222-016-9649-y. <a href="https://link.springer.com/article/10.1007/s11222-016-9649-y">Online</a>.</p>
<p>Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. <em>Statistics and Computing</em>. 27(5), 1413–1432. :10.1007/s11222-016-9696-4. <a href="https://link.springer.com/article/10.1007/s11222-016-9696-4">online</a>, <a href="https://arxiv.org/abs/1507.04544">arXiv preprint arXiv:1507.04544</a>.</p>
<p>Vehtari, A., Simpson, D., Gelman, A., Yao, Y., and Gabry, J. (2024). Pareto smoothed importance sampling. <em>Journal of Machine Learning Research</em>, 25(72):1-58. <a href="https://jmlr.org/papers/v25/19-556.html">PDF</a></p>
<p>Yao, Y., Vehtari, A., Simpson, D., and Gelman, A. (2018). Using stacking to average Bayesian predictive distributions. In Bayesian Analysis, :10.1214/17-BA1091. <a href="https://projecteuclid.org/euclid.ba/1516093227">Online</a>.</p>
</div>



<!-- code folding -->


<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
  (function () {
    var script = document.createElement("script");
    script.type = "text/javascript";
    script.src  = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
    document.getElementsByTagName("head")[0].appendChild(script);
  })();
</script>

</body>
</html>