<!DOCTYPE html>
<!-- Generated by pkgdown: do not edit by hand --><html lang="en">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Extracting and visualizing tidy residuals from Bayesian models ‚Ä¢ tidybayes</title>
<!-- favicons --><link rel="icon" type="image/png" sizes="16x16" href="../favicon-16x16.png">
<link rel="icon" type="image/png" sizes="32x32" href="../favicon-32x32.png">
<link rel="apple-touch-icon" type="image/png" sizes="180x180" href="../apple-touch-icon.png">
<link rel="apple-touch-icon" type="image/png" sizes="120x120" href="../apple-touch-icon-120x120.png">
<link rel="apple-touch-icon" type="image/png" sizes="76x76" href="../apple-touch-icon-76x76.png">
<link rel="apple-touch-icon" type="image/png" sizes="60x60" href="../apple-touch-icon-60x60.png">
<!-- jquery --><script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/3.7.1/jquery.min.js" integrity="sha512-v2CJ7UaYy4JwqLDIrZUI/4hqeoQieOmAZNXBeQyjo21dadnwR+8ZaIJVT8EE2iyI61OV8e6M8PP2/4hpQINQ/g==" crossorigin="anonymous" referrerpolicy="no-referrer"></script><!-- Bootstrap --><link href="https://cdnjs.cloudflare.com/ajax/libs/bootswatch/3.4.0/paper/bootstrap.min.css" rel="stylesheet" crossorigin="anonymous">
<script src="https://cdnjs.cloudflare.com/ajax/libs/twitter-bootstrap/3.4.1/js/bootstrap.min.js" integrity="sha256-nuL8/2cJ5NDSSwnKD8VqreErSWHtnEP9E7AySL+1ev4=" crossorigin="anonymous"></script><!-- bootstrap-toc --><link rel="stylesheet" href="../bootstrap-toc.css">
<script src="../bootstrap-toc.js"></script><!-- Font Awesome icons --><link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/all.min.css" integrity="sha256-mmgLkCYLUQbXn0B1SRqzHar6dCnv9oZFPEC1g1cwlkk=" crossorigin="anonymous">
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/v4-shims.min.css" integrity="sha256-wZjR52fzng1pJHwx4aV2AO3yyTOXrcDW7jBpJtTwVxw=" crossorigin="anonymous">
<!-- clipboard.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/2.0.6/clipboard.min.js" integrity="sha256-inc5kl9MA1hkeYUt+EC3BhlIgyp/2jDIyBLS6k3UxPI=" crossorigin="anonymous"></script><!-- headroom.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/headroom.min.js" integrity="sha256-AsUX4SJE1+yuDu5+mAVzJbuYNPHj/WroHuZ8Ir/CkE0=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/jQuery.headroom.min.js" integrity="sha256-ZX/yNShbjqsohH1k95liqY9Gd8uOiE1S4vZc+9KQ1K4=" crossorigin="anonymous"></script><!-- pkgdown --><link href="../pkgdown.css" rel="stylesheet">
<script src="../pkgdown.js"></script><link href="../extra.css" rel="stylesheet">
<meta property="og:title" content="Extracting and visualizing tidy residuals from Bayesian models">
<!-- mathjax --><script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js" integrity="sha256-nvJJv9wWKEm88qvoQl9ekL2J+k/RWIsaSScxxlsrv8k=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/config/TeX-AMS-MML_HTMLorMML.js" integrity="sha256-84DKXVJXs0/F8OTMzX4UR909+jtl4G7SPypPavF+GfA=" crossorigin="anonymous"></script><!--[if lt IE 9]>
<script src="https://oss.maxcdn.com/html5shiv/3.7.3/html5shiv.min.js"></script>
<script src="https://oss.maxcdn.com/respond/1.4.2/respond.min.js"></script>
<![endif]--><!-- Global site tag (gtag.js) - Google Analytics --><script async src="https://www.googletagmanager.com/gtag/js?id=UA-93322-5"></script><script>
  window.dataLayer = window.dataLayer || [];
  function gtag(){dataLayer.push(arguments);}
  gtag('js', new Date());

  gtag('config', 'UA-93322-5');
</script>
</head>
<body data-spy="scroll" data-target="#toc">


    <div class="container template-article">
      <header><div class="navbar navbar-default navbar-fixed-top" role="navigation">
  <div class="container">
    <div class="navbar-header">
      <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar" aria-expanded="false">
        <span class="sr-only">Toggle navigation</span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
      </button>
      <span class="navbar-brand">
        <a class="navbar-link" href="../index.html">tidybayes</a>
        <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="">3.0.7</span>
      </span>
    </div>

    <div id="navbar" class="navbar-collapse collapse">
      <ul class="nav navbar-nav">
<li>
  <a href="../articles/tidybayes.html">Get started</a>
</li>
<li>
  <a href="../reference/index.html">Reference</a>
</li>
<li class="dropdown">
  <a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" data-bs-toggle="dropdown" aria-expanded="false">
    Articles

    <span class="caret"></span>
  </a>
  <ul class="dropdown-menu" role="menu">
<li>
      <a href="../articles/tidy-brms.html">Extracting and visualizing tidy draws from brms models</a>
    </li>
    <li>
      <a href="../articles/tidy-posterior.html">Using tidybayes with the posterior package</a>
    </li>
    <li>
      <a href="../articles/tidy-rstanarm.html">Extracting and visualizing tidy draws from rstanarm models</a>
    </li>
    <li>
      <a href="../articles/tidybayes-residuals.html">Extracting and visualizing tidy residuals from Bayesian models</a>
    </li>
  </ul>
</li>
<li>
  <a href="../news/index.html">Changelog</a>
</li>
      </ul>
<ul class="nav navbar-nav navbar-right">
<li>
  <a href="https://github.com/mjskay/tidybayes/" class="external-link">
    <span class="fab fa-github fa-lg"></span>

  </a>
</li>
      </ul>
</div>
<!--/.nav-collapse -->
  </div>
<!--/.container -->
</div>
<!--/.navbar -->



      </header><div class="row">
  <div class="col-md-9 contents">
    <div class="page-header toc-ignore">
      <h1 data-toc-skip>Extracting and visualizing tidy residuals from
Bayesian models</h1>
                        <h4 data-toc-skip class="author">Matthew
Kay</h4>
            
            <h4 data-toc-skip class="date">2024-09-14</h4>
      
      <small class="dont-index">Source: <a href="https://github.com/mjskay/tidybayes/blob/HEAD/vignettes/tidybayes-residuals.Rmd" class="external-link"><code>vignettes/tidybayes-residuals.Rmd</code></a></small>
      <div class="hidden name"><code>tidybayes-residuals.Rmd</code></div>

    </div>

    
    
<style type="text/css">
.kable-table table {
  margin-left: 0;
}
img {
  border: none;
}
</style>
<div class="section level2">
<h2 id="introduction">Introduction<a class="anchor" aria-label="anchor" href="#introduction"></a>
</h2>
<p>This vignette describes how to use the <code>tidybayes</code> package
to extract <a href="https://dx.doi.org/10.18637/jss.v059.i10" class="external-link">tidy</a>
data frames of draws from residuals of Bayesian models, and also acts as
a demo for the construction of randomized quantile residuals, a generic
form of residual applicable to a wide range of models, including
censored regressions and models with discrete response variables. For a
more general introduction to <code>tidybayes</code>, see
<code><a href="../articles/tidybayes.html">vignette("tidybayes")</a></code>.</p>
</div>
<div class="section level2">
<h2 id="setup">Setup<a class="anchor" aria-label="anchor" href="#setup"></a>
</h2>
<p>The following libraries are required to run this vignette:</p>
<div class="sourceCode" id="cb1"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://dplyr.tidyverse.org" class="external-link">dplyr</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://purrr.tidyverse.org/" class="external-link">purrr</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://tidyr.tidyverse.org" class="external-link">tidyr</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://mjskay.github.io/ggdist/" class="external-link">ggdist</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://mjskay.github.io/tidybayes/">tidybayes</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://ggplot2.tidyverse.org" class="external-link">ggplot2</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://wilkelab.org/cowplot/" class="external-link">cowplot</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://mc-stan.org/rstan/" class="external-link">rstan</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://github.com/paul-buerkner/brms" class="external-link">brms</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://gganimate.com" class="external-link">gganimate</a></span><span class="op">)</span></span>
<span></span>
<span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/theme_get.html" class="external-link">theme_set</a></span><span class="op">(</span><span class="fu"><a href="http://mjskay.github.io/ggdist/reference/theme_ggdist.html" class="external-link">theme_tidybayes</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span> <span class="fu"><a href="https://wilkelab.org/cowplot/reference/panel_border.html" class="external-link">panel_border</a></span><span class="op">(</span><span class="op">)</span><span class="op">)</span></span></code></pre></div>
<p>These options help Stan run faster:</p>
<div class="sourceCode" id="cb2"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://mc-stan.org/rstan/reference/rstan_options.html" class="external-link">rstan_options</a></span><span class="op">(</span>auto_write <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/options.html" class="external-link">options</a></span><span class="op">(</span>mc.cores <span class="op">=</span> <span class="fu">parallel</span><span class="fu">::</span><span class="fu"><a href="https://rdrr.io/r/parallel/detectCores.html" class="external-link">detectCores</a></span><span class="op">(</span><span class="op">)</span><span class="op">)</span></span></code></pre></div>
</div>
<div class="section level2">
<h2 id="example-dataset">Example dataset<a class="anchor" aria-label="anchor" href="#example-dataset"></a>
</h2>
<p>We‚Äôre going to use a simple dataset to demonstrate residuals that can
be applied to discrete and censored outcomes. We‚Äôll generate data from a
Normal distribution:</p>
<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mtable><mtr><mtd columnalign="right" style="text-align: right"><msubsup><mi>y</mi><mi>i</mi><mo>*</mo></msubsup></mtd><mtd columnalign="left" style="text-align: left"><mo>‚àº</mo><mtext mathvariant="normal">Normal</mtext><mrow><mo stretchy="true" form="prefix">(</mo><mn>0.5</mn><mo>,</mo><mn>1</mn><mo stretchy="true" form="postfix">)</mo></mrow></mtd></mtr><mtr><mtd columnalign="right" style="text-align: right"><msubsup><mi>y</mi><mi>i</mi><mtext mathvariant="normal">lower</mtext></msubsup></mtd><mtd columnalign="left" style="text-align: left"><mo>=</mo><mo stretchy="false" form="prefix">‚åä</mo><msubsup><mi>y</mi><mi>i</mi><mo>*</mo></msubsup><mo stretchy="false" form="postfix">‚åã</mo></mtd></mtr><mtr><mtd columnalign="right" style="text-align: right"><msubsup><mi>y</mi><mi>i</mi><mtext mathvariant="normal">upper</mtext></msubsup></mtd><mtd columnalign="left" style="text-align: left"><mo>=</mo><mo stretchy="false" form="prefix">‚åà</mo><msubsup><mi>y</mi><mi>i</mi><mo>*</mo></msubsup><mo stretchy="false" form="postfix">‚åâ</mo></mtd></mtr></mtable><annotation encoding="application/x-tex">
\begin{align*}
y^*_i &amp;\sim \textrm{Normal}(0.5,1) \\
y^\textrm{lower}_i &amp;= \lfloor{y^*_i}\rfloor \\
y^\textrm{upper}_i &amp;= \lceil{y^*_i}\rceil
\end{align*}
</annotation></semantics></math> Or we can write the vectorized version
of the model:</p>
<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mtable><mtr><mtd columnalign="right" style="text-align: right"><msup><mi>ùê≤</mi><mo>*</mo></msup></mtd><mtd columnalign="left" style="text-align: left"><mo>‚àº</mo><mtext mathvariant="normal">Normal</mtext><mrow><mo stretchy="true" form="prefix">(</mo><mn>0.5</mn><mo>,</mo><mn>1</mn><mo stretchy="true" form="postfix">)</mo></mrow></mtd></mtr><mtr><mtd columnalign="right" style="text-align: right"><msup><mi>ùê≤</mi><mtext mathvariant="normal">lower</mtext></msup></mtd><mtd columnalign="left" style="text-align: left"><mo>=</mo><mo stretchy="false" form="prefix">‚åä</mo><msup><mi>ùê≤</mi><mo>*</mo></msup><mo stretchy="false" form="postfix">‚åã</mo></mtd></mtr><mtr><mtd columnalign="right" style="text-align: right"><msup><mi>ùê≤</mi><mtext mathvariant="normal">upper</mtext></msup></mtd><mtd columnalign="left" style="text-align: left"><mo>=</mo><mo stretchy="false" form="prefix">‚åà</mo><msup><mi>ùê≤</mi><mo>*</mo></msup><mo stretchy="false" form="postfix">‚åâ</mo></mtd></mtr></mtable><annotation encoding="application/x-tex">
\begin{align*}
\mathbf{y}^* &amp;\sim \textrm{Normal}(0.5,1) \\
\mathbf{y}^\textrm{lower} &amp;= \lfloor{\mathbf{y}^*}\rfloor \\
\mathbf{y}^\textrm{upper} &amp;= \lceil{\mathbf{y}^*}\rceil
\end{align*}
</annotation></semantics></math></p>
<p>We don‚Äôt observe each
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msubsup><mi>y</mi><mi>i</mi><mo>*</mo></msubsup><annotation encoding="application/x-tex">y^*_i</annotation></semantics></math>,
only
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msubsup><mi>y</mi><mi>i</mi><mtext mathvariant="normal">lower</mtext></msubsup><annotation encoding="application/x-tex">y^\textrm{lower}_i</annotation></semantics></math>
and
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msubsup><mi>y</mi><mi>i</mi><mtext mathvariant="normal">upper</mtext></msubsup><annotation encoding="application/x-tex">y^\textrm{upper}_i</annotation></semantics></math>,
i.e., the upper and lower bounds of an interval in which
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msubsup><mi>y</mi><mi>i</mi><mo>*</mo></msubsup><annotation encoding="application/x-tex">y^*_i</annotation></semantics></math>
lies. This is <em>interval-censored</em> data.</p>
<p>We can generate it as follows:</p>
<div class="sourceCode" id="cb3"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/Random.html" class="external-link">set.seed</a></span><span class="op">(</span><span class="fl">4118</span><span class="op">)</span></span>
<span><span class="va">n</span> <span class="op">=</span> <span class="fl">100</span></span>
<span></span>
<span><span class="va">cens_df</span> <span class="op">=</span></span>
<span>  <span class="fu"><a href="https://tibble.tidyverse.org/reference/tibble.html" class="external-link">tibble</a></span><span class="op">(</span></span>
<span>    y_star <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Normal.html" class="external-link">rnorm</a></span><span class="op">(</span><span class="va">n</span>, <span class="fl">0.5</span>, <span class="fl">1</span><span class="op">)</span>,</span>
<span>    y_lower <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/Round.html" class="external-link">floor</a></span><span class="op">(</span><span class="va">y_star</span><span class="op">)</span>,</span>
<span>    y_upper <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/Round.html" class="external-link">ceiling</a></span><span class="op">(</span><span class="va">y_star</span><span class="op">)</span>,</span>
<span>    censoring <span class="op">=</span> <span class="st">"interval"</span></span>
<span>  <span class="op">)</span></span></code></pre></div>
<p>A snapshot of the data looks like this:</p>
<div class="sourceCode" id="cb4"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/head.html" class="external-link">head</a></span><span class="op">(</span><span class="va">cens_df</span>, <span class="fl">10</span><span class="op">)</span></span></code></pre></div>
<pre><code><span><span class="co">## <span style="color: #949494;"># A tibble: 10 √ó 4</span></span></span>
<span><span class="co">##     y_star y_lower y_upper censoring</span></span>
<span><span class="co">##      <span style="color: #949494; font-style: italic;">&lt;dbl&gt;</span>   <span style="color: #949494; font-style: italic;">&lt;dbl&gt;</span>   <span style="color: #949494; font-style: italic;">&lt;dbl&gt;</span> <span style="color: #949494; font-style: italic;">&lt;chr&gt;</span>    </span></span>
<span><span class="co">## <span style="color: #BCBCBC;"> 1</span>  0.180        0       1 interval </span></span>
<span><span class="co">## <span style="color: #BCBCBC;"> 2</span>  1.05         1       2 interval </span></span>
<span><span class="co">## <span style="color: #BCBCBC;"> 3</span>  2.05         2       3 interval </span></span>
<span><span class="co">## <span style="color: #BCBCBC;"> 4</span> -<span style="color: #BB0000;">0.512</span>       -<span style="color: #BB0000;">1</span>       0 interval </span></span>
<span><span class="co">## <span style="color: #BCBCBC;"> 5</span>  0.032<span style="text-decoration: underline;">3</span>       0       1 interval </span></span>
<span><span class="co">## <span style="color: #BCBCBC;"> 6</span>  1.18         1       2 interval </span></span>
<span><span class="co">## <span style="color: #BCBCBC;"> 7</span> -<span style="color: #BB0000;">0.707</span>       -<span style="color: #BB0000;">1</span>       0 interval </span></span>
<span><span class="co">## <span style="color: #BCBCBC;"> 8</span>  0.116        0       1 interval </span></span>
<span><span class="co">## <span style="color: #BCBCBC;"> 9</span>  0.183        0       1 interval </span></span>
<span><span class="co">## <span style="color: #BCBCBC;">10</span>  0.385        0       1 interval</span></span></code></pre>
<p>This is a typical tidy format data frame: one observation per
row.</p>
<p>A graphical depiction of the censoring process might be something
like this:</p>
<div class="sourceCode" id="cb6"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">uncensored_plot</span> <span class="op">=</span> <span class="va">cens_df</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>y <span class="op">=</span> <span class="st">""</span>, x <span class="op">=</span> <span class="va">y_star</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="http://mjskay.github.io/ggdist/reference/stat_slab.html" class="external-link">stat_slab</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_jitter.html" class="external-link">geom_jitter</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>y <span class="op">=</span> <span class="fl">0.75</span>, color <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/factor.html" class="external-link">ordered</a></span><span class="op">(</span><span class="va">y_lower</span><span class="op">)</span><span class="op">)</span>, position <span class="op">=</span> <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/position_jitter.html" class="external-link">position_jitter</a></span><span class="op">(</span>height <span class="op">=</span> <span class="fl">0.2</span><span class="op">)</span>, show.legend <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/labs.html" class="external-link">ylab</a></span><span class="op">(</span><span class="cn">NULL</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/scale_continuous.html" class="external-link">scale_x_continuous</a></span><span class="op">(</span>breaks <span class="op">=</span> <span class="op">-</span><span class="fl">4</span><span class="op">:</span><span class="fl">4</span>, limits <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="op">-</span><span class="fl">4</span>, <span class="fl">4</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://wilkelab.org/cowplot/reference/background_grid.html" class="external-link">background_grid</a></span><span class="op">(</span><span class="st">"x"</span><span class="op">)</span></span>
<span></span>
<span><span class="va">censored_plot</span> <span class="op">=</span> <span class="va">cens_df</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>y <span class="op">=</span> <span class="st">""</span>, x <span class="op">=</span> <span class="op">(</span><span class="va">y_lower</span> <span class="op">+</span> <span class="va">y_upper</span><span class="op">)</span><span class="op">/</span><span class="fl">2</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_dotplot.html" class="external-link">geom_dotplot</a></span><span class="op">(</span></span>
<span>    <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>fill <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/factor.html" class="external-link">ordered</a></span><span class="op">(</span><span class="va">y_lower</span><span class="op">)</span><span class="op">)</span>,</span>
<span>    method <span class="op">=</span> <span class="st">"histodot"</span>, origin <span class="op">=</span> <span class="op">-</span><span class="fl">4</span>, binwidth <span class="op">=</span> <span class="fl">1</span>, dotsize <span class="op">=</span> <span class="fl">0.5</span>, stackratio <span class="op">=</span> <span class="fl">.8</span>, show.legend <span class="op">=</span> <span class="cn">FALSE</span>,</span>
<span>    stackgroups <span class="op">=</span> <span class="cn">TRUE</span>, binpositions <span class="op">=</span> <span class="st">"all"</span>, color <span class="op">=</span> <span class="cn">NA</span></span>
<span>  <span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_segment.html" class="external-link">geom_segment</a></span><span class="op">(</span></span>
<span>    <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>x <span class="op">=</span> <span class="va">y</span> <span class="op">+</span> <span class="fl">0.5</span>, xend <span class="op">=</span> <span class="va">y</span> <span class="op">+</span> <span class="fl">0.5</span>, y <span class="op">=</span> <span class="fl">1.75</span>, yend <span class="op">=</span> <span class="fl">1.5</span>, color <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/factor.html" class="external-link">ordered</a></span><span class="op">(</span><span class="va">y</span><span class="op">)</span><span class="op">)</span>,</span>
<span>    data <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/data.frame.html" class="external-link">data.frame</a></span><span class="op">(</span>y <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/unique.html" class="external-link">unique</a></span><span class="op">(</span><span class="va">cens_df</span><span class="op">$</span><span class="va">y_lower</span><span class="op">)</span><span class="op">)</span>, show.legend <span class="op">=</span> <span class="cn">FALSE</span>,</span>
<span>    arrow <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/grid/arrow.html" class="external-link">arrow</a></span><span class="op">(</span>type <span class="op">=</span> <span class="st">"closed"</span>, length <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/grid/unit.html" class="external-link">unit</a></span><span class="op">(</span><span class="fl">7</span>, <span class="st">"points"</span><span class="op">)</span><span class="op">)</span>, linewidth <span class="op">=</span> <span class="fl">1</span></span>
<span>  <span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/labs.html" class="external-link">ylab</a></span><span class="op">(</span><span class="cn">NULL</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/labs.html" class="external-link">xlab</a></span><span class="op">(</span><span class="st">"interval-censored y"</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/scale_continuous.html" class="external-link">scale_x_continuous</a></span><span class="op">(</span>breaks <span class="op">=</span> <span class="op">-</span><span class="fl">4</span><span class="op">:</span><span class="fl">4</span>, limits <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="op">-</span><span class="fl">4</span>, <span class="fl">4</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://wilkelab.org/cowplot/reference/background_grid.html" class="external-link">background_grid</a></span><span class="op">(</span><span class="st">"x"</span><span class="op">)</span></span>
<span></span>
<span><span class="fu"><a href="https://wilkelab.org/cowplot/reference/plot_grid.html" class="external-link">plot_grid</a></span><span class="op">(</span>align <span class="op">=</span> <span class="st">"v"</span>, ncol <span class="op">=</span> <span class="fl">1</span>, rel_heights <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="fl">1</span>, <span class="fl">2.5</span><span class="op">)</span>,</span>
<span>  <span class="va">uncensored_plot</span>,</span>
<span>  <span class="va">censored_plot</span></span>
<span><span class="op">)</span></span></code></pre></div>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-4-1.png" width="180"></p>
</div>
<div class="section level2">
<h2 id="model-for-y_star">Model for <code>y_star</code><a class="anchor" aria-label="anchor" href="#model-for-y_star"></a>
</h2>
<p>Before we fit a model to the censored data, let‚Äôs see what an ideal
model looks like: we‚Äôll fit a model to <code>y_star</code>. With
censored or discrete data in the real-world, we would not be able to fit
such a model, because we would not have observations of
<code>y_star</code>, only the intervals. But this might help us
understand what we‚Äôre looking for in a fit that‚Äôs working.</p>
<div class="sourceCode" id="cb7"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">m_ideal</span> <span class="op">=</span> <span class="fu"><a href="https://paul-buerkner.github.io/brms/reference/brm.html" class="external-link">brm</a></span><span class="op">(</span></span>
<span>  <span class="va">y_star</span> <span class="op">~</span> <span class="fl">1</span>, </span>
<span>  data <span class="op">=</span> <span class="va">cens_df</span>, </span>
<span>  family <span class="op">=</span> <span class="va">student</span>,</span>
<span>  </span>
<span>  file <span class="op">=</span> <span class="st">"models/tidybayes-residuals_m_ideal.rds"</span>  <span class="co"># cache model (can be removed)</span></span>
<span><span class="op">)</span></span></code></pre></div>
<p>The results look like this:</p>
<div class="sourceCode" id="cb8"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">m_ideal</span></span></code></pre></div>
<pre><code><span><span class="co">##  Family: student </span></span>
<span><span class="co">##   Links: mu = identity; sigma = identity; nu = identity </span></span>
<span><span class="co">## Formula: y_star ~ 1 </span></span>
<span><span class="co">##    Data: cens_df (Number of observations: 100) </span></span>
<span><span class="co">##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;</span></span>
<span><span class="co">##          total post-warmup draws = 4000</span></span>
<span><span class="co">## </span></span>
<span><span class="co">## Regression Coefficients:</span></span>
<span><span class="co">##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS</span></span>
<span><span class="co">## Intercept     0.51      0.11     0.30     0.72 1.00     2746     2546</span></span>
<span><span class="co">## </span></span>
<span><span class="co">## Further Distributional Parameters:</span></span>
<span><span class="co">##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS</span></span>
<span><span class="co">## sigma     1.02      0.09     0.85     1.21 1.00     2811     2575</span></span>
<span><span class="co">## nu       21.44     13.28     5.14    55.67 1.00     2709     2508</span></span>
<span><span class="co">## </span></span>
<span><span class="co">## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS</span></span>
<span><span class="co">## and Tail_ESS are effective sample size measures, and Rhat is the potential</span></span>
<span><span class="co">## scale reduction factor on split chains (at convergence, Rhat = 1).</span></span></code></pre>
<div class="section level3">
<h3 id="residuals">Residuals<a class="anchor" aria-label="anchor" href="#residuals"></a>
</h3>
<p><code><a href="../reference/add_predicted_draws.html">add_residual_draws()</a></code> operates much like
<code><a href="../reference/add_predicted_draws.html">add_epred_draws()</a></code> and <code><a href="../reference/add_predicted_draws.html">add_predicted_draws()</a></code>:
it gives us draws from the posterior distribution for each residual:</p>
<div class="sourceCode" id="cb10"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_residual_draws</a></span><span class="op">(</span><span class="va">m_ideal</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>x <span class="op">=</span> <span class="va">.row</span>, y <span class="op">=</span> <span class="va">.residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="http://mjskay.github.io/ggdist/reference/stat_pointinterval.html" class="external-link">stat_pointinterval</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-6-1.png" width="432"></p>
<p>This is a typical residual plot, though with the addition of the
uncertainty associated with each residual afforded by generating draws
from the distribution associated with each residual.</p>
<p>Let‚Äôs check the QQ plot:</p>
<div class="sourceCode" id="cb11"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_residual_draws</a></span><span class="op">(</span><span class="va">m_ideal</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="http://mjskay.github.io/ggdist/reference/point_interval.html" class="external-link">median_qi</a></span><span class="op">(</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">.residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq_line</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-7-1.png" width="432"></p>
<p>This is a typical way of looking at residuals. However, this type of
residual will not generalize to discrete and censored models, so let‚Äôs
consider something else‚Ä¶</p>
</div>
<div class="section level3">
<h3 id="probability-residuals-and-quantile-residuals">Probability residuals and quantile residuals<a class="anchor" aria-label="anchor" href="#probability-residuals-and-quantile-residuals"></a>
</h3>
<p>To derive a more general form of residual, we will use the posterior
predictive distribution,
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><msup><mi>ùê≤</mi><mrow><mo>*</mo><mtext mathvariant="normal">rep</mtext></mrow></msup><mo stretchy="false" form="prefix">|</mo><msup><mi>ùê≤</mi><mo>*</mo></msup></mrow><annotation encoding="application/x-tex">\mathbf{y}^{*\textrm{rep}}|\mathbf{y}^*</annotation></semantics></math>.
This gives us the model‚Äôs predictive distribution for a replication of
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msup><mi>ùê≤</mi><mo>*</mo></msup><annotation encoding="application/x-tex">\mathbf{y}^*</annotation></semantics></math>,
denoted
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msup><mi>ùê≤</mi><mrow><mo>*</mo><mtext mathvariant="normal">rep</mtext></mrow></msup><annotation encoding="application/x-tex">\mathbf{y}^{*\textrm{rep}}</annotation></semantics></math>.
In particular, for each observation, we will calculate the predicted
probability of seeing an observation less than or equal to the
observation we actually got:</p>
<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><msubsup><mi>p</mi><mi>i</mi><mrow><mo>*</mo><mtext mathvariant="normal">resid</mtext></mrow></msubsup><mo>=</mo><mi>P</mi><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>y</mi><mi>i</mi><mrow><mo>*</mo><mtext mathvariant="normal">rep</mtext></mrow></msubsup><mo>‚â§</mo><msubsup><mi>y</mi><mi>i</mi><mo>*</mo></msubsup><mo stretchy="false" form="prefix">|</mo><msup><mi>ùê≤</mi><mo>*</mo></msup><mo stretchy="true" form="postfix">)</mo></mrow></mrow><annotation encoding="application/x-tex">
p^{*\textrm{resid}}_i = P(y^{*\textrm{rep}}_i \le y^*_i|\mathbf{y}^*)
</annotation></semantics></math> This is one type of Bayesian
<strong>posterior predictive p-value</strong>; I will call it a
<strong>probability residual</strong> (<code>p_residual</code>) by
analogy to quantile residuals. If the predictive distribution is
well-calibrated, these probabilities should be uniform (in general these
will not be quite uniform due to the double use of
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>ùê≤</mi><mo>*</mo></mrow><annotation encoding="application/x-tex">\mathbf{y}*</annotation></semantics></math>;
a future iteration of this article will demonstrate the use of
leave-one-out posterior predictive distributions to address that
issue).</p>
<p>If we apply the inverse cumulative distribution function (CDF) of the
standard Normal distribution to these probability residuals, the result
should be approximately standard Normal. These are <strong>quantile
residuals</strong> (<code>z_residual</code>):</p>
<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><msubsup><mi>z</mi><mi>i</mi><mrow><mo>*</mo><mtext mathvariant="normal">resid</mtext></mrow></msubsup><mo>=</mo><msubsup><mi>F</mi><mtext mathvariant="normal">Normal</mtext><mrow><mo>‚àí</mo><mn>1</mn></mrow></msubsup><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>p</mi><mi>i</mi><mrow><mo>*</mo><mtext mathvariant="normal">resid</mtext></mrow></msubsup><mo stretchy="true" form="postfix">)</mo></mrow></mrow><annotation encoding="application/x-tex">
z^{*\textrm{resid}}_i = F_{\textrm{Normal}}^{-1}(p^{*\textrm{resid}}_i)
</annotation></semantics></math></p>
<div class="sourceCode" id="cb12"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_predicted_draws</a></span><span class="op">(</span><span class="va">m_ideal</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/summarise.html" class="external-link">summarise</a></span><span class="op">(</span></span>
<span>    p_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="va">y_star</span><span class="op">)</span>,</span>
<span>    z_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Normal.html" class="external-link">qnorm</a></span><span class="op">(</span><span class="va">p_residual</span><span class="op">)</span>,</span>
<span>    .groups <span class="op">=</span> <span class="st">"drop_last"</span></span>
<span>  <span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">z_residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_abline.html" class="external-link">geom_abline</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-8-1.png" width="432"></p>
<p>This QQ plot shows us the model fits okay.</p>
</div>
</div>
<div class="section level2">
<h2 id="basic-interval-censored-model">Basic interval-censored model<a class="anchor" aria-label="anchor" href="#basic-interval-censored-model"></a>
</h2>
<p>Assuming a situation where we do not have access to the original
(uncensored) data, we will perform a regression on the interval-censored
data. In brms, we can do this by specifying the lower
(<code>y_lower</code>) and upper (<code>y_upper</code>) bounds on each
observation, and the censoring type (<code>censoring</code>), which in
this case is <code>"interval"</code> for all observations.</p>
<div class="sourceCode" id="cb13"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">m</span> <span class="op">=</span> <span class="fu"><a href="https://paul-buerkner.github.io/brms/reference/brm.html" class="external-link">brm</a></span><span class="op">(</span></span>
<span>  <span class="va">y_lower</span> <span class="op">|</span> <span class="fu">cens</span><span class="op">(</span><span class="va">censoring</span>, <span class="va">y_upper</span><span class="op">)</span> <span class="op">~</span> <span class="fl">1</span>, </span>
<span>  data <span class="op">=</span> <span class="va">cens_df</span>,</span>
<span>  </span>
<span>  file <span class="op">=</span> <span class="st">"models/tidybayes-residuals_m.rds"</span>  <span class="co"># cache model (can be removed)</span></span>
<span><span class="op">)</span></span></code></pre></div>
<p>The results look like this:</p>
<div class="sourceCode" id="cb14"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">m</span></span></code></pre></div>
<pre><code><span><span class="co">##  Family: gaussian </span></span>
<span><span class="co">##   Links: mu = identity; sigma = identity </span></span>
<span><span class="co">## Formula: y_lower | cens(censoring, y_upper) ~ 1 </span></span>
<span><span class="co">##    Data: cens_df (Number of observations: 100) </span></span>
<span><span class="co">##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;</span></span>
<span><span class="co">##          total post-warmup draws = 4000</span></span>
<span><span class="co">## </span></span>
<span><span class="co">## Regression Coefficients:</span></span>
<span><span class="co">##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS</span></span>
<span><span class="co">## Intercept     0.56      0.11     0.33     0.78 1.00     3191     2633</span></span>
<span><span class="co">## </span></span>
<span><span class="co">## Further Distributional Parameters:</span></span>
<span><span class="co">##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS</span></span>
<span><span class="co">## sigma     1.09      0.08     0.94     1.27 1.00     2774     2513</span></span>
<span><span class="co">## </span></span>
<span><span class="co">## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS</span></span>
<span><span class="co">## and Tail_ESS are effective sample size measures, and Rhat is the potential</span></span>
<span><span class="co">## scale reduction factor on split chains (at convergence, Rhat = 1).</span></span></code></pre>
<p>Again, <code><a href="../reference/add_predicted_draws.html">add_residual_draws()</a></code> lets us look at the
residuals:</p>
<div class="sourceCode" id="cb16"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_residual_draws</a></span><span class="op">(</span><span class="va">m</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>x <span class="op">=</span> <span class="va">.row</span>, y <span class="op">=</span> <span class="va">.residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="http://mjskay.github.io/ggdist/reference/stat_pointinterval.html" class="external-link">stat_pointinterval</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<pre><code><span><span class="co">## Warning: Results may not be meaningful for censored models.</span></span></code></pre>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-10-1.png" width="432"></p>
<p>But notice the horizontal striations due to censoring. This obscures
patterns we might wish to see in the diagnostic plots. Patterns like
this often show up in censored regression, logistic regression, Poisson
regression, and other discrete regression models. Notice how
<code>brms</code> also gives us a warning message that the residuals may
not be meaningful due to censoring.</p>
<p>Let‚Äôs look at the QQ plot:</p>
<div class="sourceCode" id="cb18"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_residual_draws</a></span><span class="op">(</span><span class="va">m</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="http://mjskay.github.io/ggdist/reference/point_interval.html" class="external-link">median_qi</a></span><span class="op">(</span><span class="va">.residual</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">.residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq_line</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<pre><code><span><span class="co">## Warning: Results may not be meaningful for censored models.</span></span></code></pre>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-11-1.png" width="432"></p>
<p>Again, the striations are making interpretation difficult. But there
is hope!</p>
<div class="section level3">
<h3 id="randomized-quantile-residuals">Randomized quantile residuals<a class="anchor" aria-label="anchor" href="#randomized-quantile-residuals"></a>
</h3>
<p>We can extend the ideas of <strong>probability residuals</strong> to
make <strong>randomized probability residuals</strong>.</p>
<p>Where the probability residual for the uncensored model was something
like this:</p>
<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><msubsup><mi>p</mi><mi>i</mi><mrow><mo>*</mo><mtext mathvariant="normal">resid</mtext></mrow></msubsup><mo>=</mo><mi>P</mi><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>y</mi><mi>i</mi><mrow><mo>*</mo><mtext mathvariant="normal">rep</mtext></mrow></msubsup><mo>&lt;</mo><msubsup><mi>y</mi><mi>i</mi><mo>*</mo></msubsup><mo stretchy="false" form="prefix">|</mo><msup><mi>ùê≤</mi><mo>*</mo></msup><mo stretchy="true" form="postfix">)</mo></mrow></mrow><annotation encoding="application/x-tex">
p^{*\textrm{resid}}_i = P(y^{*\textrm{rep}}_i &lt; y^*_i|\mathbf{y}^*)
</annotation></semantics></math> We can‚Äôt calculate the residual for a
censored observation this way, as we do not know what
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msup><mi>ùê≤</mi><mo>*</mo></msup><annotation encoding="application/x-tex">\mathbf{y}^*</annotation></semantics></math>
is, only
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msup><mi>ùê≤</mi><mtext mathvariant="normal">lower</mtext></msup><annotation encoding="application/x-tex">\mathbf{y}^\textrm{lower}</annotation></semantics></math>
and
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msup><mi>ùê≤</mi><mtext mathvariant="normal">upper</mtext></msup><annotation encoding="application/x-tex">\mathbf{y}^\textrm{upper}</annotation></semantics></math>.
However, we do know that
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msubsup><mi>p</mi><mi>i</mi><mrow><mo>*</mo><mtext mathvariant="normal">resid</mtext></mrow></msubsup><annotation encoding="application/x-tex">p^{*\textrm{resid}}_i</annotation></semantics></math>
should be in the following interval:</p>
<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mtable><mtr><mtd columnalign="right" style="text-align: right"></mtd><mtd columnalign="left" style="text-align: left"></mtd><mtd columnalign="right" style="text-align: right"><msubsup><mi>p</mi><mi>i</mi><mrow><mo>*</mo><mtext mathvariant="normal">resid</mtext></mrow></msubsup></mtd><mtd columnalign="left" style="text-align: left"><mo>‚àà</mo><mo stretchy="false" form="prefix">(</mo><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">lower</mtext></msubsup><mo>,</mo><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">upper</mtext></msubsup><mo stretchy="false" form="postfix">]</mo></mtd></mtr><mtr><mtd columnalign="right" style="text-align: right"></mtd><mtd columnalign="left" style="text-align: left"><mtext mathvariant="normal">where</mtext></mtd><mtd columnalign="right" style="text-align: right"><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">lower</mtext></msubsup></mtd><mtd columnalign="left" style="text-align: left"><mo>=</mo><mi>P</mi><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>y</mi><mi>i</mi><mrow><mo>*</mo><mtext mathvariant="normal">rep</mtext></mrow></msubsup><mo>&lt;</mo><msubsup><mi>y</mi><mi>i</mi><mtext mathvariant="normal">lower</mtext></msubsup><mo stretchy="false" form="prefix">|</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">lower</mtext></msup><mo>,</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">upper</mtext></msup><mo stretchy="true" form="postfix">)</mo></mrow></mtd></mtr><mtr><mtd columnalign="right" style="text-align: right"></mtd><mtd columnalign="left" style="text-align: left"></mtd><mtd columnalign="right" style="text-align: right"><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">upper</mtext></msubsup></mtd><mtd columnalign="left" style="text-align: left"><mo>=</mo><mi>P</mi><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>y</mi><mi>i</mi><mrow><mo>*</mo><mtext mathvariant="normal">rep</mtext></mrow></msubsup><mo>‚â§</mo><msubsup><mi>y</mi><mi>i</mi><mtext mathvariant="normal">upper</mtext></msubsup><mo stretchy="false" form="prefix">|</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">lower</mtext></msup><mo>,</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">upper</mtext></msup><mo stretchy="true" form="postfix">)</mo></mrow></mtd></mtr></mtable><annotation encoding="application/x-tex">
\begin{align*}
&amp;&amp; p^{*\textrm{resid}}_i &amp;\in (p^\textrm{lower}_i, p^\textrm{upper}_i] \\
&amp;\textrm{where}&amp; p^\textrm{lower}_i &amp;= P(y^{*\textrm{rep}}_i &lt; y^\textrm{lower}_i|\mathbf{y}^\textrm{lower},\mathbf{y}^\textrm{upper}) \\
&amp;&amp;      p^\textrm{upper}_i &amp;= P(y^{*\textrm{rep}}_i \le y^\textrm{upper}_i|\mathbf{y}^\textrm{lower},\mathbf{y}^\textrm{upper})
\end{align*}
</annotation></semantics></math></p>
<p>If we define a <strong>randomized probability residual</strong>
as:</p>
<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">resid</mtext></msubsup><mo>‚àº</mo><mtext mathvariant="normal">Uniform</mtext><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">lower</mtext></msubsup><mo>,</mo><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">upper</mtext></msubsup><mo stretchy="true" form="postfix">)</mo></mrow></mrow><annotation encoding="application/x-tex">
p^{\textrm{resid}}_i \sim \textrm{Uniform}(p^\textrm{lower}_i, p^\textrm{upper}_i)
</annotation></semantics></math></p>
<p>Then these residuals will be uniformly random (see <a href="https://doi.org/10.1080/10618600.1996.10474708" class="external-link">Dunn and Smyth
1996</a>), and we can similarly define a <strong>randomized quantile
residual</strong> by passing the probability residuals through the
Normal inverse CDF:</p>
<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><msubsup><mi>z</mi><mi>i</mi><mtext mathvariant="normal">resid</mtext></msubsup><mo>=</mo><msubsup><mi>F</mi><mtext mathvariant="normal">Normal</mtext><mrow><mo>‚àí</mo><mn>1</mn></mrow></msubsup><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">resid</mtext></msubsup><mo stretchy="true" form="postfix">)</mo></mrow></mrow><annotation encoding="application/x-tex">
z^{\textrm{resid}}_i = F_{\textrm{Normal}}^{-1}(p^{\textrm{resid}}_i)
</annotation></semantics></math></p>
<p>In code, that looks like this:</p>
<div class="sourceCode" id="cb20"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_predicted_draws</a></span><span class="op">(</span><span class="va">m</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/summarise.html" class="external-link">summarise</a></span><span class="op">(</span></span>
<span>    p_lower <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="va">y_lower</span><span class="op">)</span>,</span>
<span>    p_upper <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="va">y_upper</span><span class="op">)</span>,</span>
<span>    p_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Uniform.html" class="external-link">runif</a></span><span class="op">(</span><span class="fl">1</span>, <span class="va">p_lower</span>, <span class="va">p_upper</span><span class="op">)</span>,</span>
<span>    z_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Normal.html" class="external-link">qnorm</a></span><span class="op">(</span><span class="va">p_residual</span><span class="op">)</span>,</span>
<span>    .groups <span class="op">=</span> <span class="st">"drop_last"</span></span>
<span>  <span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>x <span class="op">=</span> <span class="va">.row</span>, y <span class="op">=</span> <span class="va">z_residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_point.html" class="external-link">geom_point</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-12-1.png" width="432"></p>
<p>The striations are gone! And a QQ plot:</p>
<div class="sourceCode" id="cb21"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_predicted_draws</a></span><span class="op">(</span><span class="va">m</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/summarise.html" class="external-link">summarise</a></span><span class="op">(</span></span>
<span>    p_lower <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="va">y_lower</span><span class="op">)</span>,</span>
<span>    p_upper <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="va">y_upper</span><span class="op">)</span>,</span>
<span>    p_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Uniform.html" class="external-link">runif</a></span><span class="op">(</span><span class="fl">1</span>, <span class="va">p_lower</span>, <span class="va">p_upper</span><span class="op">)</span>,</span>
<span>    z_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Normal.html" class="external-link">qnorm</a></span><span class="op">(</span><span class="va">p_residual</span><span class="op">)</span>,</span>
<span>    .groups <span class="op">=</span> <span class="st">"drop_last"</span></span>
<span>  <span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">z_residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_abline.html" class="external-link">geom_abline</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-13-1.png" width="432"></p>
<p>Looks decent. However, our residuals are now subject to noise‚Äîwe
actually have a distribution of possible QQ plots. So, let‚Äôs check on
some of them to make sure this nice-looking one is not just a fluke.
We‚Äôll employ <a href="https://mucollective.co/project/hops-trends" class="external-link">HOPs</a> to do
that:</p>
<div class="sourceCode" id="cb22"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># NOTE: ordinarily I would use a large number of frames (k), </span></span>
<span><span class="co"># say 50 or 100, but I am keeping it small for the sake of </span></span>
<span><span class="co"># keeping these examples small</span></span>
<span><span class="va">k</span> <span class="op">=</span> <span class="fl">10</span></span>
<span></span>
<span><span class="va">p</span> <span class="op">=</span> <span class="va">cens_df</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_predicted_draws</a></span><span class="op">(</span><span class="va">m</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/summarise.html" class="external-link">summarise</a></span><span class="op">(</span></span>
<span>    p_lower <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="va">y_lower</span><span class="op">)</span>,</span>
<span>    p_upper <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="va">y_upper</span><span class="op">)</span>,</span>
<span>    p_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/list.html" class="external-link">list</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/stats/Uniform.html" class="external-link">runif</a></span><span class="op">(</span><span class="va">k</span>, <span class="va">p_lower</span>, <span class="va">p_upper</span><span class="op">)</span><span class="op">)</span>,</span>
<span>    residual_draw <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/list.html" class="external-link">list</a></span><span class="op">(</span><span class="fl">1</span><span class="op">:</span><span class="va">k</span><span class="op">)</span>,</span>
<span>    .groups <span class="op">=</span> <span class="st">"drop_last"</span></span>
<span>  <span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://tidyr.tidyverse.org/reference/unnest.html" class="external-link">unnest</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="va">p_residual</span>, <span class="va">residual_draw</span><span class="op">)</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/mutate.html" class="external-link">mutate</a></span><span class="op">(</span>z_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Normal.html" class="external-link">qnorm</a></span><span class="op">(</span><span class="va">p_residual</span><span class="op">)</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">z_residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_abline.html" class="external-link">geom_abline</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://gganimate.com/reference/transition_manual.html" class="external-link">transition_manual</a></span><span class="op">(</span><span class="va">residual_draw</span><span class="op">)</span></span>
<span></span>
<span><span class="fu"><a href="https://gganimate.com/reference/animate.html" class="external-link">animate</a></span><span class="op">(</span><span class="va">p</span>, nframes <span class="op">=</span> <span class="va">k</span>, width <span class="op">=</span> <span class="fl">384</span>, height <span class="op">=</span> <span class="fl">384</span>, units <span class="op">=</span> <span class="st">"px"</span>, res <span class="op">=</span> <span class="fl">96</span>, dev <span class="op">=</span> <span class="st">"ragg_png"</span><span class="op">)</span></span></code></pre></div>
<p><img src="tidybayes-residuals_resid_hops_1.gif"></p>
<p>Looks good.</p>
</div>
</div>
<div class="section level2">
<h2 id="what-if-the-model-does-not-fit-well">What if the model does not fit well?<a class="anchor" aria-label="anchor" href="#what-if-the-model-does-not-fit-well"></a>
</h2>
<p>Let‚Äôs instead generate the data from a t distribution, which should
give us some more extreme values:</p>
<div class="sourceCode" id="cb23"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/Random.html" class="external-link">set.seed</a></span><span class="op">(</span><span class="fl">41181</span><span class="op">)</span></span>
<span><span class="va">n</span> <span class="op">=</span> <span class="fl">100</span></span>
<span></span>
<span><span class="va">cens_df_t</span> <span class="op">=</span></span>
<span>  <span class="fu"><a href="https://tibble.tidyverse.org/reference/tibble.html" class="external-link">tibble</a></span><span class="op">(</span></span>
<span>    y <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/TDist.html" class="external-link">rt</a></span><span class="op">(</span><span class="va">n</span>, <span class="fl">3</span><span class="op">)</span> <span class="op">+</span> <span class="fl">0.5</span>,</span>
<span>    y_lower <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/Round.html" class="external-link">floor</a></span><span class="op">(</span><span class="va">y</span><span class="op">)</span>,</span>
<span>    y_upper <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/Round.html" class="external-link">ceiling</a></span><span class="op">(</span><span class="va">y</span><span class="op">)</span>,</span>
<span>    censoring <span class="op">=</span> <span class="st">"interval"</span></span>
<span>  <span class="op">)</span></span></code></pre></div>
<p>Note the presence of outliers:</p>
<div class="sourceCode" id="cb24"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">uncensored_plot</span> <span class="op">=</span> <span class="va">cens_df_t</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>y <span class="op">=</span> <span class="st">""</span>, x <span class="op">=</span> <span class="va">y</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="http://mjskay.github.io/ggdist/reference/stat_slab.html" class="external-link">stat_slab</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_jitter.html" class="external-link">geom_jitter</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>y <span class="op">=</span> <span class="fl">0.75</span>, color <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/factor.html" class="external-link">ordered</a></span><span class="op">(</span><span class="va">y_lower</span><span class="op">)</span><span class="op">)</span>, position <span class="op">=</span> <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/position_jitter.html" class="external-link">position_jitter</a></span><span class="op">(</span>height <span class="op">=</span> <span class="fl">0.2</span><span class="op">)</span>, show.legend <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/labs.html" class="external-link">ylab</a></span><span class="op">(</span><span class="cn">NULL</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/scale_continuous.html" class="external-link">scale_x_continuous</a></span><span class="op">(</span>breaks <span class="op">=</span> <span class="op">-</span><span class="fl">10</span><span class="op">:</span><span class="fl">10</span>, limits <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="op">-</span><span class="fl">10</span>, <span class="fl">10</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://wilkelab.org/cowplot/reference/background_grid.html" class="external-link">background_grid</a></span><span class="op">(</span><span class="st">"x"</span><span class="op">)</span></span>
<span></span>
<span><span class="va">censored_plot</span> <span class="op">=</span> <span class="va">cens_df_t</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>y <span class="op">=</span> <span class="st">""</span>, x <span class="op">=</span> <span class="op">(</span><span class="va">y_lower</span> <span class="op">+</span> <span class="va">y_upper</span><span class="op">)</span><span class="op">/</span><span class="fl">2</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_dotplot.html" class="external-link">geom_dotplot</a></span><span class="op">(</span></span>
<span>    <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>fill <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/factor.html" class="external-link">ordered</a></span><span class="op">(</span><span class="va">y_lower</span><span class="op">)</span><span class="op">)</span>,</span>
<span>    method <span class="op">=</span> <span class="st">"histodot"</span>, origin <span class="op">=</span> <span class="op">-</span><span class="fl">4</span>, binwidth <span class="op">=</span> <span class="fl">1</span>, dotsize <span class="op">=</span> <span class="fl">0.5</span>, stackratio <span class="op">=</span> <span class="fl">.8</span>, show.legend <span class="op">=</span> <span class="cn">FALSE</span>,</span>
<span>    stackgroups <span class="op">=</span> <span class="cn">TRUE</span>, binpositions <span class="op">=</span> <span class="st">"all"</span>, color <span class="op">=</span> <span class="cn">NA</span></span>
<span>  <span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_segment.html" class="external-link">geom_segment</a></span><span class="op">(</span></span>
<span>    <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>x <span class="op">=</span> <span class="va">y</span> <span class="op">+</span> <span class="fl">0.5</span>, xend <span class="op">=</span> <span class="va">y</span> <span class="op">+</span> <span class="fl">0.5</span>, y <span class="op">=</span> <span class="fl">1.75</span>, yend <span class="op">=</span> <span class="fl">1.5</span>, color <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/factor.html" class="external-link">ordered</a></span><span class="op">(</span><span class="va">y</span><span class="op">)</span><span class="op">)</span>,</span>
<span>    data <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/data.frame.html" class="external-link">data.frame</a></span><span class="op">(</span>y <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/unique.html" class="external-link">unique</a></span><span class="op">(</span><span class="va">cens_df_t</span><span class="op">$</span><span class="va">y_lower</span><span class="op">)</span><span class="op">)</span>, show.legend <span class="op">=</span> <span class="cn">FALSE</span>,</span>
<span>    arrow <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/grid/arrow.html" class="external-link">arrow</a></span><span class="op">(</span>type <span class="op">=</span> <span class="st">"closed"</span>, length <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/grid/unit.html" class="external-link">unit</a></span><span class="op">(</span><span class="fl">7</span>, <span class="st">"points"</span><span class="op">)</span><span class="op">)</span>, linewidth <span class="op">=</span> <span class="fl">1</span></span>
<span>  <span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/labs.html" class="external-link">ylab</a></span><span class="op">(</span><span class="cn">NULL</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/labs.html" class="external-link">xlab</a></span><span class="op">(</span><span class="st">"interval-censored y"</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/scale_continuous.html" class="external-link">scale_x_continuous</a></span><span class="op">(</span>breaks <span class="op">=</span> <span class="op">-</span><span class="fl">10</span><span class="op">:</span><span class="fl">10</span>, limits <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="op">-</span><span class="fl">10</span>, <span class="fl">10</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://wilkelab.org/cowplot/reference/background_grid.html" class="external-link">background_grid</a></span><span class="op">(</span><span class="st">"x"</span><span class="op">)</span></span>
<span></span>
<span><span class="fu"><a href="https://wilkelab.org/cowplot/reference/plot_grid.html" class="external-link">plot_grid</a></span><span class="op">(</span>align <span class="op">=</span> <span class="st">"v"</span>, ncol <span class="op">=</span> <span class="fl">1</span>, rel_heights <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="fl">1</span>, <span class="fl">2.25</span><span class="op">)</span>,</span>
<span>  <span class="va">uncensored_plot</span>,</span>
<span>  <span class="va">censored_plot</span></span>
<span><span class="op">)</span></span></code></pre></div>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-16-1.png" width="432"></p>
<div class="section level3">
<h3 id="model">Model<a class="anchor" aria-label="anchor" href="#model"></a>
</h3>
<p>Now, fitting the same model as before:</p>
<div class="sourceCode" id="cb25"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">m_t1</span> <span class="op">=</span> <span class="fu"><a href="https://paul-buerkner.github.io/brms/reference/brm.html" class="external-link">brm</a></span><span class="op">(</span></span>
<span>  <span class="va">y_lower</span> <span class="op">|</span> <span class="fu">cens</span><span class="op">(</span><span class="va">censoring</span>, <span class="va">y_upper</span><span class="op">)</span> <span class="op">~</span> <span class="fl">1</span>,</span>
<span>  data <span class="op">=</span> <span class="va">cens_df_t</span>,</span>
<span></span>
<span>  file <span class="op">=</span> <span class="st">"models/tidybayes-residuals_m_t1"</span>  <span class="co"># cache model (can be removed)</span></span>
<span><span class="op">)</span></span></code></pre></div>
<p>And checking the QQ plot:</p>
<div class="sourceCode" id="cb26"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df_t</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_residual_draws</a></span><span class="op">(</span><span class="va">m_t1</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="http://mjskay.github.io/ggdist/reference/point_interval.html" class="external-link">median_qi</a></span><span class="op">(</span><span class="va">.residual</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">.residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq_line</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<pre><code><span><span class="co">## Warning: Results may not be meaningful for censored models.</span></span></code></pre>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-17-1.png" width="432"></p>
<p>It looks like this could be problematic, but the striations make
diagnosing the problem difficult.</p>
</div>
<div class="section level3">
<h3 id="randomized-quantile-residuals-1">Randomized quantile residuals<a class="anchor" aria-label="anchor" href="#randomized-quantile-residuals-1"></a>
</h3>
<p>Let‚Äôs try randomized quantile residuals again:</p>
<div class="sourceCode" id="cb28"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df_t</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_predicted_draws</a></span><span class="op">(</span><span class="va">m_t1</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/summarise.html" class="external-link">summarise</a></span><span class="op">(</span></span>
<span>    p_lower <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="va">y_lower</span><span class="op">)</span>,</span>
<span>    p_upper <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="va">y_upper</span><span class="op">)</span>,</span>
<span>    p_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Uniform.html" class="external-link">runif</a></span><span class="op">(</span><span class="fl">1</span>, <span class="va">p_lower</span>, <span class="va">p_upper</span><span class="op">)</span>,</span>
<span>    z_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Normal.html" class="external-link">qnorm</a></span><span class="op">(</span><span class="va">p_residual</span><span class="op">)</span>,</span>
<span>    .groups <span class="op">=</span> <span class="st">"drop_last"</span></span>
<span>  <span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">z_residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_abline.html" class="external-link">geom_abline</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<pre><code><span><span class="co">## Warning: Removed 1 row containing non-finite outside the scale range (`stat_qq()`).</span></span></code></pre>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-18-1.png" width="432"></p>
<p>Does not look good for the model‚Äîthis s-shape is characteristic of
excess kurtosis (aka fat tails), which is exactly what a t distribution
should produce. So let‚Äôs try a different model.</p>
</div>
</div>
<div class="section level2">
<h2 id="robust-interval-censored-model">Robust interval-censored model<a class="anchor" aria-label="anchor" href="#robust-interval-censored-model"></a>
</h2>
<p>A common model to try to ‚Äúrobustify‚Äù linear regression is to swap out
the Normal response distribution for a Student‚Äôs t distribution, which
has fatter tails (it also happens that this is the distribution used to
generate the data above, so I am clearly cheating, but also: <em>this
had better work</em>).</p>
<div class="sourceCode" id="cb30"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">m_t2</span> <span class="op">=</span> <span class="fu"><a href="https://paul-buerkner.github.io/brms/reference/brm.html" class="external-link">brm</a></span><span class="op">(</span></span>
<span>  <span class="va">y_lower</span> <span class="op">|</span> <span class="fu">cens</span><span class="op">(</span><span class="va">censoring</span>, <span class="va">y_upper</span><span class="op">)</span> <span class="op">~</span> <span class="fl">1</span>, </span>
<span>  data <span class="op">=</span> <span class="va">cens_df_t</span>,</span>
<span>  family <span class="op">=</span> <span class="va">student</span>,</span>
<span>  </span>
<span>  file <span class="op">=</span> <span class="st">"models/tidybayes-residuals_m_t2.rds"</span>  <span class="co"># cache model (can be removed)</span></span>
<span><span class="op">)</span></span></code></pre></div>
<p>Let‚Äôs look at the QQ plot of the residuals:</p>
<div class="sourceCode" id="cb31"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df_t</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_residual_draws</a></span><span class="op">(</span><span class="va">m_t2</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="http://mjskay.github.io/ggdist/reference/point_interval.html" class="external-link">median_qi</a></span><span class="op">(</span><span class="va">.residual</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">.residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq_line</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<pre><code><span><span class="co">## Warning: Results may not be meaningful for censored models.</span></span></code></pre>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-19-1.png" width="432"></p>
<p>Odd, these don‚Äôt look great! But again, these residuals are not to be
trusted‚Ä¶</p>
<div class="section level3">
<h3 id="randomized-quantile-residuals-2">Randomized quantile residuals<a class="anchor" aria-label="anchor" href="#randomized-quantile-residuals-2"></a>
</h3>
<p>Let‚Äôs check out the randomized quantile residuals:</p>
<div class="sourceCode" id="cb33"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df_t</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_predicted_draws</a></span><span class="op">(</span><span class="va">m_t2</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/summarise.html" class="external-link">summarise</a></span><span class="op">(</span></span>
<span>    p_lower <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="va">y_lower</span><span class="op">)</span>,</span>
<span>    p_upper <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="va">y_upper</span><span class="op">)</span>,</span>
<span>    p_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Uniform.html" class="external-link">runif</a></span><span class="op">(</span><span class="fl">1</span>, <span class="va">p_lower</span>, <span class="va">p_upper</span><span class="op">)</span>,</span>
<span>    z_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Normal.html" class="external-link">qnorm</a></span><span class="op">(</span><span class="va">p_residual</span><span class="op">)</span>,</span>
<span>    .groups <span class="op">=</span> <span class="st">"drop_last"</span></span>
<span>  <span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">z_residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_abline.html" class="external-link">geom_abline</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-20-1.png" width="432"></p>
<p>Looking pretty good! As well they should, since this is the same
model that generated the data.</p>
</div>
</div>
<div class="section level2">
<h2 id="ordinal-regression">Ordinal regression<a class="anchor" aria-label="anchor" href="#ordinal-regression"></a>
</h2>
<p>But we‚Äôre not done yet‚Äîanother sensible choice of model for this data
might have been an ordinal regression. So let‚Äôs try that too.</p>
<p>First, we‚Äôll add a column that transforms the data into an ordered
factor:</p>
<div class="sourceCode" id="cb34"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df_o</span> <span class="op">=</span> <span class="va">cens_df_t</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/mutate.html" class="external-link">mutate</a></span><span class="op">(</span>y_factor <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/factor.html" class="external-link">ordered</a></span><span class="op">(</span><span class="va">y_lower</span><span class="op">)</span><span class="op">)</span></span></code></pre></div>
<p>Then we‚Äôll fit an ordinal regression model. We have to adjust the fit
parameters a bit to get it to converge:</p>
<div class="sourceCode" id="cb35"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">m_o</span> <span class="op">=</span> <span class="fu"><a href="https://paul-buerkner.github.io/brms/reference/brm.html" class="external-link">brm</a></span><span class="op">(</span></span>
<span>  <span class="va">y_factor</span> <span class="op">~</span> <span class="fl">1</span>, </span>
<span>  data <span class="op">=</span> <span class="va">cens_df_o</span>, </span>
<span>  family <span class="op">=</span> <span class="va">cumulative</span>, </span>
<span>  prior <span class="op">=</span> <span class="fu"><a href="https://paul-buerkner.github.io/brms/reference/set_prior.html" class="external-link">prior</a></span><span class="op">(</span><span class="fu"><a href="https://mc-stan.org/rstanarm/reference/priors.html" class="external-link">normal</a></span><span class="op">(</span><span class="fl">0</span>, <span class="fl">10</span><span class="op">)</span>, class <span class="op">=</span> <span class="va">Intercept</span><span class="op">)</span>, </span>
<span>  control <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/list.html" class="external-link">list</a></span><span class="op">(</span>adapt_delta <span class="op">=</span> <span class="fl">0.99</span><span class="op">)</span>,</span>
<span>  </span>
<span>  file <span class="op">=</span> <span class="st">"models/tidybayes-residuals_m_o.rds"</span>  <span class="co"># cache model (can be removed)</span></span>
<span><span class="op">)</span></span></code></pre></div>
<p>Now, let‚Äôs check out those residuals:</p>
<div class="sourceCode" id="cb36"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df_o</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_residual_draws</a></span><span class="op">(</span><span class="va">m_o</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="http://mjskay.github.io/ggdist/reference/point_interval.html" class="external-link">median_qi</a></span><span class="op">(</span><span class="va">.residual</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">.residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq_line</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<pre><code><span><span class="co">## Error: Predictive errors are not defined for ordinal or categorical models.</span></span></code></pre>
<p>Gah! <code>brms</code> will not even produce residuals for models
like these. This is very sensible, as standard residuals here are hard
to define. Fortunately, we can still construct randomized quantile
residuals.</p>
<div class="section level3">
<h3 id="randomized-quantile-residuals-for-a-discrete-distribution">Randomized quantile residuals for a discrete distribution<a class="anchor" aria-label="anchor" href="#randomized-quantile-residuals-for-a-discrete-distribution"></a>
</h3>
<p>For the interval-censored model, we had:</p>
<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mtable><mtr><mtd columnalign="right" style="text-align: right"></mtd><mtd columnalign="left" style="text-align: left"></mtd><mtd columnalign="right" style="text-align: right"><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">lower</mtext></msubsup></mtd><mtd columnalign="left" style="text-align: left"><mo>=</mo><mi>P</mi><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>y</mi><mi>i</mi><mrow><mo>*</mo><mtext mathvariant="normal">rep</mtext></mrow></msubsup><mo>&lt;</mo><msubsup><mi>y</mi><mi>i</mi><mtext mathvariant="normal">lower</mtext></msubsup><mo stretchy="false" form="prefix">|</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">lower</mtext></msup><mo>,</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">upper</mtext></msup><mo stretchy="true" form="postfix">)</mo></mrow></mtd></mtr><mtr><mtd columnalign="right" style="text-align: right"></mtd><mtd columnalign="left" style="text-align: left"></mtd><mtd columnalign="right" style="text-align: right"><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">upper</mtext></msubsup></mtd><mtd columnalign="left" style="text-align: left"><mo>=</mo><mi>P</mi><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>y</mi><mi>i</mi><mrow><mo>*</mo><mtext mathvariant="normal">rep</mtext></mrow></msubsup><mo>‚â§</mo><msubsup><mi>y</mi><mi>i</mi><mtext mathvariant="normal">upper</mtext></msubsup><mo stretchy="false" form="prefix">|</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">lower</mtext></msup><mo>,</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">upper</mtext></msup><mo stretchy="true" form="postfix">)</mo></mrow></mtd></mtr><mtr><mtd columnalign="right" style="text-align: right"></mtd><mtd columnalign="left" style="text-align: left"></mtd><mtd columnalign="right" style="text-align: right"><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">resid</mtext></msubsup></mtd><mtd columnalign="left" style="text-align: left"><mo>‚àº</mo><mtext mathvariant="normal">Uniform</mtext><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">lower</mtext></msubsup><mo>,</mo><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">upper</mtext></msubsup><mo stretchy="true" form="postfix">)</mo></mrow></mtd></mtr></mtable><annotation encoding="application/x-tex">
\begin{align*}
&amp;&amp; p^\textrm{lower}_i &amp;= P(y^{*\textrm{rep}}_i &lt; y^\textrm{lower}_i|\mathbf{y}^\textrm{lower},\mathbf{y}^\textrm{upper}) \\
&amp;&amp;      p^\textrm{upper}_i &amp;= P(y^{*\textrm{rep}}_i \le y^\textrm{upper}_i|\mathbf{y}^\textrm{lower},\mathbf{y}^\textrm{upper}) \\
&amp;&amp; p^{\textrm{resid}}_i &amp;\sim \textrm{Uniform}(p^\textrm{lower}_i, p^\textrm{upper}_i)
\end{align*}
</annotation></semantics></math></p>
<p>But that was for a model with observations of the form
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mo stretchy="false" form="prefix">(</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">lower</mtext></msup><mo>,</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">upper</mtext></msup><mo stretchy="false" form="postfix">]</mo></mrow><annotation encoding="application/x-tex">(\mathbf{y}^\textrm{lower}, \mathbf{y}^\textrm{upper}]</annotation></semantics></math>
and which can produce a posterior predictive distribution for the latent
variable
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><msup><mi>ùê≤</mi><mrow><mo>*</mo><mtext mathvariant="normal">rep</mtext></mrow></msup><mo stretchy="false" form="prefix">|</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">lower</mtext></msup><mo>,</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">upper</mtext></msup></mrow><annotation encoding="application/x-tex">\mathbf{y}^{*\textrm{rep}}|\mathbf{y}^\textrm{lower}, \mathbf{y}^\textrm{upper}</annotation></semantics></math>.
Now our observations are discrete values of
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msup><mi>ùê≤</mi><mtext mathvariant="normal">factor</mtext></msup><annotation encoding="application/x-tex">\mathbf{y}^\textrm{factor}</annotation></semantics></math>,
and our posterior predictive distribution,
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><msup><mi>ùê≤</mi><mtext mathvariant="normal">rep</mtext></msup><mo stretchy="false" form="prefix">|</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">factor</mtext></msup></mrow><annotation encoding="application/x-tex">\mathbf{y}^\textrm{rep}|\mathbf{y}^\textrm{factor}</annotation></semantics></math>,
is also discrete. So we have to tweak the definition a bit:</p>
<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mtable><mtr><mtd columnalign="right" style="text-align: right"></mtd><mtd columnalign="left" style="text-align: left"></mtd><mtd columnalign="right" style="text-align: right"><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">lower</mtext></msubsup></mtd><mtd columnalign="left" style="text-align: left"><mo>=</mo><mi>P</mi><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>y</mi><mi>i</mi><mtext mathvariant="normal">rep</mtext></msubsup><mo>&lt;</mo><msubsup><mi>y</mi><mi>i</mi><mtext mathvariant="normal">factor</mtext></msubsup><mo stretchy="false" form="prefix">|</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">factor</mtext></msup><mo stretchy="true" form="postfix">)</mo></mrow></mtd></mtr><mtr><mtd columnalign="right" style="text-align: right"></mtd><mtd columnalign="left" style="text-align: left"></mtd><mtd columnalign="right" style="text-align: right"><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">upper</mtext></msubsup></mtd><mtd columnalign="left" style="text-align: left"><mo>=</mo><mi>P</mi><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>y</mi><mi>i</mi><mtext mathvariant="normal">rep</mtext></msubsup><mo>‚â§</mo><msubsup><mi>y</mi><mi>i</mi><mtext mathvariant="normal">factor</mtext></msubsup><mo stretchy="false" form="prefix">|</mo><msup><mi>ùê≤</mi><mtext mathvariant="normal">factor</mtext></msup><mo stretchy="true" form="postfix">)</mo></mrow></mtd></mtr><mtr><mtd columnalign="right" style="text-align: right"></mtd><mtd columnalign="left" style="text-align: left"></mtd><mtd columnalign="right" style="text-align: right"><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">resid</mtext></msubsup></mtd><mtd columnalign="left" style="text-align: left"><mo>‚àº</mo><mtext mathvariant="normal">Uniform</mtext><mrow><mo stretchy="true" form="prefix">(</mo><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">lower</mtext></msubsup><mo>,</mo><msubsup><mi>p</mi><mi>i</mi><mtext mathvariant="normal">upper</mtext></msubsup><mo stretchy="true" form="postfix">)</mo></mrow></mtd></mtr></mtable><annotation encoding="application/x-tex">
\begin{align*}
&amp;&amp; p^\textrm{lower}_i &amp;= P(y^{\textrm{rep}}_i &lt; y^\textrm{factor}_i|\mathbf{y}^\textrm{factor}) \\
&amp;&amp;      p^\textrm{upper}_i &amp;= P(y^{\textrm{rep}}_i \le y^\textrm{factor}_i|\mathbf{y}^\textrm{factor}) \\
&amp;&amp; p^{\textrm{resid}}_i &amp;\sim \textrm{Uniform}(p^\textrm{lower}_i, p^\textrm{upper}_i)
\end{align*}
</annotation></semantics></math></p>
<p>Then we proceed as before:</p>
<div class="sourceCode" id="cb38"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cens_df_o</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_predicted_draws</a></span><span class="op">(</span><span class="va">m_o</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/mutate.html" class="external-link">mutate</a></span><span class="op">(</span>.prediction <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/factor.html" class="external-link">ordered</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/levels.html" class="external-link">levels</a></span><span class="op">(</span><span class="va">y_factor</span><span class="op">)</span><span class="op">[</span><span class="va">.prediction</span><span class="op">]</span>, levels <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/levels.html" class="external-link">levels</a></span><span class="op">(</span><span class="va">y_factor</span><span class="op">)</span><span class="op">)</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/summarise.html" class="external-link">summarise</a></span><span class="op">(</span></span>
<span>    p_lower <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="va">y_factor</span><span class="op">)</span>,</span>
<span>    p_upper <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="va">.prediction</span> <span class="op">&lt;=</span> <span class="va">y_factor</span><span class="op">)</span>,</span>
<span>    p_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Uniform.html" class="external-link">runif</a></span><span class="op">(</span><span class="fl">1</span>, <span class="va">p_lower</span>, <span class="va">p_upper</span><span class="op">)</span>,</span>
<span>    z_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Normal.html" class="external-link">qnorm</a></span><span class="op">(</span><span class="va">p_residual</span><span class="op">)</span>,</span>
<span>    .groups <span class="op">=</span> <span class="st">"drop_last"</span></span>
<span>  <span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">z_residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_abline.html" class="external-link">geom_abline</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-23-1.png" width="432"></p>
<p>It is not too surprising that the ordinal regression also works here,
even though it is not the data generating model. Ordinal regressions are
pretty robust.</p>
</div>
</div>
<div class="section level2">
<h2 id="a-function-for-probability-residuals">A function for probability residuals<a class="anchor" aria-label="anchor" href="#a-function-for-probability-residuals"></a>
</h2>
<p>Here‚Äôs a function that can calculate randomized probability residuals
for you. A future version of <code>tidybayes</code> will likely supply
some variant of this function, but for now, since it is experimental, I
have not folded it into the main package yet. Use at your own risk:</p>
<div class="sourceCode" id="cb39"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://rlang.r-lib.org" class="external-link">rlang</a></span><span class="op">)</span></span>
<span><span class="va">make_probability_residuals</span> <span class="op">=</span> <span class="kw">function</span><span class="op">(</span><span class="va">data</span>, <span class="va">prediction</span>, <span class="va">y</span>, <span class="va">y_upper</span> <span class="op">=</span> <span class="cn">NA</span>, <span class="va">n</span> <span class="op">=</span> <span class="fl">1</span><span class="op">)</span> <span class="op">{</span></span>
<span>  <span class="va">.prediction</span> <span class="op">=</span> <span class="fu"><a href="https://rlang.r-lib.org/reference/enquo.html" class="external-link">enquo</a></span><span class="op">(</span><span class="va">prediction</span><span class="op">)</span></span>
<span>  <span class="va">.y</span> <span class="op">=</span> <span class="fu"><a href="https://rlang.r-lib.org/reference/enquo.html" class="external-link">enquo</a></span><span class="op">(</span><span class="va">y</span><span class="op">)</span></span>
<span>  <span class="va">.y_upper</span> <span class="op">=</span> <span class="fu"><a href="https://rlang.r-lib.org/reference/enquo.html" class="external-link">enquo</a></span><span class="op">(</span><span class="va">y_upper</span><span class="op">)</span></span>
<span></span>
<span>  <span class="kw">if</span> <span class="op">(</span><span class="fu"><a href="https://rlang.r-lib.org/reference/eval_tidy.html" class="external-link">eval_tidy</a></span><span class="op">(</span><span class="fu"><a href="https://rlang.r-lib.org/reference/expr.html" class="external-link">expr</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/factor.html" class="external-link">is.factor</a></span><span class="op">(</span><span class="op">!</span><span class="op">!</span><span class="va">.prediction</span><span class="op">)</span> <span class="op">&amp;&amp;</span> <span class="op">!</span><span class="fu"><a href="https://rdrr.io/r/base/factor.html" class="external-link">is.ordered</a></span><span class="op">(</span><span class="op">!</span><span class="op">!</span><span class="va">.prediction</span><span class="op">)</span><span class="op">)</span>, <span class="va">data</span><span class="op">)</span><span class="op">)</span> <span class="op">{</span></span>
<span>    <span class="va">data</span> <span class="op">=</span> <span class="fu"><a href="https://dplyr.tidyverse.org/reference/mutate.html" class="external-link">mutate</a></span><span class="op">(</span><span class="va">data</span>, <span class="op">!</span><span class="op">!</span><span class="va">.prediction</span> <span class="op">:=</span> <span class="fu"><a href="https://rdrr.io/r/base/factor.html" class="external-link">ordered</a></span><span class="op">(</span><span class="op">!</span><span class="op">!</span><span class="va">.prediction</span>, levels <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/levels.html" class="external-link">levels</a></span><span class="op">(</span><span class="op">!</span><span class="op">!</span><span class="va">.prediction</span><span class="op">)</span><span class="op">)</span><span class="op">)</span></span>
<span>  <span class="op">}</span></span>
<span>  </span>
<span>  <span class="kw">if</span> <span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/NA.html" class="external-link">is.na</a></span><span class="op">(</span><span class="fu"><a href="https://rlang.r-lib.org/reference/set_expr.html" class="external-link">get_expr</a></span><span class="op">(</span><span class="va">.y_upper</span><span class="op">)</span><span class="op">)</span><span class="op">)</span> <span class="op">{</span></span>
<span>    <span class="co">#no y_upper provided, use y as y_upper</span></span>
<span>    <span class="va">data</span> <span class="op">=</span> <span class="fu"><a href="https://dplyr.tidyverse.org/reference/summarise.html" class="external-link">summarise</a></span><span class="op">(</span><span class="va">data</span>,</span>
<span>      .p_lower <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="op">!</span><span class="op">!</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="op">!</span><span class="op">!</span><span class="va">.y</span><span class="op">)</span>,</span>
<span>      .p_upper <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="op">!</span><span class="op">!</span><span class="va">.prediction</span> <span class="op">&lt;=</span> <span class="op">!</span><span class="op">!</span><span class="va">.y</span><span class="op">)</span>,</span>
<span>      .groups <span class="op">=</span> <span class="st">"drop_last"</span></span>
<span>    <span class="op">)</span></span>
<span>  <span class="op">}</span> <span class="kw">else</span> <span class="op">{</span></span>
<span>    <span class="co">#y_upper should be a vector, and if an entry in it is NA, use the entry from y</span></span>
<span>    <span class="va">data</span> <span class="op">=</span> <span class="fu"><a href="https://dplyr.tidyverse.org/reference/summarise.html" class="external-link">summarise</a></span><span class="op">(</span><span class="va">data</span>,</span>
<span>      .p_lower <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="op">!</span><span class="op">!</span><span class="va">.prediction</span> <span class="op">&lt;</span> <span class="op">!</span><span class="op">!</span><span class="va">.y</span><span class="op">)</span>,</span>
<span>      .p_upper <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean</a></span><span class="op">(</span><span class="op">!</span><span class="op">!</span><span class="va">.prediction</span> <span class="op">&lt;=</span> <span class="fu"><a href="https://rdrr.io/r/base/ifelse.html" class="external-link">ifelse</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/NA.html" class="external-link">is.na</a></span><span class="op">(</span><span class="op">!</span><span class="op">!</span><span class="va">.y_upper</span><span class="op">)</span>, <span class="op">!</span><span class="op">!</span><span class="va">.y</span>, <span class="op">!</span><span class="op">!</span><span class="va">.y_upper</span><span class="op">)</span><span class="op">)</span>,</span>
<span>      .groups <span class="op">=</span> <span class="st">"drop_last"</span></span>
<span>    <span class="op">)</span></span>
<span>  <span class="op">}</span></span>
<span>  </span>
<span>  <span class="va">data</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>    <span class="fu"><a href="https://dplyr.tidyverse.org/reference/mutate.html" class="external-link">mutate</a></span><span class="op">(</span></span>
<span>      .p_residual <span class="op">=</span> <span class="fu"><a href="https://purrr.tidyverse.org/reference/map2.html" class="external-link">map2</a></span><span class="op">(</span><span class="va">.p_lower</span>, <span class="va">.p_upper</span>, <span class="va">runif</span>, n <span class="op">=</span> <span class="op">!</span><span class="op">!</span><span class="va">n</span><span class="op">)</span>,</span>
<span>      .residual_draw <span class="op">=</span> <span class="fu"><a href="https://purrr.tidyverse.org/reference/map.html" class="external-link">map</a></span><span class="op">(</span><span class="va">.p_residual</span>, <span class="va">seq_along</span><span class="op">)</span></span>
<span>    <span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>    <span class="fu"><a href="https://tidyr.tidyverse.org/reference/unnest.html" class="external-link">unnest</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="va">.p_residual</span>, <span class="va">.residual_draw</span><span class="op">)</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>    <span class="fu"><a href="https://dplyr.tidyverse.org/reference/mutate.html" class="external-link">mutate</a></span><span class="op">(</span>.z_residual <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Normal.html" class="external-link">qnorm</a></span><span class="op">(</span><span class="va">.p_residual</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="op">}</span></span></code></pre></div>
</div>
<div class="section level2">
<h2 id="logistic-regression">Logistic regression<a class="anchor" aria-label="anchor" href="#logistic-regression"></a>
</h2>
<p>Consider a simple dataset with a single binary outcome:</p>
<div class="sourceCode" id="cb40"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/Random.html" class="external-link">set.seed</a></span><span class="op">(</span><span class="fl">51919</span><span class="op">)</span></span>
<span></span>
<span><span class="va">bin_df</span> <span class="op">=</span> <span class="fu"><a href="https://tibble.tidyverse.org/reference/tibble.html" class="external-link">tibble</a></span><span class="op">(</span></span>
<span>  y <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Uniform.html" class="external-link">runif</a></span><span class="op">(</span><span class="fl">100</span><span class="op">)</span> <span class="op">&gt;</span> <span class="fl">0.3</span></span>
<span><span class="op">)</span></span></code></pre></div>
<p>A model for this might be the following:</p>
<div class="sourceCode" id="cb41"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">m_bin</span> <span class="op">=</span> <span class="fu"><a href="https://paul-buerkner.github.io/brms/reference/brm.html" class="external-link">brm</a></span><span class="op">(</span></span>
<span>  <span class="va">y</span> <span class="op">~</span> <span class="fl">1</span>, </span>
<span>  data <span class="op">=</span> <span class="va">bin_df</span>, </span>
<span>  family <span class="op">=</span> <span class="va">bernoulli</span>,</span>
<span>  </span>
<span>  file <span class="op">=</span> <span class="st">"models/tidybayes-residuals_m_bin.rds"</span>  <span class="co"># cache model (can be removed)</span></span>
<span><span class="op">)</span></span></code></pre></div>
<p>Again, standard residuals are hard to interpret:</p>
<div class="sourceCode" id="cb42"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">bin_df</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_residual_draws</a></span><span class="op">(</span><span class="va">m_bin</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="http://mjskay.github.io/ggdist/reference/point_interval.html" class="external-link">median_qi</a></span><span class="op">(</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">.residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq_line</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<p><img src="tidybayes-residuals_files/figure-html/unnamed-chunk-26-1.png" width="432"></p>
<p>Randomized quantile residuals help us see that the fit is probably
reasonable. In this case (just to switch it up), we‚Äôll look at the
probability residuals instead, comparing them against a uniform
distribution, like so:</p>
<div class="sourceCode" id="cb43"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># NOTE: ordinarily I would use a large number of frames (k), </span></span>
<span><span class="co"># say 50 or 100, but I am keeping it small for the sake of </span></span>
<span><span class="co"># keeping these examples small</span></span>
<span><span class="va">k</span> <span class="op">=</span> <span class="fl">10</span></span>
<span></span>
<span><span class="va">p</span> <span class="op">=</span> <span class="va">bin_df</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="../reference/add_predicted_draws.html">add_predicted_draws</a></span><span class="op">(</span><span class="va">m_bin</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu">make_probability_residuals</span><span class="op">(</span><span class="va">.prediction</span>, <span class="va">y</span>, n <span class="op">=</span> <span class="va">k</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>sample <span class="op">=</span> <span class="va">.p_residual</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_qq.html" class="external-link">geom_qq</a></span><span class="op">(</span>distribution <span class="op">=</span> <span class="va">qunif</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_abline.html" class="external-link">geom_abline</a></span><span class="op">(</span><span class="op">)</span> <span class="op">+</span></span>
<span>  <span class="fu"><a href="https://gganimate.com/reference/transition_manual.html" class="external-link">transition_manual</a></span><span class="op">(</span><span class="va">.residual_draw</span><span class="op">)</span></span>
<span></span>
<span><span class="fu"><a href="https://gganimate.com/reference/animate.html" class="external-link">animate</a></span><span class="op">(</span><span class="va">p</span>, nframes <span class="op">=</span> <span class="va">k</span>, width <span class="op">=</span> <span class="fl">384</span>, height <span class="op">=</span> <span class="fl">384</span>, units <span class="op">=</span> <span class="st">"px"</span>, res <span class="op">=</span> <span class="fl">96</span>, dev <span class="op">=</span> <span class="st">"ragg_png"</span><span class="op">)</span></span></code></pre></div>
<p><img src="tidybayes-residuals_resid_hops_2.gif"></p>
</div>
  </div>

  <div class="col-md-3 hidden-xs hidden-sm" id="pkgdown-sidebar">

      </div>

</div>



      <footer><div class="copyright">
  <p></p>
<p>Developed by <a href="https://www.mjskay.com" class="external-link">Matthew Kay</a>.</p>
</div>

<div class="pkgdown">
  <p></p>
<p>Site built with <a href="https://pkgdown.r-lib.org/" class="external-link">pkgdown</a> 2.1.0.</p>
</div>

      </footer>
</div>






  </body>
</html>
