How to specify an extended random-effect model

Modeling
1

Hi,

Recently I have been trying to develop my basic hierarchical model into an extended one.

The basic model is a regular random-effect one like below:

data{
    int<lower=1> N;
    vector[N] coef;
    vector[N] se;
}
parameters{
    vector<lower=0>[N] coef_hat;
    real<lower=0> mu;
    real<lower=0> sigma;
}
model{
    sigma ~ cauchy( 0 , 5 );
    mu ~ normal( 2 , 5 );
    
    // measurement error model:
    coef ~ normal( coef_hat , se );
    // likelihood:
    coef_hat ~ normal( mu , sigma );
}

As you can see above, this random effect model corresponds to two layers of distributions:

coef_i = N(coefhat_i, se_i^2) and coefhat_i = N(mu , sigma^2). They can be expressed as coef_i = N(mu, se_i^2+sigma^2).

Now, imagine that the “mu” actually depends on the variance, i.e., coef_i = N(mu + beta * \sqrt{se_i^2+sigma^2}, se_i^2+sigma^2).

To model this, I tried below stan model:

data{
    int<lower=1> N;
    vector[N] coef;
    vector[N] se;
}
parameters{
    vector<lower=0>[N] coef_hat;
    real<lower=0> mu;
    real<lower=0> sigma;
    real beta
}
model{
    sigma ~ cauchy( 0 , 5 );
    mu ~ normal( 2 , 5 );
    beta ~ normal( 0 , 0.5);

    // measurement error model:
    coef ~ normal( coef_hat , se );
    // likelihood:
    coef_hat ~ normal( mu + beta*sqrt(square(se)+square(sigma)), sigma );
}

Obviously, the final line in the model block does not work, due to the incompatibility of dimision (error hint: Ill-typed arguments supplied to infix operator +.), and I am stuck here.

Any thoughts over here? Your attention and help will be highly appreciated.

Regards,
Elon

2

I’m not sure which interface to Stan you’re using, but our current parsers can parse that model if you fix the bug of not including a semicolon after the declaration of beta. I’d recommend cmdstanr or cmdstanpy to keep up with advanced in Stan.

When things don’t work, it can help to break things into a loop so you can at least get the model off the ground.

for (n in 1:N) {
  coef_hat[n] ~ normal(mu + beta * hypot(se[n], sigma), sigma);
}

Here, hypot(a, b) = sqrt(a^2 + b^2), i.e., the hypoteneuse.

3

Hi, Bob

Thanks for your clarification, and yes, after I have corrected my sloppy mistake of the missing semicolon, the model can work without loop (I am using cmdstanr).

Now, I have realized where my problem is, though I have no clue how to fix it. My actual model is below. Compared to my example model, my actual model has an important part about the imputation of “se”. Some of my observations (coef, se) lack “se” data, and so I have to impute them based on known characteristics variables X: se ~ lognormal( x * beta, sd_se);.

data{
    int<lower=0> N;
    int<lower=0> K;
    int<lower=0> N_obs;
    int<lower=0> N_mis;
    array[N_obs] int<lower = 1, upper = N_obs + N_mis> ii_obs;
    array[N_mis] int<lower = 1, upper = N_obs + N_mis> ii_mis;
    vector[N] coef;
    array[N_obs] real<lower=0> se_obs;
    matrix[N,K] x;
}
parameters{
    real<lower=0> sd_se;
    array[N_mis] real<lower=0> se_mis;
    array[N] real<lower=0> coefhat;
    real mu;
    real<lower=0> sigma;
    vector[K] beta;
    real gamma;
}
transformed parameters{
    array[N] real<lower=0> se;

    se[ii_obs] = se_obs;
    se[ii_mis] = se_mis;
}
model{   
    // priors:
    sd_se ~ normal(0 , 0.25);
    sigma ~ normal( 0 , 0.5 );
    mu ~ normal( 0.1 , 0.2 );
    beta ~ normal( 0 , 0.25);
    gamma ~ normal( 0 , 0.5);
    
    //imputation model
    se ~ lognormal( x * beta, sd_se);
    
    // measurement error model:
    coef ~ normal( coefhat , se );
    coefhat ~ normal(mu + gamma * hypot(se, sigma), sigma );
}

Due to the imputation process, I need to declare: array[N] real<lower=0> se;, which later conflicts with the final model line coefhat ~ normal(mu + gamma * hypot(se, sigma), sigma );

As the error hint says:

Ill-typed arguments supplied to infix operator *. … Instead supplied arguments of incompatible type: real, array real.

Here, the coefficient gamma is declared as real, while the se is declared as array[N] real.

Could you please spend a bit more time on my post? Thank you very much.

Warm regards,
Elon

4

Is there a reason that se needs to be an array of real values instead of a vector? If not, then it can be declared just as a vector. Alternatively, if it needs to stay an array of reals, then you could try coefhat ~ normal(mu + gamma * to_vector(hypot(se, sigma)), sigma);

I think that would solve the incompatibility error that you got

1 Like