I want to fit a varying-intercept, varying-slope multilevel model estimating a covariance matrix for the coefficients

[Andrew Gelman]

Hi, I want to fit a varying-intercept, varying-slope multilevel model estimating a covariance matrix for the coefficients. What is the recommended way to do this? In particular, I want to use the efficient parameterization, the multivariate version of non-centered parameterization. I looked in the manual but I was getting confused as there wasn’t an example for this particular case. Thanks.
A

[Rob Trangucci]

The right way to do it is on page 61 (data declaration block) and 64 (optimized code) in the manual:

data {

  int<lower=0> N; // num individuals

  int<lower=1> K; // num ind predictors

  int<lower=1> J; // num groups

  int<lower=1> L; // num group predictors

  int<lower=1,upper=J> jj[N]; // group for individual

  matrix[N,K] x; // individual predictors

  matrix[J,L] u; // group predictors

  vector[N] y; // outcomes

}

parameters {

  matrix[K,J] z;

  cholesky_factor_corr[K] L_Omega;

  vector<lower=0>[K] tau; // prior scale

  matrix[L,K] gamma; // group coeffs

  real<lower=0> sigma; // prediction error scale

}

transformed parameters {

  matrix[J,K] beta;

  beta <- u * gamma + (diag_pre_multiply(tau,L_Omega) * z)';

}

model {

  vector[N] x_beta_jj;

  // vectorized likelihood on page 62

  for (n in 1:N)

    x_beta_jj[n] <- x[n] * beta[jj[n]];

  y ~ normal(x_beta_jj, sigma);

  L_Omega ~ lkj_corr_cholesky(2);

  tau ~ cauchy(0,2.5);

  to_vector(z) ~ normal(0,1);

  to_vector(gamma) ~ normal(0,5);

}

Rob

On Fri, Apr 10, 2015 at 1:33 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

Hi, I want to fit a varying-intercept, varying-slope multilevel model estimating a covariance matrix for the coefficients.  What is the recommended way to do this?  In particular, I want to use the efficient parameterization, the multivariate version of non-centered parameterization.  I looked in the manual but I was getting confused as there wasn’t an example for this particular case.  Thanks.

A

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Andrew Gelman]

Is this the Stan 2.6 manual?  Cos I don’t see this on pages 63 and 64 there.  Is there a newer version of the manual up?

A

[Rob Trangucci]

In this version it's on page 64: http://mc-stan.org/manual.html

Rob

[Andrew Gelman]

Ahh, I didn’t realize there was an update to 2.6.2.  Thanks.

A

[Mike Lawrence]

Hm, when I try to compile this from Rstan v2.6.0, I get the error:

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 

  failed to parse Stan model 'simple_gaussian' with error message:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

no matches for function name="multiply"

    arg 0 type=row vector

    arg 1 type=row vector

available function signatures for multiply:

0.  multiply(real, real) : real

1.  multiply(vector, real) : vector

2.  multiply(row vector, real) : row vector

3.  multiply(matrix, real) : matrix

4.  multiply(row vector, vector) : real

5.  multiply(vector, row vector) : matrix

6.  multiply(matrix, vector) : vector

7.  multiply(row vector, matrix) : row vector

8.  multiply(matrix, matrix) : matrix

9.  multiply(real, vector) : vector

10.  multiply(real, row vector) : row vector

11.  multiply(real, matrix) : matrix

expression is ill formed

ERROR at line 26

 24:     // vectorized likelihood on page 62

 25:     for (n in 1:N){

 26:     x_beta_jj[n] <- x[n] * beta[jj[n]];

                                                 ^

 27:     }

[Bob Carpenter]

Stan respects standard linear algebra rather than trying
to conform vectors. So you can't multiply two row vectors.
The reason it does this is because multiplying row vectors
a and b is ambiguous --- if it's a * b', you get a scalar
but if it's a' * b, you get a matrix.

You can use .* for elementwise multiplication, transpose
one or the other to get a scalar or matrix result, or
use the dot_product() function to just take a dot product.

- Bob

[Mike Lawrence]

Ah, I simply copied-and-pasted the code Ron posted earlier. I just tried the code from the manual itself and it works fine. Seems there was a missing ' at the end of the line assigning to x_beta_jj[n]

[Mike Lawrence]

Is there a working example with data anywhere for this model?

[Richard McElreath]

Attaching a didactic script I made for my lab. It first simulates data with varying intercepts and slopes and a specified correlation. Then the first Stan model is the "centered" parameterization. Then the second Stan model is the "non-centered" parameterization that uses the Cholesky factor.

I should have written it for general number of covariates and dimensions. But hopefully extrapolating to more than two dimensions is obvious.

[Mike Lawrence]

Thanks!

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~