# Restricted cubic splines in Stan?

**URL:** <https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234>\
**Category:** Modeling\
**Created:** [May 17, 2018, 8:10pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234 "2018-05-17T20:10:01Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![africa\_stats](https://avatars.discourse-cdn.com/v4/letter/a/f475e1/32.png) [@africa\_stats](https://discourse.mc-stan.org/u/africa_stats)\
**Post date:** [May 17, 2018, 8:10pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/1 "2018-05-17T20:10:01Z")

</div>

A year ago in the Stan Google, a Marise [asked to check an implementation of restricted cubic splines](https://groups.google.com/forum/#!topic/stan-users/8jjaz1kwKcc), a method made popular (I think) by Frank Harrell. She apparently didn’t do a good job coding this and was told to learn the basics before trying to develop her own implementation, so I wasn’t helped much by the post.

Has this been implemented in Stan since then? My group is looking to transfer some of our R code, that relies on the R implementation of RCS, to Stan, if possible.

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [May 17, 2018, 8:54pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/2 "2018-05-17T20:54:16Z")

</div>

There is no function in the Stan Math library that “does” a restricted cubic spline.  
You can use pretty much any kind of spline supported by the **mgcv** package using `stan_gamm4` in **rstanarm** or `brm` in **brms**. Or there is a case study about splines at  
[http://mc-stan.org/users/documentation/case-studies/splines\_in\_stan.html](http://mc-stan.org/users/documentation/case-studies/splines_in_stan.html)  
that you could adapt to your situation.

---

<div class="post-metadata">

**Author:** ![ermeel](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ermeel/32/4180_2.png) [@ermeel](https://discourse.mc-stan.org/u/ermeel)\
**Post date:** [May 18, 2018, 7:08am UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/3 "2018-05-18T07:08:04Z")

</div>

A while ago I implemented Natural Cubic Splines (and their derivative) in Stan, as defined in [Royston & Parmar](https://www.ncbi.nlm.nih.gov/pubmed/12210632).

Haven’t tested it very thoroughly, but it could be a starting point. Any feedback is highly welcome. Below comes a snippet that should be enough to get started…

```stan
functions {
   /**
    * Return natural cubic spline basis function at a specific knot evaluated
    * at a set of x-values
    *
    * @param xs An array of reals of x-values for which the basis function should be evaluated
    * @param l A real denoting \lambda_j for the j'th basis function in Royston & Parmar
    * @param k A real denoting k_j, i.e. the position of the j'th knot
    * @param kmax A real denoting the position of the right boundary knot
    * @param kmin A real denoting the position of the left boundary knot
    * @return A vector corresponding to the values of the spline basis function evaluated at x
    */
    vector basis_function_v(real[] xs, real l, real k, real kmax, real kmin) {
        vector[size(xs)] vs;
        for(i in 1:size(xs)) 
            vs[i] = pow(max({0, xs[i] - k}),3) - l*pow(max({0, xs[i] - kmin}),3) - (1-l)*pow(max({0, xs[i] - kmax}),3);
        return vs;
    }
   /**
    * Return the derivative of the natural cubic spline basis function at a specific knot evaluated
    * at a set of x-values
    *
    * @param xs An array of reals of x-values for which the basis function should be evaluated
    * @param l A real denoting \lambda_j for the j'th basis function in Royston & Parmar
    * @param k A real denoting k_j, i.e. the position of the j'th knot
    * @param kmax A real denoting the position of the right boundary knot
    * @param kmin A real denoting the position of the left boundary knot
    * @return A vector corresponding to the values of the derivative of the spline basis function evaluated at x
    */
    vector deriv_basis_function_v(real[] xs, real l, real k, real kmax, real kmin) {
        vector[size(xs)] vs;
        for(i in 1:size(xs)) 
            vs[i] = 3*pow(max({0, xs[i] - k}),2) - 3*l*pow(max({0, xs[i] - kmin}),2) - 3*(1-l)*pow(max({0, xs[i] - kmax}),2);
        return vs;
    }
}
/ ************************************************************************************************************************ /
data {
    int<lower=1> N; // number of data points
    int<lower=0> m; // number of internal knots
    ordered[m+2] knots; // location of knots
    int<lower=1> NC; // number of covariates
    matrix[N,NC] X; // design matrix
    int<lower=0, upper=1> is_censored[N]; // delta in the paper
    vector<lower=0>[N] times; // t in the paper
}
/ ************************************************************************************************************************ /
transformed data {
    vector[N] log_times; // x in the paper
    vector[m] lambdas;
    matrix[m+1,N] basis_evals;
    matrix[m,N] deriv_basis_evals;

    log_times = log(times);
    for(i in 1:m)
        lambdas[i] = knots[m+2] - knots[1+i];
    lambdas /= (knots[m+2] - knots[1]);

    
    basis_evals[1,] = log_times'; // relative to coefficient \gamma_1
    print(basis_evals[1,1], ",", log_times[1],",", basis_evals[1,2],",", log_times[2]);
    for(j in 1:m) {
        basis_evals[j+1,] = basis_function_v(to_array_1d(log_times), 
        lambdas[j], knots[j+1], knots[m+2], knots[1])';
        
        deriv_basis_evals[j,] = deriv_basis_function_v(to_array_1d(log_times),
        lambdas[j], knots[j+1], knots[m+2], knots[1])';
    }  
}

```

---

<div class="post-metadata">

**Author:** ![harrelfe](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/harrelfe/32/2430_2.png) [@harrelfe](https://discourse.mc-stan.org/u/harrelfe)\
**Post date:** [May 18, 2018, 11:21am UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/4 "2018-05-18T11:21:07Z")

</div>

I’d like to learn what is the best medium for adding new easy-to-use spline fitting capabilities, i.e. should this become part of `brms`, `rstanarm`, or base `stan`?

I’ve liked the restricted cubic spline function because it’s easy to implement and you obtain a simple prediction equation at the end. WIth the growing use of Bayes and thinking more about entire posterior distributions and not just point estimates I wonder if the value of having a simply stated point estimate is still as high.

---

<div class="post-metadata">

**Author:** ![ermeel](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ermeel/32/4180_2.png) [@ermeel](https://discourse.mc-stan.org/u/ermeel)\
**Post date:** [May 18, 2018, 3:01pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/5 "2018-05-18T15:01:03Z")

</div>

> [@harrelfe](#):
>
> I’d like to learn what is the best medium for adding new easy-to-use spline fitting capabilities, i.e. should this become part of `brms` , `rstanarm` , or base `stan` ?

What precisely do you mean?

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [May 18, 2018, 4:03pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/6 "2018-05-18T16:03:46Z")

</div>

@harrelfe A case could be made for including spline related functions in the Stan Math C++ library, so that they could be used by all interfaces to Stan. This is somewhat involved but there plenty of people who are capable of implementing it. The main reason why I think no one has done anything on this front is that there are so many types of splines it is not obvious which ones should be implemented and with what priority.

What **brms** and **rstanarm** do is to use functions in the **mgcv** R package to create design matrices that correspond to various types of splines. Then we pass those design matrices to the Stan programs and it proceeds just like a GLM that has a non-spline design matrix. This is easier (for us) because it pushes the responsibility for the implementation onto **mgcv** , but its maintainer Simon Wood has been doing this for decades and has two editions of a textbook about it.

So, I guess for me the question is: What (if anything) to you find to be lacking about how **mgcv** does splines? If there is something, then maybe it would be easiest to add it to **mgcv**. Or is there some function in another R package that creates a design matrix for a type of spline that **brms** or **rstanarm** should have an option to call?

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [May 18, 2018, 4:09pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/7 "2018-05-18T16:09:35Z")

</div>

Seconding what Ben said but also: to get adoption outside of the R universe I’m supportive of getting (the right) splines implemented in the Stan Math C++ library. Some arguments from your end about which splines to implement (whether they’re in mgcv or not) would be helpful.

---

<div class="post-metadata">

**Author:** ![harrelfe](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/harrelfe/32/2430_2.png) [@harrelfe](https://discourse.mc-stan.org/u/harrelfe)\
**Post date:** [May 19, 2018, 2:10pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/8 "2018-05-19T14:10:40Z")

</div>

I think that implementing a variety of splines deep within `stan` is a good idea, for efficiency and precision.

I have nothing against using `mgcv` as long as I can fit a natural spline (restricted cubic spline). How does one obtain the design matrix or such a spline function from `mgcv` without using the package’s `gam` function? This is as far as I’ve gotten, but I’m not seeing how to get the natural spline—I’m getting 5 columns in the design matrix and should have 4.

> require(mgcv)  
> set.seed(1)  
> n ← 30  
> x ← exp(rnorm(n))  
> kn ← quantile(x, c(.05, .2, .5, .8, .95))
> 
> z ← s(x, bs=‘cr’, fx=TRUE, k=5, m=2)  
> X ← smooth.construct(z, data=data.frame(x), knots=list(x=kn))  
> X$X

If someone can answer a background question it would shed some light on the path forward: the easy-to-use truncated power basis for the linear tail-restricted cubic spline function has terms in the design matrix that can be very co-linear. Does this create any convergence problems in Stan? If so I’d be quicker to abandon this basis and would spend time writing back-ends that take Stan results and rephrases the regression model using the ordinary cubic spline notation, including linearly redundant terms as needed to bring about the linear tail constraints.

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [May 19, 2018, 3:01pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/9 "2018-05-19T15:01:19Z")

</div>

> [@harrelfe](#):
>
> If someone can answer a background question it would shed some light on the path forward: the easy-to-use truncated power basis for the linear tail-restricted cubic spline function has terms in the design matrix that can be very co-linear. Does this create any convergence problems in Stan?

Whether it creates problems in a simple model or not it does make sampling harder for the algorithm so if other parts of the model introduce further problems it might make the overall model fail whereas a less colinear spline implementation would work fine. Model problems like collinearity add up in a non-linear fashion so I don’t know if you’d be able to tell ahead of time in general.

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [May 19, 2018, 5:45pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/10 "2018-05-19T17:45:30Z")

</div>

Could something like the QR decomposition strategy for design matrices with collinearity problems also work in this case?

- [http://mc-stan.org/users/documentation/case-studies/qr\_regression.html](http://mc-stan.org/users/documentation/case-studies/qr_regression.html)

- [http://mc-stan.org/rstanarm/reference/QR-argument.html](http://mc-stan.org/rstanarm/reference/QR-argument.html)

---

<div class="post-metadata">

**Author:** ![harrelfe](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/harrelfe/32/2430_2.png) [@harrelfe](https://discourse.mc-stan.org/u/harrelfe)\
**Post date:** [May 19, 2018, 6:06pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/11 "2018-05-19T18:06:30Z")

</div>

I think that can definitely help. The user just has to transform back and forth, and prior specification can be tricky as mentioned in the first link you provided.

---

<div class="post-metadata">

**Author:** ![ermeel](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ermeel/32/4180_2.png) [@ermeel](https://discourse.mc-stan.org/u/ermeel)\
**Post date:** [May 19, 2018, 8:37pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/12 "2018-05-19T20:37:35Z")

</div>

> [@harrelfe](#):
>
> but I’m not seeing how to get the natural spline

For natural cubic splines (aka restricted cubic splines) I used the `ns` method from the core R splines package e.g. for a recent Royston&Parmar implementation):

> **[ns function - RDocumentation](https://www.rdocumentation.org/packages/splines/versions/3.5.0/topics/ns)**
>
> \<p\>Generate the B-spline basis matrix for a natural cubic spline.\</p\>

Would that work?

---

<div class="post-metadata">

**Author:** ![harrelfe](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/harrelfe/32/2430_2.png) [@harrelfe](https://discourse.mc-stan.org/u/harrelfe)\
**Post date:** [May 19, 2018, 9:15pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/13 "2018-05-19T21:15:49Z")

</div>

That works, although I think it is significantly slower than my truncated basis function `rcs` in the R `rms` package. But the earlier discussion expressed a desire to rely on the `mgcv` package, not the `splines` package.

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [May 19, 2018, 9:44pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/14 "2018-05-19T21:44:09Z")

</div>

Either one is fine, particularly since they are both among CRAN’s Recommended packages.

---

<div class="post-metadata">

**Author:** ![harrelfe](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/harrelfe/32/2430_2.png) [@harrelfe](https://discourse.mc-stan.org/u/harrelfe)\
**Post date:** [May 19, 2018, 9:48pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/15 "2018-05-19T21:48:07Z")

</div>

`mgcv` mentions in the help file that it has a cyclic natural spline, so I’m wondering why we can’t find a natural non-cyclic (regular) spline.

For large n (e.g., 200,000) `ns` takes about 4 times longer than `Hmisc::rcspline.eval` (`ns` takes about 0.46s).

When n=1,000,000 `ns` takes only twice as long as the simpler truncated power basis, so I don’t think computing time is as much of an issue unless computations have to be repeated hundreds of times.

I would like to find how to construct the matrix multiplication that converts between natural B-spline and truncated power basis.

---

<div class="post-metadata">

**Author:** ![mike-lawrence](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mike-lawrence/32/59_2.png) [@mike-lawrence](https://discourse.mc-stan.org/u/mike-lawrence)\
**Post date:** [May 21, 2018, 9:31pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/16 "2018-05-21T21:31:08Z")

</div>

On the topic of QR decomposition & splines, see here: [Spline fitting demo inc comparison of sparse vs non-sparse](http://discourse.mc-stan.org/t/spline-fitting-demo-inc-comparison-of-sparse-vs-non-sparse/3287/6)

---

<div class="post-metadata">

**Author:** ![harrelfe](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/harrelfe/32/2430_2.png) [@harrelfe](https://discourse.mc-stan.org/u/harrelfe)\
**Post date:** [May 22, 2018, 12:01pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/17 "2018-05-22T12:01:53Z")

</div>

If we want to use any old basis function such as the convenient truncated power basis, we can easily convert back and forth to an orthonormal basis. Don’t we want SVD instead of QR?

---

<div class="post-metadata">

**Author:** ![albertoca](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/albertoca/32/5841_2.png) [@albertoca](https://discourse.mc-stan.org/u/albertoca)\
**Post date:** [September 10, 2019, 4:05pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/18 "2019-09-10T16:05:44Z")

</div>

Good morning, in relation to the previous comments, if I use the function ns within the formula, for example fit \<- brm(days | cens(1-event) ~ ns(age,df=3) …, can someone explain how do I have to use the function ‘set\_prior’ then? I have tried to explore it with get\_prior but it doesn’t work for me? Can someone provide some example of the code and some reflection on which would be a weakly informative prior, suitable for splines?

---

<div class="post-metadata">

**Author:** ![ucfagls](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ucfagls/32/1073_2.png) [@ucfagls](https://discourse.mc-stan.org/u/ucfagls)\
**Post date:** [September 10, 2019, 5:51pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/19 "2019-09-10T17:51:24Z")

</div>

For this, I think you (a user) are better off using `mgcv::smoothCon()` to create the spline basis. The basis you get is correct, it’s just not had any identifiability constraints applied to it. You can control this and more with `smoothCon()`s arguments, in this case `absorb.cons = TRUE` will return the required 4 basis functions.

For example:

```
smoothCon(z, data = data.frame(x), knots = list(x = kn),
          absorb.cons = TRUE)[[1]]$X

```

The `fx = TRUE` isn’t needed if all you are doing is generating the design matrix; that argument will just result in the penalty, the `$S` component being set to `numeric(0)` instead of generating the penalty. ~~`m = 2` doesn’t change anything in the design matrix either.~~

Once you have this, the `PredictMat()` counterpart to `smoothCon()` can allow one to go back from whatever `smoothCon()` does (depending on the arguments). With the above, a QR decomposition is being done to the basis and is returned as an attribute, which `PredictMat()` can use.

---

<div class="post-metadata">

**Author:** ![satpal](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/satpal/32/18639_2.png) [@satpal](https://discourse.mc-stan.org/u/satpal)\
**Post date:** [October 15, 2021, 2:00pm UTC](https://discourse.mc-stan.org/t/restricted-cubic-splines-in-stan/4234/20 "2021-10-15T14:00:12Z")

</div>

Hi,

I am sorry to bump an old thread. Has there been any progress on this issue i.e., implementing the restricted cubic splines in Stan?

Thanks