# Fitting regression spline [closed]

I am reading the paper by Willemsen et al (2015), "A multivariate Bayesian model for embryonic growth", Statistics in Medicine, 34:8, 1351–1365

I have a model like $$y_{ij} = \gamma_{i2} + f((t_{ij} - \gamma_{i1}) / exp(-\gamma_{i3}) ) + \epsilon_{ij}$$ Here,$y_{ij}$= response, for ith person at jth time point.suppose i= 1,2,3,4 and j= 1,2,3

here, $$\gamma \sim MVN(0,\Sigma)$$ and $$\epsilon \sim N(0, \sigma^{2})$$

Here, f() is a spline function. I am trying to fit this model in rjags or rstan. but the problem is I have to loop through the spline function.Because, $\gamma$ is a random effect and two gamma parameters are inside the spline function.Is it possible to loop through the spline function in rjags or rstan ?

# Note :

Any kind of example for looping through the spline will be highly appreciated.

## closed as off-topic by AdamO, mdewey, jbowman, Michael R. Chernick, jldJan 25 '18 at 1:19

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – AdamO, mdewey, jbowman, Michael R. Chernick, jld
If this question can be reworded to fit the rules in the help center, please edit the question.

It somewhat depends on what sort of a spline you want, but there is a Stan case study on splines here. The essence of your Stan program would be something like:

functions {
vector f(vector t, real gamma_1, real gamma_3); // implement this spline
}
data {
int<lower=1> subjects;
int<lower=1> periods;
vector[periods] y[subjects];
vector[periods] t[subjects];
}
parameters {
vector gamma_raw[subjects]; // using non-centered parameterization
cholesky_factor_corr C:
vector<lower=0> lambda; // standard deviations for gamma
real<lower=0> sigma;
}
model {
matrix[3,3] L = diag_pre_multiply(C, lambda);
for (s in 1:subjects) {
vector gamma = L * gamma_raw[s];
target += normal_lpdf(gamma_raw[s] | 0, 1); // implies gamma is MVN
target += normal_lpdf(y[s] | gamma + f(t[s], gamma, gamma),
sigma);
}
target += lkj_corr_cholesky(C | 1); // C * C' is uniform a priori
// priors on lambda and sigma
}