Seeking Advice: Modeling Regional Time-Varying Effects in a Multilevel Framework

[Doris Rivera]

Dear all, I am working with a panel dataset for regions, and I want to predict the effect of x1 over y1 for different regions (id) and study how this is changing over time (the plan is to use this in a second model). For that I use a multilevel model with a random coefficient for x1 at the regional level (see below), in which x1 is demeaned to study within variation. After that, I use the fixed and the predicted random effect for x1 to come with the desired measure. However, that gives a time invariant prediction (fixed for each region), since none of the fixed x1 coefficient or random x1 coefficient vary with time.

So, I am asking advice on how to let this effect to change with time. These are my current ideas:

1. I could go for a multilevel model with random coefficient for x1, and repeat that model for each year. However, that would be a between regions analysis in each year (instead of a within one, right?).
2. I could use a random coefficient for x1, not only for regional level, but also for time level. Even though my years are from 2006-2022, in fact it would less since I use 5 years lagged indep. variables, and I am not sure if with these number of years (11), a random coefficient for time lelvel would remain its properties.
3. I could interact in the fixed part of the model x1*i.years and put each interaction as random coefficients for the regional level (but would be an endless estimation since I am talking about 11 coefficients in the random part). I tried this for more than a week and nothing came out.

Are any of those ideas, a plausible solution to my problem?

If not, is there any other way (not necessarily within the multilevel framework) I could obtain this effect differently for regions in each year? Below there is a small data example.

Thanks in advance.

Iteration 122: log pseudolikelihood =  -10258.05  
Iteration 123: log pseudolikelihood =  -10258.05  

Mixed-effects nbinomial regression              Number of obs     =      3,006
Overdispersion:            mean
Group variable:              id                 Number of groups  =        254

                                                Obs per group:
                                                              min =          4
                                                              avg =       11.8
                                                              max =         12

Integration method: mvaghermite                 Integration pts.  =         15

                                                Wald chi2(17)     =      67.13
Log pseudolikelihood =  -10258.05               Prob > chi2       =     0.0000
                                       (Std. Err. adjusted for 254 clusters in id)
----------------------------------------------------------------------------------
                 |               Robust
              y1 |     exp(b)   Std. Err.      z    P>|z|     [95% Conf. Interval]
-----------------+----------------------------------------------------------------
              x1 |
             L5. |   1.044958   .0216087     2.13   0.033     1.003453     1.08818
                 |
              x2 |
             L5. |   1.000673   .0130966     0.05   0.959      .975331    1.026674
                 |
              x3 |
             L5. |   .9796073   .0149627    -1.35   0.177     .9507155    1.009377
                 |
              x4 |
             L5. |   .9918395   .0124067    -0.66   0.512     .9678185    1.016457
                 |
              x5 |
             L5. |   1.003104   .0163533     0.19   0.849     .9715582    1.035673
                 |
              x6 |
             L5. |   1.016806   .0103273     1.64   0.101     .9967654     1.03725
                 |
            year |
           2006  |          1  (empty)
           2007  |          1  (empty)
           2008  |          1  (empty)
           2009  |          1  (empty)
           2010  |          1  (empty)
           2011  |   .8308796   .0507725    -3.03   0.002     .7370956    .9365962
           2012  |   .9199021   .0503215    -1.53   0.127      .826377    1.024012
           2013  |   .9671356   .0469571    -0.69   0.491     .8793447    1.063691
           2014  |   .9498539   .0480901    -1.02   0.310     .8601246    1.048944
           2015  |   .9521785   .0425776    -1.10   0.273     .8722802    1.039395
           2016  |    .945086   .0369091    -1.45   0.148     .8754448    1.020267
           2017  |   .9209528   .0349554    -2.17   0.030     .8549278    .9920769
           2018  |   .9709001   .0346993    -0.83   0.409     .9052181    1.041348
           2019  |   .9962854   .0321001    -0.12   0.908     .9353156    1.061229
           2020  |   1.015637   .0270911     0.58   0.561     .9639033    1.070147
           2021  |   1.015463   .0202364     0.77   0.441     .9765646     1.05591
           2022  |          1  (omitted)
                 |
           _cons |   17.93158    1.77702    29.13   0.000     14.76606    21.77573
-----------------+----------------------------------------------------------------
        /lnalpha |  -3.738889   .1233591                     -3.980668   -3.497109
-----------------+----------------------------------------------------------------
id               |
       var(L5.x1)|   .0274824   .0060578                      .0178414    .0423332
       var(_cons)|   2.360432   .2013258                      1.997059    2.789923
-----------------+----------------------------------------------------------------
id               |
 cov(L5.x1,_cons)|  -.0466977   .0299401    -1.56   0.119    -.1053793    .0119839
----------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation.

* Example generated by -dataex-. For more info, type help dataex
clear
input int year float(id y1 x1 x2 x3 x4 x5 x6)
2006 1  5  -1.7775646   2.723891  -1.2291875  -1.7063985  -1.5368623  2.4518554
2007 1  5  -1.2202544    1.45774           .  -1.0336447  -1.3795106   .8769412
2008 1  8   -.7455087   .3915077  -1.3092277   -.9709182  -1.1481309  -.8411472
2009 1  9   -.8487142  -1.074562           .  -1.4221123   -.9250329 -2.1297135
2010 1  9    -.807432 -1.2744807  -1.2291875  -1.0288475   -.7678502  -.9843215
2011 1  7   -.0849928  -.6080854  -.50882643   -.4270809   -.5917379 .017897163
2012 1  9    .1420595 -.07496897  -.26870614  -.09886736  -.41342854  .59059274
2013 1 11    -.642303  .12494913 -.028585836    .1738537  -.27517536 -.12527741
2014 1 12   -.1881984  -.2748877   -.1886661    .2016169    -.134725  -.6979733
2015 1  7  -.06435169  .05830998  -.10862587   .29352987   .16916154  -.2684513
2016 1  8   .05949503   .4581468   -.1886661    .5043921     .472203 .017897163
2017 1  4   .03885391  -.6080854   -.3487464   1.2704383    .6129914 -.12527741
2018 1  9    .7819343   .3915077   1.0119358   1.4209753    .7389066   .4474188
2019 1 10    .4103941 -.07496897   1.0919758     1.55902    .8878076   .8769412
2020 1 13   1.6695024  -.8746439   1.4921764  -.11110864   1.1054972  -.2684513
2021 1  9   1.6695024  -.7413643    1.812337   .25691786   1.3715247  .16107106
2022 1  6    1.607579          .           .    1.118244   1.8143405          .
2006 2 20  -1.9668324   2.536081   -1.357809  -2.0436459   -1.413372  2.6546836
2007 2 29   -.6613229   2.216574   -1.228494   -.9597974    -1.22389  1.2620628
2008 2 36   -.5426402  .21966097  -1.1638364   -.1650967  -1.0422335  -1.000947
2009 2 35     .849184  -.8187343  -1.0345212   -.9916811   -.9110399  -1.349102
2010 2 41   -.2297495  -.6589814   -.9698637   -.9121222   -.8287947  -.8268695
2011 2 52  -.05712012  -.6589814   -.3879454  -.28500873   -.7207753  -.8268695
2012 2 39   .09393056 -.09984569  -.12931518   -.3549016   -.6056876    .565752
2013 2 58  -.58579755 -.25959864  -.12931518   -.3371715   -.4664412  .39167455
2014 2 55   -.2297495  -.6589814   .06465752   -.1212194   -.2362404  -.3046362
2015 2 57   -.5642189  -.4193516   .25863037  .016521487   .11909527 -.13055874
2016 2 56   -.1758028 -.17972177   .06465752    .0959413    .4848315    .565752
2017 2 49  -.29448548  -.8187343    .4526029    .8584786    .6991033  .04351875
2018 2 59   -.9094776 -.01996883   1.2931515   1.3341718    .8478919  .04351875
2019 2 71    .3852426  -.6589814    1.228494   1.8900775   1.0198293  -.6527912
2020 2 51    .9678667 -.17972177   1.4224665  -.17276284   1.1881568  -.8268695
2021 2 64    1.658384   .4592908   1.6164395    .6946668   1.3712775  .39167455
2022 2 57   2.2625868          .           .   1.4535483   1.7182826          .
2006 3 22   -1.461945   2.505153   -1.311473    .6374172  -1.3114537   2.560008
2007 3 26   -.7937856  1.9937345  -1.2502964    1.214219  -1.2218945  2.2027972
2008 3 25   -.5662836  .20376867  -1.2502964    1.229118  -1.1317818  -.4167453
2009 3 35   -.3920703 -.17979567  -1.1891197    .6612731   -1.042213  -.4167453
2010 3 38   -.5990766 -.05194079    -1.00559    .9011381   -.9338182   -.059535
2011 3 44   -.6298202   .0119865   -.3326478   1.2566856   -.8045746   .1786053
2012 3 50   -.4720034   -.627287  -.21029447    .8686113     -.62142  -.2976753
2013 3 63   -.6195723  -.9469237   .27911794   .26373556   -.3844877  -.4167453
2014 3 54   -.6072749 -1.1387058  -.08794144  -.37933505   -.1168786  -.7739556
2015 3 59   -.6031758  -1.010851    .0955884   -.8147325   .23368683  -.4167453
2016 3 56   -.1502212  -.3076502    .5238243   -.5586359    .5688137  -.6548856
2017 3 56   .11417311  -.3076502    .4626478   -.8629529    .8004407  -.4167453
2018 3 54    .6593583  .39555085   1.3191198   -.1474514    .9431381   .2976753
2019 3 72    1.040578   .0759141   1.3191198  -.09864163    1.050073  -.4167453
2020 3 76   1.1901965  -.8190691   1.1967665  -2.2260096   1.1619887  -1.131166
2021 3 51     1.69644  .20376867    1.441473  -1.2997588   1.2593246   .1786053
2022 3 53    2.194485          .           .   -.6446698   1.5510553          .
2006 4  4  -1.5086967   2.365859   -.9710394  -1.5298297  -.10335572   .6588714
2007 4 10   -.6170927  1.7043376   -.9052062   -.5740887 -.026472345   .8635697
2008 4  9    -.636063 .013781142   -.9052062   -.3030676 -.014590368  -.8763636
2009 4 16   -.7878254   -.721243  -1.0368725  -1.3405232    -.314785 -1.3881086
2010 4 18  -.10489471 -1.3092618  -1.1027057   -.8916544   -.7942577 -2.1045513
2011 4 29   -.7688551 -.13322355  -.51020706   -.0769716  -1.1290499  -.6716653
2012 4 30 -.010043235 -.13322355   -.2468744   -.1996927  -1.3562053  -.3646186
2013 4 17   -.5791521  -.6477403   -.2468744   -.3235072   -1.381367  -.9787125
2014 4 20   -.4843006  -1.603272  -.11520804  -.23814304  -1.0028912  -.4669675
2015 4 11  .008927062 -.05972085   -.3127075   -.4089962    -.198761   .2494753
2016 4 16   -.3515086 -.13322355  -.04937492   -.3687186    .4016283   .6588714
2017 4 18    .2365706  -.4272329     .279791    .3547538     .474318   .6588714
2018 4 19  -.44636005 -.20672554   1.5306215   1.0469381    .4505541   .6588714
2019 4 21    .4452439  .16078654   1.2014556   1.2479932    .5190502   .5565225
2020 4 21     2.11463  .16078654   1.5964546  .016114943    .7056671  1.1706164
2021 4 21    1.678313   .9693134   1.7939543   1.1185893   1.2829914  1.3753147
2022 4 17   1.8111053          .           .    2.470794    2.487614          .
2006 5 31   -1.749826   2.312512  -1.2208925  -1.4052644  -1.1497936   1.350849
2007 5 24   -1.010769  1.3875076  -1.0342592   -.4284081  -1.0751946   .7930793
2008 5 32   -.5142152 -1.0791715  -1.0342592   -.3331871  -1.0174965  -1.438001
2009 5 40  -.12159124 -1.3875076    -.909837   -1.720246   -.9836946 -1.7168866
2010 5 39   -.6527883  -.7708373    -.909837   -1.270801   -.9373319  -1.298559
2011 5 42   -.6758839  -.2569458   -.4121484  -.51357704   -.8440959  -.1830182
2012 5 66  -.21397334  .15416707  -.34993735 -.020703495   -.7328357  -.4619037
2013 5 90    -.768266   .3597245   -.1633041   -.2914089   -.5656134  -.8802309
2014 5 85   -.2255211  1.1819502  -.28772616    -.160464  -.29499415   .3747514
2015 5 50  -.23706888   .4625032  -.03888193  -.27178133   .13294171   .6536369
2016 5 58   .59437007  .15416707    .1477513    .1735733   .53458154   .3747514
2017 5 63  -.39873755  -.6680586    .1477513    .8214886    .7439796  -.6013461
2018 5 83   .05162521 -.15416707   1.7030284   1.0372539    .8908288  -.1830182
2019 5 80     .744491  .05139033    1.391973   1.5860833    1.048656   .9325224
2020 5 90     1.67986 -1.1819502   1.5786064  -.09748574   1.1513892   1.490292
2021 5 82    1.772242  -.5652799    1.391973   1.1223526   1.3180498   .7930793
2022 5 64    1.726051          .           .   1.7725903   1.7806556          .
2006 6 36  -1.8859606   1.116864  -1.1116489  -1.5856335  -1.1080685   .0109135
2007 6 47   -.8505021   .6644627  -1.0483522   -.8424931  -1.0370687 -.07639617
2008 6 57   -.3454004  -1.597542  -1.1116489    -.145503   -.9727876 -1.7352732
2009 6 57 -.067594424  -2.276143   -.9850555   -1.692368   -.9274691 -1.5606548
2010 6 49  .008170831   -1.29594   -.9850555  -1.0881214   -.9066905 -1.5606548
2011 6 63  .020798374  -.6173391   -.4786815    -.295721   -.8468712 -1.2987266
2012 6 63    .1723289 -.16493784   -.4153849   -.2681947   -.7395284 -.07639617
2013 6 70  -.50955844   .4382628  .027692353  -.10241153  -.59018755   .4474602
2014 6 76 -.067594424   .4382628  -.09890115 -.020433776   -.3464224   .4474602
2015 6 86   -.3580279   .8152636   .15428586 -.064035036   .02359635   .7966972
2016 6 70   .08393609  -.4665396  -.09890115      .12803    .3840948  .09822234
2017 6 73   -.7873644  .28746194   .15428586    .6092689    .6386773  .53476906
2018 6 67   -.7242267   .5136633   1.5468142   1.1957577    .8598221   1.233244
2019 6 95    .8415887  .28746194   1.5468142   1.2844496   1.0752081  1.0586255
2020 6 92     2.73572   .7398632   1.3569238   -.2151724    1.251735   .9713157
end

[Clyde Schechter]

If I understand you correctly, you want to identify effects that vary both over region and over time. I'm assuming that in your data the variable id corresponds to the regions you refer to. In your example data, id and year jointly identify unique observations: this is simple panel data. If that is true in the data as a whole, whether you use fixed or random effects, you are, in essence, asking for effect estimates calculated from and applying to single observations! Either way, this is conceptually not possible. If you attempt it with multi-level modeling most likely the estimation will not converge because the model is probably not identifiable. But if it does converge, the coefficients will be meaningless, and probably have enormously wide confidence intervals.

Now, if you wanted to consider the possibility of a more structured time variation, such as a linear time trend, you could include c.time in your fixed effects and have a random slope for it at the region level. If you really want to stretch your data to or slightly beyond its limits, you could consider linear and quadratic terms in c.time--but given the small number of time points, made even smaller by the inclusion of lags, I think there aren't enough time points to meaningfully support that level of detail in analysis.

Another possibility would be to group some of the regions together into larger spatial entities, assuming that can be done in some meaningful way so that these larger entities abstract some kind of relevant similarity among some of the regions, and include these as another random level in the model, over region. Then you could include a random slope for time at this higher spatial level.

[Doris Rivera]

Dear Clyde Schechter, thanks for your advices. Yes, here id is for regions. So, from your explanation I understand that my first point in #1 even though if it gives some estimations (see below for year==2012) those are not reliable.

Following your advice about the linear time trend, that is interesting, since it will give a hint of how the effect of x1 evolve as time pass. But two issues here, first, how could I make this prediction I need? and second, I am having an error as initial values are not feasible (see the last code line below). In fact, I believe I cannot use factor variables in the random part. If that is the case, then I must build the interaction, right?

About your second suggestion, indeed I can aggregate regions to country level. Let me see if I follow your idea. It would be to include a random coefficient for x1 at regional level, and another random coefficient for years at country level? Some doubts here: first, is year in the fixed part modelled as a linear trend (but without being interacted with x1)? second, is year in the random coefficient for country level, also a linear trend? and most important, from this estimation, how could I get the prediction for something like the effect of x1 varying with time?
After trying what I understood from your suggestion (maybe this is wrong), I get the same error:

. menbreg y1 L5.c.x1##c.year L5.c.x2 L5.c.x3 L5.c.x4 L5.c.x5 L5.c.x6 if year<=2022 || country: c.year || id: L5.c.x1 , cov(uns)  vce(robust) eform  nolog
initial values not feasible
r(1400);

Thanks a lot for your help with these issues!!!

. menbreg y1 L5.c.x1 L5.c.x2 L5.c.x3 L5.c.x4 L5.c.x5 L5.c.x6 i.year if year==2012 || id: L5.c.x1 , cov(unstructure) vce(robust) eform intpoints(10) nolog

Mixed-effects nbinomial regression              Number of obs     =        250
Overdispersion:            mean
Group variable:              id                 Number of groups  =        250

                                                Obs per group:
                                                              min =          1
                                                              avg =        1.0
                                                              max =          1

Integration method: mvaghermite                 Integration pts.  =         10

                                                Wald chi2(6)      =      72.82
Log pseudolikelihood = -1135.1253               Prob > chi2       =     0.0000
                                       (Std. Err. adjusted for 250 clusters in id)
----------------------------------------------------------------------------------
                 |               Robust
              y1 |     exp(b)   Std. Err.      z    P>|z|     [95% Conf. Interval]
-----------------+----------------------------------------------------------------
              x1 |
             L5. |   .8721776   .1350264    -0.88   0.377     .6439111    1.181365
                 |
              x2 |
             L5. |   1.246997    .195958     1.40   0.160     .9164401    1.696784
                 |
              x3 |
             L5. |   .5777686   .1451217    -2.18   0.029     .3531444    .9452693
                 |
              x4 |
             L5. |   .7227296   .0709664    -3.31   0.001     .5962035    .8761071
                 |
              x5 |
             L5. |   .5799184   .0578084    -5.47   0.000     .4769974    .7050465
                 |
              x6 |
             L5. |   1.362799   .1664035     2.54   0.011     1.072745    1.731279
                 |
           _cons |   3.987052   1.490606     3.70   0.000     1.916109    8.296284
-----------------+----------------------------------------------------------------
        /lnalpha |  -.2055929    .135708                     -.4715757    .0603899
-----------------+----------------------------------------------------------------
id               |
       var(L5.x1)|   .3742028   1.11e-07                      .3742026     .374203
       var(_cons)|   .1690614   1.10e-08                      .1690614    .1690614
-----------------+----------------------------------------------------------------
id               |
 cov(L5.x1,_cons)|  -.2515218   3.70e-08 -6.8e+06   0.000    -.2515219   -.2515218
----------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation.

. menbreg y1 L5.c.x1##c.year L5.c.x2 L5.c.x3 L5.c.x4 L5.c.x5 L5.c.x6 if year<=2022 || id: L5.c.x1##c.year , cov(uns)  vce(robust) eform intpoints(15) nolog
initial values not feasible
r(1400);

[Clyde Schechter]

...I understand that my first point in #1 even though if it gives some estimations (see below for year==2012) those are not reliable.

Sorry, I don't think I was clear. The analysis you show in #1 does contains a random slope for x1 at the region level, but there is no region#time interaction, nor any interactions between time and any of the x variables. So it is not problematic. You can use the results of that analysis. What would be problematic is to have a random slope for a variable at the region level and also have a time-varying slope for the same variable.

menbreg y1 L5.c.x1##c.year L5.c.x2 L5.c.x3 L5.c.x4 L5.c.x5 L5.c.x6 if year<=2022 || country: c.year || id: L5.c.x1 , cov(uns) vce(robust) eform nolog

isn't quite what I had in mind. What I was thinking of (and I'm using more compact notation than you did, so distinguish that from the substantive change, which is in bold face):

menbreg y1 c.L5.x1##c. year L5.c.(x2 x3 x4 x5 x6) if year <= 2022 || country: c.L5.x1 || id: , cov(uns) vce(robust) eform nolog

In this model you avoid the problem I raised in #2 because the variations in x1 effect are at the country and year level, and country and year do not jointly identify unique observations.

Now, I don't know of any simple solutions to the initial values not feasible problem. The only thing I know of that can help with that is to simplify the model. The first thing I would do to simplify this model is remove -cov(uns)-. While it's nice not to impose restrictions on the covariance structure that you cannot provide a justification for, it also adds a large number of degrees of freedom to the calculation and can make it harder to find feasible initial values. If that alone does not resolve the problem, you may have to start stripping out some variables from your model.

[Doris Rivera]

Dear Clyde Schechter , I understand that this new suggested approach change the focus from regions to country (which is not bad for the study). However, some things are not clear to me. It is very likely I am wrong, but, if using "c.L5.x1##c. year" in the fixed part of the model, would not that imply that "|| country: c.L5.x1" is the effect of x1 for the base category (which given that year here is linear, it would be when year==0)?

I see that from the fixed part, the main effect from x1 varies with years. But from the random part, it only vary for country, not for country in different periods. So, country variation for x1 is fixed, while time variation for x1 would be the same for all countries, right? I mean, we are allowing that country 1 is in general term different from country 2 (in term of x1), but cannot say this difference is changing with time. Since time difference for x1 is the same for all countries.

This conversation is really helpful since I am realizing that maybe what I thought is not possible. This idea comes from a paper about firm's R&D from Ann Marie knott, in which she uses a 10 years period with a random coefficient for R&D in a time rolling window (one regression from 1990-2000, another regression from 1991-2001,...), and assign the fixed + random effect prediction to the last year of each regression. That way she manages to have a country variation in each year (for the last 20 years). However, she have 30 years, but in my case with this very short panel, that rolling window analysis might not be possible, since it will only leaves 3 or 4 years for the second step analysis in which I plan to use the desired prediction as indep variable. This is why I thought about using the regressions for each year as I show in #3.

On top of that, it seems that my data is not helping since I got another error when trying your solution. I knew about the difficulty of multilevel model though..

. menbreg y1 c.L5.x1##c.year L5.c.(x2 x3 x4 x5 x6) if year <= 2022 || country: c.L5.x1 || id: , cov(uns) vce(robust) eform intpoints(15)

Fitting fixed-effects model:

Iteration 0:   log likelihood = -14794.268  
Iteration 1:   log likelihood = -14147.616  
Iteration 2:   log likelihood = -14126.987  
Iteration 3:   log likelihood = -14126.833  
Iteration 4:   log likelihood = -14126.833  

Refining starting values:

Grid node 0:   log likelihood = -13024.738

Fitting full model:

Iteration 0:   log pseudolikelihood = -13024.738  (not concave)
Iteration 1:   log pseudolikelihood = -12990.091  (not concave)
Iteration 2:   log pseudolikelihood = -12922.878  (not concave)
Iteration 3:   log pseudolikelihood = -12895.784  (not concave)
Iteration 4:   log pseudolikelihood = -11052.386  
cannot compute an improvement -- discontinuous region encountered
r(430);

[Clyde Schechter]

...if using "c.L5.x1##c. year" in the fixed part of the model, would not that imply that "|| country: c.L5.x1" is the effect of x1 for the base category (which given that year here is linear, it would be when year==0)?

No. c.year is treated as a continuous variable, so there is no base category.

But you are correct that the random slope estimates corresponding to -||country: c.L5.x1- are time invariant. Given that country and year do not uniquely identify observations, you could, indeed go to - || country: c.L5.x1##c.year-, to get time-varying country-varying random slopes. (But you cannot do the analogous thing at the region level, because region and year do jointly identify unique observations.)

Whether you can get the estimation to converge is another matter altogether. You just have to try it and see what happens. Unlike simple one-level regressions, with a multi-level model there is no theorem that says that maximum likelihood estimates can be calculated by the usual algorithms. Usually they can, but the likeilhood functions of these models can be very irregular, and sometimes it just can't be done.

[Doris Rivera]

Dear Clyde Schechter, thanks a lot for your help!