Mixed effects model - cannot recover simulated value

Modeling
1

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

2

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
3

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