Do-files to replicate entire MLwiN User & MCMC Manuals

[Jamoo]

Dear All,

I've been working through some of the manual chapters with my own data.

My plot of my residuals against normality for level 1 residuals is coming out rather strangely in Stata. When I run it in MlWin using the same data it seems to work fine - and be linear. When I calculate it in Stata it truncates the std normal values at 0 (see attached).

The code I use in Stata is included below and attached the graph produced by Stata.

I'd appreciate it if someone could help me locate where I've gone wrong.

With best wishes,

Jamie

--- Stata code begins here---
do "Analysis/Repetitive Code/DeleteRandomEffects.do"
runmlwin deltabmi cons if bmichangesample==1, ///
level2(numprogid: cons, residuals(u)) level1(pid: cons,residuals(e)) nopause

***Graph normality of residuals
*level 1 graph
summarize e0
capture drop e0std e0rank e0uniform e0nscore
generate e0std = (e0 - r(mean))/r(sd)
summarize e0std
egen e0rank = rank(e0)
summarize e0rank
generate e0uniform = (e0rank - 0.5)/_N
summarize e0uniform
generate e0nscore = invnorm(e0uniform)
summarize e0nscore
scatter e0std e0nscore, ///
yline(0) xline(0) ylabel(-4(1)4) xlabel(-4(1)4) aspectratio(1) ///
note("Level 1 (participant level) residuals derived from 2 level random intercept model (Table~\ref{tab:dbmivarcomp})")

[GeorgeLeckie]

Hi Jamie,

In your runmlwin command, you use the if qualifier to restrict the estimation sample to a subset of the full sample.

However, once you pull the level-1 residuals back to Stata you perform all the post-estimation level-1 residual data manipulation steps on the full sample.

I think you need to perform all your post-estimation level-1 residual data manipulation steps on exactly the same sample as the estimation sample.

Best wishes

George

[Jamoo]

Hi George,

Thanks for the prompt reply.

[quote="GeorgeLeckie"]
I think you need to perform all your post-estimation level-1 residual data manipulation steps on exactly the same sample as the estimation sample.

I tried doing as you suggested (see code below). Unfortunately, the graph remains stubbornly truncated.

Interestingly, I tried copying the code after the model (without the if statements) to your worked example. This seemed to produce a nice linear graph in Stata . Suggests to me that maybe the problem is with my model e* estimates once they are returned to Stata (as the model produces the residuals graph fine when paused in MLwiN).

---Stata code begins---
runmlwin deltabmi cons if bmichangesample==1, ///
level2(numprogid: cons, residuals(u)) level1(pid: cons,residuals(e)) nopause

summarize e0 if bmichangesample==1
generate e0std = (e0 - r(mean))/r(sd) if bmichangesample==1
egen e0rank = rank(e0) if bmichangesample==1
generate e0uniform = (e0rank - 0.5)/_N if bmichangesample==1
generate e0nscore = invnorm(e0uniform) if bmichangesample==1
scatter e0std e0nscore if bmichangesample==1, ///
yline(0) xline(0) ylabel(-4(1)4) xlabel(-4(1)4) aspectratio(1) ///
note("Level 1 (participant level) residuals derived from 2 level random intercept model (Table~\ref{tab:dbmivarcomp})")

With best wishes,

Jamie

[GeorgeLeckie]

e0nscore is defined as the inverse normamal of the variable e0uniform and ranges from approximately -4 to 0 rather than from -4 to 4.

This suggests that e0uniform ranges from 0 to 0.5 rather than from 0 to 1.

In the line

Code: Select all

generate e0uniform = (e0rank - 0.5)/_N if bmichangesample==1

_N will refer to the full sample rather than the estimation sample. Given that e0nscore ranges from -4 to 0 rather than from -4 to 4, II suspect that the full sample is approximately twice the estiimation sample.

Try replacing this line with the following

Code: Select all

summarize e0rank
generate e0uniform = (e0rank - 0.5)/r(N) if bmichangesample==1

Best wishes

George

[Jamoo]

Dear George,

Thanks for the reply again.

Just had a chance to try the fix and can confirm that you've nailed it.

Thanks so much,

Best wishes

Jamie

[robertgrant]

Dear George and colleagues

I am wondering if you have any tips on runmlwin for multilevel factor analysis, which is not yet in the examples... or is it not there because it uses an entirely different MCMC algorithm? Do you think it is easier to run a simulation study via runmlwin or through MLwiN macros for this sort of model?

Thanks for all the examples and do-files posted so far, it's a really impressive piece of work!

Robert

[GeorgeLeckie]

Hi Rob,

Multilevel factor analysis is the one thing we have not implemented in runmlwin.

This is why there is no do-file for multilevel factor analysis in MLwiN MCMC manual.

The way to setup a multilevel factor analysis model in MLwiN is to first set up a multilevel multivariate response model.

You could therefore use runmlwin (without the nopause option) to set up such a model in MLwiN.

You would then need to use the MLwiN point and click interface to finish setting up the multilevel factor analysis model.

In general, it is much easier to use runmlwin to run a simulation study than it is to wirte a MLwiN macro for this purpose.

However, as multilevel factor analysis models are not implemented using runmlwin it is not possible to do a simulation study using runmlwin for these models.

So you will have to write a MLwiN macros I'm afraid.

Best wishes

George