Hello,
I have an unbalanced panel of country-level data: 909 observations from 180 countries. Data are available every two years over the period 2014-2024 (i.e., T=6). The data looks as follows:
I have run a two-level random intercepts model that includes two time-constant variables (income group and region_id), with standard errors clustered at the country level as shown below. This is my main model.
I wanted to test the sensitivity of my results to accounting for higher levels of clustering in the data. To this end, I have played around with the flexibility afforded by the mixed model and add another random intercept to account for additional levels of nesting. Specifically, I ran two three-level (random intercepts) models. In the first, I added a random intercepts at both the region level and the country-within-region level (I excluded the region variable from the “fixed effects” component of this model for this specification). Results are shown below.
Second, I ran model with a random intercept at both the income-group level and the country-within-income-group level. In this specification, I removed the income group variable from the "fixed" effects part of the model.
I am aware that I can’t specify vce(cluster country) in these three-level random intercept models (Stata shouts!), but I could have used “, robust”. However, I learned from the mixed helpfile that robust variances are clustered at the highest level in the multilevel model. This is no good for my data as my number of groups at the highest level is very small (3 income groups and 6 regions). So, at present, I have just used the default standard errors (as above).
Looking at the AIC and BIC statistics for my two-level and both of my three-level models, it is clear that the model with the random intercept for both income and country performs the worst; while the comparison between the two-level random intercept and three-level random intercept (|| region_id: || country: ) is ambiguous on which is better: AIC is lower for the two-level model; but BIC is lower for the three-level. I can't do LR test because Stata tells me the two-level and three-level random intercepts models are not nested.
Fundamentally, I am concerned about whether it is okay to run a three-level model at all when the number of groups at the highest level is so low. The mixed helpfile talks about specifying , reml when the number of groups is small. However, I find this discussion confusing because I don’t believe that I have a small number of groups at level 2; but I certainly do at level 3. Does a small number of groups at the highest level only still render the , reml option more appropriate than the default ML results?
Ultimately, my question has two parts. (1) Is it it even sensible to add the third level of clustering given that the number of groups is so small? (2) If it is sensible to have this kind of robustness check, is it better to report the standard ML estimates with default SEs, or is it more correct to use the, reml specification?
Thank you!
Sam
I have an unbalanced panel of country-level data: 909 observations from 180 countries. Data are available every two years over the period 2014-2024 (i.e., T=6). The data looks as follows:
Code:
* Example generated by -dataex-. For more info, type help dataex
clear
input double price_dispersion_use float TS_ce byte POWE double unem float(income region_id) int year
44.44444444444444 1 18 18.055 3 4 2014
56.666666666666664 1 18 15.418 3 4 2016
62.5 1 18 12.304 3 4 2018
60.60606060606061 1 19 11.69 3 4 2020
60 1 19 10.137 3 4 2022
50 1 19 10.25 3 4 2024
33.33333333333333 4 12 10.207 3 1 2014
35.714285714285715 6 12 10.202 3 1 2016
15 6 13 12.137 3 1 2018
50 6 13 14.057 3 1 2020
48.57142857142857 3 13 12.346 3 1 2022
42.857142857142854 6 13 11.427 3 1 2024
72.85714285714285 4 11 5.3 1 4 2014
72.85714285714285 4 11 3.3 1 4 2016
77.77777777777777 1 11 1.8 1 4 2018
68.44993141289439 1 11 2.9 1 4 2020
69.86301369863014 1 11 2.1 1 4 2022
59.09090909090908 1 11 1.4 1 4 2024
25 2 13 16.69 3 1 2020
25 2 13 14.602 3 1 2022
28.57142857142857 2 13 14.464 3 1 2024
40 2 16 7.423938916311391 1 2 2024
41.935483870967744 2 18 7.268 3 2 2014
37.75 2 18 8.085 3 2 2016
45.34920634920635 2 18 9.22 3 2 2018
18.726114649681527 2 18 11.461 3 2 2020
13.384615384615383 2 18 6.805 3 2 2022
40.055248618784525 2 18 7.876 3 2 2024
30 4 13 17.498 3 4 2014
26.666666666666668 4 14 17.617 3 4 2016
42.857142857142854 2 14 18.966 3 4 2018
47.5 1 18 18.175 3 4 2020
48.23529411764706 1 19 13.379 3 4 2022
35 1 19 13.329 3 4 2024
78.93318965517241 1 19 6.078 1 6 2014
73.84341637010677 1 19 5.711 1 6 2016
82.34126984126985 1 19 5.3 1 6 2018
71.02189781021899 4 18 6.456 1 6 2020
68.45524542829644 4 18 3.728 1 6 2022
60.58098915241773 4 18 4.072 1 6 2024
80.61224489795919 3 13 5.674 1 4 2014
80 3 15 6.064 1 4 2016
80 3 15 4.933 1 4 2018
82.45614035087719 3 17 5.201 1 4 2020
68.35820895522387 3 16 4.992 1 4 2022
83.07692307692308 3 16 5.439 1 4 2024
24 4 12 4.91 3 4 2014
56.25 1 13 5 3 4 2016
23.076923076923077 4 14 4.9 3 4 2018
47.05882352941177 4 14 7.24 3 4 2020
55.55555555555556 4 14 5.65 3 4 2022
48.88888888888889 1 14 5.594 3 4 2024
48.658536585365916 1 10 13.8 1 2 2014
40.22346368715088 1 10 12.7 1 2 2016
61.08949416342412 1 10 12.027 1 2 2020
81.76100628930817 1 10 8.463 1 2 2024
40 2 15 1.2 1 3 2018
34.78260869565218 2 15 1.781 1 3 2020
58.333333333333336 2 15 1.326 1 3 2022
58.333333333333336 2 15 1.102 1 3 2024
15.789473684210526 5 13 4.416 3 5 2014
22.727272727272727 5 15 4.35 3 5 2016
33.33333333333333 5 15 4.407 3 5 2018
32.142857142857146 5 15 5.436 3 5 2020
25 5 15 4.593 3 5 2022
30.864197530864196 5 15 4.68 3 5 2024
79.32850559578671 1 13 12.17 1 2 2014
81.25 1 13 8.247 1 2 2016
45.23433385992628 1 16 8.322 1 2 2018
78.84615384615384 1 16 8.365 1 2 2022
79.98999499749875 1 15 7.529 1 2 2024
35.625 4 14 5.902 3 4 2014
31.914893617021278 4 15 5.844 3 4 2016
30.645161290322577 4 15 4.763 3 4 2018
25.71428571428572 4 15 4.049 3 4 2020
23.958333333333332 4 15 3.574 3 4 2022
47.26027397260275 4 15 3.361 3 4 2024
80.82901554404145 3 14 8.523 1 4 2014
81.64556962025317 3 15 7.83 1 4 2016
83.33333333333334 3 15 5.941 1 4 2018
85.29411764705883 3 15 5.545 1 4 2020
72.5 3 15 5.57 1 4 2022
65.21739130434783 3 15 5.488 1 4 2024
41.66666666666667 1 10 7 3 2 2016
40 1 10 7.896 3 2 2018
50 1 10 10.784 3 2 2020
50.391644908616186 1 10 8.763 3 2 2022
55.55555555555556 1 10 7 3 2 2024
20 2 15 1.784 3 1 2016
22.5 2 15 1.41 3 1 2018
47.368421052631575 2 15 1.502 3 1 2020
50 2 17 1.722 3 1 2024
35.714285714285715 2 16 2.021 3 2 2014
24 5 16 3.498 3 2 2016
40 4 14 3.519 3 2 2018
50 4 18 3.552 3 2 2022
36.40776699029126 4 16 3.091 3 2 2024
67.44186046511628 3 10 27.517 3 4 2014
60.416666666666664 3 10 25.408 3 4 2016
68.96551724137932 3 10 18.4 3 4 2018
end
label values TS_ce TS_ce_l
label def TS_ce_l 1 "1. specific uniform", modify
label def TS_ce_l 2 "2. advalorem uniform", modify
label def TS_ce_l 3 "3. mixed uniform", modify
label def TS_ce_l 4 "4. specific_tiered", modify
label def TS_ce_l 5 "5. advalorem tiered", modify
label def TS_ce_l 6 "6. mixed tiered", modify
label values income income_l
label def income_l 1 "1. High", modify
label def income_l 3 "3. Middle", modify
label values region_id region_id_l
label def region_id_l 1 "1. AFR", modify
label def region_id_l 2 "2. AMR", modify
label def region_id_l 3 "3. EMR", modify
label def region_id_l 4 "4. EUR", modify
label def region_id_l 5 "5. SEA", modify
label def region_id_l 6 "6. WPR", modify
Code:
. mixed price_dispersion_use i.TS_ce unem POWE i.income i.region_id i.year || country: , vce(cluster country)
Performing EM optimization ...
Performing gradient-based optimization:
Iteration 0: Log pseudolikelihood = -3684.8293
Iteration 1: Log pseudolikelihood = -3684.8293
Computing standard errors ...
Mixed-effects regression Number of obs = 909
Group variable: country Number of groups = 180
Obs per group:
min = 1
avg = 5.0
max = 6
Wald chi2(19) = 485.11
Log pseudolikelihood = -3684.8293 Prob > chi2 = 0.0000
(Std. err. adjusted for 180 clusters in country)
---------------------------------------------------------------------------------------
| Robust
price_dispersion_use | Coefficient std. err. z P>|z| [95% conf. interval]
----------------------+----------------------------------------------------------------
TS_ce |
2. advalorem uniform | -12.79406 3.177566 -4.03 0.000 -19.02198 -6.56615
3. mixed uniform | -.1719755 2.991038 -0.06 0.954 -6.034301 5.69035
4. specific_tiered | -12.01681 3.243756 -3.70 0.000 -18.37445 -5.659162
5. advalorem tiered | -14.82506 4.008941 -3.70 0.000 -22.68244 -6.967678
6. mixed tiered | -12.34974 4.721909 -2.62 0.009 -21.60451 -3.094967
|
unem | -.2082926 .1604439 -1.30 0.194 -.5227567 .1061716
POWE | .3630399 .3414757 1.06 0.288 -.3062401 1.03232
|
income |
2. Low | -21.74503 3.60037 -6.04 0.000 -28.80162 -14.68843
3. Middle | -15.07393 2.614105 -5.77 0.000 -20.19749 -9.950383
|
region_id |
2. AMR | 7.216272 3.616945 2.00 0.046 .1271903 14.30535
3. EMR | -3.9681 4.507714 -0.88 0.379 -12.80306 4.866857
4. EUR | 9.939152 4.166454 2.39 0.017 1.773052 18.10525
5. SEA | -6.913486 5.509264 -1.25 0.210 -17.71144 3.884473
6. WPR | 6.218809 4.986658 1.25 0.212 -3.55486 15.99248
|
year |
2016 | 1.136412 1.125618 1.01 0.313 -1.069759 3.342582
2018 | 1.815113 1.429828 1.27 0.204 -.9872978 4.617523
2020 | 2.493011 1.442856 1.73 0.084 -.3349358 5.320958
2022 | 2.409131 1.480869 1.63 0.104 -.4933187 5.311581
2024 | 3.506218 1.70807 2.05 0.040 .1584622 6.853975
|
_cons | 60.35922 5.927152 10.18 0.000 48.74222 71.97623
---------------------------------------------------------------------------------------
------------------------------------------------------------------------------
| Robust
Random-effects parameters | Estimate std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
country: Identity |
var(_cons) | 146.1695 23.78864 106.2493 201.0885
-----------------------------+------------------------------------------------
var(Residual) | 135.2595 12.23187 113.2899 161.4894
------------------------------------------------------------------------------
.
end of do-file
. estat ic
Akaike's information criterion and Bayesian information criterion
-----------------------------------------------------------------------------
Model | N ll(null) ll(model) df AIC BIC
-------------+---------------------------------------------------------------
. | 909 . -3684.829 22 7413.659 7519.53
-----------------------------------------------------------------------------
Note: BIC uses N = number of observations. See [R] IC note.
. estat icc
Residual intraclass correlation
------------------------------------------------------------------------------
Level | ICC Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
country | .5193833 .0510441 .4198939 .6173588
------------------------------------------------------------------------------
Code:
. mixed price_dispersion_use i.TS_ce unem POWE i.income i.year || region_id: || country:
Performing EM optimization ...
Performing gradient-based optimization:
Iteration 0: Log likelihood = -3692.3964
Iteration 1: Log likelihood = -3692.3964
Computing standard errors ...
Mixed-effects ML regression Number of obs = 909
Grouping information
-------------------------------------------------------------
| No. of Observations per group
Group variable | groups Minimum Average Maximum
----------------+--------------------------------------------
region_id | 6 34 151.5 302
country | 180 1 5.0 6
-------------------------------------------------------------
Wald chi2(14) = 175.87
Log likelihood = -3692.3964 Prob > chi2 = 0.0000
---------------------------------------------------------------------------------------
price_dispersion_use | Coefficient Std. err. z P>|z| [95% conf. interval]
----------------------+----------------------------------------------------------------
TS_ce |
2. advalorem uniform | -13.26848 2.33598 -5.68 0.000 -17.84692 -8.69004
3. mixed uniform | .0820774 2.364169 0.03 0.972 -4.551609 4.715763
4. specific_tiered | -12.63683 2.414591 -5.23 0.000 -17.36934 -7.904314
5. advalorem tiered | -15.75917 4.300878 -3.66 0.000 -24.18874 -7.329606
6. mixed tiered | -12.74558 3.013626 -4.23 0.000 -18.65218 -6.83898
|
unem | -.2302968 .1583598 -1.45 0.146 -.5406762 .0800826
POWE | .3434487 .2773258 1.24 0.216 -.2000998 .8869972
|
income |
2. Low | -23.69504 4.213635 -5.62 0.000 -31.95362 -15.43647
3. Middle | -15.98282 2.564233 -6.23 0.000 -21.00862 -10.95701
|
year |
2016 | 1.107993 1.366686 0.81 0.418 -1.570663 3.78665
2018 | 1.779171 1.420145 1.25 0.210 -1.004262 4.562603
2020 | 2.426264 1.424773 1.70 0.089 -.3662399 5.218767
2022 | 2.326682 1.457251 1.60 0.110 -.5294783 5.182841
2024 | 3.404323 1.47028 2.32 0.021 .522627 6.286019
|
_cons | 64.83817 5.434994 11.93 0.000 54.18578 75.49057
---------------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
region_id: Identity |
var(_cons) | 19.51945 18.35721 3.08993 123.3066
-----------------------------+------------------------------------------------
country: Identity |
var(_cons) | 153.1753 19.94595 118.672 197.7103
-----------------------------+------------------------------------------------
var(Residual) | 135.2647 7.100675 122.0396 149.9229
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 361.60 Prob > chi2 = 0.0000
Note: LR test is conservative and provided only for reference.
. estat ic
Akaike's information criterion and Bayesian information criterion
-----------------------------------------------------------------------------
Model | N ll(null) ll(model) df AIC BIC
-------------+---------------------------------------------------------------
. | 909 . -3692.396 18 7420.793 7507.415
-----------------------------------------------------------------------------
Note: BIC uses N = number of observations. See [R] IC note.
. estat icc
Residual intraclass correlation
------------------------------------------------------------------------------
Level | ICC Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
region_id | .0633832 .0563688 .0104144 .3032105
country|region_id | .5607712 .0399014 .4817034 .6368698
------------------------------------------------------------------------------
.
end of do-file
Code:
. mixed price_dispersion_use i.TS_ce unem POWE i.region_id i.year || income: || id:
Performing EM optimization ...
Performing gradient-based optimization:
Iteration 0: Log likelihood = -3690.9924
Iteration 1: Log likelihood = -3690.9924
Computing standard errors ...
Mixed-effects ML regression Number of obs = 909
Grouping information
-------------------------------------------------------------
| No. of Observations per group
Group variable | groups Minimum Average Maximum
----------------+--------------------------------------------
income | 3 98 303.0 498
id | 180 1 5.0 6
-------------------------------------------------------------
Wald chi2(17) = 141.36
Log likelihood = -3690.9924 Prob > chi2 = 0.0000
---------------------------------------------------------------------------------------
price_dispersion_use | Coefficient Std. err. z P>|z| [95% conf. interval]
----------------------+----------------------------------------------------------------
TS_ce |
2. advalorem uniform | -12.84855 2.354166 -5.46 0.000 -17.46263 -8.234469
3. mixed uniform | -.1261898 2.394097 -0.05 0.958 -4.818533 4.566154
4. specific_tiered | -12.12209 2.432416 -4.98 0.000 -16.88953 -7.35464
5. advalorem tiered | -14.97557 4.313191 -3.47 0.001 -23.42927 -6.521872
6. mixed tiered | -12.39842 3.032808 -4.09 0.000 -18.34262 -6.454229
|
unem | -.2092789 .1600007 -1.31 0.191 -.5228745 .1043167
POWE | .3769979 .27681 1.36 0.173 -.1655398 .9195356
|
region_id |
2. AMR | 8.073306 3.468474 2.33 0.020 1.275223 14.87139
3. EMR | -3.270375 4.061883 -0.81 0.421 -11.23152 4.690769
4. EUR | 11.1213 3.61721 3.07 0.002 4.0317 18.2109
5. SEA | -6.254756 5.694607 -1.10 0.272 -17.41598 4.906468
6. WPR | 7.118405 3.968743 1.79 0.073 -.6601881 14.897
|
year |
2016 | 1.127195 1.366599 0.82 0.409 -1.551289 3.805679
2018 | 1.795766 1.42039 1.26 0.206 -.9881475 4.57968
2020 | 2.475048 1.424955 1.74 0.082 -.3178118 5.267907
2022 | 2.38783 1.45853 1.64 0.102 -.470836 5.246497
2024 | 3.484633 1.47202 2.37 0.018 .5995274 6.369738
|
_cons | 47.50363 6.870172 6.91 0.000 34.03834 60.96892
---------------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
income: Identity |
var(_cons) | 70.4047 64.12905 11.81064 419.6911
-----------------------------+------------------------------------------------
id: Identity |
var(_cons) | 149.3076 19.26522 115.9447 192.2706
-----------------------------+------------------------------------------------
var(Residual) | 135.2595 7.097993 122.0393 149.912
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 435.46 Prob > chi2 = 0.0000
Note: LR test is conservative and provided only for reference.
. estat ic
Akaike's information criterion and Bayesian information criterion
-----------------------------------------------------------------------------
Model | N ll(null) ll(model) df AIC BIC
-------------+---------------------------------------------------------------
. | 909 . -3690.992 21 7423.985 7525.044
-----------------------------------------------------------------------------
Note: BIC uses N = number of observations. See [R] IC note.
. estat icc
Residual intraclass correlation
------------------------------------------------------------------------------
Level | ICC Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
income | .1983388 .1454561 .0395546 .5977974
id|income | .618957 .0728953 .469873 .7485514
------------------------------------------------------------------------------
.
end of do-file
Looking at the AIC and BIC statistics for my two-level and both of my three-level models, it is clear that the model with the random intercept for both income and country performs the worst; while the comparison between the two-level random intercept and three-level random intercept (|| region_id: || country: ) is ambiguous on which is better: AIC is lower for the two-level model; but BIC is lower for the three-level. I can't do LR test because Stata tells me the two-level and three-level random intercepts models are not nested.
Fundamentally, I am concerned about whether it is okay to run a three-level model at all when the number of groups at the highest level is so low. The mixed helpfile talks about specifying , reml when the number of groups is small. However, I find this discussion confusing because I don’t believe that I have a small number of groups at level 2; but I certainly do at level 3. Does a small number of groups at the highest level only still render the , reml option more appropriate than the default ML results?
Ultimately, my question has two parts. (1) Is it it even sensible to add the third level of clustering given that the number of groups is so small? (2) If it is sensible to have this kind of robustness check, is it better to report the standard ML estimates with default SEs, or is it more correct to use the, reml specification?
Thank you!
Sam
Comment