RStan Random Effects Model doesn't match with my simulated data from R

Hey guys,
First of all, i appreciate it very much being part of this group.

I’m currently working at a R/RStan project where i want to build a linear model with splines and random effects. So far, i tried my model on simulated data. The simulation was made in R.
So basically, the splines are working quite stable and good, i pretty much borrowed this part from Kharratzadeh case study (https://github.com/milkha/Splines_in_Stan , thank you at this point).
Unfortunately, the random effects i get from the stan-model do not match with my simulated random effects from R. I am getting results, but not the results i want to have… the random effects i get from my stan model should be as similar as possible as my random effects i generated in R.
I’m not really sure where my mistake is, so maybe someone of you could help me.
First:
How i simulated my data in R:

#Splines
X1 <- seq(from=-5, to=5, length=150) # generating inputs
B <- t(bs(X1, knots=seq(-5,5,1), degree=3, intercept = TRUE)) # creating the B-splines

num_data <- length(X1) 
num_basis <- nrow(B)

a0 <- 0.2 # intercept
a <- rnorm(num_basis, 0, 1) # coefficients of B-splines

#Random Effects
set.seed(123)
sigma_u0 <- 1
u_0j     <- rnorm(num_basis, mean = 0, sd = sigma_u0) #"Random intercepts"
beta_0j  <- rep(u_0j, length=num_data) # Varying intercepts 
cluster  <- rep(1:num_basis, length=num_data) # Cluster indicator

Y_true <- beta_0j + as.vector(a0*X1 + a%*%B)   # generating the output

So… i’m not 100% sure if this Y_true is well simulated for a mixed linear model with random effects and Splines, but it looks quite good to me. Any suggestions …?

#combine
dat      <- list(Ni      = length(Y_true2),
                 Nj      = num_basis,
                 cluster = cluster,
                 Y       = Y_true2,
                 X       = X1,
                 B       = B
)

sooooo and now to my Stan model. So far i just do normal b-splines without any smoothing splines or sth.

data {
 // Define variables in data
 int<lower=0> Ni;            // Number of observations (an integer) 
 int<lower=0> Nj;            // Number of clusters (an integer) 
 int<lower=1> cluster[Ni];   // Cluster IDs 
 vector[Ni] Y;               // output plus noise, model vector?
 vector[Ni] X;               // X as in intervalls
 matrix[Nj, Ni] B;           // splines-matrix

 }
 
  
parameters { 
  row_vector[Nj] a_raw;          //rows of the b-splines 
  real a0;                       //intercept
  real<lower=0> sigma;           //std.dev.
  real<lower=0> tau;             //coeff of b-splines
  real u_0j[Nj];                 // Random effects
  real<lower=0> sigma_u0;        // Standard deviation of random effects

} 
 
transformed parameters  {
  row_vector[Nj] a;
  vector[Ni] Y_hat;
  real mu[Ni];                // Individual mean
  
  a = a_raw*tau;  
  for (i in 1:Ni) {
   mu[i] = u_0j[cluster[i]];  // Individual mean 
 }
   Y_hat =  a0*X + to_vector(a*B);           //Creating Splines

 
 }
 
 model { 
  a_raw ~ normal(0, 1); 
  tau ~ normal(0, 1); 
  sigma ~ normal(0, 1); 
  sigma_u0 ~ normal(0, 1); //try
  u_0j ~ normal(0, sigma_u0);    // Random effects distribution
  for (i in 1:Ni)
  Y ~ normal(mu[Ni] + Y_hat, sigma);     //model

} 

If someone has any ideas for getting a solution for my problem, also how to optimize my model or making it more efficient, i would be forever grateful.

Thanks in advance,
Annika

I don’t really know much about splines, but it looked like the transposed splines matrix had columns of zeros. I switched to the splines2 package and it gave me a warning about knots=seq(-5, 5, 1) putting knots on the boundaries. If I moved the knot positions in that warning went away. I don’t have specific advice, but I think this is something work investigating. Any parameter multiplied by zero isn’t gonna be recoverable.

The a0 term is labeled intercept. Adding it to Y_true like a0 * X1 makes it a linear term. Is this a mistake? Also, the spline term should be able to fit a linear term. I’d expect an non-identifiability here between the linear term and the splines.

It might be useful to just try fitting the splines themselves first (something like y ~ normal(a * B, sigma)).

How is it that you’re diagnosing that the model is not working, by the way?

1 Like

Hello bbbales2,
thank you very much for your answer!

I didn’t know the splines2 package so far but i’ll definitely use with this one now.
That’s a very good point, thank you, i’ll definitely work on this one…!

So i want to build a linear regression model, and in this i want to implement splines and random effects. I thought this might be the right way, bringing the splines into a linear term… but I know that there are some problems discussed about intercept and splines. I’m reading some papers right now about this topic and maybe I have to leave it out in the end.
Do you think this possible non-identifiability here between the linear term and the splines could have influence on my random effects?

I simulated all my data, the random effects, the splines, the output and so on.
So I know my real random effects. If i run my stan model it should reproduce as precise as possible the random effects i simulated. But the random effects i get from my stan model are completely different to the simulated ones. So that’s my point, why i think my model is not working properly…
The splines are working really good so far. They match with my simulated point cloud, so they do what i want.
But i need to change my stan model, that i get the random effects i simulated in R, and i don’t know how to do this… Maybe there’s also a mistake in the simulation which i don’t see.

Thanks again,
Annika

I don’t know if there’s a difference or if one is better. I just didn’t understand the intercept argument or the zeros and tried a different package. The zeros were a problem, I think.

The intercept I still don’t understand. I think the way to interpret that is that the intercept term determines the value of the spline at zero.

Like, y = a x + b, the intercept term (b) determines what y is when x == 0. The difference here might be that the spline intercept might specify only what happens at x == 0, and not what happens everywhere else (in the linear regression example b adds on everywhere). I’d recommend investigating this. I don’t know for sure what’s going on. That’s just my guess.

Yeah. Get rid of the linear term for now. Even if in an idealized world there’d be no effect, any sort of big non-identifiability in the model can make things run way slower than they should.

The way to get this working is to break this into pieces. Get a model with just splines working (figure out for sure if you want a spline intercept or not). Get a model with just the group level random effects (or whatever sort of effects these are) working. Once they’re working individually, put em’ together.

1 Like

Hey bbbales2, thanks again for your time in helping me on my problem.

I just had a second thought about this one. As I told you i borrowed this splines code pretty much from the case study from Milad Kharratzadeh. Could it be, that a0 is not an intercept, it is a slope? So it’s just wrong labeled? That might explain why there’s no prior for a0 in the model part.
If I’m right, and it’s a slope, would you still leave this part out?

It’s conceivable you do a regression with a linear term + a non-parametric thing.

The issue is that if the non-parametric thing can explain the linear term, you might not get back the coefficients you expect (cause the data can be explained with linear + non-parameteric or only the non-parametric thing). It’s my experience that this sorta thing is hard to get right.

I think I’ve seen splines stuff where there’d be some sorta parametric term + a non-parametric thing that was guaranteed to not be able to fit that first thing, but you gotta do this manually.

I had a look at the pdf just now. I think get rid of that linear term still. But you can try generating from and fitting to just that line + spline model. I could be wrong and it works fine.

1 Like