Using suest after margins with a Tobit model

[Katy Tabero]

Hello,

I would like to know if there is a way to use suest after using the margins command following a Tobit regression.

I have been using regression coefficients in a test calculation for my project; this involves running two regressions (the same data and independent variables, but differing dependent variables) and referencing coefficients in each as part of a larger calculation. Here is a simplified example of the code below:

// Regress

// Create and Store coefficients

regress team_correct endoimposed endonotimposed exoimposed exonotimposed, noconstant
est store reg_correct

regress team_tetris endoimposed endonotimposed exoimposed exonotimposed, noconstant
est store reg_tetris

suest reg_correct reg_tetris, vce(cluster session)

// Run Tests 

test ([reg_correct_mean]endoimposed*3*0.2)+([reg_tetris_mean]endoimposed/60*0.16)=([reg_correct_mean]exoimposed*3*0.2)+([reg_tetris_mean]exoimposed/60*0.16)

However, I have since tried running the same test with a Tobit model (see code below), and receive the error: ape_correct was estimated with a nonstandard vce (delta).

// Create and Store coefficients

tobit team_correct endoimposed endonotimposed exoimposed exonotimposed, ll(0) noconstant
margins, predict(ystar(0,.)) dydx(*) post
est store ape_correct

tobit team_tetris endoimposed endonotimposed exoimposed exonotimposed, ll(0) noconstant
margins, predict(ystar(0,.)) dydx(*) post
est store ape_tetris

suest ape_correct ape_tetris


// Run Tests 

test (0.2*[ape_correct_team_correct]endoimposed*3)+([ape_tetris_team_tetris]endoimposed/60*0.16) = (0.2*[ape_correct_team_correct]exoimposed*3)+([ape_tetris_team_tetris]exoimposed/60*0.16)

Is there a way that I can access the APEs in the same way as the regression coefficients? I have tried several variations of the same code and commands, but it either causes the same error or saves the APEs as a scalar (in which case the test does not produce results).

Please let me know if you have any ideas, or if I am in the wrong area entirely.

Best wishes,

Katy

(Stata 17)

[Andrew Musau]

You cannot combine margins results using suest as margins outputs a nonstandard VCE. This is a tough one as the tobit includes an ancillary parameter ("var(e.depvar)"), and a stacked model will thus not generate the correct results. It would be possible to replace the coefficient matrix and the variance matrix from glm and combine using suest, but you would be responsible for any flaws inherent in this method.

[Katy Tabero]

Thank you for your quick response. I have taken a look at the documentation for the glm command and there doesn't seem to be an option for censored data - do you know if there is a way around this if I was to try to try the method you suggest?

[Andrew Musau]

You do not do it using glm, but you replace its results with those from margins. The caveat is that you are responsible for any issues with this method. Below, I use erepost from SSC.

ssc install erepost, replace

webuse mroz87, clear
set seed 03302022
gen outcome2= runiformint(1, 4920)
tobit whrs75 nwinc wedyrs wexper c.wexper#c.wexper wifeage kl6 k618, ll(0)
margins, predict(ystar(0,.)) dydx(*) post
mat b= e(b)
mat V= e(V)
qui glm whrs75 nwinc wedyrs wexper wifeage kl6 k618, nocons
erepost b= b, rename
erepost V= V, rename
est store ape_correct
tobit outcome2 nwinc wedyrs wexper c.wexper#c.wexper wifeage kl6 k618, ll(0)
margins, predict(ystar(0,.)) dydx(*) post
mat b= e(b)
mat V= e(V)
qui glm outcome2 nwinc wedyrs wexper wifeage kl6 k618, nocons
erepost b= b, rename
erepost V= V, rename
est store ape_tetris
suest ape_correct ape_tetris

Res.:

. tobit whrs75 nwinc wedyrs wexper c.wexper#c.wexper wifeage kl6 k618, ll(0)

Refining starting values:

Grid node 0:   log likelihood = -3961.1577

Fitting full model:

Iteration 0:   log likelihood = -3961.1577  
Iteration 1:   log likelihood = -3836.8928  
Iteration 2:   log likelihood = -3819.2637  
Iteration 3:   log likelihood = -3819.0948  
Iteration 4:   log likelihood = -3819.0946  

Tobit regression                                Number of obs     =        753
                                                   Uncensored     =        428
Limits: lower = 0                                  Left-censored  =        325
        upper = +inf                               Right-censored =          0

                                                LR chi2(7)        =     271.59
                                                Prob > chi2       =     0.0000
Log likelihood = -3819.0946                     Pseudo R2         =     0.0343

-----------------------------------------------------------------------------------
           whrs75 |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
------------------+----------------------------------------------------------------
            nwinc |  -8.814227   4.459089    -1.98   0.048    -17.56808   -.0603708
           wedyrs |   80.64541   21.58318     3.74   0.000     38.27441    123.0164
           wexper |    131.564   17.27935     7.61   0.000     97.64211     165.486
                  |
c.wexper#c.wexper |  -1.864153   .5376606    -3.47   0.001    -2.919661   -.8086455
                  |
          wifeage |  -54.40491   7.418483    -7.33   0.000     -68.9685   -39.84133
              kl6 |  -894.0202   111.8777    -7.99   0.000    -1113.653   -674.3875
             k618 |  -16.21805    38.6413    -0.42   0.675    -92.07668    59.64057
            _cons |   965.3068   446.4351     2.16   0.031     88.88827    1841.725
------------------+----------------------------------------------------------------
     var(e.whrs75)|    1258927   93304.48                       1088458     1456093
-----------------------------------------------------------------------------------

. 
. margins, predict(ystar(0,.)) dydx(*) post

Average marginal effects                        Number of obs     =        753
Model VCE    : OIM

Expression   : E(whrs75*|whrs75>0), predict(ystar(0,.))
dy/dx w.r.t. : nwinc wedyrs wexper wifeage kl6 k618

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       nwinc |  -5.188619   2.621409    -1.98   0.048    -10.32649   -.0507516
      wedyrs |   47.47306    12.6214     3.76   0.000     22.73558    72.21054
      wexper |    48.7931   3.587269    13.60   0.000     41.76218    55.82402
     wifeage |  -32.02622    4.29211    -7.46   0.000     -40.4386   -23.61384
         kl6 |  -526.2776   64.70619    -8.13   0.000    -653.0994   -399.4558
        k618 |  -9.546987   22.75224    -0.42   0.675    -54.14056    35.04659
------------------------------------------------------------------------------


. 
. tobit outcome2 nwinc wedyrs wexper c.wexper#c.wexper wifeage kl6 k618, ll(0)

Refining starting values:

Grid node 0:   log likelihood = -6550.7516

Fitting full model:

Iteration 0:   log likelihood = -6550.7516  
Iteration 1:   log likelihood = -6550.7516  

Tobit regression                                Number of obs     =        753
                                                   Uncensored     =        753
Limits: lower = 0                                  Left-censored  =          0
        upper = +inf                               Right-censored =          0

                                                LR chi2(7)        =       5.70
                                                Prob > chi2       =     0.5753
Log likelihood = -6550.7516                     Pseudo R2         =     0.0004

-----------------------------------------------------------------------------------
         outcome2 |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
------------------+----------------------------------------------------------------
            nwinc |    -2.8083   4.923527    -0.57   0.569    -12.47392    6.857317
           wedyrs |   23.37747   25.07164     0.93   0.351    -25.84191    72.59684
           wexper |   25.51334   19.28185     1.32   0.186    -12.33979    63.36648
                  |
c.wexper#c.wexper |  -.8997776   .6281175    -1.43   0.152    -2.132866    .3333107
                  |
          wifeage |     7.6236   8.445607     0.90   0.367    -8.956385    24.20359
              kl6 |   159.3098   113.8887     1.40   0.162    -64.27071    382.8903
             k618 |   47.07904   44.85406     1.05   0.294    -40.97618    135.1343
            _cons |   1768.561   524.0627     3.37   0.001     739.7478    2797.374
------------------+----------------------------------------------------------------
   var(e.outcome2)|    2107889   108633.9                       1905058     2332315
-----------------------------------------------------------------------------------

. 
. margins, predict(ystar(0,.)) dydx(*) post

Average marginal effects                        Number of obs     =        753
Model VCE    : OIM

Expression   : E(outcome2*|outcome2>0), predict(ystar(0,.))
dy/dx w.r.t. : nwinc wedyrs wexper wifeage kl6 k618

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       nwinc |   -2.69342   4.721821    -0.57   0.568    -11.94802     6.56118
      wedyrs |   22.42115   24.04194     0.93   0.351    -24.70018    69.54248
      wexper |   6.134151   8.223431     0.75   0.456    -9.983477    22.25178
     wifeage |   7.311737   8.098861     0.90   0.367    -8.561739    23.18521
         kl6 |   152.7928    109.191     1.40   0.162    -61.21766    366.8033
        k618 |   45.15315   43.01024     1.05   0.294    -39.14537    129.4517
------------------------------------------------------------------------------



. 
. suest ape_correct ape_tetris

Simultaneous results for ape_correct, ape_tetris

                                                Number of obs     =        753

------------------------------------------------------------------------------
             |               Robust
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
ape_correct  |
       nwinc |  -5.188619   6.718937    -0.77   0.440    -18.35749    7.980255
      wedyrs |   47.47306   138.7741     0.34   0.732    -224.5192    319.4653
      wexper |    48.7931   6.426354     7.59   0.000     36.19767    61.38852
     wifeage |  -32.02622   58.63383    -0.55   0.585    -146.9464    82.89398
         kl6 |  -526.2776   364.1324    -1.45   0.148    -1239.964    187.4088
        k618 |  -9.546987   178.3852    -0.05   0.957    -359.1756    340.0817
-------------+----------------------------------------------------------------
ape_tetris   |
       nwinc |   -2.69342   9.643416    -0.28   0.780    -21.59417    16.20733
      wedyrs |   22.42115   225.0471     0.10   0.921     -418.663    463.5053
      wexper |   6.134151   15.78436     0.39   0.698    -24.80262    37.07092
     wifeage |   7.311737   96.06982     0.08   0.939    -180.9816    195.6051
         kl6 |   152.7928   550.7591     0.28   0.781    -926.6751    1232.261
        k618 |   45.15315   284.6732     0.16   0.874    -512.7961    603.1024
------------------------------------------------------------------------------

.

[Katy Tabero]

I understand what you mean now. Thank you for demonstrating the code, it is very helpful. I am certainly responsible for any issues using this method as I understand it is not using Stata as intended to some degree.

Given this, I have tested it on my data and it does produce sensible-looking results and allow me to use the APE's in my tests without the previous error occurring. However, unfortunately, the standard errors are not the same as those produced by the margins code, which makes me think that the test results are unreliable, even though I know the APE values themselves are correct. As a quick test, I ran a simple 'test _b[endoimposed] = 10' after the margins result and then 'test [ape_correct]endoimposed = 10' after the suest command, and I do get different results. I would think these should be the same if the method was working as intended.

[Andrew Musau]

Here is a stacked model where you are constraining the residual variance to be the same across both groups. This restriction means that the coefficients will not exactly match with the separate regressions. You need to generate the coefficients manually (not as interactions) to correctly run margins on the stacked results.

webuse mroz87, clear
set seed 03302022
gen outcome2= runiformint(1, 4920)
expand 2, g(new)
gen outcome= cond(new, outcome2, whrs75)
gen cons=1
tobit outcome i.new#(c.nwinc c.wedyrs c.wexper c.wexper#c.wexper c.wifeage c.kl6 c.k618), nocons ll(0)

Res.:

. tobit outcome i.new#(c.nwinc c.wedyrs c.wexper c.wexper#c.wexper c.wifeage c.kl6 c.k618), nocons ll(0)

Refining starting values:

Grid node 0:   log likelihood = -10492.499

Fitting full model:

Iteration 0:   log likelihood = -10492.499  
Iteration 1:   log likelihood = -10397.712  
Iteration 2:   log likelihood = -10392.972  
Iteration 3:   log likelihood = -10392.939  
Iteration 4:   log likelihood = -10392.939  

Tobit regression                                Number of obs     =      1,506
                                                   Uncensored     =      1,181
Limits: lower = 0                                  Left-censored  =        325
        upper = +inf                               Right-censored =          0

                                                F(  14,   1492)   =     204.89
                                                Prob > F          =     0.0000
Log likelihood = -10392.939                     Pseudo R2         =     0.0675

 ( 1)  [outcome]_cons = 0
---------------------------------------------------------------------------------------
              outcome |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
----------------------+----------------------------------------------------------------
          new#c.nwinc |
                   0  |  -10.81984   5.414174    -2.00   0.046    -21.44004    -.199635
                   1  |   -4.47991   4.597772    -0.97   0.330    -13.49869    4.538873
                      |
         new#c.wedyrs |
                   0  |   117.7736   20.81981     5.66   0.000     76.93436    158.6128
                   1  |   75.32679   18.57433     4.06   0.000     38.89221    111.7614
                      |
         new#c.wexper |
                   0  |   153.5471   20.41576     7.52   0.000     113.5004    193.5937
                   1  |   39.26225   17.68933     2.22   0.027     4.563663    73.96084
                      |
new#c.wexper#c.wexper |
                   0  |  -2.381474   .6322801    -3.77   0.000    -3.621727   -1.141222
                   1  |  -1.418335   .5716371    -2.48   0.013    -2.539633   -.2970377
                      |
        new#c.wifeage |
                   0  |   -48.7055   5.757885    -8.46   0.000    -59.99991   -37.41109
                   1  |     30.512   4.723714     6.46   0.000     21.24617    39.77782
                      |
            new#c.kl6 |
                   0  |  -926.4693   129.9171    -7.13   0.000    -1181.309   -671.6298
                   1  |   286.6867   100.8546     2.84   0.005     88.85488    484.5186
                      |
           new#c.k618 |
                   0  |   21.60309   42.36491     0.51   0.610    -61.49802    104.7042
                   1  |   112.2164   38.00266     2.95   0.003     37.67209    186.7607
                      |
                _cons |          0  (omitted)
----------------------+----------------------------------------------------------------
        var(e.outcome)|    1856983   78158.97                       1709828     2016803
---------------------------------------------------------------------------------------

[FernandoRios]

Hi Katy
Some additional food for thought here.
They way Suest works (based on my understanding), is that standard errors are re-estimated, making use score functions associated with the estimated models.
When you estimate marginal effects, however, there is no score function to be created. Just a set of variance and covarainces. I think the Trick that Andrew suggests allows you to bypass this problem with suest, but because there is no proper score, the re estimation of the variance-covariance matrix is incorrect, which gives you the discrepancies your report, and that Andrew shows in his example.

I don't think trying to apply suest to margins is the correct approach. Instead, you should be thinking about applying suest (or other simultaneous regression) to the tobit model (or try my own take with "ssc install mbitobit").

Once the simultaneous tobit is estimated, you can get the marginal effects that account for correlation.

If you use suest, you will have to do additional work using the function "expression", since you will need to construct the parameter of interest (the expression equivalent to ystar(0,.)). If you use mbitobit, I think I added many options in predict that may suit your needs.
Hope this helps
F
Example, following Andrew's setup:

webuse mroz87, clear
set seed 03302022
gen outcome2= runiformint(1, 4920)
mbitobit (whrs75 =  nwinc wedyrs wexper c.wexper#c.wexper wifeage kl6 k618) ///
(outcome2 =  nwinc wedyrs wexper c.wexper#c.wexper wifeage kl6 k618)
margins, dydx(*) predict(ys1) predict(ys2)

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
nwinc        |
    _predict |
          1  |  -5.186877   2.621095    -1.98   0.048    -10.32413   -.0496256
          2  |  -2.693418   4.721821    -0.57   0.568    -11.94802    6.561181
-------------+----------------------------------------------------------------
wedyrs       |
    _predict |
          1  |   47.42546   12.62255     3.76   0.000     22.68571     72.1652
          2  |    22.4211   24.04194     0.93   0.351    -24.70023    69.54242
-------------+----------------------------------------------------------------
wexper       |
    _predict |
          1  |   48.80014   3.587215    13.60   0.000     41.76933    55.83095
          2  |   6.134165    8.22343     0.75   0.456    -9.983462    22.25179
-------------+----------------------------------------------------------------
wifeage      |
    _predict |
          1  |  -32.01075   4.292556    -7.46   0.000      -40.424   -23.59749
          2  |   7.311757    8.09886     0.90   0.367    -8.561718    23.18523
-------------+----------------------------------------------------------------
kl6          |
    _predict |
          1  |  -525.8472   64.72207    -8.12   0.000    -652.7001   -398.9943
          2  |   152.7932    109.191     1.40   0.162    -61.21721    366.8037
-------------+----------------------------------------------------------------
k618         |
    _predict |
          1  |  -9.519593   22.75379    -0.42   0.676    -54.11621    35.07702
          2  |   45.15318   43.01024     1.05   0.294    -39.14533    129.4517
------------------------------------------------------------------------------

[Andrew Musau]

Thanks FernandoRios for the suggestion to use mbitobit from SSC.

[Katy Tabero]

Thank you both, these suggestions are definitely helping. Hopefully, 2 quick last questions. First, when using mbitobit, the results vary slightly from those when using tobit (this can be seen in slight variations between your example results in posts #4 and #7, and I also experienced this when using my data) - do you know why this is? Second, is there a no constant option when using mbitobit (or a way that I could implement this)? I cannot find anything about it in the documentation, but it would be useful in my case.

[FernandoRios]

Hi Katy
The reason they vary is because they are not estimated independent from each other. Instead, they are estimated simultaneously. This affects the objective function slightly.
It is kind of the same reason why a two step heckman and a ML heckman have different results, even in terms of the probit itself (selection model).

And, No, I didn't add a no-constant option. Since I didn't need that when I wrote the command
F

[Katy Tabero]

Hi Fernando,

Thank you for the explanation, that makes sense. I will keep looking for methods that will allow no constant as an option as I need a coefficient for each variable to use in my 'test' calculations, but this has been very useful even so.

[Andrew Musau]

Note that I forgot to include the constant that I generated in #6, otherwise at least one set of coefficients should match exactly.

webuse mroz87, clear
set seed 03302022
gen outcome2= runiformint(1, 4920)
tobit whrs75 c.nwinc c.wedyrs c.wexper c.wexper#c.wexper c.wifeage c.kl6 c.k618, ll(0)
tobit outcome2 c.nwinc c.wedyrs c.wexper c.wexper#c.wexper c.wifeage c.kl6 c.k618, ll(0)
expand 2, g(new)
gen outcome= cond(new, outcome2, whrs75)
gen cons=1
tobit outcome i.new#(c.nwinc c.wedyrs c.wexper c.wexper#c.wexper c.wifeage c.kl6 c.k618 c.cons), nocons ll(0) vce(robust)

Res.:

. tobit whrs75 c.nwinc c.wedyrs c.wexper c.wexper#c.wexper c.wifeage c.kl6 c.k618, ll(0)

Refining starting values:

Grid node 0:   log likelihood = -3961.1577

Fitting full model:

Iteration 0:   log likelihood = -3961.1577  
Iteration 1:   log likelihood = -3836.8928  
Iteration 2:   log likelihood = -3819.2637  
Iteration 3:   log likelihood = -3819.0948  
Iteration 4:   log likelihood = -3819.0946  

Tobit regression                                Number of obs     =        753
                                                   Uncensored     =        428
Limits: lower = 0                                  Left-censored  =        325
        upper = +inf                               Right-censored =          0

                                                LR chi2(7)        =     271.59
                                                Prob > chi2       =     0.0000
Log likelihood = -3819.0946                     Pseudo R2         =     0.0343

-----------------------------------------------------------------------------------
           whrs75 |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
------------------+----------------------------------------------------------------
            nwinc |  -8.814227   4.459089    -1.98   0.048    -17.56808   -.0603708
           wedyrs |   80.64541   21.58318     3.74   0.000     38.27441    123.0164
           wexper |    131.564   17.27935     7.61   0.000     97.64211     165.486
                  |
c.wexper#c.wexper |  -1.864153   .5376606    -3.47   0.001    -2.919661   -.8086455
                  |
          wifeage |  -54.40491   7.418483    -7.33   0.000     -68.9685   -39.84133
              kl6 |  -894.0202   111.8777    -7.99   0.000    -1113.653   -674.3875
             k618 |  -16.21805    38.6413    -0.42   0.675    -92.07668    59.64057
            _cons |   965.3068   446.4351     2.16   0.031     88.88827    1841.725
------------------+----------------------------------------------------------------
     var(e.whrs75)|    1258927   93304.48                       1088458     1456093
-----------------------------------------------------------------------------------

. 
. tobit outcome2 c.nwinc c.wedyrs c.wexper c.wexper#c.wexper c.wifeage c.kl6 c.k618, ll(0)

Refining starting values:

Grid node 0:   log likelihood = -6550.7516

Fitting full model:

Iteration 0:   log likelihood = -6550.7516  
Iteration 1:   log likelihood = -6550.7516  

Tobit regression                                Number of obs     =        753
                                                   Uncensored     =        753
Limits: lower = 0                                  Left-censored  =          0
        upper = +inf                               Right-censored =          0

                                                LR chi2(7)        =       5.70
                                                Prob > chi2       =     0.5753
Log likelihood = -6550.7516                     Pseudo R2         =     0.0004

-----------------------------------------------------------------------------------
         outcome2 |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
------------------+----------------------------------------------------------------
            nwinc |    -2.8083   4.923527    -0.57   0.569    -12.47392    6.857317
           wedyrs |   23.37747   25.07164     0.93   0.351    -25.84191    72.59684
           wexper |   25.51334   19.28185     1.32   0.186    -12.33979    63.36648
                  |
c.wexper#c.wexper |  -.8997776   .6281175    -1.43   0.152    -2.132866    .3333107
                  |
          wifeage |     7.6236   8.445607     0.90   0.367    -8.956385    24.20359
              kl6 |   159.3098   113.8887     1.40   0.162    -64.27071    382.8903
             k618 |   47.07904   44.85406     1.05   0.294    -40.97618    135.1343
            _cons |   1768.561   524.0627     3.37   0.001     739.7478    2797.374
------------------+----------------------------------------------------------------
   var(e.outcome2)|    2107889   108633.9                       1905058     2332315
-----------------------------------------------------------------------------------

. 
. expand 2, g(new)
(753 observations created)

. 
. gen outcome= cond(new, outcome2, whrs75)

. 
. gen cons=1

. 
. tobit outcome i.new#(c.nwinc c.wedyrs c.wexper c.wexper#c.wexper c.wifeage c.kl6 c.k618 c.cons), nocons ll(0) vce(robust)

Refining starting values:

Grid node 0:   log likelihood = -10484.914

Fitting full model:

Iteration 0:   log pseudolikelihood = -10484.914  
Iteration 1:   log pseudolikelihood = -10390.013  
Iteration 2:   log pseudolikelihood = -10385.218  
Iteration 3:   log pseudolikelihood = -10385.186  
Iteration 4:   log pseudolikelihood = -10385.186  

Tobit regression                                Number of obs     =      1,506
                                                   Uncensored     =      1,181
Limits: lower = 0                                  Left-censored  =        325
        upper = +inf                               Right-censored =          0

                                                F(  16,   1490)   =     172.03
                                                Prob > F          =     0.0000
Log pseudolikelihood = -10385.186               Pseudo R2         =     0.0682

 ( 1)  [outcome]_cons = 0
---------------------------------------------------------------------------------------
                      |               Robust
              outcome |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
----------------------+----------------------------------------------------------------
          new#c.nwinc |
                   0  |  -9.898906   5.016097    -1.97   0.049    -19.73827   -.0595439
                   1  |    -2.8083   4.995507    -0.56   0.574    -12.60727    6.990673
                      |
         new#c.wedyrs |
                   0  |   92.22536   23.68473     3.89   0.000     45.76641    138.6843
                   1  |   23.37747   24.24086     0.96   0.335    -24.17236     70.9273
                      |
         new#c.wexper |
                   0  |   145.1355   20.42069     7.11   0.000     105.0792    185.1919
                   1  |   25.51334   18.89915     1.35   0.177    -11.55842    62.58511
                      |
new#c.wexper#c.wexper |
                   0  |  -2.091404   .6268803    -3.34   0.001    -3.321066   -.8617426
                   1  |  -.8997776   .6033226    -1.49   0.136     -2.08323    .2836742
                      |
        new#c.wifeage |
                   0  |  -59.32494   7.826249    -7.58   0.000    -74.67657    -43.9733
                   1  |     7.6236   8.594299     0.89   0.375     -9.23461    24.48181
                      |
            new#c.kl6 |
                   0  |  -981.8592   129.8309    -7.56   0.000     -1236.53   -727.1885
                   1  |   159.3098    110.641     1.44   0.150    -57.71899    376.3385
                      |
           new#c.k618 |
                   0  |  -10.78523   42.96339    -0.25   0.802    -95.06037    73.48992
                   1  |   47.07904   44.20634     1.06   0.287    -39.63424    133.7923
                      |
           new#c.cons |
                   0  |   856.3279   486.9199     1.76   0.079    -98.79335    1811.449
                   1  |   1768.561   536.1366     3.30   0.001     716.8983    2820.224
                      |
                _cons |          0  (omitted)
----------------------+----------------------------------------------------------------
        var(e.outcome)|    1828479   57590.13                       1718932     1945008
---------------------------------------------------------------------------------------

[Katy Tabero]

Hi Andrew,

Thank you for pointing this out, your new results look good. Unfortunately, as all 4 of my variables are treatment dummies, this method removes a variable each time for collinearity when I edit it for my data/models (despite adding the no constant option). As you can see below, the results do not match with either of the individual tobit results even closely:

. tobit team_correct endoimposed endonotimposed exoimposed exonotimposed, ll(0) nocons

Refining starting values:

Grid node 0:   log likelihood = -599.69996

Fitting full model:

Iteration 0:   log likelihood = -599.69996  
Iteration 1:   log likelihood = -591.04115  
Iteration 2:   log likelihood = -590.73055  
Iteration 3:   log likelihood = -590.72922  
Iteration 4:   log likelihood = -590.72922  

Tobit regression                                    Number of obs     =    184
                                                           Uncensored =    140
Limits: Lower =    0                                    Left-censored =     44
        Upper = +inf                                   Right-censored =      0

                                                    F(4, 180)         =  30.45
                                                    Prob > F          = 0.0000
Log likelihood = -590.72922                         Pseudo R2         = 0.0613

 ( 1)  [team_correct]_cons = 0
-------------------------------------------------------------------------------------
       team_correct | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
--------------------+----------------------------------------------------------------
        endoimposed |   15.97486   2.055118     7.77   0.000     11.91964    20.03008
     endonotimposed |   6.700489   1.700729     3.94   0.000     3.344559    10.05642
         exoimposed |    11.3178    1.93295     5.86   0.000     7.503645    15.13196
      exonotimposed |   6.537463   1.823243     3.59   0.000     2.939783    10.13514
              _cons |          0  (omitted)
--------------------+----------------------------------------------------------------
 var(e.team_correct)|   149.1048   18.84752                      116.1897    191.3444
-------------------------------------------------------------------------------------

. tobit team_tetris endoimposed endonotimposed exoimposed exonotimposed, ll(0) nocons

Refining starting values:

Grid node 0:   log likelihood = -1574.4946

Fitting full model:

Iteration 0:   log likelihood = -1574.4946  
Iteration 1:   log likelihood = -1573.8712  
Iteration 2:   log likelihood = -1573.8678  
Iteration 3:   log likelihood = -1573.8678  

Tobit regression                                    Number of obs     =    184
                                                           Uncensored =    172
Limits: Lower =    0                                    Left-censored =     12
        Upper = +inf                                   Right-censored =      0

                                                    F(4, 180)         =  49.76
                                                    Prob > F          = 0.0000
Log likelihood = -1573.8678                         Pseudo R2         = 0.0390

 ( 1)  [team_tetris]_cons = 0
------------------------------------------------------------------------------------
       team_tetris | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------------+----------------------------------------------------------------
       endoimposed |   1080.047   360.5997     3.00   0.003     368.5005    1791.593
    endonotimposed |   2884.201    282.266    10.22   0.000     2327.225    3441.177
        exoimposed |   1520.906   333.3076     4.56   0.000     863.2129    2178.599
     exonotimposed |   2427.749   301.1811     8.06   0.000     1833.449    3022.049
             _cons |          0  (omitted)
-------------------+----------------------------------------------------------------
 var(e.team_tetris)|    4501757     491219                       3629716     5583306
------------------------------------------------------------------------------------

. expand 2, g(new)
(184 observations created)

. gen outcome= cond(new, team_correct, team_tetris)

. gen cons=1

. tobit outcome i.new#(endoimposed endonotimposed exoimposed exonotimposed c.cons), nocons ll(0) vce(robust)
note: 1.new#1.endonotimposed omitted because of collinearity.
note: 1.new#1.exoimposed omitted because of collinearity.
note: 1.new#0.exonotimposed omitted because of collinearity.
note: 1.new#1.exonotimposed omitted because of collinearity.
note: 0.new#c.cons omitted because of collinearity.
note: 1.new#c.cons omitted because of collinearity.

Refining starting values:

Grid node 0:   log likelihood = -2785.3011

Fitting full model:

Iteration 0:   log pseudolikelihood = -2785.3011  
Iteration 1:   log pseudolikelihood = -2777.9594  
Iteration 2:   log pseudolikelihood = -2777.7602  
Iteration 3:   log pseudolikelihood = -2777.7597  
Iteration 4:   log pseudolikelihood = -2777.7597  

Tobit regression                                    Number of obs     =    368
                                                           Uncensored =    312
Limits: Lower =    0                                    Left-censored =     56
        Upper = +inf                                   Right-censored =      0

                                                    F(8, 360)         =  27.61
                                                    Prob > F          = 0.0000
Log pseudolikelihood = -2777.7597                   Pseudo R2         = 0.0427

 ( 1)  [outcome]_cons = 0
-------------------------------------------------------------------------------------------
                          |               Robust
                  outcome | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
--------------------------+----------------------------------------------------------------
          new#endoimposed |
  Original observation#1  |   1141.228   321.8277     3.55   0.000     508.3299    1774.127
Duplicated observation#0  |   -106.633   169.5589    -0.63   0.530    -440.0834    226.8175
Duplicated observation#1  |   246.8417   271.1162     0.91   0.363    -286.3288    780.0122
                          |
       new#endonotimposed |
  Original observation#1  |   2903.422   326.7339     8.89   0.000     2260.876    3545.969
Duplicated observation#0  |   -40.3819   153.8199    -0.26   0.793    -342.8804    262.1166
Duplicated observation#1  |          0  (omitted)
                          |
           new#exoimposed |
  Original observation#1  |   1541.231   240.7531     6.40   0.000     1067.772    2014.691
Duplicated observation#0  |  -334.0119   136.6501    -2.44   0.015    -602.7446   -65.27919
Duplicated observation#1  |          0  (omitted)
                          |
        new#exonotimposed |
  Original observation#1  |   2443.039   300.9255     8.12   0.000     1851.246    3034.831
Duplicated observation#0  |          0  (omitted)
Duplicated observation#1  |          0  (omitted)
                          |
               new#c.cons |
    Original observation  |          0  (omitted)
  Duplicated observation  |          0  (omitted)
                          |
                    _cons |          0  (omitted)
--------------------------+----------------------------------------------------------------
            var(e.outcome)|    2454255     239521                       2025663     2973530
-------------------------------------------------------------------------------------------

[Andrew Musau]

If you do not add the -c.- prefix, Stata will treat the variables as categorical, which is not what you want. Also, you don't need to create and interact a constant variable if you are specifying the -nocons- option in the separate regressions.

tobit team_correct endoimposed endonotimposed exoimposed exonotimposed, ll(0) nocons
tobit team_tetris endoimposed endonotimposed exoimposed exonotimposed, ll(0) nocons
expand 2, g(new)
gen outcome= cond(new, team_tetris, team_correct)
tobit outcome i.new#(c.endoimposed c.endonotimposed c.exoimposed c.exonotimposed), nocons ll(0) vce(robust)

[Katy Tabero]

Thank you for this correction and sorry for the delay in response - it has been busy lately!

The results display correctly in the table now (without the omitted variables), however, both sets of results differ from the individual Tobit regressions (as you mentioned they might), so I still cannot be sure that the end results of my tests are correct.

I think unless I can see that both Tobit regressions have correctly run together, showing the same results as when run individually, I cannot avoid my issue sadly. I am not skilled enough in the mechanics of Stata to trust the results when they don't match visually, and as you have said, I am responsible for ensuring the method and results are correct.