How to specify a SURE (or bivariate) Poisson model using -gsem- ?

[Marco Ercolani]

Dear Statalisters,

I need to estimate non-linear SURE models (Poisson, negative binomial, ordered probit, zero-inflated) in StataNow 18.5 but am stuck on how to do so. Stata’s -gsem- command seems the best option because it allows one to specify various non-linear SUREs and it reports the joint log-likelihood value. Within -gsem-, I seem incapable of specifying the cross-equation covariance when the disturbances are non-Gaussian. (I’m moderately well-versed in estimating linear SURE models and aware of the semantic difference between SURE and multivariate models.)

For the sake of argument, assume I wish to estimate a Poisson SURE model. This allows me to demonstrate my problem using the user-contributed -bivpoisson- command and its data. The SURE example in -bivpoisson- works well (btw, it does not report the joint log-likelihood):

use "https://github.com/zhangyl334/bivpoisson/raw/main/Health_Data.dta", clear
bivpoisson (ofp = privins black numchron) (ofnp = privins black numchron age)
...
Iteration 9:  f(p) = -832.69538 
                                                           Number of obs = 207

------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
ofp          |
     privins |   .3997619   .1830324     2.18   0.029     .0410251    .7584988
       black |  -.1335777   .1905021    -0.70   0.483    -.5069549    .2397995
    numchron |   .2380122    .053071     4.48   0.000      .133995    .3420294
       _cons |   .6682984   .1939622     3.45   0.001     .2881395    1.048457
-------------+----------------------------------------------------------------
ofnp         |
     privins |   1.305625   .4458126     2.93   0.003     .4318481    2.179401
       black |  -2.151162    .919045    -2.34   0.019    -3.952457   -.3498668
    numchron |   .2358258   .1392374     1.69   0.090    -.0370744    .5087261
         age |  -.0809188   .3125795    -0.26   0.796    -.6935634    .5317259
       _cons |  -2.271814   2.292566    -0.99   0.322    -6.765161    2.221533
-------------+----------------------------------------------------------------
sigmasq1     |
       _cons |   .8514477    .130599     6.52   0.000     .5954784    1.107417
-------------+----------------------------------------------------------------
sigmasq2     |
       _cons |    3.47855   .6043017     5.76   0.000      2.29414     4.66296
-------------+----------------------------------------------------------------
sigma12      |
       _cons |   .4178387   .2111369     1.98   0.048      .004018    .8316593
------------------------------------------------------------------------------

If I ignore the cross-equation covariance, -gsem- can easily estimate the same two-equation system:

gsem (ofp <- privins black numchron, family(poisson) link(log)) ///
(ofnp <- privins black numchron age, family(poisson) link(log))
...
Iteration 5:  Log likelihood = -1269.1883  

Generalized structural equation model                      Number of obs = 207

Response: ofp    
Family:   Poisson
Link:     Log    

Response: ofnp   
Family:   Poisson
Link:     Log    

Log likelihood = -1269.1883

------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
ofp          |
     privins |   .3544521   .0806143     4.40   0.000      .196451    .5124532
       black |   .2887277   .0998223     2.89   0.004     .0930796    .4843758
    numchron |   .1585311   .0197342     8.03   0.000     .1198528    .1972094
       _cons |   1.171029    .083035    14.10   0.000     1.008283    1.333774
-------------+----------------------------------------------------------------
ofnp         |
     privins |   1.283982   .2688524     4.78   0.000     .7570411    1.810923
       black |  -1.665092   .5855751    -2.84   0.004    -2.812798   -.5173859
    numchron |   .1740241    .041588     4.18   0.000     .0925131    .2555352
         age |  -.0513212   .1111205    -0.46   0.644    -.2691134    .1664709
       _cons |  -.7977588    .843642    -0.95   0.344    -2.451267    .8557492
------------------------------------------------------------------------------

However, I’m stuck on how to use -gsem- to estimate the same system with cross-equation covariance. The following code does not work because of the constrained constants but it might include something useful:

gsem (ofp <- _cons@0 privins black numchron Const1, family(poisson) link(log)) ///
(ofnp <- _cons@0 privins black numchron age Const2, family(poisson) link(log)), ///
vsquish covstruct(_lexogenous, diagonal) latent(Const1 Const2) cov(Const1*Const2) nocapslatent
...
Iteration 9:  Log likelihood = -835.80108  

Generalized structural equation model                      Number of obs = 207

Response: ofp    
Family:   Poisson
Link:     Log    

Response: ofnp   
Family:   Poisson
Link:     Log    

Log likelihood = -835.80108

 ( 1)  [ofp]_cons = 0
 ( 2)  [ofp]Const1 = 1
 ( 3)  [ofnp]_cons = 0
 ( 4)  [ofnp]Const2 = 1
------------------------------------------------------------------------------------
                   | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------------+----------------------------------------------------------------
ofp                |
           privins |   .9009955   .1267399     7.11   0.000     .6525899    1.149401
             black |   .3153675   .2670429     1.18   0.238     -.208027     .838762
          numchron |   .3150931   .0543793     5.79   0.000     .2085117    .4216746
            Const1 |          1  (constrained)
             _cons |          0  (omitted)
-------------------+----------------------------------------------------------------
ofnp               |
           privins |   1.619999   .5899405     2.75   0.006     .4637373    2.776262
             black |  -2.312566   1.181007    -1.96   0.050    -4.627297    .0021655
          numchron |   .2799475   .1482755     1.89   0.059    -.0106672    .5705621
               age |  -.4640784   .0925598    -5.01   0.000    -.6454924   -.2826644
            Const2 |          1  (constrained)
             _cons |          0  (omitted)
-------------------+----------------------------------------------------------------
        var(Const1)|   .9947875   .1421176                      .7518405    1.316239
        var(Const2)|   4.505426   1.070492                      2.828068    7.177644
-------------------+----------------------------------------------------------------
 cov(Const1,Const2)|   .5386765    .260694     2.07   0.039     .0277256    1.049627
------------------------------------------------------------------------------------

How can I use the -gsem- command to get close to replicating the -bivpoisson- estimate above? I feel I'm nearly there, but just can't work out the exact -gsem- options.

[Marco Ercolani]

It’s been a month, and there have been no suggestion, so it’s either too hard or too easy for most. The best I could come up with was the following, which gives slightly different regression results compared to -bivpoisson-.

gsem ///
(ofp  <- privins black numchron     Const1       , family(poisson) link(log)) ///
(ofnp <- privins black numchron age Const2       , family(poisson) link(log)) ///
, vsquish covstruct(_lexogenous, diagonal) latent(Const1 Const2) cov(Const1*Const2) nocapslatent
…
Iteration 9:  Log likelihood = -830.94026 
 
Generalized structural equation model                      Number of obs = 207
 
Response: ofp   
Family:   Poisson
Link:     Log   
 
Response: ofnp  
Family:   Poisson
Link:     Log   
 
Log likelihood = -830.94026
 
 ( 1)  [ofp]Const1 = 1
 ( 2)  [ofnp]Const2 = 1
------------------------------------------------------------------------------------
                   | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------------+----------------------------------------------------------------
ofp                |
           privins |   .4475027   .1923388     2.33   0.020     .0705256    .8244799
             black |   .0235707   .2688634     0.09   0.930    -.5033919    .5505333
          numchron |   .2351314   .0572258     4.11   0.000     .1229709     .347292
            Const1 |          1  (constrained)
             _cons |   .6206109   .2033017     3.05   0.002     .2221468    1.019075
-------------------+----------------------------------------------------------------
ofnp               |
           privins |   1.582423   .6057297     2.61   0.009     .3952148    2.769632
             black |  -2.329756   1.184758    -1.97   0.049    -4.651839   -.0076717
          numchron |   .2571534   .1468214     1.75   0.080    -.0306113    .5449181
               age |  -.0554717   .3289951    -0.17   0.866    -.7002903    .5893469
            Const2 |          1  (constrained)
             _cons |  -2.937929   2.514607    -1.17   0.243    -7.866468    1.990611
-------------------+----------------------------------------------------------------
        var(Const1)|   .8503109   .1313055                      .6282533    1.150855
        var(Const2)|   4.512871   1.068897                      2.836883    7.179008
-------------------+----------------------------------------------------------------
 cov(Const1,Const2)|   .4843369    .236391     2.05   0.040     .0210191    .9476548
------------------------------------------------------------------------------------

The “1” parameters ‘seem’ irrelevant and changing them only affects the variances and covariance. Thus, -gsem- either gives better results than -bivpoisson-, or incorrect results.

[Erik Ruzek]

I am not sure this can be done in gsem. Check out this article by Hu & Hardin (2016) in Stata Journal for their bivcnto program. It estimates a bivariate poisson (or nbreg) model that utilizes different a copula function approach for the joint probabilities. Note that certain models report a lambda coefficient, the sign of which indicates the direction of the correlation (positive or negative). Others report the correlation directly. All models estimated by this program report the final loglikelihood of the model.