Mixed effects model with more than one dependent variable

[Kjell Weyde]

I am going to analyse data where there are several dependent variables. These are different kinds of results from different cognitive tests on children. In addition, children are nested within different Schools, and therefore, I thought about using a mixed model with school as random effect. However, how can I account for multiple dep vars? Should I do one test for each dep var (there are about 20), or is there any method in Stata to analyse all the dep vars at once (ie, like mvreg, but for hierarchical data)?

Best,
Kjell Weyde

[Clyde Schechter]

You can, in principle, do this with -gsem-. See -help sem-, click the link to Examples, and then select the link to example 39g to see how to set up a three-level model in -gsem-. The example given is for negative binomial regression, but it will work equally well with simple regression, and incorporating multiple dependent variables is just a matter of adding them to the list of dependent variables.

That said, the advantage of doing this simultaneous estimation approach for the different outcomes (as opposed to doing separate -mixed- analyses) will only be realized by also specifying a covariance structure among the outcomes. With 20 dependent variables, this will entail estimating 20*19/2 = 190 covariance parameters. This will only be possible with a pretty large data set. Moreover, with that many added parameters it will take a long time to converge, and it might not converge at all.

The syntax will look something like this:

gsem (outcome1 outcome2 /*etc*/ <- predictor1 predictor2 /*etc*/ ///
    M1[school] M2[school>child]), cov(e.outcome1*e.outcome2) /*etc.*/

Good luck!

[Weiwen Ng]

Another way that you can do this in principle is move to a random effects framework, either within sem/gsem or within a Bayesian framework. You could assume that the effects for each outcome arise from some distribution - in your case, that is probably a crossed random effect, if you're in still in the frequentist framework. It's similar to a meta-analysis in principle. I've seen this done in a presentation, and I can grasp (if barely) the theoretical rationale for this, but I haven't seen this done in practice in my field. Also, you have to put all variables on the same scale (as far as I know) for this approach to work. It is perfectly acceptable to just run 20 regressions for each dependent variable.

[Maoliang Ling]

Hi everyone, I met a similar problem where there are 5 ordinal dependent variables measued by 5-point Likert-scale. I want to estimate the influences of a set of indepdendent variables on each of them (with -ologit-). However, these dependent variables are found to be highly correlated with each other. So I think it would be better to model them simultaneously and take into account their correlations. However, it seems that -gsem- cannot cope with this situation as it cannot specify the covariance structure between ordinal outcomes. Is there another way to solve this problem? Thank you very much!

[Joseph Coveney]

Maybe you can do something like the following. Begin at the "Begin here" comment; the first part of the do-file just creates a fictional dataset for illustration.

.ÿ
.ÿversionÿ16.0

.ÿ
.ÿclearÿ*

.ÿ
.ÿsetÿseedÿ`=strreverse("1531457")'

.ÿ
.ÿtempnameÿCorr

.ÿmatrixÿdefineÿ`Corr'ÿ=ÿJ(4,ÿ4,ÿ0.75)ÿ+ÿI(4)ÿ*ÿ0.25

.ÿquietlyÿdrawnormÿlat1ÿlat2ÿlat3ÿlat4,ÿdoubleÿcorr(`Corr')ÿn(250)

.ÿgenerateÿintÿpidÿ=ÿ_n

.ÿ
.ÿlocalÿcut_list

.ÿforvaluesÿcutÿ=ÿ0.2(0.2)0.8ÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿlocalÿcut_listÿ`cut_list'ÿ`=invnormal(`cut')'
ÿÿ3.ÿ}

.ÿforvaluesÿiÿ=ÿ1/4ÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿgrologitÿlat`i',ÿgenerate(out`i')ÿcuts(`cut_list')
ÿÿ3.ÿ}

.ÿ
.ÿforvaluesÿiÿ=ÿ1/2ÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿgenerateÿdoubleÿpr`i'ÿ=ÿruniform(-0.5,ÿ0.5)
ÿÿ3.ÿ}

.ÿ
.ÿ*
.ÿ*ÿBeginÿhere
.ÿ*
.ÿforvaluesÿiÿ=ÿ1/4ÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿquietlyÿologitÿout`i'ÿc.(pr1ÿpr2)
ÿÿ3.ÿÿÿÿÿÿÿÿÿestimatesÿstoreÿModel`i'
ÿÿ4.ÿ}

.ÿ
.ÿsuestÿModel?,ÿvce(clusterÿpid)

SimultaneousÿresultsÿforÿModel1,ÿModel2,ÿModel3,ÿModel4

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿÿÿÿ250

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ(Std.ÿErr.ÿadjustedÿforÿ250ÿclustersÿinÿpid)
------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿRobust
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿCoef.ÿÿÿStd.ÿErr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
Model1_out1ÿÿ|
ÿÿÿÿÿÿÿÿÿpr1ÿ|ÿÿÿ.2154766ÿÿÿ.4274275ÿÿÿÿÿ0.50ÿÿÿ0.614ÿÿÿÿÿ-.622266ÿÿÿÿ1.053219
ÿÿÿÿÿÿÿÿÿpr2ÿ|ÿÿÿ.3256259ÿÿÿ.4159809ÿÿÿÿÿ0.78ÿÿÿ0.434ÿÿÿÿ-.4896816ÿÿÿÿ1.140933
-------------+----------------------------------------------------------------
/Model1ÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿcut1ÿ|ÿÿÿÿ-.65183ÿÿÿ.1356932ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ-.9177837ÿÿÿ-.3858763
ÿÿÿÿÿÿÿÿcut2ÿ|ÿÿ-.2279872ÿÿÿ.1301684ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ-.4831125ÿÿÿÿ.0271381
ÿÿÿÿÿÿÿÿcut3ÿ|ÿÿÿ.3861968ÿÿÿ.1317042ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.1280612ÿÿÿÿ.6443323
ÿÿÿÿÿÿÿÿcut4ÿ|ÿÿÿ.7496668ÿÿÿ.1380844ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.4790264ÿÿÿÿ1.020307
-------------+----------------------------------------------------------------
Model2_out2ÿÿ|
ÿÿÿÿÿÿÿÿÿpr1ÿ|ÿÿÿ.4864051ÿÿÿ.4341661ÿÿÿÿÿ1.12ÿÿÿ0.263ÿÿÿÿ-.3645449ÿÿÿÿ1.337355
ÿÿÿÿÿÿÿÿÿpr2ÿ|ÿÿÿ.6244087ÿÿÿ.3965368ÿÿÿÿÿ1.57ÿÿÿ0.115ÿÿÿÿ-.1527893ÿÿÿÿ1.401607
-------------+----------------------------------------------------------------
/Model2ÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿcut1ÿ|ÿÿ-.7340508ÿÿÿ.1362721ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ-1.001139ÿÿÿ-.4669625
ÿÿÿÿÿÿÿÿcut2ÿ|ÿÿ-.2147664ÿÿÿÿ.130228ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ-.4700085ÿÿÿÿ.0404758
ÿÿÿÿÿÿÿÿcut3ÿ|ÿÿÿ.2892741ÿÿÿ.1298267ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.0348185ÿÿÿÿ.5437297
ÿÿÿÿÿÿÿÿcut4ÿ|ÿÿÿ.6991941ÿÿÿÿ.137494ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.4297107ÿÿÿÿ.9686775
-------------+----------------------------------------------------------------
Model3_out3ÿÿ|
ÿÿÿÿÿÿÿÿÿpr1ÿ|ÿÿÿ.2463621ÿÿÿ.4150305ÿÿÿÿÿ0.59ÿÿÿ0.553ÿÿÿÿ-.5670827ÿÿÿÿ1.059807
ÿÿÿÿÿÿÿÿÿpr2ÿ|ÿÿÿ.1213054ÿÿÿ.3695268ÿÿÿÿÿ0.33ÿÿÿ0.743ÿÿÿÿ-.6029538ÿÿÿÿ.8455646
-------------+----------------------------------------------------------------
/Model3ÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿcut1ÿ|ÿÿ-.5513776ÿÿÿ.1317963ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ-.8096937ÿÿÿ-.2930616
ÿÿÿÿÿÿÿÿcut2ÿ|ÿÿ-.0569693ÿÿÿ.1267136ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ-.3053233ÿÿÿÿ.1913847
ÿÿÿÿÿÿÿÿcut3ÿ|ÿÿÿ.4310226ÿÿÿ.1290643ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.1780612ÿÿÿÿ.6839841
ÿÿÿÿÿÿÿÿcut4ÿ|ÿÿÿ.7271841ÿÿÿ.1355828ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.4614467ÿÿÿÿ.9929215
-------------+----------------------------------------------------------------
Model4_out4ÿÿ|
ÿÿÿÿÿÿÿÿÿpr1ÿ|ÿÿ-.0358479ÿÿÿÿ.436053ÿÿÿÿ-0.08ÿÿÿ0.934ÿÿÿÿÿ-.890496ÿÿÿÿ.8188003
ÿÿÿÿÿÿÿÿÿpr2ÿ|ÿÿÿ.5275623ÿÿÿ.4046977ÿÿÿÿÿ1.30ÿÿÿ0.192ÿÿÿÿ-.2656306ÿÿÿÿ1.320755
-------------+----------------------------------------------------------------
/Model4ÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿcut1ÿ|ÿÿ-.4178492ÿÿÿÿ.130821ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ-.6742536ÿÿÿ-.1614448
ÿÿÿÿÿÿÿÿcut2ÿ|ÿÿÿ.0072064ÿÿÿ.1286196ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ-.2448833ÿÿÿÿ.2592961
ÿÿÿÿÿÿÿÿcut3ÿ|ÿÿÿÿ.347783ÿÿÿ.1309459ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.0911336ÿÿÿÿ.6044323
ÿÿÿÿÿÿÿÿcut4ÿ|ÿÿÿ.7268718ÿÿÿ.1374533ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.4574683ÿÿÿÿ.9962753
------------------------------------------------------------------------------

.ÿtestÿ[Model1_out1]pr1ÿ[Model2_out2]pr1ÿ[Model3_out3]pr1ÿ[Model4_out4]pr1ÿ

ÿ(ÿ1)ÿÿ[Model1_out1]pr1ÿ=ÿ0
ÿ(ÿ2)ÿÿ[Model2_out2]pr1ÿ=ÿ0
ÿ(ÿ3)ÿÿ[Model3_out3]pr1ÿ=ÿ0
ÿ(ÿ4)ÿÿ[Model4_out4]pr1ÿ=ÿ0

ÿÿÿÿÿÿÿÿÿÿÿchi2(ÿÿ4)ÿ=ÿÿÿÿ1.41
ÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿ=ÿÿÿÿ0.8421

.ÿ
.ÿexit

endÿofÿdo-file


.

[Joseph Coveney]

The -vce(cluster )- option only matters if you have multiple observations per cluster, which isn't the case in my example. But it didn't harm anything, either, yielding the conventional sandwich estimates as would be the case with not specifying any -vce()- option with -suest-.

[Maoliang Ling]

Thank you for the help, Joseph. It works!