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