File: loo2-elpd.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 (557 lines) | stat: -rw-r--r-- 34,030 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="Bruno Nicenboim" />

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

<title>Holdout validation and K-fold cross-validation of Stan programs with 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">Holdout validation and K-fold cross-validation of Stan programs with the loo package</h1>
<h4 class="author">Bruno Nicenboim</h4>
<h4 class="date">2025-12-22</h4>


<div id="TOC">
<ul>
<li><a href="#introduction">Introduction</a></li>
<li><a href="#example-eradication-of-roaches-using-holdout-validation-approach">Example: Eradication of Roaches using holdout validation approach</a>
<ul>
<li><a href="#coding-the-stan-model">Coding the Stan model</a></li>
<li><a href="#setup">Setup</a></li>
</ul></li>
<li><a href="#holdout-validation">Holdout validation</a>
<ul>
<li><a href="#splitting-the-data-between-train-and-test">Splitting the data between train and test</a></li>
<li><a href="#fitting-the-model-with-rstan">Fitting the model with RStan</a></li>
<li><a href="#computing-holdout-elpd">Computing holdout elpd:</a></li>
</ul></li>
<li><a href="#k-fold-cross-validation">K-fold cross validation</a>
<ul>
<li><a href="#splitting-the-data-in-folds">Splitting the data in folds</a></li>
<li><a href="#fitting-and-extracting-the-log-pointwise-predictive-densities-for-each-fold">Fitting and extracting the log pointwise predictive densities for each fold</a></li>
<li><a href="#computing-k-fold-elpd">Computing K-fold elpd:</a></li>
</ul></li>
<li><a href="#references">References</a></li>
</ul>
</div>

<!--
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{Holdout validation and K-fold cross-validation of Stan programs with the loo package}
-->
<div id="introduction" class="section level1">
<h1>Introduction</h1>
<p>This vignette demonstrates how to do holdout validation and K-fold cross-validation with <strong>loo</strong> for a Stan program.</p>
</div>
<div id="example-eradication-of-roaches-using-holdout-validation-approach" class="section level1">
<h1>Example: Eradication of Roaches using holdout validation approach</h1>
<p>This vignette uses the same example as in the vignettes <a href="http://mc-stan.org/loo/articles/loo2-example.html"><em>Using the loo package (version &gt;= 2.0.0)</em></a> and <a href="https://mc-stan.org/loo/articles/loo2-moment-matching.html"><em>Avoiding model refits in leave-one-out cross-validation with moment matching</em></a>.</p>
<div id="coding-the-stan-model" class="section level2">
<h2>Coding the Stan model</h2>
<p>Here is the Stan code for fitting a Poisson regression model:</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="co"># Note: some syntax used in this Stan program requires RStan &gt;= 2.26 (or CmdStanR)</span></span>
<span id="cb1-2"><a href="#cb1-2" aria-hidden="true" tabindex="-1"></a><span class="co"># To use an older version of RStan change the line declaring `y` to: int y[N];</span></span>
<span id="cb1-3"><a href="#cb1-3" aria-hidden="true" tabindex="-1"></a>stancode <span class="ot">&lt;-</span> <span class="st">&quot;</span></span>
<span id="cb1-4"><a href="#cb1-4" aria-hidden="true" tabindex="-1"></a><span class="st">data {</span></span>
<span id="cb1-5"><a href="#cb1-5" aria-hidden="true" tabindex="-1"></a><span class="st">  int&lt;lower=1&gt; K;</span></span>
<span id="cb1-6"><a href="#cb1-6" aria-hidden="true" tabindex="-1"></a><span class="st">  int&lt;lower=1&gt; N;</span></span>
<span id="cb1-7"><a href="#cb1-7" aria-hidden="true" tabindex="-1"></a><span class="st">  matrix[N,K] x;</span></span>
<span id="cb1-8"><a href="#cb1-8" aria-hidden="true" tabindex="-1"></a><span class="st">  array[N] int y;</span></span>
<span id="cb1-9"><a href="#cb1-9" aria-hidden="true" tabindex="-1"></a><span class="st">  vector[N] offset_; // offset is reserved keyword in Stan so use offset_</span></span>
<span id="cb1-10"><a href="#cb1-10" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb1-11"><a href="#cb1-11" aria-hidden="true" tabindex="-1"></a><span class="st">  real beta_prior_scale;</span></span>
<span id="cb1-12"><a href="#cb1-12" aria-hidden="true" tabindex="-1"></a><span class="st">  real alpha_prior_scale;</span></span>
<span id="cb1-13"><a href="#cb1-13" aria-hidden="true" tabindex="-1"></a><span class="st">}</span></span>
<span id="cb1-14"><a href="#cb1-14" aria-hidden="true" tabindex="-1"></a><span class="st">parameters {</span></span>
<span id="cb1-15"><a href="#cb1-15" aria-hidden="true" tabindex="-1"></a><span class="st">  vector[K] beta;</span></span>
<span id="cb1-16"><a href="#cb1-16" aria-hidden="true" tabindex="-1"></a><span class="st">  real intercept;</span></span>
<span id="cb1-17"><a href="#cb1-17" aria-hidden="true" tabindex="-1"></a><span class="st">}</span></span>
<span id="cb1-18"><a href="#cb1-18" aria-hidden="true" tabindex="-1"></a><span class="st">model {</span></span>
<span id="cb1-19"><a href="#cb1-19" aria-hidden="true" tabindex="-1"></a><span class="st">  y ~ poisson(exp(x * beta + intercept + offset_));</span></span>
<span id="cb1-20"><a href="#cb1-20" aria-hidden="true" tabindex="-1"></a><span class="st">  beta ~ normal(0,beta_prior_scale);</span></span>
<span id="cb1-21"><a href="#cb1-21" aria-hidden="true" tabindex="-1"></a><span class="st">  intercept ~ normal(0,alpha_prior_scale);</span></span>
<span id="cb1-22"><a href="#cb1-22" aria-hidden="true" tabindex="-1"></a><span class="st">}</span></span>
<span id="cb1-23"><a href="#cb1-23" aria-hidden="true" tabindex="-1"></a><span class="st">generated quantities {</span></span>
<span id="cb1-24"><a href="#cb1-24" aria-hidden="true" tabindex="-1"></a><span class="st">  vector[N] log_lik;</span></span>
<span id="cb1-25"><a href="#cb1-25" aria-hidden="true" tabindex="-1"></a><span class="st">  for (n in 1:N)</span></span>
<span id="cb1-26"><a href="#cb1-26" aria-hidden="true" tabindex="-1"></a><span class="st">    log_lik[n] = poisson_lpmf(y[n] | exp(x[n] * beta + intercept + offset_[n]));</span></span>
<span id="cb1-27"><a href="#cb1-27" aria-hidden="true" tabindex="-1"></a><span class="st">}</span></span>
<span id="cb1-28"><a href="#cb1-28" aria-hidden="true" tabindex="-1"></a><span class="st">&quot;</span></span></code></pre></div>
<p>Following the usual approach recommended in <a href="http://mc-stan.org/loo/articles/loo2-with-rstan.html"><em>Writing Stan programs for use with the loo package</em></a>, we compute the log-likelihood for each observation in the <code>generated quantities</code> block of the Stan program.</p>
</div>
<div id="setup" class="section level2">
<h2>Setup</h2>
<p>In addition to <strong>loo</strong>, we load the <strong>rstan</strong> package for fitting the model. We will also need the <strong>rstanarm</strong> package for the data.</p>
<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">library</span>(<span class="st">&quot;rstan&quot;</span>)</span></code></pre></div>
<pre><code>Warning: package &#39;StanHeaders&#39; was built under R version 4.4.3</code></pre>
<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><span class="fu">library</span>(<span class="st">&quot;loo&quot;</span>)</span>
<span id="cb4-2"><a href="#cb4-2" aria-hidden="true" tabindex="-1"></a>seed <span class="ot">&lt;-</span> <span class="dv">9547</span></span>
<span id="cb4-3"><a href="#cb4-3" aria-hidden="true" tabindex="-1"></a><span class="fu">set.seed</span>(seed)</span></code></pre></div>
</div>
</div>
<div id="holdout-validation" class="section level1">
<h1>Holdout validation</h1>
<p>For this approach, the model is first fit to the “train” data and then is evaluated on the held-out “test” data.</p>
<div id="splitting-the-data-between-train-and-test" class="section level2">
<h2>Splitting the data between train and test</h2>
<p>The data is divided between train (80% of the data) and test (20%):</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><span class="co"># Prepare data</span></span>
<span id="cb5-2"><a href="#cb5-2" aria-hidden="true" tabindex="-1"></a><span class="fu">data</span>(roaches, <span class="at">package =</span> <span class="st">&quot;rstanarm&quot;</span>)</span>
<span id="cb5-3"><a href="#cb5-3" aria-hidden="true" tabindex="-1"></a>roaches<span class="sc">$</span>roach1 <span class="ot">&lt;-</span> <span class="fu">sqrt</span>(roaches<span class="sc">$</span>roach1)</span>
<span id="cb5-4"><a href="#cb5-4" aria-hidden="true" tabindex="-1"></a>roaches<span class="sc">$</span>offset <span class="ot">&lt;-</span> <span class="fu">log</span>(roaches[,<span class="st">&quot;exposure2&quot;</span>])</span>
<span id="cb5-5"><a href="#cb5-5" aria-hidden="true" tabindex="-1"></a><span class="co"># 20% of the data goes to the test set:</span></span>
<span id="cb5-6"><a href="#cb5-6" aria-hidden="true" tabindex="-1"></a>roaches<span class="sc">$</span>test <span class="ot">&lt;-</span> <span class="dv">0</span></span>
<span id="cb5-7"><a href="#cb5-7" aria-hidden="true" tabindex="-1"></a>roaches<span class="sc">$</span>test[<span class="fu">sample</span>(.<span class="dv">2</span> <span class="sc">*</span> <span class="fu">seq_len</span>(<span class="fu">nrow</span>(roaches)))] <span class="ot">&lt;-</span> <span class="dv">1</span></span>
<span id="cb5-8"><a href="#cb5-8" aria-hidden="true" tabindex="-1"></a><span class="co"># data to &quot;train&quot; the model</span></span>
<span id="cb5-9"><a href="#cb5-9" aria-hidden="true" tabindex="-1"></a>data_train <span class="ot">&lt;-</span> <span class="fu">list</span>(<span class="at">y =</span> roaches<span class="sc">$</span>y[roaches<span class="sc">$</span>test <span class="sc">==</span> <span class="dv">0</span>],</span>
<span id="cb5-10"><a href="#cb5-10" aria-hidden="true" tabindex="-1"></a>                   <span class="at">x =</span> <span class="fu">as.matrix</span>(roaches[roaches<span class="sc">$</span>test <span class="sc">==</span> <span class="dv">0</span>,</span>
<span id="cb5-11"><a href="#cb5-11" aria-hidden="true" tabindex="-1"></a>                                         <span class="fu">c</span>(<span class="st">&quot;roach1&quot;</span>, <span class="st">&quot;treatment&quot;</span>, <span class="st">&quot;senior&quot;</span>)]),</span>
<span id="cb5-12"><a href="#cb5-12" aria-hidden="true" tabindex="-1"></a>                   <span class="at">N =</span> <span class="fu">nrow</span>(roaches[roaches<span class="sc">$</span>test <span class="sc">==</span> <span class="dv">0</span>,]),</span>
<span id="cb5-13"><a href="#cb5-13" aria-hidden="true" tabindex="-1"></a>                   <span class="at">K =</span> <span class="dv">3</span>,</span>
<span id="cb5-14"><a href="#cb5-14" aria-hidden="true" tabindex="-1"></a>                   <span class="at">offset_ =</span> roaches<span class="sc">$</span>offset[roaches<span class="sc">$</span>test <span class="sc">==</span> <span class="dv">0</span>],</span>
<span id="cb5-15"><a href="#cb5-15" aria-hidden="true" tabindex="-1"></a>                   <span class="at">beta_prior_scale =</span> <span class="fl">2.5</span>,</span>
<span id="cb5-16"><a href="#cb5-16" aria-hidden="true" tabindex="-1"></a>                   <span class="at">alpha_prior_scale =</span> <span class="fl">5.0</span></span>
<span id="cb5-17"><a href="#cb5-17" aria-hidden="true" tabindex="-1"></a>                   )</span>
<span id="cb5-18"><a href="#cb5-18" aria-hidden="true" tabindex="-1"></a><span class="co"># data to &quot;test&quot; the model</span></span>
<span id="cb5-19"><a href="#cb5-19" aria-hidden="true" tabindex="-1"></a>data_test <span class="ot">&lt;-</span> <span class="fu">list</span>(<span class="at">y =</span> roaches<span class="sc">$</span>y[roaches<span class="sc">$</span>test <span class="sc">==</span> <span class="dv">1</span>],</span>
<span id="cb5-20"><a href="#cb5-20" aria-hidden="true" tabindex="-1"></a>                   <span class="at">x =</span> <span class="fu">as.matrix</span>(roaches[roaches<span class="sc">$</span>test <span class="sc">==</span> <span class="dv">1</span>,</span>
<span id="cb5-21"><a href="#cb5-21" aria-hidden="true" tabindex="-1"></a>                                         <span class="fu">c</span>(<span class="st">&quot;roach1&quot;</span>, <span class="st">&quot;treatment&quot;</span>, <span class="st">&quot;senior&quot;</span>)]),</span>
<span id="cb5-22"><a href="#cb5-22" aria-hidden="true" tabindex="-1"></a>                   <span class="at">N =</span> <span class="fu">nrow</span>(roaches[roaches<span class="sc">$</span>test <span class="sc">==</span> <span class="dv">1</span>,]),</span>
<span id="cb5-23"><a href="#cb5-23" aria-hidden="true" tabindex="-1"></a>                   <span class="at">K =</span> <span class="dv">3</span>,</span>
<span id="cb5-24"><a href="#cb5-24" aria-hidden="true" tabindex="-1"></a>                   <span class="at">offset_ =</span> roaches<span class="sc">$</span>offset[roaches<span class="sc">$</span>test <span class="sc">==</span> <span class="dv">1</span>],</span>
<span id="cb5-25"><a href="#cb5-25" aria-hidden="true" tabindex="-1"></a>                   <span class="at">beta_prior_scale =</span> <span class="fl">2.5</span>,</span>
<span id="cb5-26"><a href="#cb5-26" aria-hidden="true" tabindex="-1"></a>                   <span class="at">alpha_prior_scale =</span> <span class="fl">5.0</span></span>
<span id="cb5-27"><a href="#cb5-27" aria-hidden="true" tabindex="-1"></a>                   )</span></code></pre></div>
</div>
<div id="fitting-the-model-with-rstan" class="section level2">
<h2>Fitting the model with RStan</h2>
<p>Next we fit the model to the “test” data in Stan using the <strong>rstan</strong> package:</p>
<div class="sourceCode" id="cb6"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb6-1"><a href="#cb6-1" aria-hidden="true" tabindex="-1"></a><span class="co"># Compile</span></span>
<span id="cb6-2"><a href="#cb6-2" aria-hidden="true" tabindex="-1"></a>stanmodel <span class="ot">&lt;-</span> <span class="fu">stan_model</span>(<span class="at">model_code =</span> stancode)</span>
<span id="cb6-3"><a href="#cb6-3" aria-hidden="true" tabindex="-1"></a><span class="co"># Fit model</span></span>
<span id="cb6-4"><a href="#cb6-4" aria-hidden="true" tabindex="-1"></a>fit <span class="ot">&lt;-</span> <span class="fu">sampling</span>(stanmodel, <span class="at">data =</span> data_train, <span class="at">seed =</span> seed, <span class="at">refresh =</span> <span class="dv">0</span>)</span></code></pre></div>
<p>We recompute the generated quantities using the posterior draws conditional on the training data, but we now pass in the held-out data to get the log predictive densities for the test data. Because we are using independent data, the log predictive density coincides with the log likelihood of the test data.</p>
<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>gen_test <span class="ot">&lt;-</span> <span class="fu">gqs</span>(stanmodel, <span class="at">draws =</span> <span class="fu">as.matrix</span>(fit), <span class="at">data=</span> data_test)</span>
<span id="cb7-2"><a href="#cb7-2" aria-hidden="true" tabindex="-1"></a>log_pd <span class="ot">&lt;-</span> <span class="fu">extract_log_lik</span>(gen_test)</span></code></pre></div>
</div>
<div id="computing-holdout-elpd" class="section level2">
<h2>Computing holdout elpd:</h2>
<p>Now we evaluate the predictive performance of the model on the test data using <code>elpd()</code>.</p>
<div class="sourceCode" id="cb8"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb8-1"><a href="#cb8-1" aria-hidden="true" tabindex="-1"></a>(elpd_holdout <span class="ot">&lt;-</span> <span class="fu">elpd</span>(log_pd))</span></code></pre></div>
<pre><code>
Computed from 4000 by 52 log-likelihood matrix using the generic elpd function

     Estimate    SE
elpd  -1734.5 288.3
ic     3468.9 576.6</code></pre>
<p>When one wants to compare different models, the function <code>loo_compare()</code> can be used to assess the difference in performance.</p>
</div>
</div>
<div id="k-fold-cross-validation" class="section level1">
<h1>K-fold cross validation</h1>
<p>For this approach the data is divided into folds, and each time one fold is tested while the rest of the data is used to fit the model (see Vehtari et al., 2017).</p>
<div id="splitting-the-data-in-folds" class="section level2">
<h2>Splitting the data in folds</h2>
<p>We use the data that is already pre-processed and we divide it in 10 random folds using <code>kfold_split_random</code></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"># Prepare data</span></span>
<span id="cb10-2"><a href="#cb10-2" aria-hidden="true" tabindex="-1"></a>roaches<span class="sc">$</span>fold <span class="ot">&lt;-</span> <span class="fu">kfold_split_random</span>(<span class="at">K =</span> <span class="dv">10</span>, <span class="at">N =</span> <span class="fu">nrow</span>(roaches))</span></code></pre></div>
</div>
<div id="fitting-and-extracting-the-log-pointwise-predictive-densities-for-each-fold" class="section level2">
<h2>Fitting and extracting the log pointwise predictive densities for each fold</h2>
<p>We now loop over the 10 folds. In each fold we do the following. First, we fit the model to all the observations except the ones belonging to the left-out fold. Second, we compute the log pointwise predictive densities for the left-out fold. Last, we store the predictive density for the observations of the left-out fold in a matrix. The output of this loop is a matrix of the log pointwise predictive densities of all the observations.</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="co"># Prepare a matrix with the number of post-warmup iterations by number of observations:</span></span>
<span id="cb11-2"><a href="#cb11-2" aria-hidden="true" tabindex="-1"></a>log_pd_kfold <span class="ot">&lt;-</span> <span class="fu">matrix</span>(<span class="at">nrow =</span> <span class="dv">4000</span>, <span class="at">ncol =</span> <span class="fu">nrow</span>(roaches))</span>
<span id="cb11-3"><a href="#cb11-3" aria-hidden="true" tabindex="-1"></a><span class="co"># Loop over the folds</span></span>
<span id="cb11-4"><a href="#cb11-4" aria-hidden="true" tabindex="-1"></a><span class="cf">for</span>(k <span class="cf">in</span> <span class="dv">1</span><span class="sc">:</span><span class="dv">10</span>){</span>
<span id="cb11-5"><a href="#cb11-5" aria-hidden="true" tabindex="-1"></a>  data_train <span class="ot">&lt;-</span> <span class="fu">list</span>(<span class="at">y =</span> roaches<span class="sc">$</span>y[roaches<span class="sc">$</span>fold <span class="sc">!=</span> k],</span>
<span id="cb11-6"><a href="#cb11-6" aria-hidden="true" tabindex="-1"></a>                   <span class="at">x =</span> <span class="fu">as.matrix</span>(roaches[roaches<span class="sc">$</span>fold <span class="sc">!=</span> k,</span>
<span id="cb11-7"><a href="#cb11-7" aria-hidden="true" tabindex="-1"></a>                                         <span class="fu">c</span>(<span class="st">&quot;roach1&quot;</span>, <span class="st">&quot;treatment&quot;</span>, <span class="st">&quot;senior&quot;</span>)]),</span>
<span id="cb11-8"><a href="#cb11-8" aria-hidden="true" tabindex="-1"></a>                   <span class="at">N =</span> <span class="fu">nrow</span>(roaches[roaches<span class="sc">$</span>fold <span class="sc">!=</span> k,]),</span>
<span id="cb11-9"><a href="#cb11-9" aria-hidden="true" tabindex="-1"></a>                   <span class="at">K =</span> <span class="dv">3</span>,</span>
<span id="cb11-10"><a href="#cb11-10" aria-hidden="true" tabindex="-1"></a>                   <span class="at">offset_ =</span> roaches<span class="sc">$</span>offset[roaches<span class="sc">$</span>fold <span class="sc">!=</span> k],</span>
<span id="cb11-11"><a href="#cb11-11" aria-hidden="true" tabindex="-1"></a>                   <span class="at">beta_prior_scale =</span> <span class="fl">2.5</span>,</span>
<span id="cb11-12"><a href="#cb11-12" aria-hidden="true" tabindex="-1"></a>                   <span class="at">alpha_prior_scale =</span> <span class="fl">5.0</span></span>
<span id="cb11-13"><a href="#cb11-13" aria-hidden="true" tabindex="-1"></a>                   )</span>
<span id="cb11-14"><a href="#cb11-14" aria-hidden="true" tabindex="-1"></a>  data_test <span class="ot">&lt;-</span> <span class="fu">list</span>(<span class="at">y =</span> roaches<span class="sc">$</span>y[roaches<span class="sc">$</span>fold <span class="sc">==</span> k],</span>
<span id="cb11-15"><a href="#cb11-15" aria-hidden="true" tabindex="-1"></a>                   <span class="at">x =</span> <span class="fu">as.matrix</span>(roaches[roaches<span class="sc">$</span>fold <span class="sc">==</span> k,</span>
<span id="cb11-16"><a href="#cb11-16" aria-hidden="true" tabindex="-1"></a>                                         <span class="fu">c</span>(<span class="st">&quot;roach1&quot;</span>, <span class="st">&quot;treatment&quot;</span>, <span class="st">&quot;senior&quot;</span>)]),</span>
<span id="cb11-17"><a href="#cb11-17" aria-hidden="true" tabindex="-1"></a>                   <span class="at">N =</span> <span class="fu">nrow</span>(roaches[roaches<span class="sc">$</span>fold <span class="sc">==</span> k,]),</span>
<span id="cb11-18"><a href="#cb11-18" aria-hidden="true" tabindex="-1"></a>                   <span class="at">K =</span> <span class="dv">3</span>,</span>
<span id="cb11-19"><a href="#cb11-19" aria-hidden="true" tabindex="-1"></a>                   <span class="at">offset_ =</span> roaches<span class="sc">$</span>offset[roaches<span class="sc">$</span>fold <span class="sc">==</span> k],</span>
<span id="cb11-20"><a href="#cb11-20" aria-hidden="true" tabindex="-1"></a>                   <span class="at">beta_prior_scale =</span> <span class="fl">2.5</span>,</span>
<span id="cb11-21"><a href="#cb11-21" aria-hidden="true" tabindex="-1"></a>                   <span class="at">alpha_prior_scale =</span> <span class="fl">5.0</span></span>
<span id="cb11-22"><a href="#cb11-22" aria-hidden="true" tabindex="-1"></a>                   )</span>
<span id="cb11-23"><a href="#cb11-23" aria-hidden="true" tabindex="-1"></a>  fit <span class="ot">&lt;-</span> <span class="fu">sampling</span>(stanmodel, <span class="at">data =</span> data_train, <span class="at">seed =</span> seed, <span class="at">refresh =</span> <span class="dv">0</span>)</span>
<span id="cb11-24"><a href="#cb11-24" aria-hidden="true" tabindex="-1"></a>  gen_test <span class="ot">&lt;-</span> <span class="fu">gqs</span>(stanmodel, <span class="at">draws =</span> <span class="fu">as.matrix</span>(fit), <span class="at">data=</span> data_test)</span>
<span id="cb11-25"><a href="#cb11-25" aria-hidden="true" tabindex="-1"></a>  log_pd_kfold[, roaches<span class="sc">$</span>fold <span class="sc">==</span> k] <span class="ot">&lt;-</span> <span class="fu">extract_log_lik</span>(gen_test)</span>
<span id="cb11-26"><a href="#cb11-26" aria-hidden="true" tabindex="-1"></a>}</span></code></pre></div>
</div>
<div id="computing-k-fold-elpd" class="section level2">
<h2>Computing K-fold elpd:</h2>
<p>Now we evaluate the predictive performance of the model on the 10 folds using <code>elpd()</code>.</p>
<div class="sourceCode" id="cb12"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb12-1"><a href="#cb12-1" aria-hidden="true" tabindex="-1"></a>(elpd_kfold <span class="ot">&lt;-</span> <span class="fu">elpd</span>(log_pd_kfold))</span></code></pre></div>
<pre><code>
Computed from 4000 by 262 log-likelihood matrix using the generic elpd function

     Estimate     SE
elpd  -5558.1  729.0
ic    11116.2 1457.9</code></pre>
<p>If one wants to compare several models (with <code>loo_compare</code>), one should use the same folds for all the different models.</p>
</div>
</div>
<div id="references" class="section level1">
<h1>References</h1>
<p>Gelman, A., and Hill, J. (2007). <em>Data Analysis Using Regression and Multilevel Hierarchical Models.</em> Cambridge University Press.</p>
<p>Stan Development Team (2020) <em>RStan: the R interface to Stan, Version 2.21.1</em> <a href="https://mc-stan.org" class="uri">https://mc-stan.org</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. Links: <a href="https://link.springer.com/article/10.1007/s11222-016-9696-4">published</a> | <a href="https://arxiv.org/abs/1507.04544">arXiv preprint</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>