Manually plotting coefficients after suest command

[Vera Mihailovna]

Hello,

I cannot plot the coefficients by using coefplot after I run suest command as suest command does not store the coefficients. Is there a way of manually plotting effect sizes with confidence intervals in Stata?

[Andrew Musau]

I think that if you followed the FAQ Advice and presented a reproducible example of your problem, you would have made progress by now. See FAQ Advice #12 for details on this and on the need to identify the provenance of community-contributed commands. That said, I do not see what extra thing you get from suest that you do not get from the estimation commands. If it is that suest adjusts a standard VCE for clustering, so can the estimation commands.

[Vera Mihailovna]

Hi Andrew, let me copy past from my previous question:

webuse auto.dta reg price weight length trunk if foreign == 1 est sto m1 reg price weight length trunk if foreign == 0 est sto m2 suest m1 m2

So I would like to obtain and plot the coefficients produced after the suest command

[Andrew Musau]

suest leaves behind results in e(b) and e(V), so theoretically, you can construct the confidence intervals and use a matrix to plot the results. But as I said in #2, the only additional element to suest is that adjusts a standard VCE for clustering, which can be done with the commands. By default, regress reports t-statistics whilst suest reports z-statistics. If you are fussy about this, use glm in place of regress which outputs z-statistics by default. coefplot is from the Stata Journal, as you are asked to explain (FAQ Advice #12). You can change the coefficient labels in the final graph, of course.

sysuse auto, clear
eststo m1: qui reg price weight length trunk if foreign 
eststo m2: qui reg price weight length trunk if !foreign 
suest m1 m2
gen cons=1
eststo m3: glm price i.foreign#(c.weight c.length c.trunk c.cons), nocons robust
coefplot m3, drop(*.foreign#c.cons) baselevels scheme(s1mono)

Res.:

. suest m1 m2

Simultaneous results for m1, m2

                                                Number of obs     =         74

------------------------------------------------------------------------------
             |               Robust
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
m1_mean      |
      weight |   4.508842   1.266932     3.56   0.000     2.025701    6.991983
      length |   22.94693   35.98976     0.64   0.524    -47.59169    93.48556
       trunk |   112.7666   83.84186     1.34   0.179    -51.56045    277.0936
       _cons |  -9211.551   3998.236    -2.30   0.021    -17047.95   -1375.152
-------------+----------------------------------------------------------------
m1_lnvar     |
       _cons |   14.27296   .2011034    70.97   0.000     13.87881    14.66712
-------------+----------------------------------------------------------------
m2_mean      |
      weight |   6.290123   1.876977     3.35   0.001     2.611315    9.968931
      length |  -132.7666   70.44623    -1.88   0.059    -270.8386    5.305508
       trunk |   54.14187   97.99361     0.55   0.581    -137.9221    246.2058
       _cons |   10448.88   7463.625     1.40   0.162    -4179.552    25077.32
-------------+----------------------------------------------------------------
m2_lnvar     |
       _cons |   15.34627   .2138415    71.76   0.000     14.92715    15.76539
------------------------------------------------------------------------------

. 
. gen cons=1

. 
. eststo m3: glm price i.foreign#(c.weight c.length c.trunk c.cons), nocons robust

Iteration 0:   log pseudolikelihood = -661.26069  

Generalized linear models                         Number of obs   =         74
Optimization     : ML                             Residual df     =         66
                                                  Scale parameter =    3792151
Deviance         =  250281937.1                   (1/df) Deviance =    3792151
Pearson          =  250281937.1                   (1/df) Pearson  =    3792151

Variance function: V(u) = 1                       [Gaussian]
Link function    : g(u) = u                       [Identity]

                                                  AIC             =   18.08813
Log pseudolikelihood = -661.2606908               BIC             =   2.50e+08

----------------------------------------------------------------------------------
                 |               Robust
           price |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-----------------+----------------------------------------------------------------
foreign#c.weight |
       Domestic  |   6.290123   1.876977     3.35   0.001     2.611315    9.968931
        Foreign  |   4.508842   1.266932     3.56   0.000     2.025701    6.991983
                 |
foreign#c.length |
       Domestic  |  -132.7666   70.44623    -1.88   0.059    -270.8386    5.305508
        Foreign  |   22.94693   35.98976     0.64   0.524    -47.59169    93.48556
                 |
 foreign#c.trunk |
       Domestic  |   54.14187   97.99361     0.55   0.581    -137.9221    246.2058
        Foreign  |   112.7666   83.84186     1.34   0.179    -51.56045    277.0936
                 |
  foreign#c.cons |
       Domestic  |   10448.88   7463.625     1.40   0.162    -4179.552    25077.32
        Foreign  |  -9211.551   3998.236    -2.30   0.021    -17047.95   -1375.152
----------------------------------------------------------------------------------

[Vera Mihailovna]

Thank you very much for your kind reply, but indeed, this was just a toy example. Not all of the variables exist in two regressions that I run with my data. Basically, I would like to get a comparison of the coefficient estimates across models in which some of the variables are different.

[Andrew Musau]

That's why we insist that your example should be representative of your problem. Looking closer at the coefplot options, all you need are e(b) and e(V), which you get from suest. Multiple equations pose a challenge, but you can overcome this by simply renaming the column names to include the column equation names and uniquely renaming the column equations.

sysuse auto, clear
eststo m1: qui reg price weight length trunk if foreign
eststo m2: qui reg price weight length trunk if !foreign
suest m1 m2
mat b= e(b)
mat V=e(V)
local cnames: colnames b
local eqnames: coleq b
local mycolnames
forval i=1/`=wordcount("`cnames'")'{
    local mycolnames "`mycolnames' `=word("`eqnames'",`i')'`=word("`cnames'",`i')'"
}
local mycoleqnames
forval i=1/`=wordcount("`eqnames'")'{
    local mycoleqnames "`mycoleqnames' `=word("`eqnames'",1)'"
}
mat colnames b= `mycolnames'
mat coleq b= `mycoleqnames'
coefplot mat(b), v(V) drop(*_cons) order(*weight *length *trunk) scheme(s1mono)

Res.:

[Vera Mihailovna]

Dear Andrew thanks a tonne. I thought a toy example would do the job for my question but obviously I was wrong. My apologies for this.
Just a follow-up question. If I would like to plot it only for a single variable, let's say m1_meanweight and m2_meanweight, how should I proceed in this case?

[Andrew Musau]

Last line:

coefplot mat(b), v(V) keep(*weight) scheme(s1mono)

[Vera Mihailovna]

Sorry again for not being clear, but what if in different regressions, I have different number of variables and I might not find a corresponding variable in the other regression.
So for example:

 
 sysuse auto, clear gen exhange_rate= runiform(0, 1) if foreign gen size = runiform(0, 1000) if foreign eststo m1: qui reg price weight length size exhange_rate if foreign eststo m2: qui reg price weight length trunk if !foreign suest m1 m2

What if I have size and inflation variables in the first regression and not in the second? In this case, should I just replace missing variables with 1s or 0s? e.g.

 replace  exchange_rate = 1 if missing(exchange_rate)

[Andrew Musau]

Your code in #9 is all in one line.

sysuse auto, clear
gen exhange_rate= runiform(0, 1) if foreign
gen size = runiform(0, 1000) if foreign
eststo m1: qui reg price weight length size exhange_rate if foreign
eststo m2: qui reg price weight length trunk if !foreign
suest m1 m2

What if I have size and inflation variables in the first regression and not in the second? In this case, should I just replace missing variables with 1s or 0s? e.g.

If I understand you correctly, your question no longer relates to graphing. You have a variable that is observed for a particular subsample but not for another. If you are estimating the samples separately and then using suest, then replacing these values with a constant does nothing as they will be collinear with the intercept and will be dropped. In this case, there is no comparison to make as you do not have information for one group. If pooling the samples, there is some literature relating to censored regressors and suggestions on how to handle them. Unfortunately, this literature never got traction, which signals that the suggested methods are not widely accepted. I would advise that you start a new thread that specifically deals with your question.

[Vera Mihailovna]

Hi Andrew, actually it still relates to the graphing. In my original dateset, I have two samples one for woman the other for man. I am interested in identifying how changing dosage of a drug impacts the mortality rate differently for man and woman. There are different peculiarities that only impact woman (e.g. previous pregnancy) or man (e.g. andropause), therefore I would like to run two separate regressions and then compare the main variable of interest with suest command. Actually, I knew that they would be collinear and won't be estimated but I thought it could actually help with multiple equations problem. If we get back to our toy example, I just want to plot, let's say "weight" variable without necessarily renaming the others.

[Andrew Musau]

Maybe it's me on a Friday evening, but I am not getting it. If we run the code in #10, the results are:

. suest m1 m2

Simultaneous results for m1, m2

                                                Number of obs     =         74

------------------------------------------------------------------------------
             |               Robust
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
m1_mean      |
      weight |   4.047623   1.179271     3.43   0.001     1.736295    6.358951
      length |    49.6821   35.38175     1.40   0.160    -19.66486    119.0291
        size |   .8006065   .7976607     1.00   0.316    -.7627797    2.363993
exhange_rate |   2086.069   752.9904     2.77   0.006     610.2346    3561.903
       _cons |  -12857.83   3627.147    -3.54   0.000     -19966.9   -5748.748
-------------+----------------------------------------------------------------
m1_lnvar     |
       _cons |   14.15744   .1954989    72.42   0.000     13.77427    14.54061
-------------+----------------------------------------------------------------
m2_mean      |
      weight |   6.290123   1.876977     3.35   0.001     2.611315    9.968931
      length |  -132.7666   70.44623    -1.88   0.059    -270.8386    5.305508
       trunk |   54.14187   97.99361     0.55   0.581    -137.9221    246.2058
       _cons |   10448.88   7463.625     1.40   0.162    -4179.552    25077.32
-------------+----------------------------------------------------------------
m2_lnvar     |
       _cons |   15.34627   .2138415    71.76   0.000     14.92715    15.76539
------------------------------------------------------------------------------

We have 7 coefficients in total. "weight" and "length" are present in both samples. Appending the highlighted code in #6, I can plot all 7 coefficients.

Now, you can drop some of the coefficients using the -drop()- or -keep()- options, order them using the -order()- option and even label a specific coefficient, e.g., to "Foreign: Weight" instead of "m1_meanweight".

If we get back to our toy example, I just want to plot, let's say "weight" variable without necessarily renaming the others.

There are 2 weight coefficients. You can plot either or both. You can label them in the axis however you want. By default, as I explained, I need to rename the columns in the matrix so that I am able to plot all the coefficients at once as coefplot (by default) will only plot coefficients corresponding to the first equation. As your final objective is to plot a coefficient, whether or not you rename the matrix columns is irrelevant. What am I not getting?

[Vera Mihailovna]

Dear Andrew,

No, the problem was with my code. Sorry and thanks a lot

[Erdem Yilmaz]

Dear Andrew,

I have a similar problem and saw this thread. Andrew Musau .

Thank you very much for your reply. I get the following error:

mat colnames b= `mycolnames'
m1__t612b:  operator invalid

How could I solve that ?

[Andrew Musau]

Actually, my workaround is not needed. coefplot can handle multiple equations. I will illustrate using the example above.

sysuse auto, clear
eststo m1: qui reg price weight length trunk if foreign
eststo m2: qui reg price weight length trunk if !foreign
eststo m3: suest m1 m2
coefplot m3, keep(*:) drop(_cons) eqlabel(Foreign Domestic) scheme(s1color)

Res.: