Question about three-level hierarchical meta-analysis

Modeling
1

Dear all,

I have a question on how to fit a random effects meta-analysis with three levels. I have observations (effect sizes), some of which are nested in papers, and for which I would like to estimate an overall distribution of effect sizes. An additional difficulty is that only some of the papers have several effect sizes, while the majority of papers report only one. Here is the model with only effect sizes and the overall distribution, which works fine:

data {
int<lower=1> N;
real g[N]; //this is the effect size
real se[N]; //this is the standard error
}
parameters {
real theta[N];//effect size parameters for each observation
real mu;
real<lower=0> tau;
}
model {
g ~ normal(theta, se);
theta ~ normal(mu, tau);
mu ~ normal(0,10);
tau ~ cauchy(0,5);
}

I have now tried for some time to model that the parameters theta could be themselves distributed with mean mu_p and standard deviation tau_p (defined at the paper level), with mu_p ~ normal(mu,tau), but I cannot seem to make this work. I have also played around with loops instead of the vector notation. While I at least got the program to work in that case, it is extremely slow and still does not give me what I want. What happens is that in some range I get the same effect size mu_p, while for some papers the estimation seems to work. Here is my code for that:

data { 
  int<lower=1> N;
  real g[N];
  real se[N];
  int<lower=1,upper=N> study[N];
  int<lower=1,upper=143> paper[N];
  int<lower=0,upper=1> unique[N];//dummy indicating that the is only 1 effect in a paper
} 
parameters {
  vector[N] theta;
  vector[143] mu_p;
  real mu;"
  vector<lower=0>[143] tau_p;
  real<lower=0> tau;
} 
model {
 for (n in 1:N){
 for (k in 1:143){
    g[n] ~ normal(theta[study[n]],se[n]);
if (unique[n]==1) {
 theta[study[n]] ~ normal(mu_p[paper[k]],se[n]);
} else {
 theta[study[n]] ~ normal(mu_p[paper[k]],tau_p[paper[k]]);
  }
mu_p[paper[k]] ~ normal(mu,tau);
 }
}
mu ~ normal(0,10);
tau ~ cauchy(0,5);
}

Anybody out there who can help me with this?

Cheers,

Ferdinand

2

The main problem is in your program in you modeling block. For all N studies, you are looping over all 143 papers. That is a lot of extra and unnecessary work for Stan.

I think you want something like this without loops.


// expected effect size per paper
mu_p ~ normal(mu, tau);
tau_p ~ ... ;
// If you want different tau_p per study, a hierarchical structure just like for the mu,
// is going to help I think. 
// Or you can just make tau_p a scalar. Does it make sense to have the
// heterogeneity different per paper? Maybe start with a scalar and build up
// your stan program if that works.

// expected effect size per study. tau_p captures the heterogeneity between 
// studies in one paper.
theta ~ normal(mu_p[paper], tau_p[paper]);

// measurement error
g ~ normal(theta, se)

Given the a large g effect size is .8, I also think that your priors for mu and tau are too wide. Without knowing what you are studying, mu ~ normal(0, 1) seems more sensible to me.

3

Hi Stijn,

thanks so much for your help, this solved my problem. It works like a charm
now. Below I post the complete working model in case anybody else wants to
take a look at this.

Cheers!

Ferdinand

4

You didn’t include the program. You can just paste it into the message between triple back ticks:

```
***your code here***
```
5

Apologies, not sure why this did not come along. Here comes another attempt:

data {
  int<lower=1> N;//number of effect sizes
  int<lower=1> K;//number of papers
  vector[N] g; //effect  size
  vector[N] se; //SE of effect  size
  int<lower=1,upper=K> paper[N]; //index of paper
}
parameters {
 real theta[N];
 real mu_p[K];
 real<lower=0> tau_p[K];
  real mu;
  real<lower=0> tau;
}
model {
    g ~ normal(theta,se);
    theta ~ normal(mu_p[paper],tau_p[paper]);
    mu_p ~ normal(mu,tau);
    tau_p ~ cauchy(0,2);
mu ~ normal(0,1);
tau ~ cauchy(0,2);
}
6

As always, I recommend consistent indentation—without it, programs are very hard to scan. And whatever you do, don’t use tabs (no matter what Richard on Silicon Valley says—it’s a TV show!).

You probably are going to want a non-centered parameterization for mu_p (see the manual).

7

Hi Bob,

thanks for your tips. I have been mulling over them for some time now (as well as developing the model somewhat further). However, now it seems like I am stuck again, on two points: 1) reparameterizing this further to make the sampling more efficient; and 2) building random slopes into the model.

Here is my current (working) code for the random intercepts model:

    data { 
  int<lower=1> N;
  int<lower=1> G;
  vector[N] g;
  vector[N] se;
  int<lower=1,upper=G> group[N];
} 
parameters {
  vector[N] theta;
  vector[G] mu;
  real K;
  real<lower=0> sigma_int;
  vector<lower=0>[G] tau;
  real<lower=0> df;
} 
model {
  vector[N] intercept;
  intercept = K + mu[group];
  g ~ normal(theta , se);
  theta ~ student_t(df , intercept , tau[group]);
  K ~ normal(0 , 2);
  mu ~ normal(0 , sigma_int);
  tau ~ cauchy(0 , 1);
  sigma_int ~ cauchy(0 , 1);
  df ~ normal(10 , 100);
}

  1. As far as reparametrization is concerned, I tried to make the sampling more efficient (trying to follow your advise [here]) by rewriting the model as

    intercept = K + mu[group]*sigma_int;

However, that yields much worse sampling results in this context. I must be making an obvious mistake here, but I just cannot figure out what it is (it persists if I use a normal instead of the student_t distribution or play around with the priors).

  1. It should be possible to add a random slope in a straightforward manner, yet I could not figure out how to do this in the current context of the model. Lets say I want to add a regression parameter for the independent variable x, and I want this parameter to differ by group. I should be able to write something like:

    model {
    vector[N] intercept;
    vector[N] slope;
    intercept = K + mu[group];
    slope = S + b[group];
    g ~ normal(theta , se);
    theta ~ student_t(df , intercept + slope*x , tau[group]);
    K ~ normal(0,2);
    S ~ normal(0,2);
    mu ~ normal(0 , sigma_int);
    b ~ normal(0 , sigma_sl);
    tau ~ cauchy(0 , 1);
    sigma_int ~ cauchy(0 , 1);
    sigma_sl ~ cauchy(0 , 1);
    df ~ normal(10 , 100);
    }

If I run this, I get back an error that “the expression is ill formed”. I can of course see why that is the case, but I could not yet figure out a way around this issue–any pointers you could give me would be much appreciated!

Cheers,

Ferdinand

8

It helps if you send the entire program that’s causing an error and the exact output from Stan that you’re getting. Otherwise I have to paste all of this together and run through Stan.

You are probably also going to want to use a non-centered parameterization for mu

9

Apologies. I attach the two models I have troubles with, plus the data. When I run the random_slopes model, I get the following error:

--- Translating Stan model to C++ code ---

bin/stanc testmodels/random_slope.stan --o=testmodels/random_slope.hpp
Model name=random_slope_model
Input file=testmodels/random_slope.stan
Output file=testmodels/random_slope.hpp

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

No matches for:

vector * vector

Available argument signatures for operator*:

real * real
vector * real
row vector * real
matrix * real
row vector * vector
vector * row vector
matrix * vector
row vector * matrix
matrix * matrix
real * vector
real * row vector
real * matrix

No matches for:

vector + ill formed

Available argument signatures for operator+:

int + int
real + real
vector + vector
row vector + row vector
matrix + matrix
vector + real
row vector + real
matrix + real
real + vector
real + row vector
real + matrix
+int
+real
+vector
+row vector
+matrix

expression is ill formed
error in ‘testmodels/random_slope.stan’ at line 25, column 48

23: slope = S + b[group];
24: g ~ normal(theta,se);
25: theta ~ student_t(nr,intercept + slope*year_z,tau[group]);
^
26: K ~ normal(0,2);

make: *** [testmodels/random_slope.hpp] Error 253
##############################

Output from sampling

##############################

/bin/bash: ././testmodels/random_slope: No such file or directory

File /Users/nando/cmdstan/output*.csv not found
No valid input files, exiting.

I do of course realize that the vector*vector multiplication is what causes me grief here, but I just cannot figure a way out of it without creating some other error about incompatibility of dimensions.

In the random_intercept_improved model I have tried to move towards a non-centered parametrization as you suggest, but I am not at all sure that this is what it is supposed to look like, or that it is sufficient. The trace plots in the attached models look somewhat better than in earlier attempts, but still far from what I think they should look like, which again makes me suspect that I am doing something wrong here.

Cheers,

Ferdinand

random_slope.stan (669 Bytes)

random_intercept_improved.stan (624 Bytes)

data_example.csv (8.5 KB)