## Comparing variance estimates from two MCMC models

Welcome to the forum for runmlwin users. Feel free to post your question about runmlwin here. The Centre for Multilevel Modelling take no responsibility for the accuracy of these posts, we are unable to monitor them closely. Do go ahead and post your question and thank you in advance if you find the time to post any answers!

Go to runmlwin: Running MLwiN from within Stata >> http://www.bristol.ac.uk/cmm/software/runmlwin/
slavutich89
Posts: 5
Joined: Sat Apr 04, 2020 11:44 am

### Comparing variance estimates from two MCMC models

I'd like to compare variance components from two models estimated with MCMC. Model 1 is an unconditional model and Model 2 includes Level 1 predictors. I need to calculate the relative percentage change in Level 1 variance between two models, such that this change also contains a standard error and p-value (I saw it in one paper using MLwiN and would like to replicate it).

Using one of the MLwiN examples, here is what I intend doing (my model differs from the example below, but the logic is the same):

Code: Select all

``````use "http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial.dta", clear
runmlwin normexam cons, level2(school: cons) level1(student: cons) nopause
runmlwin normexam cons, level2(school: cons) level1(student: cons) mcmc(on) initsprevious nopause
estimates store model1
runmlwin normexam cons girl standlrt, level2(school: cons) level1(student: cons) nopause
runmlwin normexam cons girl standlrt, level2(school: cons) level1(student: cons) mcmc(on) initsprevious nopause
estimates store model2``````
After estimating these two models, I'd like to run something like:

Code: Select all

``nlcom ([model1][RP1]var(cons)-[model2][RP1]var(cons))/([model1][RP1]var(cons)+[model1][RP2]var(cons))``
Normally, it could be done by first running "suest model1 model2", but it is not possible to use with runmlwin or other multilevel models.

Is it possible to simply store MCMC chains from both estimations (mcmcsum, getchains), manually merge them and then do the following?

Code: Select all

``````gen dif=(RP1_var_cons_model1-RP1_var_cons_model2)/(RP1_var_cons_model1+RP2_var_cons_model1)
mcmcsum dif, variables``````
ChrisCharlton
Posts: 1250
Joined: Mon Oct 19, 2009 10:34 am

### Re: Comparing variance estimates from two MCMC models

I don't think it is possible to do what they want as the percentage change in the level-1 variance involves two different models fitted to the sample data. So while you could mechanically do what is suggested, I don't think it would be statistically valid. That is, you could save a chain of 1000 values for the level-1 variance from each of your two models as two separate columns of data and then calculate a third variable as the percentage change in the level-1 variance, then report the mean and SD of this as the point estimate and SE. The problem is that you have assumed the sampling co-variability of the level-1 variance parameter process the two models is zero when it would be positive.

One approach which might be more fruitful would be to take a bootstrapping approach and estimate both models on each of the B bootstrapped samples and calculate the percentage change in the level-1 variance on each of the B bootstrapped samples and go about it that way.
slavutich89
Posts: 5
Joined: Sat Apr 04, 2020 11:44 am

### Re: Comparing variance estimates from two MCMC models

Dear Chris,

Would it be possible for you to share a Stata code how I can use a bootstrapping approach in my MCMC example from the previous post?

I know it could have been potentially done through comparing scalars in a similar way that is described here:
https://www.statalist.org/forums/forum/ ... ll-be-save
I have tried bootstrapping a difference in scalars as in this Statalist post by simply storing scalars for RP1:var(cons) from each models. But I cannot use "return scalar" after "runmlwin" because Stata returns the following error: "non r-class program may not set r()"
ChrisCharlton
Posts: 1250
Joined: Mon Oct 19, 2009 10:34 am

### Re: Comparing variance estimates from two MCMC models

runmlwin is an estimation command, so you need need to retrieve the returned results using ereturn rather than return (see https://www.stata.com/help.cgi?ereturn). You can find the model estimates in the matrix e(b), and their covariance matrix in e(V).

I am not aware of any example code for the bootstrapping approach and runmlwin, but will check with George and post again if we do.
ChrisCharlton
Posts: 1250
Joined: Mon Oct 19, 2009 10:34 am

### Re: Comparing variance estimates from two MCMC models

Unfortunately we don't have any examples, however George gave the following outline of how this would work:

You would need to set up a bootstrapping routine in Stata, which would call a program that you define. In this you would fit each model and store the quantities of interest as scalars. The bootstrapping routine would then monitor these scalars.

If you look in the bootstrapping topic in the Stata manuals (https://www.stata.com/manuals13/rbootstrap.pdf) you should find examples of how to do this. For your purposes you will need to do a clustered bootstrap.
slavutich89
Posts: 5
Joined: Sat Apr 04, 2020 11:44 am

### Re: Comparing variance estimates from two MCMC models

Dear Chris,

Thank you for your help with this. Greatly appreciated.
slavutich89
Posts: 5
Joined: Sat Apr 04, 2020 11:44 am

### Re: Comparing variance estimates from two MCMC models

Dear Chris,

I would like to follow up on this post, if it is fine with you. I have managed to write a bootstrap program in Stata and it worked fine for a 2-level model, as the one above.

I now have a question regarding having more levels in the model and using the same produce as you and George have recommended. If I have firms (Level 2) observed over several years (Level 1) and which are nested within subnational regions (Level 3) and these regions are nested within countries (Level 4), I am wondering whether bootstrap clustering should be done at the aggregate (country in this case) level? I only have 10 countries. Or perhaps clustering at the firm level would be sufficient? Thank you.
ChrisCharlton
Posts: 1250
Joined: Mon Oct 19, 2009 10:34 am

### Re: Comparing variance estimates from two MCMC models

Sorry, I'm afraid that we don't know the answer to this. George thought that the following paper, along with the papers that it references may be of some help:

Austin, P. C., & Leckie, G. (2020). Bootstrapped inference for variance parameters, measures of heterogeneity and random effects in multilevel logistic regression models. Journal of Statistical Computation and Simulation, 1-25.
LINDAOKEEFFE2020
Posts: 1
Joined: Mon Nov 02, 2020 12:40 pm

### Re: Comparing variance estimates from two MCMC models

Hi All,

I am examining change over time in SBP for females and males (2 level model with individual level random effects and occasion level for intercept and slopes) - models are sex-stratified. I want to formally compare estimates from the two models producing an estimate of the difference in SBP between females and males at different points on the trajectory. I have stored the estimation results for both the female and male model but can't seem to formally test whether (for example) _cons for model 1 is different from _cons for model 2. Slavutich89 I wondered if you could share your bootstrapping code as I think this may work for me?

Thanks,
Linda
slavutich89
Posts: 5
Joined: Sat Apr 04, 2020 11:44 am

### Re: Comparing variance estimates from two MCMC models

Hello Linda,

You could look how they do bootstrapping for a difference of two coefficients from two estimations:
https://www.statalist.org/forums/forum/ ... difference

For clustered bootstrap (panel data), see this link:
https://www.statalist.org/forums/forum/ ... s-together

Just use runmlwin commands instead of commands in these examples on Statalist. You will need to see how the estimate you need is called in the matrix produced by runmlwin. You can use matrix list e(b) after each estimation.

The explanation why clustered bootstrap should be done this way is here:
https://www.stata.com/support/faqs/stat ... anel-data/

Hope it helps.