Why do I get different results with mixed and meglm commands?

[Susanne Urban]

Dear statalists members
I have a dependent variable that is the logarithm of minutes of commuting time from home to school, and I want to do a multilevel regression with random intercept at municipal and neighborhood level. Reading the stata manual https://www.stata.com/manuals13/me.pdf, and previous questions on this forum, it seems as if I should get the same result from using mixed or meglm command, but I dont! I would really appreciate if someone could explain to me wich one I should shoose and why. Below I show an example using ethnicity as a categorical dependent variable, Kommun=municipality, desoomrade=neighborhood (wich are smaller areas within municipalities).
Best wishes
Susanne


1. Results from meglm:
meglm LogCommutingTimeMinutes i.ethnicity|| Kommun :|| desoomrade :

Integration method: mvaghermite Integration pts. = 7

Wald chi2(4) = 139.44
Log likelihood = -26444.087 Prob > chi2 = 0.0000
------------------------------------------------------------------------------------------------
LogCommutingTimeMinutes | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------------------+----------------------------------------------------------------
ethnicity |
2. gn1st | .1769682 .0270593 6.54 0.000 .123933 .2300034
3. gs1st | .1529214 .0180929 8.45 0.000 .1174601 .1883828
4. gn2nd | .060759 .0276835 2.19 0.028 .0065004 .1150175
5. gs2nd | .0978085 .013486 7.25 0.000 .0713763 .1242407
|
_cons | 3.416905 .0062768 544.37 0.000 3.404603 3.429207
-------------------------------+----------------------------------------------------------------
Kommun |
var(_cons)| .1288427 .0491893 .0609664 .2722887
-------------------------------+----------------------------------------------------------------
Kommun>desoomrade |
var(_cons)| .1641032 .0147931 .1375262 .1958163
-------------------------------+----------------------------------------------------------------
var(e.LogCommutingTimeMinutes)| .366027 .0038894 .3584827 .3737301
------------------------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 0.00 Prob > chi2 = 1.0000



Results from mixed:
mixed LogCommutingTimeMinutes i.ethnicity|| Kommun :|| desoomrade :
Mixed-effects ML regression Number of obs = 18,068



Wald chi2(4) = 20.90
Log likelihood = -15305.266 Prob > chi2 = 0.0003

-----------------------------------------------------------------------------------------
LogCommutingTimeMinutes | Coef. Std. Err. z P>|z| [95% Conf. Interval]
------------------------+----------------------------------------------------------------
ethnicity |
2. gn1st | .1095635 .0248939 4.40 0.000 .0607723 .1583547
3. gs1st | .0309414 .0178491 1.73 0.083 -.0040422 .0659251
4. gn2nd | .0064226 .0248707 0.26 0.796 -.0423231 .0551684
5. gs2nd | .0066133 .0139866 0.47 0.636 -.0208 .0340266
|
_cons | 3.568386 .0569499 62.66 0.000 3.456766 3.680006
-----------------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
Kommun: Identity |
var(_cons) | .0774919 .0238052 .0424395 .1414955
-----------------------------+------------------------------------------------
desoomrade: Identity |
var(_cons) | .1263648 .0060637 .1150219 .1388264
-----------------------------+------------------------------------------------
var(Residual) | .2774865 .0030279 .2716149 .283485
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 5912.76 Prob > chi2 = 0.0000

[Joseph Coveney]

I would really appreciate if someone could explain to me wich one I should shoose and why.

Neither, at least not just yet. You're doing something wrong and you don't know which if either is correct. The log-likelihoods aren't even close—it's a guess, but it looks as if you're fitting the same model to two different subsets of the dataset.

There's a lot that you're not showing us. Take a look at the output below, paying attention to the grouping information and other information in the regression output table headers.

.ÿ
.ÿversionÿ16.1

.ÿ
.ÿclearÿ*

.ÿ
.ÿversionÿ16.1

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

.ÿ
.ÿquietlyÿsetÿobsÿ50

.ÿgenerateÿbyteÿcidÿ=ÿ_n

.ÿgenerateÿdoubleÿcid_uÿ=ÿrnormal()

.ÿ
.ÿquietlyÿexpandÿ50

.ÿbysortÿcid:ÿgenerateÿbyteÿdidÿ=ÿ_n

.ÿgenerateÿdoubleÿdid_uÿ=ÿrnormal()

.ÿ
.ÿquietlyÿexpandÿ10

.ÿgenerateÿbyteÿracÿ=ÿruniformint(1,ÿ4)

.ÿ
.ÿgenerateÿdoubleÿoutÿ=ÿcid_uÿ+ÿdid_uÿ+ÿ0ÿ*ÿracÿ+ÿrnormal()

.ÿ
.ÿmeglmÿoutÿi.racÿ||ÿcid:ÿ||ÿdid:ÿ,ÿfamily(gaussian)ÿlink(identity)ÿnolrtestÿnolog

Mixed-effectsÿGLMÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿ25,000
Family:ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿGaussian
Link:ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿidentity

-------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿNo.ÿofÿÿÿÿÿÿÿObservationsÿperÿGroup
ÿGroupÿVariableÿ|ÿÿÿÿÿGroupsÿÿÿÿMinimumÿÿÿÿAverageÿÿÿÿMaximum
----------------+--------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿcidÿ|ÿÿÿÿÿÿÿÿÿ50ÿÿÿÿÿÿÿÿ500ÿÿÿÿÿÿ500.0ÿÿÿÿÿÿÿÿ500
ÿÿÿÿÿÿÿÿÿÿÿÿdidÿ|ÿÿÿÿÿÿ2,500ÿÿÿÿÿÿÿÿÿ10ÿÿÿÿÿÿÿ10.0ÿÿÿÿÿÿÿÿÿ10
-------------------------------------------------------------

Integrationÿmethod:ÿmvaghermiteÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿIntegrationÿpts.ÿÿ=ÿÿÿÿÿÿÿÿÿÿ7

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿWaldÿchi2(3)ÿÿÿÿÿÿ=ÿÿÿÿÿÿÿ5.44
Logÿlikelihoodÿ=ÿ-38714.868ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿÿÿÿÿÿÿ=ÿÿÿÿÿ0.1423
------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿoutÿ|ÿÿÿÿÿÿCoef.ÿÿÿStd.ÿErr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿracÿ|
ÿÿÿÿÿÿÿÿÿÿ2ÿÿ|ÿÿ-.0265771ÿÿÿ.0189788ÿÿÿÿ-1.40ÿÿÿ0.161ÿÿÿÿ-.0637749ÿÿÿÿ.0106207
ÿÿÿÿÿÿÿÿÿÿ3ÿÿ|ÿÿ-.0288326ÿÿÿ.0188656ÿÿÿÿ-1.53ÿÿÿ0.126ÿÿÿÿ-.0658085ÿÿÿÿ.0081432
ÿÿÿÿÿÿÿÿÿÿ4ÿÿ|ÿÿ-.0430679ÿÿÿ.0188684ÿÿÿÿ-2.28ÿÿÿ0.022ÿÿÿÿ-.0800492ÿÿÿ-.0060866
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿ_consÿ|ÿÿÿ.1622936ÿÿÿ.1357167ÿÿÿÿÿ1.20ÿÿÿ0.232ÿÿÿÿ-.1037063ÿÿÿÿ.4282934
-------------+----------------------------------------------------------------
cidÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿvar(_cons)|ÿÿÿ.8932119ÿÿÿ.1828568ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.5979956ÿÿÿÿÿ1.33417
-------------+----------------------------------------------------------------
cid>didÿÿÿÿÿÿ|
ÿÿÿvar(_cons)|ÿÿÿ.9515111ÿÿÿ.0301113ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.8942871ÿÿÿÿ1.012397
-------------+----------------------------------------------------------------
ÿÿÿvar(e.out)|ÿÿÿ1.018295ÿÿÿ.0096006ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.9996505ÿÿÿÿ1.037286
------------------------------------------------------------------------------

.ÿgenerateÿbyteÿmeglm_usedÿ=ÿe(sample)

.ÿ
.ÿmixedÿoutÿi.racÿ||ÿcid:ÿ||ÿdid:ÿ,ÿmlÿnolrtestÿnolog

Mixed-effectsÿMLÿregressionÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿ25,000

-------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿNo.ÿofÿÿÿÿÿÿÿObservationsÿperÿGroup
ÿGroupÿVariableÿ|ÿÿÿÿÿGroupsÿÿÿÿMinimumÿÿÿÿAverageÿÿÿÿMaximum
----------------+--------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿcidÿ|ÿÿÿÿÿÿÿÿÿ50ÿÿÿÿÿÿÿÿ500ÿÿÿÿÿÿ500.0ÿÿÿÿÿÿÿÿ500
ÿÿÿÿÿÿÿÿÿÿÿÿdidÿ|ÿÿÿÿÿÿ2,500ÿÿÿÿÿÿÿÿÿ10ÿÿÿÿÿÿÿ10.0ÿÿÿÿÿÿÿÿÿ10
-------------------------------------------------------------

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿWaldÿchi2(3)ÿÿÿÿÿÿ=ÿÿÿÿÿÿÿ5.44
Logÿlikelihoodÿ=ÿ-38714.868ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿÿÿÿÿÿÿ=ÿÿÿÿÿ0.1423

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿoutÿ|ÿÿÿÿÿÿCoef.ÿÿÿStd.ÿErr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿracÿ|
ÿÿÿÿÿÿÿÿÿÿ2ÿÿ|ÿÿ-.0265771ÿÿÿ.0189788ÿÿÿÿ-1.40ÿÿÿ0.161ÿÿÿÿ-.0637749ÿÿÿÿ.0106207
ÿÿÿÿÿÿÿÿÿÿ3ÿÿ|ÿÿ-.0288326ÿÿÿ.0188654ÿÿÿÿ-1.53ÿÿÿ0.126ÿÿÿÿ-.0658082ÿÿÿÿ.0081429
ÿÿÿÿÿÿÿÿÿÿ4ÿÿ|ÿÿ-.0430679ÿÿÿ.0188683ÿÿÿÿ-2.28ÿÿÿ0.022ÿÿÿÿ-.0800491ÿÿÿ-.0060866
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿ_consÿ|ÿÿÿ.1622936ÿÿÿ.1357169ÿÿÿÿÿ1.20ÿÿÿ0.232ÿÿÿÿ-.1037066ÿÿÿÿ.4282937
------------------------------------------------------------------------------

------------------------------------------------------------------------------
ÿÿRandom-effectsÿParametersÿÿ|ÿÿÿEstimateÿÿÿStd.ÿErr.ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-----------------------------+------------------------------------------------
cid:ÿIdentityÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿvar(_cons)ÿ|ÿÿÿ.8932144ÿÿÿ.1828574ÿÿÿÿÿÿ.5979972ÿÿÿÿ1.334173
-----------------------------+------------------------------------------------
did:ÿIdentityÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿvar(_cons)ÿ|ÿÿÿ.9515128ÿÿÿ.0301113ÿÿÿÿÿÿ.8942887ÿÿÿÿ1.012399
-----------------------------+------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿvar(Residual)ÿ|ÿÿÿ1.018294ÿÿÿ.0096006ÿÿÿÿÿÿ.9996504ÿÿÿÿ1.037286
------------------------------------------------------------------------------

.ÿgenerateÿbyteÿmixed_usedÿ=ÿe(sample)

.ÿ
.ÿassertÿmeglm_usedÿ==ÿmixed_used

.ÿ
.ÿexit

endÿofÿdo-file


.

[Susanne Urban]

Thank you Joseph for the response! I run the models again following the specified commands that you used, but with the same result. The number of observations is exactly the same. I am wondering if it matters that the "desoomrade" (neighbourhood) is a string variable containing both numbers and letters. Kommun (Municipality) is integer and contains only numbers. The data set is also not sorted.

Here is the output (sorry I dont know how to paste text from the log file so neatly as you did):

. meglm LogCommutingTimeMinutes i.ethnicity || Kommun :|| desoomrade : , family (gaussian) link (identity) nolrtest nolog

Mixed-effects GLM Number of obs = 18,068
Family: Gaussian
Link: identity

-------------------------------------------------------------
| No. of Observations per Group
Group Variable | Groups Minimum Average Maximum
----------------+--------------------------------------------
Kommun | 26 79 694.9 6,298
desoomrade | 1,283 1 14.1 43
-------------------------------------------------------------

Integration method: mvaghermite Integration pts. = 7

Wald chi2(4) = 139.44
Log likelihood = -26444.087 Prob > chi2 = 0.0000
------------------------------------------------------------------------------------------------
LogCommutingTimeMinutes | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------------------+----------------------------------------------------------------
ethnicity |
2. gn1st | .1769682 .0270593 6.54 0.000 .123933 .2300034
3. gs1st | .1529214 .0180929 8.45 0.000 .1174601 .1883828
4. gn2nd | .060759 .0276835 2.19 0.028 .0065004 .1150175
5. gs2nd | .0978085 .013486 7.25 0.000 .0713763 .1242407
|
_cons | 3.416905 .0062768 544.37 0.000 3.404603 3.429207
-------------------------------+----------------------------------------------------------------
Kommun |
var(_cons)| .1288427 .0491893 .0609664 .2722887
-------------------------------+----------------------------------------------------------------
Kommun>desoomrade |
var(_cons)| .1641032 .0147931 .1375262 .1958163
-------------------------------+----------------------------------------------------------------
var(e.LogCommutingTimeMinutes)| .366027 .0038894 .3584827 .3737301
------------------------------------------------------------------------------------------------

. generate byte meglm_used =e(sample)

. mixed LogCommutingTimeMinutes i.ethnicity || Kommun :|| desoomrade : , ml nolrtest nolog

Mixed-effects ML regression Number of obs = 18,068

-------------------------------------------------------------
| No. of Observations per Group
Group Variable | Groups Minimum Average Maximum
----------------+--------------------------------------------
Kommun | 26 79 694.9 6,298
desoomrade | 1,283 1 14.1 43
-------------------------------------------------------------

Wald chi2(4) = 20.90
Log likelihood = -15305.266 Prob > chi2 = 0.0003

-----------------------------------------------------------------------------------------
LogCommutingTimeMinutes | Coef. Std. Err. z P>|z| [95% Conf. Interval]
------------------------+----------------------------------------------------------------
ethnicity |
2. gn1st | .1095635 .0248939 4.40 0.000 .0607723 .1583547
3. gs1st | .0309414 .0178491 1.73 0.083 -.0040422 .0659251
4. gn2nd | .0064226 .0248707 0.26 0.796 -.0423231 .0551684
5. gs2nd | .0066133 .0139866 0.47 0.636 -.0208 .0340266
|
_cons | 3.568386 .0569499 62.66 0.000 3.456766 3.680006
-----------------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
Kommun: Identity |
var(_cons) | .0774919 .0238052 .0424395 .1414955
-----------------------------+------------------------------------------------
desoomrade: Identity |
var(_cons) | .1263648 .0060637 .1150219 .1388264
-----------------------------+------------------------------------------------
var(Residual) | .2774865 .0030279 .2716149 .283485
------------------------------------------------------------------------------

[Susanne Urban]

Hello again
I transformed the neighbourhoodvariable to a numeric variable. This made the result from the meglm command change sligtly, but the result from mixed is still the same. Now I see that that the data set are sorted after municipality and neigbourhood. Now I dont know how to proceed.

See output below.


encode desoomrade , generate(DESONMBR)

. meglm LogCommutingTimeMinutes i.ethnicity || Kommun :|| DESONMBR : , family (gaussian) link (identity) nolrtest nolog

Mixed-effects GLM Number of obs = 18,068
Family: Gaussian
Link: identity

-------------------------------------------------------------
| No. of Observations per Group
Group Variable | Groups Minimum Average Maximum
----------------+--------------------------------------------
Kommun | 26 79 694.9 6,298
DESONMBR | 1,283 1 14.1 43
-------------------------------------------------------------

Integration method: mvaghermite Integration pts. = 7

Wald chi2(4) = 138.40
Log likelihood = -24838.226 Prob > chi2 = 0.0000
-----------------------------------------------------------------------------------------------
LogCommutingTimeMinutes | Coef. Std. Err. z P>|z| [95% Conf. Interval]
------------------------------+----------------------------------------------------------------
|
ethnicity |
2. gn1st | .1780269 .0264245 6.74 0.000 .1262359 .2298179
3. gs1st | .1516614 .0177609 8.54 0.000 .1168506 .1864721
4. gn2nd | .0454636 .027009 1.68 0.092 -.007473 .0984003
5. gs2nd | .0938808 .0132683 7.08 0.000 .0678755 .1198861
|
_cons | 3.437338 .0066181 519.39 0.000 3.424367 3.450309
------------------------------+----------------------------------------------------------------
Kommun |
var(_cons)| .1099987 .0389791 .0549232 .2203024
------------------------------+----------------------------------------------------------------
Kommun>DESONMBR |
var(_cons)| .1621076 .0125959 .139208 .1887742
------------------------------+----------------------------------------------------------------
var(e.LogCommutingTimeMinutes)| .3461477 .0036911 .3389883 .3534583
-----------------------------------------------------------------------------------------------

. generate byte meglm_used =e(sample)

.
.
.
. mixed LogCommutingTimeMinutes i.ethnicity || Kommun :|| DESONMBR : , ml nolrtest nolog

Mixed-effects ML regression Number of obs = 18,068

-------------------------------------------------------------
| No. of Observations per Group
Group Variable | Groups Minimum Average Maximum
----------------+--------------------------------------------
Kommun | 26 79 694.9 6,298
DESONMBR | 1,283 1 14.1 43
-------------------------------------------------------------

Wald chi2(4) = 20.90
Log likelihood = -15305.266 Prob > chi2 = 0.0003

-----------------------------------------------------------------------------------------
LogCommutingTimeMinutes | Coef. Std. Err. z P>|z| [95% Conf. Interval]
------------------------+----------------------------------------------------------------
ethnicity |
2. gn1st | .1095635 .0248939 4.40 0.000 .0607723 .1583547
3. gs1st | .0309414 .0178491 1.73 0.083 -.0040422 .0659251
4. gn2nd | .0064226 .0248707 0.26 0.796 -.0423231 .0551684
5. gs2nd | .0066133 .0139866 0.47 0.636 -.0208 .0340266
|
_cons | 3.568386 .0569499 62.66 0.000 3.456766 3.680006
-----------------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
Kommun: Identity |
var(_cons) | .0774919 .0238052 .0424395 .1414955
-----------------------------+------------------------------------------------
DESONMBR: Identity |
var(_cons) | .1263648 .0060637 .1150219 .1388264
-----------------------------+------------------------------------------------
var(Residual) | .2774865 .0030279 .2716149 .283485
------------------------------------------------------------------------------

[Joseph Coveney]

You might have an ugly likelihood surface with multiple modes; meglm found one local maximum and mixed, starting with expectation maximization, homed in on another. You could try using one model's fitted parameter estimates as starting values for the other and vice versa. You could try profiling the likelihood over a range of parameter values, but you've got a lot of parameters and a feasible grid might be too coarse. Or you could chance it and go with the fit that has the greater likelihood in the hope that there is a global maximum and that that's it.

[Susanne Urban]

Thank you. Since I get the same result with meglm and mixed when I only use neighbourhood as a level, I think that will skip using municipality as a level, and only add municipalities as explaining variables.