<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/simultaneous-intervals-for-smooths-revisited/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.fromthebottomoftheheap.net/2016/12/15/simultaneous-interval-revisited/"> From the Bottom of the Heap - R</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>
Eighteen months ago I <a target="_blank" rel="nofollow" href="http://www.fromthebottomoftheheap.net/2014/06/16/simultaneous-confidence-intervals-for-derivatives/">wrote a post</a> in which I described the use of simulation from the posterior distribution of a fitted GAM to derive simultaneous confidence intervals for the derivatives of a penalised spline. It was a nice post that attracted some interest. It was also wrong. I have no idea what I was thinking when I thought the intervals described in that post were simultaneous. Here I hope to rectify that past mistake.
</p>
<p>
I’ll tackle the issue of simultaneous intervals for the derivatives of penalised spline in a follow-up post. Here, I demonstrate one way to compute a simultaneous interval for a penalised spline in a fitted GAM. As example data, I’ll use the strontium isotope data set included in the <strong>SemiPar</strong> package, and which is extensively analyzed in the monograph <em>Semiparametric Regression</em> <span data-cites="Ruppert2003-pt">(Ruppert, Wand, and Carroll 2003)</span>. First, load the packages we’ll need as well as the data, which is data set <code>fossil</code>. If you don’t have <strong>SemiPar</strong> installed, install it using <code>install.packages(“SemiPar”)</code> before proceeding
</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>&quot;mgcv&quot;</span><span>)</span>
<span>library</span><span>(</span><span>&quot;ggplot2&quot;</span><span>)</span>
<span>theme_set</span><span>(</span><span>theme_bw</span><span>())</span>
<span>data</span><span>(</span><span>fossil</span><span>,</span> <span>package</span> <span>=</span> <span>&quot;SemiPar&quot;</span><span>)</span></code></pre>
</figure>
<p>
The <code>fossil</code> data set includes two variables and is a time series of strontium isotope measurements on samples from a sediment core. The data are shown below using <code>ggplot()</code>
</p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>(</span><span>fossil</span><span>,</span> <span>aes</span><span>(</span><span>x</span> <span>=</span> <span>age</span><span>,</span> <span>y</span> <span>=</span> <span>strontium.ratio</span><span>))</span> <span>+</span>
    <span>geom_point</span><span>()</span></code></pre>
</figure>
<figure>
<img alt="The strontium isotope example data used in the post" src="https://i0.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/simultaneous-intervals-revisited-plot-fossil-data-1.png?w=456" data-recalc-dims="1">
<figcaption>
The strontium isotope example data used in the post<br>
</figcaption>
</figure>
<p>
The aim of the analysis of these data is to model how the measured strontium isotope ratio changed through time, using a GAM to estimate the clearly non-linear change in the response. I won’t cover how the GAM is fitted and what all the options are here, but a reasonable GAM for these data is fitted using <strong>mgcv</strong> and <code>gam()</code>
</p>
<figure>
<pre><code data-lang="r"><span>m</span> <span>&lt;-</span> <span>gam</span><span>(</span><span>strontium.ratio</span> <span>~</span> <span>s</span><span>(</span><span>age</span><span>,</span> <span>k</span> <span>=</span> <span>20</span><span>),</span> <span>data</span> <span>=</span> <span>fossil</span><span>,</span> <span>method</span> <span>=</span> <span>&quot;REML&quot;</span><span>)</span></code></pre>
</figure>
<p>
The essentially arbitrary default for <code>k</code>, the basis dimension of the spline, is changed to <code>20</code> as there is a modest amount of non-linearity in the strontium isotopes ratio time series. By using <code>method = “REML”</code>, the penalised spline model is expressed as a linear mixed model with the wiggly bits of the spline treated as random effects, and is estimated using restricted maximum likelihood; <code>method = “ML”</code> would also work here.
</p>
<p>
The fitted model uses ~12 effective degrees of freedom (which wouldn’t have been achievable with the default of <code>k = 10</code>!)
</p>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>m</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Family: gaussian 
Link function: identity 

Formula:
strontium.ratio ~ s(age, k = 20)

Parametric coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 7.074e-01  2.435e-06  290527   &lt;2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
         edf Ref.df     F p-value    
s(age) 11.52  13.88 62.07  &lt;2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.891   Deviance explained = 90.3%
-REML = -932.05  Scale est. = 6.2839e-10  n = 106</code></pre>
</figure>
<p>
The fitted spline captures the main variation in strontium isotope ratio values; the output from <code>plot.gam()</code> is shown below
</p>
<figure>
<pre><code data-lang="r"><span>plot</span><span>(</span><span>m</span><span>,</span> <span>shade</span> <span>=</span> <span>TRUE</span><span>,</span> <span>seWithMean</span> <span>=</span> <span>TRUE</span><span>,</span> <span>residuals</span> <span>=</span> <span>TRUE</span><span>,</span> <span>pch</span> <span>=</span> <span>16</span><span>,</span> <span>cex</span> <span>=</span> <span>0.8</span><span>)</span></code></pre>
</figure>
<figure>
<img alt="The fitted penalised spline with approximate 95% point-wise confidence interval, as produced with plot.gam()" src="https://i0.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/simultaneous-intervals-revisited-gam-plot-1.png?w=456" data-recalc-dims="1">
<figcaption>
The fitted penalised spline with approximate 95% point-wise confidence interval, as produced with <code>plot.gam()</code><br>
</figcaption>
</figure>
<p>
The confidence interval shown around the fitted spline is a 95% Bayesian credible interval. For reasons that don’t need to concern us right now, this interval has a surprising frequentist interpretation as a 95% <em>“across the function”</em> interval <span data-cites="Nychka1988-rz Marra2012-bq">(Nychka 1988; Marra and Wood 2012)</span>; under repeated resampling from the population 95% of such confidence intervals will contain the true function. Such “across the function” intervals are quite intuitive, but, as we’ll see shortly, they don’t reflect the uncertainty in the fitted function; far fewer than 95% of splines drawn from the posterior distribution of the fitted GAM would lie within the confidence interval shown in the plot above.
</p>
<p>
How to compute a simultaneous interval for a spline is a well studied problem and a number of solutions have been proposed in the literature. Here I follow <span data-cites="Ruppert2003-pt">Ruppert, Wand, and Carroll (2003)</span> and use a simulation-based approach to generate a simultaneous interval. We proceed by considering a simultaneous confidence interval for a function <span>(f(x))</span> at a set of <span>(M)</span> locations in <span>(x)</span>; we’ll refer to these locations, following the notation of <span data-cites="Ruppert2003-pt">Ruppert, Wand, and Carroll (2003)</span> by
</p>
<p>
<span>[ = (g_1, g_2, , g_M) ]</span>
</p>
<p>
The true function over <span>()</span>, <span>()</span>, is defined as the vector of evaluations of <span>(f)</span> at each of the <span>(M)</span> locations
</p>
<p>
<span>[
<p>]
</p></span>
<p>
and the corresponding estimate of the true function given by the fitted GAM as <span>()</span>. The difference between the true function and our unbiased estimator is given by
</p>
<p>
<span>[
<p>]
</p></span>
<p>
where <span>()</span> is the evaluation of the basis functions at the locations <span>()</span>, and the thing in square brackets is the bias in the estimated model coefficients, which we assume to be mean 0 and follows, approximately, a multivariate normal distribution with mean vector <span>()</span> and covariance matrix <span>()</span>
</p>
<p>
<span>[
<p> N (, ) ]
</p></span>
<p>
Having got those definitions out of the way, the 100(1 – <span>()</span>)% simultaneous confidence interval is
</p>
<p>
<span>[
<p>]
</p></span>
<p>
where <span>(m_{1 – })</span> is the 1 – <span>()</span> quantile of the random variable
</p>
<p>
<span>[ <em>{x } | | </em>{1 M} | | ]</span>
</p>
<p>
Yep, that was <em>exactly</em> my reaction when I first read this section of <span data-cites="Ruppert2003-pt">Ruppert, Wand, and Carroll (2003)</span>!
</p>
<p>
Let’s deal with the left-hand side of the equation first. The <span>()</span> refers to the <em>supremum</em> or the <em>least upper bound</em>; this is the least value of <span>()</span>, the set of all values of which we observed subset <span>(x)</span>, that is <em>greater</em> than all of the values in the subset. Often this is the maximum value of the subset. This is what is indicated by the right-hand side of the equation; we want the maximum (absolute) value of the ratio over all values in <span>()</span>.
</p>
<p>
The fractions in both sides of the equation correspond to the standardized deviation between the true function and the model estimate, and we consider the <em>maximum absolute</em> standardized deviation. We don’t usually know the distribution of the maximum absolute standardized deviation but we need this to access its quantiles. However, we can closely approximate the distribution via simulation. The difference here is that rather than simulating from the posterior of the model as we have done in earlier posts on this blog, this time we simulate from the multivariate normal distribution with mean vector <span>()</span> and covariance matrix <span>()</span>, the Bayesian covariance matrix of the fitted model. For each simulation we find the maximum absolute standardized deviation of the fitted function from the true function over the grid of <span>(x)</span> values we are considering. Then we collect all these maxima, sort them and either take the 1 – <span>()</span> probability quantile of the maxima, or the maximum with rank <span>((1 – ) / N )</span>.
</p>
<p>
OK, that’s enough of words and crazy equations. Implementing this in R is going to be easier than those equations might suggest. I’ll run through the code we need line by line. First we define a simple function to generate random values from a multivariate normal: this is in the manual for <strong>mgcv</strong> and saves us loading another package just for this:
</p>
<figure>
<pre><code data-lang="r"><span>rmvn</span> <span>&lt;-</span> <span>function</span><span>(</span><span>n</span><span>,</span> <span>mu</span><span>,</span> <span>sig</span><span>)</span> <span>{</span> <span>## MVN random deviates
</span>    <span>L</span> <span>&lt;-</span> <span>mroot</span><span>(</span><span>sig</span><span>)</span>
    <span>m</span> <span>&lt;-</span> <span>ncol</span><span>(</span><span>L</span><span>)</span>
    <span>t</span><span>(</span><span>mu</span> <span>+</span> <span>L</span> <span>%*%</span> <span>matrix</span><span>(</span><span>rnorm</span><span>(</span><span>m</span><span>*</span><span>n</span><span>),</span> <span>m</span><span>,</span> <span>n</span><span>))</span>
<span>}</span></code></pre>
</figure>
<p>
Next we extract a few things that we need from the fitted GAM
</p>
<figure>
<pre><code data-lang="r"><span>Vb</span> <span>&lt;-</span> <span>vcov</span><span>(</span><span>m</span><span>)</span>
<span>newd</span> <span>&lt;-</span> <span>with</span><span>(</span><span>fossil</span><span>,</span> <span>data.frame</span><span>(</span><span>age</span> <span>=</span> <span>seq</span><span>(</span><span>min</span><span>(</span><span>age</span><span>),</span> <span>max</span><span>(</span><span>age</span><span>),</span> <span>length</span> <span>=</span> <span>200</span><span>)))</span>
<span>pred</span> <span>&lt;-</span> <span>predict</span><span>(</span><span>m</span><span>,</span> <span>newd</span><span>,</span> <span>se.fit</span> <span>=</span> <span>TRUE</span><span>)</span>
<span>se.fit</span> <span>&lt;-</span> <span>pred</span><span>$</span><span>se.fit</span></code></pre>
</figure>
<p>
The first is the Bayesian covariance matrix of the model coefficients, <span>()</span>. This <span>()</span> is conditional upon the smoothing parameter(s). If you want a version that adjusts for the smoothing parameters being estimated rather than known values, add <code>unconditional = TRUE</code> to the <code>vcov()</code> call. Second, we define our grid of <span>(x)</span> values over which we want a confidence band. Then we generate predictions and standard errors from the model for the grid of values. The last line just extracts out the standard errors of the fitted values for use later.
</p>
<p>
Now we are ready to generate simulations of the maximum absolute standardized deviation of the fitted model from the true model. We set the pseudo-random seed to make the results reproducible and specify the number of simulations to generate.
</p>
<figure>
<pre><code data-lang="r"><span>set.seed</span><span>(</span><span>42</span><span>)</span>
<span>N</span> <span>&lt;-</span> <span>10000</span></code></pre>
</figure>
<p>
Next, we want <code>N</code> draws from <span>(
<p>), which is approximately distributed multivariate normal with mean vector <span>()</span> and covariance matrix <code>Vb</code>
</p></span>
<figure>
<pre><code data-lang="r"><span>BUdiff</span> <span>&lt;-</span> <span>rmvn</span><span>(</span><span>N</span><span>,</span> <span>mu</span> <span>=</span> <span>rep</span><span>(</span><span>0</span><span>,</span> <span>nrow</span><span>(</span><span>Vb</span><span>)),</span> <span>sig</span> <span>=</span> <span>Vb</span><span>)</span></code></pre>
</figure>
<p>
Now we calculate <span>((x) – f(x))</span>, which is given by <span>(
<p>) evaluated at the grid of <span>(x)</span> values
</p></span>
<figure>
<pre><code data-lang="r"><span>Cg</span> <span>&lt;-</span> <span>predict</span><span>(</span><span>m</span><span>,</span> <span>newd</span><span>,</span> <span>type</span> <span>=</span> <span>&quot;lpmatrix&quot;</span><span>)</span>
<span>simDev</span> <span>&lt;-</span> <span>Cg</span> <span>%*%</span> <span>t</span><span>(</span><span>BUdiff</span><span>)</span></code></pre>
</figure>
<p>
The first line evaluates the basis function at <span>()</span> and the second line computes the deviations between the fitted and true parameters. Then we find the absolute values of the standardized deviations from the true model. Here we do this in a single step for all simulations using <code>sweep()</code>
</p>
<figure>
<pre><code data-lang="r"><span>absDev</span> <span>&lt;-</span> <span>abs</span><span>(</span><span>sweep</span><span>(</span><span>simDev</span><span>,</span> <span>1</span><span>,</span> <span>se.fit</span><span>,</span> <span>FUN</span> <span>=</span> <span>&quot;/&quot;</span><span>))</span></code></pre>
</figure>
<p>
The maximum of the absolute standardized deviations at the grid of <span>(x)</span> values for each simulation is computed via an <code>apply()</code> call
</p>
<figure>
<pre><code data-lang="r"><span>masd</span> <span>&lt;-</span> <span>apply</span><span>(</span><span>absDev</span><span>,</span> <span>2L</span><span>,</span> <span>max</span><span>)</span></code></pre>
</figure>
<p>
The last step is to find the critical value used to scale the standard errors to yield the simultaneous interval; here we calculate the critical value for a 95% simultaneous confidence interval/band
</p>
<figure>
<pre><code data-lang="r"><span>crit</span> <span>&lt;-</span> <span>quantile</span><span>(</span><span>masd</span><span>,</span> <span>prob</span> <span>=</span> <span>0.95</span><span>,</span> <span>type</span> <span>=</span> <span>8</span><span>)</span></code></pre>
</figure>
<p>
The critical value estimated above is 3.205. Intervals generated using this value will be 1.6 times larger than the point-wise interval shown above.
</p>
<p>
Now that we have the critical value, we can calculate the simultaneous confidence interval. In the code block below I first add the grid of values (<code>newd</code>) to the fitted values and standard errors at those new values and then augment this with upper and lower limits for a 95% simultaneous confidence interval (<code>uprS</code> and <code>lwrS</code>), as well as the usual 95% point-wise intervals for comparison (<code>uprP</code> and <code>lwrP</code>). Then I plot the two intervals:
</p>
<figure>
<pre><code data-lang="r"><span>pred</span> <span>&lt;-</span> <span>transform</span><span>(</span><span>cbind</span><span>(</span><span>data.frame</span><span>(</span><span>pred</span><span>),</span> <span>newd</span><span>),</span>
                  <span>uprP</span> <span>=</span> <span>fit</span> <span>+</span> <span>(</span><span>2</span> <span>*</span> <span>se.fit</span><span>),</span>
                  <span>lwrP</span> <span>=</span> <span>fit</span> <span>-</span> <span>(</span><span>2</span> <span>*</span> <span>se.fit</span><span>),</span>
                  <span>uprS</span> <span>=</span> <span>fit</span> <span>+</span> <span>(</span><span>crit</span> <span>*</span> <span>se.fit</span><span>),</span>
                  <span>lwrS</span> <span>=</span> <span>fit</span> <span>-</span> <span>(</span><span>crit</span> <span>*</span> <span>se.fit</span><span>))</span>
<span>ggplot</span><span>(</span><span>pred</span><span>,</span> <span>aes</span><span>(</span><span>x</span> <span>=</span> <span>age</span><span>))</span> <span>+</span>
    <span>geom_ribbon</span><span>(</span><span>aes</span><span>(</span><span>ymin</span> <span>=</span> <span>lwrS</span><span>,</span> <span>ymax</span> <span>=</span> <span>uprS</span><span>),</span> <span>alpha</span> <span>=</span> <span>0.2</span><span>,</span> <span>fill</span> <span>=</span> <span>&quot;red&quot;</span><span>)</span> <span>+</span>
    <span>geom_ribbon</span><span>(</span><span>aes</span><span>(</span><span>ymin</span> <span>=</span> <span>lwrP</span><span>,</span> <span>ymax</span> <span>=</span> <span>uprP</span><span>),</span> <span>alpha</span> <span>=</span> <span>0.2</span><span>,</span> <span>fill</span> <span>=</span> <span>&quot;red&quot;</span><span>)</span> <span>+</span>
    <span>labs</span><span>(</span><span>y</span> <span>=</span> <span>&quot;Strontium isotope ratio&quot;</span><span>,</span>
         <span>x</span> <span>=</span> <span>&quot;Age [Ma BP]&quot;</span><span>)</span></code></pre>
</figure>
<figure>
<img alt="Comparison of point-wise and simultaneous 95% confidence intervals for the fitted GAM" src="https://i0.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/simultaneous-intervals-revisited-confidence-intervals-1.png?w=456" data-recalc-dims="1">
<figcaption>
Comparison of point-wise and simultaneous 95% confidence intervals for the fitted GAM<br>
</figcaption>
</figure>
<p>
Finally, I’m going to look at the coverage properties of the interval we just created, which is something I should have done in the older post as it would have shown, as we’ll see, that the old interval I wrote about wasn’t even close to having the correct coverage properties.
</p>
<p>
Start by drawing a large sample from the posterior distribution of the fitted model. Note that this time, we’re simulating from a multivariate normal with mean vector given by the estimated model coefficients
</p>
<figure>
<pre><code data-lang="r"><span>sims</span> <span>&lt;-</span> <span>rmvn</span><span>(</span><span>N</span><span>,</span> <span>mu</span> <span>=</span> <span>coef</span><span>(</span><span>m</span><span>),</span> <span>sig</span> <span>=</span> <span>Vb</span><span>)</span>
<span>fits</span> <span>&lt;-</span> <span>Cg</span> <span>%*%</span> <span>t</span><span>(</span><span>sims</span><span>)</span></code></pre>
</figure>
<p>
<code>fits</code> now contains N = 10<sup>4</sup> draws from the model posterior. Before we look at how many of the 10<sup>4</sup> samples from the posterior are entirely contained within the simultaneous interval, choose 30 at random and stack them in so-called tidy form for use with <code>ggplot()</code>
</p>
<figure>
<pre><code data-lang="r"><span>nrnd</span> <span>&lt;-</span> <span>30</span>
<span>rnd</span> <span>&lt;-</span> <span>sample</span><span>(</span><span>N</span><span>,</span> <span>nrnd</span><span>)</span>
<span>stackFits</span> <span>&lt;-</span> <span>stack</span><span>(</span><span>as.data.frame</span><span>(</span><span>fits</span><span>[,</span> <span>rnd</span><span>]))</span>
<span>stackFits</span> <span>&lt;-</span> <span>transform</span><span>(</span><span>stackFits</span><span>,</span> <span>age</span> <span>=</span> <span>rep</span><span>(</span><span>newd</span><span>$</span><span>age</span><span>,</span> <span>length</span><span>(</span><span>rnd</span><span>)))</span></code></pre>
</figure>
<p>
What we’ve done in this post can be summarized in the figure below
</p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>(</span><span>pred</span><span>,</span> <span>aes</span><span>(</span><span>x</span> <span>=</span> <span>age</span><span>,</span> <span>y</span> <span>=</span> <span>fit</span><span>))</span> <span>+</span>
    <span>geom_ribbon</span><span>(</span><span>aes</span><span>(</span><span>ymin</span> <span>=</span> <span>lwrS</span><span>,</span> <span>ymax</span> <span>=</span> <span>uprS</span><span>),</span> <span>alpha</span> <span>=</span> <span>0.2</span><span>,</span> <span>fill</span> <span>=</span> <span>&quot;red&quot;</span><span>)</span> <span>+</span>
    <span>geom_ribbon</span><span>(</span><span>aes</span><span>(</span><span>ymin</span> <span>=</span> <span>lwrP</span><span>,</span> <span>ymax</span> <span>=</span> <span>uprP</span><span>),</span> <span>alpha</span> <span>=</span> <span>0.2</span><span>,</span> <span>fill</span> <span>=</span> <span>&quot;red&quot;</span><span>)</span> <span>+</span>
    <span>geom_path</span><span>(</span><span>lwd</span> <span>=</span> <span>2</span><span>)</span> <span>+</span>
    <span>geom_path</span><span>(</span><span>data</span> <span>=</span> <span>stackFits</span><span>,</span> <span>mapping</span> <span>=</span> <span>aes</span><span>(</span><span>y</span> <span>=</span> <span>values</span><span>,</span> <span>x</span> <span>=</span> <span>age</span><span>,</span> <span>group</span> <span>=</span> <span>ind</span><span>),</span>
              <span>alpha</span> <span>=</span> <span>0.4</span><span>,</span> <span>colour</span> <span>=</span> <span>&quot;grey20&quot;</span><span>)</span> <span>+</span>
    <span>labs</span><span>(</span><span>y</span> <span>=</span> <span>&quot;Strontium isotope ratio&quot;</span><span>,</span>
         <span>x</span> <span>=</span> <span>&quot;Age [Ma BP]&quot;</span><span>,</span>
         <span>title</span> <span>=</span> <span>&quot;Point-wise &amp; Simultaneous 95% confidence intervals for fitted GAM&quot;</span><span>,</span>
         <span>subtitle</span> <span>=</span> <span>sprintf</span><span>(</span><span>&quot;Each line is one of %i draws from the Bayesian posterior distribution of the model&quot;</span><span>,</span> <span>nrnd</span><span>))</span></code></pre>
</figure>
<figure>
<img alt="Summary plot showing 30 random draws from the model posterior and approximate 95% simultaneous and point-wise confidence intervals for the the fitted GAM" src="https://i1.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/simultaneous-intervals-revisited-plot-intervals-and-posterior-draws-1.png?w=456" data-recalc-dims="1">
<figcaption>
Summary plot showing 30 random draws from the model posterior and approximate 95% simultaneous and point-wise confidence intervals for the the fitted GAM<br>
</figcaption>
</figure>
<p>
It shows the fitted model and the 95% simultaneous and point-wise confidence intervals, and is augmented with 30 draws from the posterior distribution of the GAM. As you can see, many of the lines lie outside the point-wise confidence interval. The situation is quite different with the simultaneous interval; only a couple of the posterior draws go outside of the 95% simultaneous interval, which is what we’d expect for a 95% interval. So that’s encouraging!
</p>
<p>
As a final check we’ll look at the proportion of all the posterior simulations that lie entirely within the simultaneous interval. To facilitate this we create a little wrapper function, <code>inCI()</code>, which returns <code>TRUE</code> if all the evaluation points <span>()</span> lie within the stated interval and <code>FALSE</code> otherwise. This is then applied to each posterior simulation (column of <code>fits</code>) and we do this for the simultaneous intervals and the point-wise version. The final two lines work out what proportion of the posterior simulations lie within the two confidence intervals.
</p>
<figure>
<pre><code data-lang="r"><span>inCI</span> <span>&lt;-</span> <span>function</span><span>(</span><span>x</span><span>,</span> <span>upr</span><span>,</span> <span>lwr</span><span>)</span> <span>{</span>
    <span>all</span><span>(</span><span>x</span> <span>>=</span> <span>lwr</span> <span>&amp;</span> <span>x</span> <span>&lt;=</span> <span>upr</span><span>)</span>
<span>}</span>

<span>fitsInPCI</span> <span>&lt;-</span> <span>apply</span><span>(</span><span>fits</span><span>,</span> <span>2L</span><span>,</span> <span>inCI</span><span>,</span> <span>upr</span> <span>=</span> <span>pred</span><span>$</span><span>uprP</span><span>,</span> <span>lwr</span> <span>=</span> <span>pred</span><span>$</span><span>lwrP</span><span>)</span>
<span>fitsInSCI</span> <span>&lt;-</span> <span>apply</span><span>(</span><span>fits</span><span>,</span> <span>2L</span><span>,</span> <span>inCI</span><span>,</span> <span>upr</span> <span>=</span> <span>pred</span><span>$</span><span>uprS</span><span>,</span> <span>lwr</span> <span>=</span> <span>pred</span><span>$</span><span>lwrS</span><span>)</span>

<span>sum</span><span>(</span><span>fitsInPCI</span><span>)</span> <span>/</span> <span>length</span><span>(</span><span>fitsInPCI</span><span>)</span>      <span># Point-wise
</span><span>sum</span><span>(</span><span>fitsInSCI</span><span>)</span> <span>/</span> <span>length</span><span>(</span><span>fitsInSCI</span><span>)</span>      <span>#</span> <span>Simultaneous</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">[1] 0.3028
[1] 0.9526</code></pre>
</figure>
<p>
As you can see, the point-wise confidence interval includes just a small proportion of the posterior simulations, but the simultaneous interval contains approximately the right number of simulations for a 95% interval.
</p>
<p>
So how bad are the intervals I created in the old post? They should as bad as the 95% point-wise interval, and they are
</p>
<figure>
<pre><code data-lang="r"><span>oldCI</span> <span>&lt;-</span> <span>apply</span><span>(</span><span>fits</span><span>,</span> <span>1L</span><span>,</span> <span>quantile</span><span>,</span> <span>probs</span> <span>=</span> <span>c</span><span>(</span><span>0.025</span><span>,</span> <span>0.975</span><span>))</span>
<span>pred</span> <span>&lt;-</span> <span>transform</span><span>(</span><span>pred</span><span>,</span> <span>lwrOld</span> <span>=</span> <span>oldCI</span><span>[</span><span>1</span><span>,</span> <span>],</span> <span>uprOld</span> <span>=</span> <span>oldCI</span><span>[</span><span>2</span><span>,</span> <span>])</span>
<span>fitsInOldCI</span> <span>&lt;-</span> <span>apply</span><span>(</span><span>fits</span><span>,</span> <span>2L</span><span>,</span> <span>inCI</span><span>,</span> <span>upr</span> <span>=</span> <span>pred</span><span>$</span><span>uprOld</span><span>,</span> <span>lwr</span> <span>=</span> <span>pred</span><span>$</span><span>lwrOld</span><span>)</span>
<span>sum</span><span>(</span><span>fitsInOldCI</span><span>)</span> <span>/</span> <span>length</span><span>(</span><span>fitsInOldCI</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">[1] 0.2655</code></pre>
</figure>
<p>
So, there we have it — a proper 95% simultaneous confidence interval for a penalised spline. Now I just need to go back to that old post and strike out all reference to <em>simultaneous</em>…
</p>
<div>
<h3>
References<br>
</h3>
<div>
<p>
Marra, Giampiero, and Simon N Wood. 2012. “Coverage Properties of Confidence Intervals for Generalized Additive Model Components.” <em>Scandinavian Journal of Statistics, Theory and Applications</em> 39 (1). Blackwell Publishing Ltd: 53–74. doi:<a target="_blank" rel="nofollow" href="http://doi.org/10.1111/j.1467-9469.2011.00760.x">http://doi.org/10.1111/j.1467-9469.2011.00760.x</a>.
</p>
</div>
<div>
<p>
Nychka, Douglas. 1988. “Bayesian Confidence Intervals for Smoothing Splines.” <em>Journal of the American Statistical Association</em> 83 (404): 1134–43. doi:<a target="_blank" rel="nofollow" href="http://doi.org/10.1080/01621459.1988.10478711">http://doi.org/10.1080/01621459.1988.10478711</a>.
</p>
</div>
<div>
<p>
Ruppert, David, M P Wand, and R J Carroll. 2003. <em>Semiparametric Regression</em>. Cambridge University Press.
</p>
</div>
</div>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.fromthebottomoftheheap.net/2016/12/15/simultaneous-interval-revisited/"> From the Bottom of the Heap - R</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/JY7z88OZjWc"></a></div>