gsem and collect table

[Rich Goldstein]

I am using -gsem- to estimate a mediation model where the outcome variable is ordinal; let's use the auto data as an example and there are two
mediators; let's say that the direct effect is the variable "foreign" and the two mediators are price and mpg (note that there is no implication
here that this makes substantive sense); while I can get basic results in a comparative table for several models (say, a model with no mediator
and two models each with one mediator and a final model with both mediators), I cannot figure out how to (1) show the indirect and total effects in my
table (calculated using -nlcom- as shown in the documentation) and (2) I would like to suppress the showing of the 4 cutpoints (or possibly just
move them; the issue here is that if the model with no mediation is shown first, the cutpoints will be automatically shown above the mediators
in the table that has comparative results); (3) for the models with rep78 as the outcome, showing the exponentiated coefficients would be preferred
if possible (results of -estat eform-)

following is the do file and below that is that output from running the do file

sysuse auto, clear

collect clear
collect create Models1
    
    qui collect _r_b _r_p,name(Models1) tag(model[foreign]): gsem (i.foreign -> rep78, family(ordinal) link(logit)), nocapslatent
    qui collect AIC=r(S)[1,5] BIC=r(S)[1,6], name(Models1) tag(model[foreign]): estat ic

    qui collect _r_b _r_p,name(Models1) tag(model[foreignp]): gsem (i.foreign -> price, ) ///
(i.foreign -> rep78 , family(ordinal) link(logit)) ///    
(price -> rep78 , family(ordinal) link(logit)), nocapslatent    
    qui collect AIC=r(S)[1,5] BIC=r(S)[1,6], name(Models1) tag(model[foreignp]): estat ic
    
    qui collect _r_b _r_p, name(Models1) tag(model[foreignm]): gsem (i.foreign -> mpg, ) ///
(i.foreign -> rep78 , family(ordinal) link(logit)) ///
(mpg -> rep78, family(ordinal) link(logit)), nocapslatent
    qui collect AIC=r(S)[1,5] BIC=r(S)[1,6], name(Models1) tag(model[foreignm]): estat ic
    
    qui collect _r_b, name(Models1) tag(model[foreignmp]): gsem (i.foreign -> mpg, ) ///
(i.foreign -> price, ) (i.foreign -> rep78, family(ordinal) link(logit)) /// 
(mpg -> rep78, family(ordinal) link(logit)) ///
(price -> rep78, family(ordinal) link(logit)), nocapslatent
    qui collect AIC=r(S)[1,5] BIC=r(S)[1,6], name(Models1) tag(model[foreignmp]): estat ic
    
collect style showbase off
collect style cell, nformat(%5.3f)
collect style header result, level(hide)
collect style row stack, nobinder
collect style header result[N AIC BIC], level(label)
collect label levels result N "N", modify
*collect label levels result r2_p "Pseudo R2", modify
collect style cell result[N], warn margin(top, width(72) )
*noi collect label list result

noi collect layout (colname#result result[N AIC BIC]) (model)

*example code for indirect and total effects
noi di "indirect effect for price: "
noi nlcom _b[rep78:price]*_b[price:1.foreign], eform
noi di "indirect effect for mpg: "
noi nlcom _b[rep78:mpg]*_b[mpg:1.foreign], eform
noi di "Total effect: "
noi nlcom _b[rep78:1.foreign]+(_b[rep78:price]*_b[price:1.foreign]) ///
+ (_b[rep78:mpg]*_b[mpg:1.foreign]), eform

output from above:

Collection: Models1
      Rows: colname#result result[N AIC BIC]
   Columns: model
   Table 1: 18 x 4

---------------------------------------------------
              | foreign foreignp foreignm foreignmp
--------------+------------------------------------
Car origin    |                                    
  Foreign     |   2.982                            
              |   0.000                            
cut1          |  -3.158   -3.080   -1.885    -0.712
cut2          |  -1.363   -1.283   -0.092     1.091
cut3          |   1.232    1.313    2.525     3.719
cut4          |   3.246    3.327    4.581     5.801
Price         |            0.000              0.000
              |            0.871              0.312
var(e.price)  |          8.6e+06            8.6e+06
Mileage (mpg) |                     0.067     0.098
              |                     0.174     0.095
var(e.mpg)    |                    27.910    27.910
Intercept     |         6072.423   19.827          
              |            0.000    0.000          
N             |  69.000   74.000   74.000    74.000
AIC           | 168.058 1567.282  630.525  2028.731
BIC           | 179.229 1588.018  651.262  2058.684
---------------------------------------------------
r; t=0.15 9:54:22

. 
. *example code for indirect and total effects
. noi di "indirect effect for price: "
indirect effect for price: 
r; t=0.00 9:54:22

. noi nlcom _b[rep78:price]*_b[price:1.foreign], eform

       _nl_1: _b[rep78:price]*_b[price:1.foreign]

------------------------------------------------------------------------------
             |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _nl_1 |   1.029483   .0771883     0.39   0.698     .8887874     1.19245
------------------------------------------------------------------------------
r; t=0.15 9:54:22

. noi di "indirect effect for mpg: "
indirect effect for mpg: 
r; t=0.00 9:54:22

. noi nlcom _b[rep78:mpg]*_b[mpg:1.foreign], eform

       _nl_1: _b[rep78:mpg]*_b[mpg:1.foreign]

------------------------------------------------------------------------------
             |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _nl_1 |   1.626355   .5196207     1.52   0.128     .8694737    3.042107
------------------------------------------------------------------------------
r; t=0.09 9:54:22

. noi di "Total effect: "
Total effect: 
r; t=0.00 9:54:22

. noi nlcom _b[rep78:1.foreign]+(_b[rep78:price]*_b[price:1.foreign]) ///
> + (_b[rep78:mpg]*_b[mpg:1.foreign]), eform

       _nl_1: _b[rep78:1.foreign]+(_b[rep78:price]*_b[price:1.foreign]) + (_b[r
> ep78:mpg]*_b[mpg:1.foreign])

------------------------------------------------------------------------------
             |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _nl_1 |   19.51854   12.52709     4.63   0.000     5.548072    68.66772

[Jeff Pitblado (StataCorp)]

Hi Rich, thanks for the working example.

For (1), it is a matter of using collect get after nlcom to pull the values of interest. I'm assuming you are interested in estimated effect and it's p-value. Add custom levels for coleq and colname in the tags() option to keep your layout from getting too busy.

For (2), I suggest you use dimension coleq with colname in the layout. You will then have control over which equations show up, and the order in which they show up on the table. You can hide the levels of coleq from the row header later if they are redundant.

For (3), some coding is required. gsem uses a not-documented command gsem_display for reporting it's table, as does estat eform. This gsem_display command behaves much like _coef_table, allows you to select which equations to show, and supports option eform. So I would write a program that uses gsem_display to exponentiate the requested equations and stack the resulting r(table) matrix on the r(table) from the not-exponentiated coefficient table. Then post this stacked matrix to r(table) and use collect get r() as you would after commands like estat eform and margins.

Here is my do-file based on your example.

clear all

sysuse auto, clear

* code up my custom coefficient table, since we want some equations to
* be -eform-'d
program MyCoefTable, rclass
    version 17
    syntax [namelist(name=eform)]

    local eqlist : coleq e(b)
    local eqlist : list uniq eqlist
    local eqlist : list eqlist - eform

    tempname T
    if "`eform'" != "" {
        gsem_display, eqselect(`eform') eform noheader nodvheader
        matrix `T' = r(table)
    }
    if "`eqlist'" != "" {
        gsem_display, eqselect(`eqlist') noheader nodvheader
        matrix `T' = nullmat(`T'), r(table)
    }
    return matrix table `T'
    return scalar level = r(N)
    return scalar N = e(N)
end

collect clear
* once you create a collection, it becomes current and you do not need
* the -name()- option unless you are working with more than one
* collection (simultaneously) or if the -collect- command was created
* by the dialog
collect create Models1

* we have a custom table, so no need to -collect- from -gsem- directly
gsem (i.foreign -> rep78, family(ordinal) link(logit)), nocapslatent
* custom table
MyCoefTable rep78
collect get r() , tag(model[foreign])
* collect model fit statistics
collect AIC=r(S)[1,5] BIC=r(S)[1,6],  tag(model[foreign]): estat ic
* collect direct/total effect
nlcom _b[rep78:1.foreign], eform
collect get ///
    _r_b=(r(table)["b",1]) ///
    _r_p=(r(table)["pvalue",1]) ///
    , tag(model[foreign] colname[Total] coleq[Effects])

* some style choices
collect style cell, nformat(%5.3f)
collect style header result, level(hide)
collect style header result[N AIC BIC], level(label)
collect style row stack, nobinder
collect style cell result[N], warn margin(top, width(72) )

* results of interest
collect style autolevels result _r_b _r_p, clear
* precog the order of equations, but we could have done this later
collect style autolevels coleq rep78 price mpg /rep78 / Effects, clear

* arrange table items for the first model; add dimension coleq to the
* layout since some levels of colname appear in more than one equation
collect layout (coleq#colname#result result[N AIC BIC]) (model)

gsem (i.foreign -> price, ) ///
    (i.foreign -> rep78 , family(ordinal) link(logit)) ///
    (price -> rep78 , family(ordinal) link(logit)), nocapslatent
MyCoefTable rep78
collect get r() , tag(model[foreignp])
collect AIC=r(S)[1,5] BIC=r(S)[1,6],  tag(model[foreignp]): estat ic
* collect indirect effect, via price
nlcom _b[rep78:price]*_b[price:1.foreign], eform
collect get ///
    _r_b=(r(table)["b",1]) ///
    _r_p=(r(table)["pvalue",1]) ///
    , tag(model[foreignp] colname[price] coleq[Effects])
* collect total effect
nlcom _b[rep78:1.foreign] + _b[rep78:price]*_b[price:1.foreign], eform
collect get ///
    _r_b=(r(table)["b",1]) ///
    _r_p=(r(table)["pvalue",1]) ///
    , tag(model[foreignp] colname[Total] coleq[Effects])

* verify arrangement still works, some rows may not be in the order we
* want, but we can fix that later
collect preview

gsem (i.foreign -> mpg, ) ///
    (i.foreign -> rep78 , family(ordinal) link(logit)) ///
    (mpg -> rep78, family(ordinal) link(logit)), nocapslatent
MyCoefTable rep78
collect get r() , tag(model[foreignm])
collect AIC=r(S)[1,5] BIC=r(S)[1,6],  tag(model[foreignm]): estat ic
* collect indirect effect, via mpg
nlcom _b[rep78:mpg]*_b[mpg:1.foreign], eform
collect get ///
    _r_b=(r(table)["b",1]) ///
    _r_p=(r(table)["pvalue",1]) ///
    , tag(model[foreignm] colname[mpg] coleq[Effects])
* collect total effect
nlcom _b[rep78:1.foreign] + _b[rep78:mpg]*_b[mpg:1.foreign], eform
collect get ///
    _r_b=(r(table)["b",1]) ///
    _r_p=(r(table)["pvalue",1]) ///
    , tag(model[foreignm] colname[Total] coleq[Effects])

* verify the results we want are still getting reported
collect preview

gsem (i.foreign -> mpg, ) ///
    (i.foreign -> price, ) ///
    (i.foreign -> rep78, family(ordinal) link(logit)) ///
    (mpg -> rep78, family(ordinal) link(logit)) ///
    (price -> rep78, family(ordinal) link(logit)), nocapslatent
MyCoefTable rep78
collect get r() , tag(model[foreignmp])
collect AIC=r(S)[1,5] BIC=r(S)[1,6],  tag(model[foreignmp]): estat ic
* collect indirect effect, via price
nlcom _b[rep78:price]*_b[price:1.foreign], eform
collect get ///
    _r_b=(r(table)["b",1]) ///
    _r_p=(r(table)["pvalue",1]) ///
    , tag(model[foreignmp] colname[price] coleq[Effects])
* collect indirect effect, via mpg
nlcom _b[rep78:mpg]*_b[mpg:1.foreign], eform
collect get ///
    _r_b=(r(table)["b",1]) ///
    _r_p=(r(table)["pvalue",1]) ///
    , tag(model[foreignmp] colname[mpg] coleq[Effects])
* collect total effect
nlcom _b[rep78:1.foreign] ///
    + _b[rep78:price]*_b[price:1.foreign] ///
    + _b[rep78:mpg]*_b[mpg:1.foreign] ///
    , eform
collect get ///
    _r_b=(r(table)["b",1]) ///
    _r_p=(r(table)["pvalue",1]) ///
    , tag(model[foreignmp] colname[Total] coleq[Effects])

* verify the results we want are still getting reported
collect preview

* clean up order of levels in -colname-
collect style autolevels ///
    colname ///
    1.foreign price mpg _cons ///
    cut1 cut2 cut3 cut4 ///
    var(e.price) ///
    var(e.mpg) ///
    Total ///
    , clear
collect preview

* get rid of the ".b" p-values
collect recode result _r_p = junk, ///
    fortags(colname[cut1 cut2 cut3 cut4 var(e.price) var(e.mpg)])
collect preview

Here is the final table.

--------------------------------------------------------
                   | foreign foreignp foreignm foreignmp
-------------------+------------------------------------
Repair record 1978 |
  Car origin       |
    Foreign        |  19.718   19.782   13.451    11.658
                   |   0.000    0.000    0.000     0.000
  Price            |            1.000              1.000
                   |            0.871              0.312
  Mileage (mpg)    |                     1.070     1.103
                   |                     0.174     0.095
Price              |
  Car origin       |
    Foreign        |          312.259            312.259
                   |            0.675              0.675
  Intercept        |         6072.423           6072.423
                   |            0.000              0.000
Mileage (mpg)      |
  Car origin       |
    Foreign        |                     4.946     4.946
                   |                     0.000     0.000
  Intercept        |                    19.827    19.827
                   |                     0.000     0.000
/rep78             |
  cut1             |  -3.158   -3.080   -1.885    -0.712
  cut2             |  -1.363   -1.283   -0.092     1.091
  cut3             |   1.232    1.313    2.525     3.719
  cut4             |   3.246    3.327    4.581     5.801
/                  |
  var(e.price)     |          8.6e+06            8.6e+06
  var(e.mpg)       |                    27.910    27.910
Effects            |
  Price            |            1.004              1.029
                   |            0.879              0.698
  Mileage (mpg)    |                     1.395     1.626
                   |                     0.202     0.128
  Total            |  19.718   19.860   18.762    19.519
                   |   0.000    0.000    0.000     0.000
N                  |  69.000   74.000   74.000    74.000
AIC                | 168.058 1567.282  630.525  2028.731
BIC                | 179.229 1588.018  651.262  2058.684
--------------------------------------------------------

[Rich Goldstein]

Jeff Pitblado (StataCorp) thank you very much - it will take me a while to be sure I fully understand so I can transport to my real data but I am very impressed!