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, jld Jan 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[3] gamma_raw[subjects]; // using non-centered parameterization
  cholesky_factor_corr[3] C:
  vector<lower=0>[3] 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[3] gamma = L * gamma_raw[s];
    target += normal_lpdf(gamma_raw[s] | 0, 1); // implies gamma is MVN
    target += normal_lpdf(y[s] | gamma[2] + f(t[s], gamma[1], gamma[3]),
  target += lkj_corr_cholesky(C | 1); // C * C' is uniform a priori
  // priors on lambda and sigma

Not the answer you're looking for? Browse other questions tagged or ask your own question.