A while ago I implemented Natural Cubic Splines (and their derivative) in Stan, as defined in Royston & Parmar.

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…

```
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])';
}
}
```