## Using gsem to combine estimation results

**gsem** is a very flexible command that allows us to fit very sophisticated models. However, it is also useful in situations that involve simple models.

For example, when we want to compare parameters among two or more models, we usually use **suest**, which combines the estimation results under one parameter vector and creates a simultaneous covariance matrix of the robust type. This covariance estimate is described in the *Methods and formulas* of **[R] suest** as the robust variance from a “stacked model”. Actually, **gsem** can estimate these kinds of “stacked models”, even if the estimation samples are not the same and eventually overlap. By using the option **vce(robust)**, we can replicate the results from **suest** if the models are available for **gsem**. In addition, **gsem** allows us to combine results from some estimation commands that are not supported by **suest**, like models including random effects.

### Example: Comparing parameters from two models

Let’s consider the **childweight** dataset, described in **[ME] mixed**. Consider the following models, where weights of boys and girls are modeled using the age and the age-squared:

. webuse childweight, clear (Weight data on Asian children) . regress weight age c.age#c.age if girl == 0, noheader ------------------------------------------------------------------------------ weight | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- age | 7.985022 .6343855 12.59 0.000 6.725942 9.244101 | c.age#c.age | -1.74346 .2374504 -7.34 0.000 -2.214733 -1.272187 | _cons | 3.684363 .3217223 11.45 0.000 3.045833 4.322893 ------------------------------------------------------------------------------ . regress weight age c.age#c.age if girl == 1, noheader ------------------------------------------------------------------------------ weight | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- age | 7.008066 .5164687 13.57 0.000 5.982746 8.033386 | c.age#c.age | -1.450582 .1930318 -7.51 0.000 -1.833798 -1.067365 | _cons | 3.480933 .2616616 13.30 0.000 2.961469 4.000397 ------------------------------------------------------------------------------

To test whether birthweights are the same for the two groups, we need to test whether the intercepts in the two regressions are the same. Using **suest**, we would proceed as follows:

. quietly regress weight age c.age#c.age if girl == 0, noheader . estimates store boys . quietly regress weight age c.age#c.age if girl == 1, noheader . estimates store girls . suest boys girls Simultaneous results for boys, girls Number of obs = 198 ------------------------------------------------------------------------------ | Robust | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- boys_mean | age | 7.985022 .4678417 17.07 0.000 7.068069 8.901975 | c.age#c.age | -1.74346 .2034352 -8.57 0.000 -2.142186 -1.344734 | _cons | 3.684363 .1719028 21.43 0.000 3.34744 4.021286 -------------+---------------------------------------------------------------- boys_lnvar | _cons | .4770289 .1870822 2.55 0.011 .1103546 .8437032 -------------+---------------------------------------------------------------- girls_mean | age | 7.008066 .4166916 16.82 0.000 6.191365 7.824766 | c.age#c.age | -1.450582 .1695722 -8.55 0.000 -1.782937 -1.118226 | _cons | 3.480933 .1556014 22.37 0.000 3.17596 3.785906 -------------+---------------------------------------------------------------- girls_lnvar | _cons | .0097127 .1351769 0.07 0.943 -.2552292 .2746545 ------------------------------------------------------------------------------

Invoking an estimation command with the option **coeflegend** will give us a legend we can use to refer to the parameters when we use postestimation commands like **test**.

. suest, coeflegend Simultaneous results for boys, girls Number of obs = 198 ------------------------------------------------------------------------------ | Coef. Legend -------------+---------------------------------------------------------------- boys_mean | age | 7.985022 _b[boys_mean:age] | c.age#c.age | -1.74346 _b[boys_mean:c.age#c.age] | _cons | 3.684363 _b[boys_mean:_cons] -------------+---------------------------------------------------------------- boys_lnvar | _cons | .4770289 _b[boys_lnvar:_cons] -------------+---------------------------------------------------------------- girls_mean | age | 7.008066 _b[girls_mean:age] | c.age#c.age | -1.450582 _b[girls_mean:c.age#c.age] | _cons | 3.480933 _b[girls_mean:_cons] -------------+---------------------------------------------------------------- girls_lnvar | _cons | .0097127 _b[girls_lnvar:_cons] ------------------------------------------------------------------------------ . test _b[boys_mean:_cons] = _b[girls_mean:_cons] ( 1) [boys_mean]_cons - [girls_mean]_cons = 0 chi2( 1) = 0.77 Prob > chi2 = 0.3803

We find no evidence that the intercepts are different.

Now, let’s replicate those results by using the **gsem** command. We generate the variable **weightboy**, a copy of weight for boys and missing otherwise, and the variable **weightgirl**, a copy of weight for girls and missing otherwise.

. quietly generate weightboy = weight if girl == 0 . quietly generate weightgirl = weight if girl == 1 . gsem (weightboy <- age c.age#c.age) (weightgirl <- age c.age#c.age), /// > nolog vce(robust) Generalized structural equation model Number of obs = 198 Log pseudolikelihood = -302.2308 ------------------------------------------------------------------------------- | Robust | Coef. Std. Err. z P>|z| [95% Conf. Interval] -----------------+------------------------------------------------------------- weightboy <- | age | 7.985022 .4678417 17.07 0.000 7.068069 8.901975 | c.age#c.age | -1.74346 .2034352 -8.57 0.000 -2.142186 -1.344734 | _cons | 3.684363 .1719028 21.43 0.000 3.34744 4.021286 -----------------+------------------------------------------------------------- weightgirl <- | age | 7.008066 .4166916 16.82 0.000 6.191365 7.824766 | c.age#c.age | -1.450582 .1695722 -8.55 0.000 -1.782937 -1.118226 | _cons | 3.480933 .1556014 22.37 0.000 3.17596 3.785906 -----------------+------------------------------------------------------------- var(e.weightboy)| 1.562942 .3014028 1.071012 2.280821 var(e.weightgirl)| .978849 .1364603 .7448187 1.286414 ------------------------------------------------------------------------------- . gsem, coeflegend Generalized structural equation model Number of obs = 198 Log pseudolikelihood = -302.2308 ------------------------------------------------------------------------------- | Coef. Legend -----------------+------------------------------------------------------------- weightboy <- | age | 7.985022 _b[weightboy:age] | c.age#c.age | -1.74346 _b[weightboy:c.age#c.age] | _cons | 3.684363 _b[weightboy:_cons] -----------------+------------------------------------------------------------- weightgirl <- | age | 7.008066 _b[weightgirl:age] | c.age#c.age | -1.450582 _b[weightgirl:c.age#c.age] | _cons | 3.480933 _b[weightgirl:_cons] -----------------+------------------------------------------------------------- var(e.weightboy)| 1.562942 _b[var(e.weightboy):_cons] var(e.weightgirl)| .978849 _b[var(e.weightgirl):_cons] ------------------------------------------------------------------------------- . test _b[weightgirl:_cons]= _b[weightboy:_cons] ( 1) - [weightboy]_cons + [weightgirl]_cons = 0 chi2( 1) = 0.77 Prob > chi2 = 0.3803

**gsem** allowed us to fit models on different subsets simultaneously. By default, the model is assumed to be a linear regression, but several links and families are available; for example, you can combine two Poisson models or a multinomial logistic model with a regular logistic model. See **[SEM] sem and gsem** for details.

Here, I use the **vce(robust)** option to replicate the results for **suest**. However, when estimation samples don’t overlap, results from both estimations are assumed to be independent, and thus the option **vce(robust)** is not needed. When performing the estimation without the **vce(robust)** option, the joint covariance matrix will contain two blocks with the covariances from the original models and 0s outside those blocks.

### An example with random effects

The **childweight** dataset contains repeated measures, and it is, in the documentation, analyzed used the **mixed** command, which allows us to account for the intra-individual correlation via random effects.

Now, let’s use the techniques described above to combine results from two random-effects models. Here are the two separate models:

. mixed weight age c.age#c.age if girl == 0 || id:, nolog Mixed-effects ML regression Number of obs = 100 Group variable: id Number of groups = 34 Obs per group: min = 1 avg = 2.9 max = 5 Wald chi2(2) = 1070.28 Log likelihood = -149.05479 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ weight | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- age | 8.328882 .4601093 18.10 0.000 7.427084 9.230679 | c.age#c.age | -1.859798 .1722784 -10.80 0.000 -2.197458 -1.522139 | _cons | 3.525929 .2723617 12.95 0.000 2.99211 4.059749 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ id: Identity | var(_cons) | .7607779 .2439115 .4058409 1.426133 -----------------------------+------------------------------------------------ var(Residual) | .7225673 .1236759 .5166365 1.010582 ------------------------------------------------------------------------------ LR test vs. linear regression: chibar2(01) = 30.34 Prob >= chibar2 = 0.0000 . mixed weight age c.age#c.age if girl == 1 || id:, nolog Mixed-effects ML regression Number of obs = 98 Group variable: id Number of groups = 34 Obs per group: min = 1 avg = 2.9 max = 5 Wald chi2(2) = 2141.72 Log likelihood = -114.3008 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ weight | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- age | 7.273082 .3167266 22.96 0.000 6.652309 7.893854 | c.age#c.age | -1.538309 .118958 -12.93 0.000 -1.771462 -1.305156 | _cons | 3.354834 .2111793 15.89 0.000 2.94093 3.768738 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ id: Identity | var(_cons) | .6925554 .1967582 .396848 1.208606 -----------------------------+------------------------------------------------ var(Residual) | .3034231 .0535359 .2147152 .4287799 ------------------------------------------------------------------------------ LR test vs. linear regression: chibar2(01) = 47.42 Prob >= chibar2 = 0.0000

Random effects can be included in a **gsem** model by incorporating latent variables at the group level; these are the latent variables **M1[id]** and **M2[id]** below. By default, **gsem** will try to estimate a covariance when it sees two latent variables at the same level. This can be easily solved by restricting this covariance term to 0. Option **vce(robust)** should be used whenever we want to produce the mechanism used by **suest**.

. gsem (weightboy <- age c.age#c.age M1[id]) /// > (weightgirl <- age c.age#c.age M2[id]), /// > cov(M1[id]*M2[id]@0) vce(robust) nolog Generalized structural equation model Number of obs = 198 Log pseudolikelihood = -263.35559 ( 1) [weightboy]M1[id] = 1 ( 2) [weightgirl]M2[id] = 1 (Std. Err. adjusted for clustering on id) ------------------------------------------------------------------------------- | Robust | Coef. Std. Err. z P>|z| [95% Conf. Interval] -----------------+------------------------------------------------------------- weightboy <- | age | 8.328882 .4211157 19.78 0.000 7.50351 9.154253 | c.age#c.age | -1.859798 .1591742 -11.68 0.000 -2.171774 -1.547823 | M1[id] | 1 (constrained) | _cons | 3.525929 .1526964 23.09 0.000 3.22665 3.825209 -----------------+------------------------------------------------------------- weightgirl <- | age | 7.273082 .3067378 23.71 0.000 6.671887 7.874277 | c.age#c.age | -1.538309 .120155 -12.80 0.000 -1.773808 -1.30281 | M2[id] | 1 (constrained) | _cons | 3.354834 .1482248 22.63 0.000 3.064319 3.64535 -----------------+------------------------------------------------------------- var(M1[id])| .7607774 .2255575 .4254915 1.360268 var(M2[id])| .6925553 .1850283 .4102429 1.169144 -----------------+------------------------------------------------------------- var(e.weightboy)| .7225674 .1645983 .4623572 1.129221 var(e.weightgirl)| .3034231 .0667975 .1970877 .4671298 -------------------------------------------------------------------------------

Above, we have the joint output from the two models, which would allow us to perform tests among parameters in both models. Notice that option **vce(robust)** implies that standard errors will be clustered on the groups determined by **id**.

**gsem**, when called with the **vce(robust)** option, will complain if there are inconsistencies among the groups in the models (for example, if the random effects in both models were crossed).

### Checking that you are fitting the same model

In the previous model, **gsem**‘s default covariance structure included a term that wasn’t in the original two models, so we needed to include an additional restriction. This can be easy to spot in a simple model, but if you don’t want to rely just on a visual inspection, you can write a small loop to make sure that all the estimates in the joint model are actually also in the original models.

Let’s see an example with random effects, this time with overlapping data.

. *fit first model and save the estimates . gsem (weightboy <- age c.age#c.age M1[id]), nolog Generalized structural equation model Number of obs = 100 Log likelihood = -149.05479 ( 1) [weightboy]M1[id] = 1 ------------------------------------------------------------------------------- | Coef. Std. Err. z P>|z| [95% Conf. Interval] ----------------+-------------------------------------------------------------- weightboy <- | age | 8.328882 .4609841 18.07 0.000 7.425369 9.232394 | c.age#c.age | -1.859798 .1725233 -10.78 0.000 -2.197938 -1.521659 | M1[id] | 1 (constrained) | _cons | 3.525929 .2726322 12.93 0.000 2.99158 4.060279 ----------------+-------------------------------------------------------------- var(M1[id])| .7607774 .2439114 .4058407 1.426132 ----------------+-------------------------------------------------------------- var(e.weightboy)| .7225674 .1236759 .5166366 1.010582 ------------------------------------------------------------------------------- . mat b1 = e(b) . *fit second model and save the estimates . gsem (weight <- age M2[id]), nolog Generalized structural equation model Number of obs = 198 Log likelihood = -348.32402 ( 1) [weight]M2[id] = 1 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- weight <- | age | 3.389281 .1152211 29.42 0.000 3.163452 3.615111 | M2[id] | 1 (constrained) | _cons | 5.156913 .1803059 28.60 0.000 4.80352 5.510306 -------------+---------------------------------------------------------------- var(M2[id])| .6076662 .2040674 .3146395 1.173591 -------------+---------------------------------------------------------------- var(e.weight)| 1.524052 .1866496 1.198819 1.937518 ------------------------------------------------------------------------------ . mat b2 = e(b) . *stack estimates from first and second models . mat stacked = b1, b2 . *estimate joint model and save results . gsem (weightboy <- age c.age#c.age M1[id]) /// > (weight <- age M2[id]), cov(M1[id]*M2[id]@0) vce(robust) nolog Generalized structural equation model Number of obs = 198 Log pseudolikelihood = -497.37881 ( 1) [weightboy]M1[id] = 1 ( 2) [weight]M2[id] = 1 (Std. Err. adjusted for clustering on id) ------------------------------------------------------------------------------- | Robust | Coef. Std. Err. z P>|z| [95% Conf. Interval] ----------------+-------------------------------------------------------------- weightboy <- | age | 8.328882 .4211157 19.78 0.000 7.50351 9.154253 | c.age#c.age | -1.859798 .1591742 -11.68 0.000 -2.171774 -1.547823 | M1[id] | 1 (constrained) | _cons | 3.525929 .1526964 23.09 0.000 3.22665 3.825209 ----------------+-------------------------------------------------------------- weight <- | age | 3.389281 .1157835 29.27 0.000 3.16235 3.616213 | M2[id] | 1 (constrained) | _cons | 5.156913 .1345701 38.32 0.000 4.89316 5.420665 ----------------+-------------------------------------------------------------- var(M1[id])| .7607774 .2255575 .4254915 1.360268 var(M2[id])| .6076662 .1974 .3214791 1.148623 ----------------+-------------------------------------------------------------- var(e.weightboy)| .7225674 .1645983 .4623572 1.129221 var(e.weight)| 1.524052 .1705637 1.223877 1.897849 ------------------------------------------------------------------------------- . mat b = e(b) . *verify that estimates from the joint model are the same as . *from models 1 and 2 . local stripes : colfullnames(b) . foreach l of local stripes{ 2. matrix r1 = b[1,"`l'"] 3. matrix r2 = stacked[1,"`l'"] 4. assert reldif(el(r1,1,1), el(r2,1,1))<1e-5 5. }

The loop above verifies that all the labels in the second model correspond to estimates in the first and that the estimates are actually the same. If you omit the restriction for the variance in the joint model, the **assert** command will produce an error.

### Technical note

As documented in **[U] 20.21.2 Correlated errors: Cluster-robust standard errors**, the formula for the robust estimator of the variance is

\[

V_{robust} = \hat V(\sum_{j=1}^N u'_ju_j) \hat V

\]

where \(N\) is the number of observations, \(\hat V\) is the conventional estimator of the variance, and for each observation \(j\), \(u_j\) is a row vector (with as many columns as parameters), which represents the contribution of this observation to the gradient. (If we stack the rows \(u_j\), the columns of this matrix are the scores.)

When we apply **suest**, the matrix \(\hat V\) is constructed as the stacked block-diagonal conventional variance estimates from the original submodels; this is the variance you will see if you apply **gsem** to the joint model without the **vce(robust)** option. The \(u_j\) values used by **suest** are now the values from both estimations, so we have as many \(u_j\) values as the sum of observations in the two original models and each row contains as many columns as the total number of parameters in both models. This is the exact operation that **gsem, vce(robust)** does.

When random effects are present, standard errors will be clustered on groups. Instead of observation-level contributions to the gradient, we would use cluster-level contributions. This means that observations in the two models would need to be clustered in a consistent manner; observations that are common to the two estimations would need to be in the same cluster in the two estimations.