# Mixed effects model - cannot recover simulated value

Hi Stan community,

I am trying to model symptom trajectories of simulated data from N = 10 schizophrenia patients who were measured repeatedly over the course of 12 weeks.
I read a very interesting paper by Sorensen & Vasishth (2016) and was trying to replicate their varying intercept, varying slope mixed effects model in stan as follows:

``````data {
int<lower=1> N;                  // number of subjects
int<lower=0> t;                  // number of time points (predictor)
real Y[N * t];                   // psychopathology outcome
}

parameters {
vector[2] beta;                  // intercept and slope
real<lower=0> sigma_e;           // error sd
vector<lower=0>[2] sigma_u;      // subj sd
cholesky_factor_corr[2] L_u;
matrix[2, N] z_u;                // variance-covariance matrix for slopeXinter
}

transformed parameters{
matrix[2, N] u;
u = diag_pre_multiply(sigma_u, L_u) * z_u;	    // varying slopes and intercepts
}

model {
real mu;
int count;
//priors
L_u ~ lkj_corr_cholesky(2.0);       // definition of the prior
to_vector(z_u) ~ normal(0, 1);      // places normal distribution on z_u
//likelihood
count = 0;
for (i in 1:N){
// Posterior parameter distribution of the mean
for (j in 0:(t-1)){
count = count + 1;
mu = beta[1] + u[1, i] + (beta[2] + u[2, i]) * j;
// likelihood
Y[count] ~ lognormal(mu, sigma_e);
}
}
}
``````

The data we simulated for this purpose was defined as

``````set.seed(08202018)

# define variables
N        <- 10  # subjects
int      <- 85  # intercept
sd_int   <- 17
slope    <- -1.5 # slope
sd_slope <- 1
t        <- 12  # weeks

beta0 <- rnorm(N, int, sd_int)
beta1 <- rnorm(N, slope, sd_slope)

df <- data.frame(Y = NA,
week=rep(0:t, N),
ID = rep(1:N, each=t+1))

# linear regression to estimate the symptom trajectories
count <- 0
for (i in 1:N) {
for(j in 0:t){
count <- count + 1
df\$Y[count] <- beta0[i] + beta1[i] * j + rnorm(1, 0, 1)
}
}
``````

When we now run stan with

``````stanDat<-list(subj = as.integer(df\$ID),
Y = df\$Y,
t = t+1,
N = N)

## run stan
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

fit <- stan(file = "longi.stan", data = stanDat,
iter = 2000, chains = 4,
control = list(max_treedepth = 10, adapt_delta = .8))
``````

… we find that not only is the estimated intercept much lower than the simulated intercept, but also do the traceplots show no values for L_u, the Cholesky factor used to manipulate the variance-covariance matrix.

Therefore, I would be very greatful for your opinion on whether or not

1. the model as we set it up makes sense,
2. it is a reasonable choice to have the parameter L_u defined as “cholesky_factor_corr[2] L_u;”, and
3. you see any obvious reason why the beta (for the intercept) that we get with the print function is so way off.

It would be fantastic to get some thoughts and (to me very valuable) advice on this!
Thanks!
Stephanie

Is it intentional that the generated data uses a normal distribution:

``````df\$Y[count] <- beta0[i] + beta1[i] * j + rnorm(1, 0, 1)
``````

And the model to fit a lognormal?

``````Y[count] ~ lognormal(mu, sigma_e);
``````

Could the issue come from this?

1 Like

Thanks for the helpful hint, Ben!
That was indeed not intended and changing the model accordingly resulted in a way more reasonable beta 0!

1 Like