### Stata Code Fragment Compute SRMR from SEM using Mata

Version info: Code for this page was tested in Stata 12.1.

This code fragment page shows an example using Mata to write a function that calculates the SRMR by comparing the expected covariance from a saturated model to that of the model of interest. Thus it generalizes to models with non complete data.

mata
/* mata function to calculate the SRMR given
the expected covariance matrix under the
saturated model and model of interest
*/
real scalar srmr(string scalar saturated, ///
string scalar model)
{
sat = st_matrix(saturated) /*read in matrices */
mod = st_matrix(model)
res = corr(sat) :- corr(mod) /*covar to cor and diff*/
n = cols(res) /*n manifest vars*/
res = lowertriangle(res) /*extract lower triangle including diag*/
std = sqrt((2 * sum(res :* res))/(n * (n + 1))) /*the SRMR*/
return(std)
}
/* function to calculate covariance of data matrix using
an N rather than N - 1 matrix to match Stata's sem
*/
real matrix cov(string scalar varlist)
{
X = st_data(., varlist)
n = rows(X)
one = J(n, 1, 1)
X = X - one * one' * X :* (1/n)
Sigma = X'X :* (1/n)
return(Sigma)
}
end

/*load a built in Stata dataset*/
sysuse auto, clear

/*Our first goal is to demonstrate that
A: under a properly parameterized saturated model
with complete data, $E(\hat{\Sigma} | \theta) = \Sigma$
B: ignoring the first moments, the above mata function,
srmr, is equivalent to Stata's internal calculations
*/

/* saturated model */
quietly sem (<- price mpg weight), nomeans
quietly estat framework, fitted
mat satSigma = r(Sigma) /*store model expected covariance matrix*/

/*note that these two covariance matrices are identical */
mat list satSigma /*from the model*/

symmetric satSigma[3,3]
observed:   observed:   observed:
price         mpg      weight
observed:price   8581964.8
observed:mpg   -7888.225   33.019722
observed:weight     1217990  -3580.3798   595867.28

mata cov("price mpg weight") /*from the data*/

[symmetric]
1              2              3
+----------------------------------------------+
1 |   8581964.812                                |
2 |  -7888.224982    33.01972243                 |
3 |   1217990.004   -3580.379839    595867.2754  |
+----------------------------------------------+

/*fit the model of interest and extract the
expected covariance matrix.  This is a dull
model, but it serves its pedogogical purpose*/
quietly sem (price mpg <- weight), nomeans
quietly estat framework, fitted
mat mSigma = r(Sigma)

/*because we have complete data, we can compare
Stata's SRMR estimate with that of the srmr function*/
quietly estat gof, stats(res)
display r(srmr)

.01381637

mata srmr("satSigma", "mSigma")

.0138163701

/*This method generalizes to non complete data
however, the 'sample' covariance matrix, is not longer observed
but estimated, and note we do not residuals separately for
different missing data patterns*/

/*add some missingness (how much is irrelevant)*/
replace price = . in 1/10

/*fit a saturated model on all data using maximum likelihood
again we use estat framework, fitted to get the expected
covariance matrix under the model which becomes our reference
for srmr*/
quietly sem (<- price mpg weight), method(mlmv)
quietly estat framework, fitted
mat satSigma = r(Sigma)

/*fit model of interest on all data using maximum likelihood*/
quietly sem (price mpg <- weight), method(mlmv)
quietly estat framework, fitted
mat mSigma = r(Sigma)

/*calculate SRMR, in this case the standardized root
mean residual from a saturated model (not necessarily
reality)*/
mata srmr("satSigma", "mSigma")
.0154135804


The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.