Three-level Multivariate Response model

[Jascalla]

Hi All,
Thanks ChrisCharlton for your earlier comment.
Here is an update with the syntax inserted.
I am very glad to join this forum.

I am currently running multivariate response models with three continuous outcome variables.
I was able to successfully fit single and two-level multivariate response models on the three outcome variables using runmlwin from within stata.

I am trying to extend the model to three-level but mlwin keep crashing anytime the model is running and as a result, I decided to run the three-level model on single outcome variable using runmlwin and that worked well.
I read the literature several times but I can't find anything on the three-level analysis for multivariate response models but only two-level.

Please, does it mean the 'runmlwin' can't fit the three-level multivariate response models but can do for only two-level multivariate response models? or am I missing out on something?
Below are the codes for this problem:

Sorting the data according to model hierarchy with sdist, hhid and childid been level 3, 2 and 1 respectively

Code: Select all

sort sdist hhid childid 

Code: Select all

. runmlwin ///
(waz cons hw1new m5new1 m5new2 b0new b11 b4 h11new h22new m18new1 m18new2 hw57new m15new, eq(1)) ///
(whz cons hw1new m5new1 m5new2 b0new b11 b4 h11new h22new m18new1 m18new2 hw57new m15new, eq(2)) ///
(haz cons hw1new m5new1 m5new2 b0new b11 b4 h11new h22new m18new1 m18new2 hw57new m15new, eq(3)), ///
level3(sdist: (cons, eq(1)) (cons, eq(2)) (cons, eq(3))) ///
level2(hhid: (cons, eq(1)) (cons, eq(2)) (cons, eq(3))) ///
level1(childid: (cons, eq(1)) (cons, eq(2)) (cons, eq(3))) maxiterations(200) nopause

The Error message is: "MLwiN has stopped working. A problem caused the program to stop working correctly. Windows will close the program and notify you if a solution is available".

Any help is highly appreciated.

Many thanks in advance.

[ChrisCharlton]

MLwiN should be able to fit this type of model. Would it be possible for you to provide example syntax and output demonstrating the problem that you are getting?

[Jascalla]

Dear All,

I have identified the problem to my query and I can confirm that the model is now running smoothly.
The problem is not with the syntax as I used same to achieve my desired results now.

This is the problem: I have very few missing observations in the dataset but the model ran successfully for the two-level multivariate model. When I introduced the third level, the MLwiN report an error which clearly did not talk about model specification error or wrong syntax but just keep crashing with the error "MLwiN has stopped working. A problem caused the program to stop working correctly. Windows will close the program and notify you if a solution is available".

Upon several readings with no success, I have decided to run the three-level model on 'fictitious data' which I created myself and is complete. To my surprise, the three-level model now ran correctly on the fictitious dataset with same syntax so I suspect those few missing observations in my dataset.

As a result, I dropped those few missing observations in my original dataset using "drop if missing(waz, whz, haz, hw1new, m5new, ....)" and the three-level model now ran correctly.

So what it means to me is that missing observations could create problem for some analysis at the higher levels of hierarchy like level 3 but not at level 2 or 1 in multivariate analysis as it is in my case. I hope this comments can help others in the near future who might face similar situations.

I remember sometime ago when I was running multivariate mixed effect model using MCMC and I keep getting errors due to missingness until GeorgeLeckie asked me to drop them and I ran the model successfully afterwards but I can't believe it can happen in this case too as listwise deletion is done automatically except in MCMC.

Finally, please, I have extracted my level-2 and level-1 residuals for each equation (I have 3 equations) and I used them to do the q-q plot to investigate multivariate normality for each outcome variable but seems like I am just working on the marginal distributions. However, I am unable to perform multivariate normality test (graphical or empirical) for the 3 outcomes jointly. How can I perform the joint multivariate normality test (graphical or empirical) for my fitted model with runmlwin?

Please, any form of help is highly appreciated.

Many thanks.

-J

[GeorgeLeckie]

Hi Jascalla,

Well done on exploring and identifying the source of the MLwiN crash. It would be great if you could email us offline the do-file and dataset to replicate the crash. We will then look into putting in a fix.

In terms of your second query concerning testing for multivariate normality, I am not sure. Of course you could plot a scatterplot matrix at each level as a rather informal inspection, but I am not sure what would be best formal approach. Let me turn the question around, how would you test for multivariate normality graphically or empirically had you instead fitted your model by -mixed-? If you can find the answer to this question then you can probably carry out the corresponding steps on the -runmlwin- residuals.

Best wishes

George

[Jascalla]

Hi George,

Good to hear from you.

For the first part, I will discuss that in email.

For the second part, there is a a package in Stata called 'multnorm' (see http://www.stata.com/stb/stb2/srd2/multnorm.hlp) which is a graphical procedure for examining multivariate normality diagnostics from a standard linear regression for at least two variables. The method plots Mahalanobis distance against its chi-square distribution and approximate straight line about 45 degrees in the x-y plane suggest multivariate normality. There are other empirical tests for multivariate normality like 'mvtest' in Stata which is not based on any fitted model residuals.

Currently, after obtaining my output from runmlwin, this is what I am doing with the residuals I retrieved from level 2 and level 1:

Fit level 2 model

Code: Select all

. runmlwin ///
(waz cons hw1new m5new1 m5new2 b0new b11 b4 h11new h22new m18new1 m18new2 hw57new m15new, eq(1)) ///
(whz cons hw1new m5new1 m5new2 b0new b11 b4 h11new h22new m18new1 m18new2 hw57new m15new, eq(2)) ///
(haz cons hw1new m5new1 m5new2 b0new b11 b4 h11new h22new m18new1 m18new2 hw57new m15new, eq(3)) ///
, ///
level2(hhid: (cons, eq(1)) (cons, eq(2)) (cons, eq(3)), residuals(u)) ///
level1(childid: (cons, eq(1)) (cons, eq(2)) (cons, eq(3)), residuals(e)) ///
maxiterations(200) nopause

Normal probability plot for checking level 2 normality of residuals for first equation with confidence bands

Code: Select all

. pnorm u0 if pickone_household==1, scale(1.96) aspectratio(1)

Normal probability plot for checking level 2 normality of residuals for second equation with confidence bands

Code: Select all

. pnorm u1 if pickone_household==1, scale(1.96) aspectratio(1)

Normal probability plot for checking level 2 normality of residuals for third equation with confidence bands

Code: Select all

. pnorm u2 if pickone_household==1, scale(1.96) aspectratio(1)

Also, I repeat same for level 1 residuals.

However, this procedure seems like a kind of an informal way of checking multivariate normality as you also rightly put it, although is the same approach I normally use when I fit multilevel models on only one outcome variable (univariate normality) using 'xtmixed' but which is for univariate normality.

From the literature, it was obvious that overall multivariate normality check is difficult to examine formally, especially for mixed effect models after fitting such models hence the need to use the above procedures to informally inspect.

So my query was to help me know if any such formal multivariate normality check is available and thank you George for your comments.

Many thanks.
-J

[GeorgeLeckie]

Hi J,

The multivariate normality assumption concerns the residuals rather than the response variable. So even if you have one response variable, as soon as you introduce a random slope you will be assuming that the random intercept and random slope are bivariate normally distributed. I am not aware of researchers generally going beyond the univariate normality plots you have described. All I can suggest is looking at a scatter plot matrix of the three sets of residuals at each level to see whether there are unusual outliers.

Best wishes

George

[Jascalla]

Hi George,
Many thanks for all these.
I am so happy about this part of your comments "I am not aware of researchers generally going beyond the univariate normality plots you have described". This is because I have been wondering for sometime now as to how I can do overall multivariate normality plots, that is going beyond the univariate normality plots I usually do.
You will soon hear from me on the other issue offline about the dataset and the do-files for the fix.
Thank you.
-J