Simultaneous confidence intervals for derivatives of splines in GAMs
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.
Last time out I looked at one of the complications of time series modelling with smoothers; you have a nonlinear trend which may be statistically significant but it may not be increasing or decreasing everywhere. How do we identify where in the series the data are changing? In that post I explained how we can use the first derivatives of the model splines for this purpose, and used the method of finite differences to estimate them. To assess statistical significance of the derivative (the rate of change) I relied upon asymptotic normality and the usual pointwise confidence interval. That interval is fine if looking at just one point on the spline (not of much practical use), but when considering more points at once we have a multiple comparisons issue. Instead, a simultaneous interval is required, and for that we need to revisit a technique I blogged about a few years ago; posterior simulation from the fitted GAM.
To get a headstart on this, I’ll reuse the model we fitted to the CET time series from the previous post. Just copy and paste the code below into your R session
<span class="c1">## Load the CET data and process as per other blog post</span>
tmpf <span class="o"><</span> tempfile<span class="p">()</span>
download.file<span class="p">(</span><span class="s">"https://gist.github.com/gavinsimpson/b52f6d375f57d539818b/raw/2978362d97ee5cc9e7696d2f36f94762554eefdf/loadprocesscetmonthly.R"</span><span class="p">,</span>
tmpf<span class="p">,</span> method <span class="o">=</span> <span class="s">"wget"</span><span class="p">)</span>
source<span class="p">(</span>tmpf<span class="p">)</span>
<span class="c1">## Load mgcv and fit the model</span>
require<span class="p">(</span><span class="s">"mgcv"</span><span class="p">)</span>
ctrl <span class="o"><</span> list<span class="p">(</span>niterEM <span class="o">=</span> <span class="m">0</span><span class="p">,</span> msVerbose <span class="o">=</span> <span class="kc">TRUE</span><span class="p">,</span> optimMethod<span class="o">=</span><span class="s">"LBFGSB"</span><span class="p">)</span>
m2 <span class="o"><</span> gamm<span class="p">(</span>Temperature <span class="o">~</span> s<span class="p">(</span>nMonth<span class="p">,</span> bs <span class="o">=</span> <span class="s">"cc"</span><span class="p">,</span> k <span class="o">=</span> <span class="m">12</span><span class="p">)</span> <span class="o">+</span> s<span class="p">(</span>Time<span class="p">,</span> k <span class="o">=</span> <span class="m">20</span><span class="p">),</span>
data <span class="o">=</span> cet<span class="p">,</span> correlation <span class="o">=</span> corARMA<span class="p">(</span>form <span class="o">=</span> <span class="o">~</span> <span class="m">1</span><span class="o"></span>Year<span class="p">,</span> p <span class="o">=</span> <span class="m">2</span><span class="p">),</span>
control <span class="o">=</span> ctrl<span class="p">)</span>
<span class="c1">## prediction data</span>
want <span class="o"><</span> seq<span class="p">(</span><span class="m">1</span><span class="p">,</span> nrow<span class="p">(</span>cet<span class="p">),</span> length.out <span class="o">=</span> <span class="m">200</span><span class="p">)</span>
pdat <span class="o"><</span> with<span class="p">(</span>cet<span class="p">,</span>
data.frame<span class="p">(</span>Time <span class="o">=</span> Time<span class="p">[</span>want<span class="p">],</span> Date <span class="o">=</span> Date<span class="p">[</span>want<span class="p">],</span>
nMonth <span class="o">=</span> nMonth<span class="p">[</span>want<span class="p">]))</span>
Here, I’ll use a version of the Deriv()
function used in the last post modified to do the posterior simulation; derivSimulCI()
. Let’s load that too
<span class="c1">## download the derivatives gist</span>
tmpf <span class="o"><</span> tempfile<span class="p">()</span>
download.file<span class="p">(</span><span class="s">"https://gist.githubusercontent.com/gavinsimpson/ca18c9c789ef5237dbc6/raw/295fc5cf7366c831ab166efaee42093a80622fa8/derivSimulCI.R"</span><span class="p">,</span>
tmpf<span class="p">,</span> method <span class="o">=</span> <span class="s">"wget"</span><span class="p">)</span>
source<span class="p">(</span>tmpf<span class="p">)</span>
Posterior simulation
The sorts of GAMs fitted by mgcv::gam()
are, if we assume normally distributed errors, really just a linear regression. Instead of being a linear model in the original data however, the linear model is fitted using the basis functions as the covariates^{1}. As with any other linear model, we get back from it the point estimate, the ( _j ), and their standard errors. Consider the simple linear regression of x on y. Such a model has two terms
 the constant term (the model intercept), and
 the effect on y of a unit change in x.
In fitting the model we get a point estimate for each term, plus their standard errors in the form of the variancecovariance (VCOV) matrix of the terms^{2}. Taken together, the point estimates of the model terms and the VCOV describe a multivariate normal distribution. In the case of the simple linear regression, this is a bivariate normal. Note that the point estimates are known as the mean vector of the multivariate normal; each point estimate is the mean, or expectation, of a single random normal variable whose variance is given by the standard error of the point estimate.
Computers are good at simulating data and you’ll most likely be familiar with rnorm()
to generate random, normally distributed values from a distribution with mean 0 and unit standard deviation. Well, simulating from a multivariate normal is just as simple^{3}, as long as you have the mean vector and the variance covariance matrix of the parameters.
Returning to the simple linear regression case, let’s do a little simulation from a known model and look at the multivariate normal distribution of the model parameters.
<span class="o">></span> set.seed<span class="p">(</span><span class="m">1</span><span class="p">)</span>
<span class="o">></span> N <span class="o"><</span> <span class="m">100</span>
<span class="o">></span> dat <span class="o"><</span> data.frame<span class="p">(</span>x <span class="o">=</span> runif<span class="p">(</span>N<span class="p">,</span> min <span class="o">=</span> <span class="m">1</span><span class="p">,</span> max <span class="o">=</span> <span class="m">20</span><span class="p">))</span>
<span class="o">></span> dat <span class="o"><</span> transform<span class="p">(</span>dat<span class="p">,</span> y <span class="o">=</span> <span class="m">3</span> <span class="o">+</span> <span class="p">(</span><span class="m">1.45</span> <span class="o">*</span> x<span class="p">)</span> <span class="o">+</span> rnorm<span class="p">(</span>N<span class="p">,</span> mean <span class="o">=</span> <span class="m">2</span><span class="p">,</span> sd <span class="o">=</span> <span class="m">3</span><span class="p">))</span>
<span class="o">></span> <span class="c1">## sort dat on x to make things easier later</span>
<span class="o">></span> dat <span class="o"><</span> dat<span class="p">[</span>order<span class="p">(</span>dat<span class="o">$</span>x<span class="p">),</span> <span class="p">]</span>
<span class="o">></span> mod <span class="o"><</span> lm<span class="p">(</span>y <span class="o">~</span> x<span class="p">,</span> data <span class="o">=</span> dat<span class="p">)</span>
The mean vector for the multivariate normal is just the set of model coefficients for mod
, which are extracted using the coef()
function, and the vcov()
function is used to extract the VCOV of the fitted model.
<span class="o">></span> coef<span class="p">(</span>mod<span class="p">)</span>
(Intercept) x
4.413 1.499
<span class="o">></span> <span class="p">(</span>vc <span class="o"><</span> vcov<span class="p">(</span>mod<span class="p">))</span>
(Intercept) x
(Intercept) 0.44563 0.033760
x 0.03376 0.003115
Remember, the standard error is the square root of the diagonal elements of the VCOV
<span class="o">></span> coef<span class="p">(</span>summary<span class="p">(</span>mod<span class="p">))[,</span> <span class="s">"Std. Error"</span><span class="p">]</span>
(Intercept) x
0.66756 0.05581
<span class="o">></span> sqrt<span class="p">(</span>diag<span class="p">(</span>vc<span class="p">))</span>
(Intercept) x
0.66756 0.05581
The multivariate normal distribution is not part of the base R distributions set. Several implementations are available in a range of packages, but here I’ll use the one in the MASS package which ships with all versions of R. To draw a nice plot, I’ll simulate a large number of values but we’ll just show the first few below
<span class="o">></span> require<span class="p">(</span><span class="s">"MASS"</span><span class="p">)</span>
<span class="o">></span> set.seed<span class="p">(</span><span class="m">10</span><span class="p">)</span>
<span class="o">></span> nsim <span class="o"><</span> <span class="m">5000</span>
<span class="o">></span> sim <span class="o"><</span> mvrnorm<span class="p">(</span>nsim<span class="p">,</span> mu <span class="o">=</span> coef<span class="p">(</span>mod<span class="p">),</span> Sigma <span class="o">=</span> vc<span class="p">)</span>
<span class="o">></span> head<span class="p">(</span>sim<span class="p">)</span>
(Intercept) x
[1,] 4.398 1.477
[2,] 4.536 1.496
[3,] 5.328 1.427
[4,] 4.811 1.446
[5,] 4.216 1.510
[6,] 4.153 1.528
Each row of sim
contains a pair of values, one intercept and one ( _x ), from the implied multivariate normal. The models implied by each row are all consistent with the fitted model. To visualize the multivariate normal for mod
I’ll use a bivariate kernel density estimate to estimate the density of points over a grid of simulated intercept and slope values
<span class="o">></span> kde <span class="o"><</span> kde2d<span class="p">(</span>sim<span class="p">[,</span><span class="m">1</span><span class="p">],</span> sim<span class="p">[,</span><span class="m">2</span><span class="p">],</span> n <span class="o">=</span> <span class="m">75</span><span class="p">)</span>
<span class="o">></span> plot<span class="p">(</span>sim<span class="p">,</span> pch <span class="o">=</span> <span class="m">19</span><span class="p">,</span> cex <span class="o">=</span> <span class="m">0.5</span><span class="p">,</span> col <span class="o">=</span> <span class="s">"darkgrey"</span><span class="p">)</span>
<span class="o">></span> contour<span class="p">(</span>kde<span class="o">$</span>x<span class="p">,</span> kde<span class="o">$</span>y<span class="p">,</span> kde<span class="o">$</span>z<span class="p">,</span> add <span class="o">=</span> <span class="kc">TRUE</span><span class="p">,</span> col <span class="o">=</span> <span class="s">"red"</span><span class="p">,</span> lwd <span class="o">=</span> <span class="m">2</span><span class="p">,</span> drawlabels <span class="o">=</span> <span class="kc">FALSE</span><span class="p">)</span>
The large spread in the points (from top left to bottom right) is illustrative of greater uncertainty in the intercept term than in ( _x ).
As I said earlier, each point on the plot represents a valid model consistent with the estimates we achieved for the sample of data used to fit the model. If we were to multiple the second column of sim
with the observed data and add on the first column of sim
, we’d obtain fitted values for the observed x
values for 5000 simulations from the fitted model as shown in the plot below
<span class="o">></span> plot<span class="p">(</span>y <span class="o">~</span> x<span class="p">,</span> data <span class="o">=</span> dat<span class="p">)</span>
<span class="o">></span> set.seed<span class="p">(</span><span class="m">42</span><span class="p">)</span>
<span class="o">></span> take <span class="o"><</span> sample<span class="p">(</span>nrow<span class="p">(</span>sim<span class="p">),</span> <span class="m">50</span><span class="p">)</span> <span class="c1">## take 50 simulations at random</span>
<span class="o">></span> fits <span class="o"><</span> cbind<span class="p">(</span><span class="m">1</span><span class="p">,</span> dat<span class="o">$</span>x<span class="p">)</span> <span class="o">%*%</span> t<span class="p">(</span>sim<span class="p">[</span>take<span class="p">,</span> <span class="p">])</span>
<span class="o">></span> matlines<span class="p">(</span>dat<span class="o">$</span>x<span class="p">,</span> fits<span class="p">,</span> col <span class="o">=</span> <span class="s">"#A9A9A97D"</span><span class="p">,</span> lty <span class="o">=</span> <span class="s">"solid"</span><span class="p">,</span> lwd <span class="o">=</span> <span class="m">2</span><span class="p">)</span>
<span class="o">></span> abline<span class="p">(</span>mod<span class="p">,</span> col <span class="o">=</span> <span class="s">"red"</span><span class="p">,</span> lwd <span class="o">=</span> <span class="m">1</span><span class="p">)</span>
<span class="o">></span> matlines<span class="p">(</span>dat<span class="o">$</span>x<span class="p">,</span> predict<span class="p">(</span>mod<span class="p">,</span> interval <span class="o">=</span> <span class="s">"confidence"</span><span class="p">)[,</span><span class="m">1</span><span class="p">],</span> col <span class="o">=</span> <span class="s">"red"</span><span class="p">,</span> lty <span class="o">=</span> <span class="s">"dashed"</span><span class="p">)</span>
The grey lines show the model fits for a random sample of 50 pairs of coefficients from the set of simulated values.
Posterior simulation for additive models
You’ll be pleased to know that there is very little difference (non really) between what I just went through above for a simple linear regression and what is required to simulate from the posterior distribution of a GAM. However, instead of dealing with two or just a few regression coefficients, we now have to concern ourselves with the potentially larger number of coefficients corresponding to the basis functions that combine to form the fitted splines. The only practical difference is that instead of multiplying each simulation by the observed data^{4} with mgcv we generate the linear predictor matrix for the observations and multiply that by the model coefficients to get simulations. If you’ve read the previous post you should be somewhat familiar with the lpmatrix
now.
Before we get to posterior simulations for the derivatives of the CET additive model fitted earlier, let’s look at some simulations for the trend term in that model, m2
. If you look back at an earlier code block, I created a grid of 200 points over the range of the data which we’ll use to evaluate properties of the fitted model. This is in object pdat
. First we generate the linear predictor matrix using predict()
and grab the model coefficients and the variance covariance matrix of the coefficients
<span class="o">></span> lp <span class="o"><</span> predict<span class="p">(</span>m2<span class="o">$</span>gam<span class="p">,</span> newdata <span class="o">=</span> pdat<span class="p">,</span> type <span class="o">=</span> <span class="s">"lpmatrix"</span><span class="p">)</span>
<span class="o">></span> coefs <span class="o"><</span> coef<span class="p">(</span>m2<span class="o">$</span>gam<span class="p">)</span>
<span class="o">></span> vc <span class="o"><</span> vcov<span class="p">(</span>m2<span class="o">$</span>gam<span class="p">)</span>
Next, generate a small sample from the posterior of the model, just for the purposes of illustration; we’ll generate far larger samples later when we estimate a confidence interval on the derivatives of the trend spline.
<span class="o">></span> set.seed<span class="p">(</span><span class="m">35</span><span class="p">)</span>
<span class="o">></span> sim <span class="o"><</span> mvrnorm<span class="p">(</span><span class="m">25</span><span class="p">,</span> mu <span class="o">=</span> coefs<span class="p">,</span> Sigma <span class="o">=</span> vc<span class="p">)</span>
The linear predictor matrix, lp
, has a column for every basis function, plus the constant term, in the model, but because the model is additive we can ignore the columns relating to the nMonth
spline and the constant term and just work with the coefficients and columns of lp
that pertain to the trend spline. Let’s identify those
<span class="o">></span> want <span class="o"><</span> grep<span class="p">(</span><span class="s">"Time"</span><span class="p">,</span> colnames<span class="p">(</span>lp<span class="p">))</span>
Again, a simple bit of matrix multiplication gets us fitted values for the trend spline only
<span class="o">></span> fits <span class="o"><</span> lp<span class="p">[,</span> want<span class="p">]</span> <span class="o">%*%</span> t<span class="p">(</span>sim<span class="p">[,</span> want<span class="p">])</span>
<span class="o">></span> dim<span class="p">(</span>fits<span class="p">)</span> <span class="c1">## 25 columns, 1 per simulation, 200 rows, 1 per evaln point</span>
[1] 200 25
We can now draw out each of these posterior simulations as follows
<span class="o">></span> ylims <span class="o"><</span> range<span class="p">(</span>fits<span class="p">)</span>
<span class="o">></span> plot<span class="p">(</span>Temperature <span class="o">~</span> Date<span class="p">,</span> data <span class="o">=</span> cet<span class="p">,</span> pch <span class="o">=</span> <span class="m">19</span><span class="p">,</span> ylim <span class="o">=</span> ylims<span class="p">,</span> type <span class="o">=</span> <span class="s">"n"</span><span class="p">)</span>
<span class="o">></span> matlines<span class="p">(</span>pdat<span class="o">$</span>Date<span class="p">,</span> fits<span class="p">,</span> col <span class="o">=</span> <span class="s">"black"</span><span class="p">,</span> lty <span class="o">=</span> <span class="s">"solid"</span><span class="p">)</span>
Posterior simulation for the first derivatives of a spline
As we saw in the previous post, the linear predictor matrix can be used to generate finite differencesbased estimates of the derivatives of a spline in a GAM fitted by mgcv. And as we just went through, we can combine posterior simulations with the linear predictor matrix. The main steps in the process of computing the finite differences and doing the posterior simulation are
X0 <span class="o"><</span> predict<span class="p">(</span>mod<span class="p">,</span> newDF<span class="p">,</span> type <span class="o">=</span> <span class="s">"lpmatrix"</span><span class="p">)</span>
newDF <span class="o"><</span> newDF <span class="o">+</span> eps
X1 <span class="o"><</span> predict<span class="p">(</span>mod<span class="p">,</span> newDF<span class="p">,</span> type <span class="o">=</span> <span class="s">"lpmatrix"</span><span class="p">)</span>
Xp <span class="o"><</span> <span class="p">(</span>X1 <span class="o"></span> X0<span class="p">)</span> <span class="o">/</span> eps
where two linear predictor matrices are created, offset from one another by a small amount eps
, and differenced to get the slope of the spline, and
<span class="kr">for</span><span class="p">(</span>i <span class="kr">in</span> seq_len<span class="p">(</span>nt<span class="p">))</span> <span class="p">{</span>
Xi <span class="o"><</span> Xp <span class="o">*</span> <span class="m">0</span>
want <span class="o"><</span> grep<span class="p">(</span>t.labs<span class="p">[</span>i<span class="p">],</span> colnames<span class="p">(</span>X1<span class="p">))</span>
Xi<span class="p">[,</span> want<span class="p">]</span> <span class="o"><</span> Xp<span class="p">[,</span> want<span class="p">]</span>
df <span class="o"><</span> Xi <span class="o">%*%</span> t<span class="p">(</span>simu<span class="p">[,</span> want<span class="p">])</span> <span class="c1"># derivatives</span>
<span class="p">}</span>
which loops over the terms in the model, selects the relevant columns from the differenced predictor matrix, and computes the derivatives by a matrix multiplication with the set of posterior simulations. simu
is the matrix of random draws from the posterior, multivariate normal distribution of the fitted model’s parameters. Note that the code in derivSimulCI()
is slightly different to this, but it does the same thing.
To cut to the chase then, here is the code required to generate posterior simulations for the first derivatives of the spline terms in an additive model
<span class="o">></span> fd <span class="o"><</span> derivSimulCI<span class="p">(</span>m2<span class="p">,</span> samples <span class="o">=</span> <span class="m">10000</span><span class="p">)</span>
fd
is a list, the first n terms of which relate the the n terms in the model. Here n = 2. The names of the first two components are the names of the terms referenced in the model formula used to fit the model
<span class="o">></span> str<span class="p">(</span>fd<span class="p">,</span> max <span class="o">=</span> <span class="m">1</span><span class="p">)</span>
List of 5
$ nMonth :List of 2
$ Time :List of 2
$ gamModel:List of 31
.. attr(*, "class")= chr "gam"
$ eps : num 1e07
$ eval : num [1:200, 1:2] 1 1.06 1.11 1.17 1.22 ...
.. attr(*, "dimnames")=List of 2
 attr(*, "class")= chr "derivSimulCI"
As I haven’t yet written a confint()
method, we’ll need to compute the confidence interval by hand, which is no bad thing of course! We do this by by taking two extreme quantiles of the distribution of the 10,000 posterior simulations we generated for the first derivative at each of the 200 points we wanted to evaluate the derivative. One of the reasons I did 10,000 simulations is that for a 95% confidence interval we only need sort the simulated derivatives in ascending order and extract the 250th and the 9750th of these ordered values. In practice we’ll let the quantile()
function do the hard work
<span class="o">></span> CI <span class="o"><</span> lapply<span class="p">(</span>fd<span class="p">[</span><span class="m">1</span><span class="o">:</span><span class="m">2</span><span class="p">],</span>
<span class="o">+</span> <span class="kr">function</span><span class="p">(</span>x<span class="p">)</span> apply<span class="p">(</span>x<span class="o">$</span>simulations<span class="p">,</span> <span class="m">1</span><span class="p">,</span> quantile<span class="p">,</span>
<span class="o">+</span> probs <span class="o">=</span> c<span class="p">(</span><span class="m">0.025</span><span class="p">,</span> <span class="m">0.975</span><span class="p">)))</span>
CI
is now a list with two components, each of which contains a matrix with two rows (the two probability quantiles we asked for) and 200 columns (the number of locations at which the first derivative was evaluated).
There is a plot()
method, which by default produces plots of all the terms in the model and includes the simultaneous confidence interval as well
<span class="o">></span> plot<span class="p">(</span>fd<span class="p">,</span> sizer <span class="o">=</span> <span class="kc">TRUE</span><span class="p">)</span>
Wrapping up
derivSimulCI()
computes the actual derivative as well as the derivatives for each simulation. Rather than rely upon the plot()
method we could draw our own plot with the confidence interval. To extract the derivative of the fitted spline use
<span class="o">></span> fit.fd <span class="o"><</span> fd<span class="p">[[</span><span class="m">2</span><span class="p">]]</span><span class="o">$</span>deriv
and then produce a plot with the actual derivative, the 95% simultaneous confidence interval, and 20 of the derivatives for the posterior simulations, we can use
<span class="o">></span> set.seed<span class="p">(</span><span class="m">76</span><span class="p">)</span>
<span class="o">></span> take <span class="o"><</span> sample<span class="p">(</span>nrow<span class="p">(</span>fd<span class="p">[[</span><span class="m">2</span><span class="p">]]</span><span class="o">$</span>simulations<span class="p">),</span> <span class="m">20</span><span class="p">)</span>
<span class="o">></span> plot<span class="p">(</span>pdat<span class="o">$</span>Date<span class="p">,</span> fit.fd<span class="p">,</span> type <span class="o">=</span> <span class="s">"l"</span><span class="p">,</span> ylim <span class="o">=</span> range<span class="p">(</span>CI<span class="p">[[</span><span class="m">2</span><span class="p">]]),</span> lwd <span class="o">=</span> <span class="m">2</span><span class="p">)</span>
<span class="o">></span> matlines<span class="p">(</span>pdat<span class="o">$</span>Date<span class="p">,</span> t<span class="p">(</span>CI<span class="p">[[</span><span class="m">2</span><span class="p">]]),</span> lty <span class="o">=</span> <span class="s">"dashed"</span><span class="p">,</span> col <span class="o">=</span> <span class="s">"red"</span><span class="p">)</span>
<span class="o">></span> matlines<span class="p">(</span>pdat<span class="o">$</span>Date<span class="p">,</span> fd<span class="p">[[</span><span class="m">2</span><span class="p">]]</span><span class="o">$</span>simulations<span class="p">[,</span> take<span class="p">],</span> lty <span class="o">=</span> <span class="s">"solid"</span><span class="p">,</span>
<span class="o">+</span> col <span class="o">=</span> <span class="s">"grey"</span><span class="p">)</span>

It is a little bit more complex than this, of course. If you allow
gam()
to select the degree of smoothness then you need to fit a penalized regression. Plus, the time series models fitted to the CET data aren’t fitted viagam()
but viagamm()
, where we are using the observation that a penalized regression can be expressed as a linear mixed model, with random effects being used to represent some of the penalty terms. If you specify the degree of smoothing to use, these complications go away.↩ 
The (squares of the) standard errors are on the diagonal of the VCOV, with the relationship between pairs of parameters being contained in the offdiagonal elements.↩

In practice. I suspect it is not quite so simple if one had to sit down and implement it….↩

Or a set of new values at which you want to evaluate the confidence interval.↩
Rbloggers.com offers daily email updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/datascience job.
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.