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