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/
Post Reply
slavutich89
Posts: 3
Joined: Sat Apr 04, 2020 11:44 am

Comparing variance estimates from two MCMC models

Post by slavutich89 »

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: 1165
Joined: Mon Oct 19, 2009 10:34 am

Re: Comparing variance estimates from two MCMC models

Post by ChrisCharlton »

I asked George about this and this was his response:
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: 3
Joined: Sat Apr 04, 2020 11:44 am

Re: Comparing variance estimates from two MCMC models

Post by slavutich89 »

Dear Chris,

Thank you for your answer.

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: 1165
Joined: Mon Oct 19, 2009 10:34 am

Re: Comparing variance estimates from two MCMC models

Post by ChrisCharlton »

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: 1165
Joined: Mon Oct 19, 2009 10:34 am

Re: Comparing variance estimates from two MCMC models

Post by ChrisCharlton »

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: 3
Joined: Sat Apr 04, 2020 11:44 am

Re: Comparing variance estimates from two MCMC models

Post by slavutich89 »

Dear Chris,

Thank you for your help with this. Greatly appreciated.

Post Reply