Hierarchical model of antibody levels

[Miracle]

Hi,
I’m working on a hierarchical model which has measurement-specific antibody levels for different individuals. There are two sets of measurements collected at two different periods of time for all the individuals. So we expect that there should be measurement error bias because of some changes that might have occurred in the mode of measurement in the second period as compared to the first time. The traditional approach would be to estimate the parameters for each set of measurements independently, using the data for each. But I need to use a hierarchical model for this problem. I generated synthetic data with 50 measurements and 100 individuals where the first 30 measurements have a different measurement error (with zero mean) from the last 20 measurements (which has a non-zero mean). And I’m looking at how the hierarchical model for this problem would look in order to properly estimate the parameters of the model from the data. I have attached my code for this. Not all parameters were well estimated when I tried and I guess it could be a problem of how I formulated the model. I have attached the code and the simulated data as well as the plot of the simulated data. I would appreciate your timely assistance and guidance on how to solve this.

True values are:
beta=-10
muu=200
sigmau=1

mug_2=15
sigmag_1=4
sigmag_2=0.5

BR, Miracle

Diff_measurement_errors.csv (82.6 KB)

Diff_measurement_errors_2.csv (82.6 KB)

Case_2_hierarchy_R.R (920 Bytes)

Case2_hierarchy.stan (982 Bytes)

[martinmodrak]

Looking at a slightly extended pairs plot:

pairs(fit, pars = c("beta", "muu","sigmau", "sigmag_1","mug_2", "U[1]", "lp__"), las = 1)

It appears there are two distinct modes, but one has substantially lower log probability, so it should probalby not be relevant for the posterior, but some chains get stuck in this mode. The high probability mode recovers your parameters pretty well. The other mode can be characterized as:

parameter         mean          sd        2.5%         25%         50%         75%       97.5%
  beta      -9.8596005 0.084921209 -10.0204176  -9.9153187  -9.8676944  -9.8043096  -9.6767373
  muu        7.6930757 2.003437201   4.7650798   5.7910531   7.9143447   9.5193770  11.1417022
  sigmau     1.0952892 0.067344331   1.0125373   1.0430084   1.0752789   1.1386108   1.2552557
  mug_2    206.5754087 2.113732342 203.0213513 204.5585816 206.4127521 208.4012331 210.0167018
  sigmag_1 144.5988477 2.775781604 139.6795985 142.7552799 144.4016960 146.6372281 149.9807145
  sigmag_2   0.5118676 0.009525023   0.4955096   0.5058612   0.5110425   0.5168931   0.5383117

So large muu and small mug_2, sigmag_1 are switched for small muu and large mug_2 and sigmag_1. So it appears that the second mode corresponds to a very poor fit in the first part which catches up in the second part.

I do not immediately see why any local change from the second mode results in a decrease in log-posterior, so that’s something you’ll need to investigate. It is plausible that if you have prior knowledge that mug_2 and sigmag_1 cannot be so big, adding this as a prior could make the second mode go away.
Additionally, if one is sure the second mode is superfluous, one can likely avoid it completely by a good initialization. But the best way would be to understand why the second mode appears and modify the model to completely avoid it.

A good strategy would be to simplify the model and see when the problems disappear (either by removing some aspects of the model altogether, or taking some parameters as known and passing them as data instead of estimating them).

If I replace the gumbel with normal in the model, the problems are still there, so it is likely not caused by some special property of the gumbel distribution…

Best of luck in further investigations.

[Miracle]

Hi, thank you for your comment. Following your suggestions, I tried to investigate more.
First, I gave upper bounds to sigmag_1 and mug_2, and the results improved only in the sense that sampled parameter values got really close to the true parameters. But the pairs plot looks really weird, don’t understand what’s going on there (see with_upper_bound.png). Next, I fixed mug_2 (with no upper bound for sigmag_1), and other parameters seem to be well estimated except for sigmau which I think the chain didn’t mix very well (nevertheless, the true value of sigmau was captured in the posterior samples…see fixed_mug_2.png). Finally, I fixed sigmag_1 (with no upper bound for mug_2), parameters were also well estimated and the same poor mixing occured again for sigmau. However, the pairs plot for this case shows two different colors and I don’t quite understand the reason for that (see fixed_sigma_1.png). It’s not so clear to me right now, what could be the challenge with getting all parameters properly estimated. Your comments are welcomed.

[martinmodrak]

@Miracle

However, the pairs plot for this case shows two different colors and I don’t quite understand the reason for that

If you mean the yellow dots, then those indicate iterations where the sampler reached maximum treedepth which means that it is exploring very slowly.

I took a second look and realized that reparametrizing U as non-centered could help the model, i.e. having:

parameters {
 vector[N] U_raw;
 ...
}


transformed parameters {
  vector[N] U = U_raw * sigmau + muu;
  ...
}

model {
        U_raw ~ std_normal();
        ...
}

And indeed the model seems much better behaved in this setting as the U_raw are now uncorrelated.

Additionally, it seems like you could marginalize out U completely - this is essentially an autoregressive (AR) model and could AFAIK be treated as multivariate normal with carefully constructed covariance matrix, which could further improve the speed and general behaviour of them model.

[Miracle]

Hi Martin,

I’m happy with the results after reparameterizing U as non-centered as you suggested. The results greatly improved and we can work with that. In the future when I have a
similar problem, I will consider using the multivariate normal distribution as you also suggested. Thank you so much for your assistance.