I am currently learning Stan (MCMC C++ engine with wrappers in python and R) and implemented a linear mixed model

$y_{i,j} = \beta_0 + \mathbf{x}_{i,j}^T\beta + \alpha_i + \epsilon_{i,j},\ 1\leq i\leq I,\ 1\leq j\leq J$

where $\epsilon_{i,j}\sim\mathcal{N}(0,\sigma_{\epsilon}^2),\ \beta_k\sim\mathcal{N}(0,100),\ \alpha_l\sim\mathcal{N}(0,\sigma_r^2),\ \sigma_r^2\sim\text{Half-Cauchy}(0,25)$ and $I=4500,\ J=20$

I set $\beta_0=1,\beta_2=2,\beta_3=3$ and $\alpha_l$ where initialized to uniform random numbers between $1$ and $4$. So there is a random intercept effect per group.

plot of responses

After the sampler is done

fit = pystan.stan(model_code=stan_code, data=model_data, iter=10000, chains=3, thin=10, warmup=6000, n_jobs=3)

the estimates for $\beta_1,\beta_2$ are spot on but $\beta_0$ and the random effects are totally off.

glmm fit

Disregarding stan, is there something about the mixed model that I misunderstand?

Here is the stan code (using pystan)

model_data = {}
model_data['N'] = len(Y)
model_data['D_fixed'] = X_fixed_effect.shape[1]
model_data['D_random'] = X_random_effect.shape[1]
model_data['Y'] = Y
model_data['X_fixed'] = X_fixed_effect
model_data['X_random'] = X_random_effect
model_data['prior_mean_beta_fixed'] = np.zeros_like(beta_fixed)
model_data['prior_cov_beta_fixed'] = 100 * np.eye(len(beta_fixed))
model_data['prior_mean_beta_random'] = np.zeros_like(beta_random)
model_data['prior_cov_beta_random'] = np.eye(len(beta_random))

stan_code = """
data {
    int<lower=0> N;
    int<lower=0> D_fixed;
    int<lower=0> D_random;
    matrix[N,D_fixed]  X_fixed;
    matrix[N,D_random] X_random;
    vector[N] Y;

    vector[D_fixed]           prior_mean_beta_fixed;
    matrix[D_fixed,D_fixed]   prior_cov_beta_fixed;

    vector[D_random]            prior_mean_beta_random;
    matrix[D_random,D_random]   prior_cov_beta_random;
parameters {
    vector[D_fixed] beta_fixed;
    vector[D_random] beta_random;
    real<lower=0> sigma_eps;
    real<lower=0> sigma_random;
transformed parameters {

model {
    sigma_eps ~ cauchy(0,25);
    beta_fixed  ~ multi_normal(prior_mean_beta_fixed, prior_cov_beta_fixed);
    beta_random ~ multi_normal(prior_mean_beta_random, sigma_random * prior_cov_beta_random);

    Y ~ normal(X_fixed * beta_fixed + X_random * beta_random, sigma_eps);

fit = pystan.stan(model_code=stan_code, data=model_data, iter=10000, chains=3, thin=10, warmup=6000, n_jobs=3)

Thanks for the help

  • $\begingroup$ You need to give us more details. Why the results are off? Notice that effective sample size is extremely small, so you probably need more iterations. Did you check any convergence diagnostics? $\endgroup$ – Tim Jul 9 '18 at 11:13
  • $\begingroup$ I expected the posterior mean of $\hat{\beta}_0$ to be close to $1$, as this is the true value. I am confused why the effective sample size is so small. I was following the implementation by Hilbe in "Bayesian Models for Astrophysical Data" (ch. 8.2.4). He also had 4500 observations, 20 groups, 10k iterations, 3 jobs and 6000 warumup samples. Still he gets a much larger effective sample size for $\beta_0$ and all the random effect coefficients (~700). $\endgroup$ – tenticon Jul 9 '18 at 11:50
  • $\begingroup$ So you simulated the data? Are you sure that the authors you refer to used exactly the same model and software? If not, then there is no reason to expect that the same setting should work. $\endgroup$ – Tim Jul 9 '18 at 12:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.