## Calculating power using Monte Carlo simulations, part 5: Structural equation models

In our last four posts in this series, we showed you how to calculate power for a *t* test using Monte Carlo simulations, how to integrate your simulations into Stata’s **power** command, and how to do this for linear and logistic regression models and multilevel models. In today’s post, I’m going to show you how to estimate power for structural equation models (SEM) using simulations.

Our goal is to write a program that will calculate power for a given SEM at different sample sizes. We’ll follow the same general procedure as the previous two posts, but the way we’ll go about simulating data is a bit different. Rather than individually simulating each variable for our specified model, we’ll be simulating all our variables simultaneously from a given covariance matrix. Means for each of the variables can also be used to simulate the data if your SEM has a mean structure, such as in group analysis or growth curve analysis.

There are three ways you can obtain a covariance matrix to simulate SEM data:

- Using a covariance matrix published in a paper or other source.
- Using the reticular action model (RAM) to derive the model-implied covariance matrix using expected parameter estimates.
- Extracting the model-implied covariance matrix after performing an
**sem**analysis in Stata either on your own pilot study data or on another data source.

The RAM method will be demonstrated at the end of this post. Extracting model-implied covariances after using the **sem** command will be demonstrated below. Whichever method you choose to obtain a covariance matrix to simulate from, the remainder of our procedure will be the same as in the previous two posts. We will work through the following steps:

- Obtain or derive a covariance matrix (and means, if applicable) that corresponds with your hypothesized model under the alternative hypothesis.
- Simulate a single dataset and fit the model.
- Write a program to create the datasets, fit the models, and use
**simulate**to test the program. - Write a program called
**power_cmd_simsem**that allows you to run your simulations with**power**. - Optional: Write a program called
**power_cmd_simsem_init**so that you can see convergence rates for your model at different sample sizes.

We are planning a new study to evaluate the interaction effect between age and sex on health. We can use the NHANES dataset to obtain a reasonable covariance matrix from which we can simulate new data. We will define health as a latent variable measured by systolic blood pressure (**bpsystol**), diastolic blood pressure (**bpdiast**), serum cholesterol (**tcresult**), and serum triglyercides (**tgresult**). Specifically, we would like to fit the following model:

**Figure 1: Path diagram of the hypothesized model**

**Step 1: Obtain or derive a covariance matrix (and means, if applicable) that corresponds with your hypothesized model**

We are going to use the model-implied covariance matrix from fitting the above model to the NHANES data. First, we need to load the dataset, create the interaction variable, and then fit our model.

. webuse nhanes2 . gen ageXfemale = age*female . sem (Health -> bpsystol bpdiast tcresult tgresult) > (age female ageXfemale -> Health), nolog (5301 observations with missing values excluded) Endogenous variables Measurement: bpsystol bpdiast tcresult tgresult Latent: Health Exogenous variables Observed: age female ageXfemale Structural equation model Number of obs = 5,050 Estimation method: ml Log likelihood = -140667.36 ( 1) [bpsystol]Health = 1 ------------------------------------------------------------------------------ | OIM | Coefficient std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- Structural | Health | age | .4371085 .0237304 18.42 0.000 .3905977 .4836193 female | -20.89331 1.67081 -12.50 0.000 -24.16804 -17.61859 ageXfemale | .3395419 .0328812 10.33 0.000 .2750961 .4039878 -------------+---------------------------------------------------------------- Measurement | bpsystol | Health | 1 (constrained) _cons | 111.2516 1.207464 92.14 0.000 108.885 113.6182 -----------+---------------------------------------------------------------- bpdiast | Health | .4404182 .0096514 45.63 0.000 .4215018 .4593346 _cons | 72.88132 .5554644 131.21 0.000 71.79263 73.97001 -----------+---------------------------------------------------------------- tcresult | Health | .6630559 .0373398 17.76 0.000 .5898712 .7362406 _cons | 203.6104 1.217816 167.19 0.000 201.2235 205.9973 -----------+---------------------------------------------------------------- tgresult | Health | 1.122149 .0726778 15.44 0.000 .9797029 1.264594 _cons | 123.03 2.275853 54.06 0.000 118.5694 127.4906 -------------+---------------------------------------------------------------- var(e.bpsy~l)| 67.62701 8.855679 52.31863 87.41462 var(e.bpdi~t)| 80.1877 2.164596 76.05544 84.54446 var(e.tcre~t)| 2083.505 42.75414 2001.372 2169.01 var(e.tgre~t)| 8724.877 176.5499 8385.617 9077.862 var(e.Health)| 340.7899 11.08864 319.7351 363.2312 ------------------------------------------------------------------------------ LR test of model vs. saturated: chi2(11) = 1542.14 Prob > chi2 = 0.0000

Age, sex, and their interaction are all significant predictors of Health when performed on 5,050 observations, but perhaps you don’t have the resources to collect a sample that large. How large does your sample need to be to have 80% power to detect an interaction between age and sex? To simulate data for our power analysis, we can use the model-implied covariance reported with the **fitted** option of the **estat framework** command. This will provide a lot of output, so we’ll run it **quietly** and just return what we’re interested in: the fitted covariance and the fitted means. Our hypothesized SEM doesn’t include mean structure, so including the means for our power analysis is unnecessary. However, I will include them in all the steps below so that our code can be generalized to other types of SEMs.

. quietly estat framework, fitted . matrix list r(mu) r(mu)[1,8] Observed: Observed: Observed: Observed: Latent: Observed: bpsystol bpdiast tcresult tgresult Health age mu 129.84614 81.070693 215.9396 143.89584 18.594545 47.941386 Observed: Observed: female ageXfemale mu .51663366 24.836832 . matrix list r(Sigma) symmetric r(Sigma)[8,8] Observed: Observed: Observed: Observed: bpsystol bpdiast tcresult tgresult Observed:bpsystol 532.05333 Observed:bpdiast 204.54181 170.27163 Observed:tcresult 307.94061 135.62265 2287.6872 Observed:tgresult 521.15538 229.52632 345.55515 9309.6905 Latent:Health 464.42632 204.54181 307.94061 521.15538 Observed:age 179.49835 79.05434 119.01744 201.42383 Observed:female -1.1112203 -.48940167 -.73680119 -1.2469544 Observed:ageXfemale 64.672663 28.483019 42.881591 72.572343 Latent: Observed: Observed: Observed: Health age female ageXfemale Latent:Health 464.42632 Observed:age 179.49835 293.25241 Observed:female -1.1112203 .06869774 .24972332 Observed:ageXfemale 64.672663 155.35796 12.005288 729.20189

Notice that both our covariance and our means include observed as well as latent variables. We just want the information that corresponds with the observed variables, so we’ll need to extract just those portions of these matrices and save them to a new matrix. We can specify a range of rows or columns using **..** between the starting and ending row or column. When concatenating matrices together, commas are used to signify a new column and backslashes (**\**) are used to signify a new row.

. matrix mu = r(mu)[1,1..4],r(mu)[1,6..8] . matrix Sigma = r(Sigma)[1..4,1..4],r(Sigma)[1..4,6..8]\r(Sigma)[6..8,1..4], > r(Sigma)[6..8,6..8] . matrix list mu mu[1,7] Observed: Observed: Observed: Observed: Observed: Observed: bpsystol bpdiast tcresult tgresult age female mu 129.84614 81.070693 215.9396 143.89584 47.941386 .51663366 Observed: ageXfemale mu 24.836832 . matrix list Sigma symmetric Sigma[7,7] Observed: Observed: Observed: Observed: bpsystol bpdiast tcresult tgresult Observed:bpsystol 532.05333 Observed:bpdiast 204.54181 170.27163 Observed:tcresult 307.94061 135.62265 2287.6872 Observed:tgresult 521.15538 229.52632 345.55515 9309.6905 Observed:age 179.49835 79.05434 119.01744 201.42383 Observed:female -1.1112203 -.48940167 -.73680119 -1.2469544 Observed:ageXfemale 64.672663 28.483019 42.881591 72.572343 Observed: Observed: Observed: age female ageXfemale Observed:age 293.25241 Observed:female .06869774 .24972332 Observed:ageXfemale 155.35796 12.005288 729.20189

Now we have the covariance matrix and means we need to simulate our data!

**Step 2: Simulate a single dataset assuming the alternative hypothesis, and fit the model**

Next we create a simulated dataset from our covariance matrix (and means) using the **drawnorm** command. **drawnorm** simulates a variable or set of variables based on sample size, means, and covariance. Here we’ll use a sample size of 200.

. set seed 12345 . clear . drawnorm bpsystol bpdiast tcresult tgresult age female ageXfemale, > n(200) means(mu) cov(Sigma) (obs 200)

We can then use **sem** to fit the hypothesized model to our simulated data.

. sem (Health -> bpsystol bpdiast tcresult tgresult) > (age female ageXfemale -> Health), nolog Endogenous variables Measurement: bpsystol bpdiast tcresult tgresult Latent: Health Exogenous variables Observed: age female ageXfemale Structural equation model Number of obs = 200 Estimation method: ml Log likelihood = -5599.7049 ( 1) [bpsystol]Health = 1 ------------------------------------------------------------------------------ | OIM | Coefficient std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- Structural | Health | age | .5016981 .1243038 4.04 0.000 .2580671 .7453291 female | -18.12394 8.529735 -2.12 0.034 -34.84191 -1.405964 ageXfemale | .2329537 .1611524 1.45 0.148 -.0828992 .5488066 -------------+---------------------------------------------------------------- Measurement | bpsystol | Health | 1 (constrained) _cons | 111.2897 6.191827 17.97 0.000 99.15397 123.4255 -----------+---------------------------------------------------------------- bpdiast | Health | .4507904 .0627156 7.19 0.000 .3278701 .5737107 _cons | 73.00124 3.087596 23.64 0.000 66.94967 79.05282 -----------+---------------------------------------------------------------- tcresult | Health | .5630107 .183409 3.07 0.002 .2035357 .9224857 _cons | 207.841 6.011418 34.57 0.000 196.0588 219.6231 -----------+---------------------------------------------------------------- tgresult | Health | .8698218 .3926678 2.22 0.027 .1002072 1.639436 _cons | 126.6153 11.5315 10.98 0.000 104.0139 149.2166 -------------+---------------------------------------------------------------- var(e.bpsy~l)| 137.2281 53.00586 64.36608 292.5695 var(e.bpdi~t)| 75.17043 12.98479 53.58117 105.4586 var(e.tcre~t)| 2235.489 226.4431 1832.949 2726.432 var(e.tgre~t)| 8917.533 902.1103 7313.681 10873.1 var(e.Health)| 303.1287 58.65825 207.4491 442.9376 ------------------------------------------------------------------------------ LR test of model vs. saturated: chi2(11) = 13.10 Prob > chi2 = 0.2866

To test the null hypothesis that the interaction term equals zero, it will be easier to conduct a Wald test using the **test** command with **sem** than using the **lrtest** command used in previous posts. In our case, we could have extracted the *p*-value from the table above, but this step can be adapted to test any hypothesis of interest.

. test ageXfemale ( 1) [Health]ageXfemale = 0 chi2( 1) = 2.09 Prob > chi2 = 0.1483

The *p*-value for our test is 0.148, so we would not be able to reject the null hypothesis that the interaction term equals zero.

As discussed in the previous post, the model won’t always converge when fit to simulated datasets. Like **mixed**, **sem** will store a value of 1 to **e(converged)** if the model converges and 0 otherwise. We will store the value of **e(converged)** from the model to a local macro named **conv** to keep track of the convergence of our models. We’ll also create a local variable, **reject**, that will keep track of whether the estimated interaction term is significant. The value that is stored in **reject** is determined by the conditional function **cond()**. If the model converged, then **r(p)<0.05** will be evaluated. A value of 1 will be stored to **reject** if the Wald test *p*-value, **r(p)**, is less than the alpha level, here specified as 0.05, and 0 otherwise. If the model failed to converge, then a missing value will be stored in **reject**.

. sem (Health -> bpsystol bpdiast tcresult tgresult) > (age female ageXfemale -> Health) . local conv = e(converged) . test ageXfemale . local reject = cond(`conv'==1, r(p)<.05, .)

**Step 3: Write a program to create the datasets, fit the models, and use simulate to test the program**

Next let’s write a program that creates datasets under the alternative hypothesis, fits **sem** models, tests the null hypothesis of interest, and uses **simulate** to run many iterations of the program. The primary difference between this program and the programs we have written in previous posts is that this program needs to accept matrices as input arguments. This can be a little tricky. What we need to do is specify the type of input as a string rather than a matrix. Then within the program we will define our matrices by their name, passed to the program as a string, that is, **mat C = `cov’**. Finally, we can use these matrices to simulate our data with **drawnorm**, fit our model, and test the null hypothesis. We’ve also added **capture** in front of the **sem** command to capture errors in case the model doesn’t converge. The code block below contains the syntax for this program, called **simsem**.

capture program drop simsem program simsem, rclass version 17 // PARSE INPUT syntax, n(integer) /// cov(string) /// [ means(string) /// alpha(real 0.05) ] // COMPUTE POWER quietly { drop _all mat C =`cov' mat m = `means' drawnorm bpsystol bpdiast tcresult tgresult age female ageXfemale, /// n(`n') means(m) cov(C) capture sem (Health -> bpsystol bpdiast tcresult tgresult) /// (age female ageXfemale -> Health) local conv = e(converged) test ageXfemale local reject = cond(`conv'==1, r(p)<`alpha', .) } // RETURN RESULTS return scalar reject = `reject' return scalar conv = `conv' end

We then use **simulate** to run **simsem** 10 times using our covariance matrix and means for a sample size of 200.

. simulate reject=r(reject) converged=r(conv), reps(10) seed(12345): > simsem, n(200) means(mu) cov(Sigma) Command: simsem, n(200) means(mu) cov(Sigma) reject: r(reject) converged: r(conv) Simulations (10) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 ...x...x..

**simulate** saved two variables: **reject** and **converged**. The mean of **converged** is the convergence rate. The mean of **reject** is the power to test the null hypothesis that the age by sex interaction term equals zero.

. summarize reject Variable | Obs Mean Std. dev. Min Max -------------+--------------------------------------------------------- reject | 8 .5 .5345225 0 1 . summarize converged Variable | Obs Mean Std. dev. Min Max -------------+--------------------------------------------------------- converged | 10 .8 .421637 0 1

With a sample size of 200, the model converged in 8 of the 10 repetitions (80%) and had a power of 50% to detect the interaction effect on Health.

**Step 4: Write a program called power_cmd_simsem, which allows you to run your simulations with power**

We could stop with our quick simulation if we were interested only in a specific set of assumptions. But it’s easy to write an additional program named **power_cmd_simsem** that will allow us to use Stata’s **power** command to create tables and graphs for a range of sample sizes. We just need to include the input syntax as we did in the **simsem** command, **simulate** with the **simsem** command, and return the results.

capture program drop power_cmd_simsem program power_cmd_simsem, rclass version 17 // PARSE INPUT syntax, n(integer) /// cov(string) /// [ alpha(real 0.05) /// means(string) /// reps(integer 100) ] // COMPUTE POWER quietly simulate reject=r(reject) converged=r(conv), reps(`reps'): /// simsem, n(`n') means(`means') cov(`cov') alpha(`alpha') // RETURN RESULTS return scalar N = `n' return scalar alpha = `alpha' quietly summarize reject return scalar power = r(mean) quietly summarize converged return scalar conv_rate = r(mean) end

**Step 5: (Optional) Write a program called power_cmd_simsem_init so that you can see convergence rates for your model at different sample sizes.**

Convergence is often an issue in SEM simulation. If your model is having trouble converging at smaller sample sizes, it may be useful to add a convergence rate column to the power output table.

capture program drop power_cmd_simsem_init program power_cmd_simsem_init, sclass version 17 sreturn clear // ADD COLUMNS TO THE OUTPUT TABLE sreturn local pss_colnames "conv_rate" end

**Using power simsem**

At last, we can use **power simsem** to simulate power and convergence rate for a range of sample sizes. The example below simulates power for sample sizes ranging from 200 to 500, using 100 repetitions each.

. power simsem, n(200(100)500) means(mu) cov(Sigma) reps(100) > table(N conv_rate power) graph xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx Estimated power Two-sided test +---------------------------+ | N conv_rate power | |---------------------------| | 200 .86 .5 | | 300 .92 .7283 | | 400 .93 .8602 | | 500 .97 .8866 | +---------------------------+

**Figure 2: Estimated power for a structural equation model**

The table and graph above indicate that 80% power is achieved around a sample size of 350.

The procedure demonstrated here can be used to perform a power analysis on any SEM. You just need to define your model of interest and simulate data based on a covariance matrix. Including the previous posts in this series, we have now given examples of how you can use Stata to perform power analysis by simulation for a variety of models. You can use similar methods for virtually any power analysis you need to perform.

**Deriving covariances based on SEM path coefficients using the RAM method**

As discussed at the beginning of this post, you can use the RAM to directly specify population parameters or use results from a publication to obtain a covariance matrix. The RAM method uses a set of matrices to represent the model, which can be combined using matrix algebra to derive the corresponding covariance matrix. Three matrices are required: \(\mathbf{A}\), \(\mathbf{S}\), and \(\mathbf{F}\). These contain the path coefficients, factor loadings, variances, and covariances of the model you would like to simulate data from. Additionally, a matrix of means, \(\mathbf{M}\), is required if the SEM involves mean structure.

To demonstrate how to use model results to build these matrices, we will use the results from our data analysis on the NHANES datset. We will use the **noxconditional** option in the **sem** command because we need the estimated variances or covariances of the exogenous variables in order to create our matrices.

. webuse nhanes2, clear . gen ageXfemale = age*female . sem (Health -> bpsystol bpdiast tcresult tgresult) > (age female ageXfemale -> Health), noxconditional nolog (5301 observations with missing values excluded) Endogenous variables Measurement: bpsystol bpdiast tcresult tgresult Latent: Health Exogenous variables Observed: age female ageXfemale Structural equation model Number of obs = 5,050 Estimation method: ml Log likelihood = -140667.36 ( 1) [bpsystol]Health = 1 ------------------------------------------------------------------------------ | OIM | Coefficient std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- Structural | Health | age | .4371086 .0237305 18.42 0.000 .3905977 .4836194 female | -20.89331 1.67081 -12.50 0.000 -24.16804 -17.61858 ageXfemale | .3395419 .0328812 10.33 0.000 .275096 .4039878 -------------+---------------------------------------------------------------- Measurement | bpsystol | Health | 1 (constrained) _cons | 111.2516 1.207464 92.14 0.000 108.885 113.6182 -----------+---------------------------------------------------------------- bpdiast | Health | .4404182 .0096514 45.63 0.000 .4215018 .4593346 _cons | 72.88132 .5554643 131.21 0.000 71.79263 73.97001 -----------+---------------------------------------------------------------- tcresult | Health | .6630558 .0373398 17.76 0.000 .5898711 .7362405 _cons | 203.6104 1.217816 167.19 0.000 201.2235 205.9973 -----------+---------------------------------------------------------------- tgresult | Health | 1.122149 .0726778 15.44 0.000 .9797027 1.264594 _cons | 123.03 2.275853 54.06 0.000 118.5694 127.4906 -------------+---------------------------------------------------------------- mean(age)| 47.94139 .2409767 198.95 0.000 47.46908 48.41369 mean(female)| .5166337 .0070321 73.47 0.000 .502851 .5304163 mean(ageXf~e)| 24.83683 .3799953 65.36 0.000 24.09205 25.58161 -------------+---------------------------------------------------------------- var(e.bpsy~l)| 67.62698 8.85568 52.3186 87.41459 var(e.bpdi~t)| 80.1877 2.164596 76.05545 84.54447 var(e.tcre~t)| 2083.505 42.75414 2001.372 2169.01 var(e.tgre~t)| 8724.877 176.5499 8385.617 9077.862 var(e.Health)| 340.79 11.08864 319.7351 363.2313 var(age)| 293.2524 5.835941 282.0344 304.9166 var(female)| .2497233 .0049697 .2401704 .2596562 var(ageXfe~e)| 729.2019 14.51166 701.3071 758.2062 -------------+---------------------------------------------------------------- cov(age,| female)| .0686977 .1204256 0.57 0.568 -.167332 .3047275 cov(age,| ageXfemale)| 155.358 6.864694 22.63 0.000 141.9034 168.8125 cov(female,| ageXfemale)| 12.00529 .2541636 47.23 0.000 11.50714 12.50344 ------------------------------------------------------------------------------ LR test of model vs. saturated: chi2(11) = 1542.14 Prob > chi2 = 0.0000

The first matrix we will create is the \(\mathbf{A}\) matrix. The \(\mathbf{A}\) matrix is an asymmetric matrix that stores the single-headed path coefficients in the model. The rows and columns of our matrix represent each of the observed and latent variables. The cells represent paths from the column variable to the row variable. I will start by creating a matrix of 0s of the size that I need using **J(***r***,***c***,***z***)**, an \(r \times c\) matrix containing elements *z*. Then I will substitute blocks of cells with the paths from the model into the \(\mathbf{A}\) matrix. We only need to specify the upper left cell of the cells we are substituting. Finally, I will add row and column names to this matrix to make it easier to interpret.

. matrix A = J(8,8,0) . matrix A[1,8] = (1.00\0.44\0.66\1.12) . matrix A[8,5] = (0.44,-20.89,0.34) . matrix rownames A = bpsystol bpdiast tcresult tgresult age female > ageXfemale Health . matrix colnames A = bpsystol bpdiast tcresult tgresult age female > ageXfemale Health . matrix list A A[8,8] bpsystol bpdiast tcresult tgresult age bpsystol 0 0 0 0 0 bpdiast 0 0 0 0 0 tcresult 0 0 0 0 0 tgresult 0 0 0 0 0 age 0 0 0 0 0 female 0 0 0 0 0 ageXfemale 0 0 0 0 0 Health 0 0 0 0 .44 female ageXfemale Health bpsystol 0 0 1 bpdiast 0 0 .44 tcresult 0 0 .66 tgresult 0 0 1.12 age 0 0 0 female 0 0 0 ageXfemale 0 0 0 Health -20.89 .34 0

The \(\mathbf{S}\) matrix is a symmetric matrix that stores double-headed path coefficients, variances, and covariances of our model. I start by creating a diagnal matrix of all the estimated variances from the model, and then I add the covariances of the exogenous variables.

. matrix S = diag((67.63,80.19,2083.51,8724.88,0,0,0,340.79)) . matrix S[5,5] = (293.25,0.07,155.36\0.07,0.25,12.01\155.36,12.01,729.20) . matrix rownames S = bpsystol bpdiast tcresult tgresult age female > ageXfemale Health . matrix colnames S = bpsystol bpdiast tcresult tgresult age female > ageXfemale Health . matrix list S symmetric S[8,8] bpsystol bpdiast tcresult tgresult age bpsystol 67.63 bpdiast 0 80.19 tcresult 0 0 2083.51 tgresult 0 0 0 8724.88 age 0 0 0 0 293.25 female 0 0 0 0 .07 ageXfemale 0 0 0 0 155.36 Health 0 0 0 0 0 female ageXfemale Health female .25 ageXfemale 12.01 729.2 Health 0 0 340.79

The \(\mathbf{F}\) matrix is a rectangular matrix to distinguish observed and latent variables. The rows represent each observed variable, and the columns represent each observed and latent variable. Each row gets a `1′ in the cell of that variable’s corresponding column, creating a diagonal submatrix. A diagonal matrix of 1s is also called the identity matrix. We can specify an identity matrix with **I(***n***)** where *n* is the number of rows/columns. I add a column vector of 0s with **J()**.

. matrix F = I(7),J(7,1,0) . matrix rownames F = bpsystol bpdiast tcresult tgresult age female > ageXfemale . matrix colnames F = bpsystol bpdiast tcresult tgresult age female > ageXfemale Health . matrix list F F[7,8] bpsystol bpdiast tcresult tgresult age bpsystol 1 0 0 0 0 bpdiast 0 1 0 0 0 tcresult 0 0 1 0 0 tgresult 0 0 0 1 0 age 0 0 0 0 1 female 0 0 0 0 0 ageXfemale 0 0 0 0 0 female ageXfemale Health bpsystol 0 0 0 bpdiast 0 0 0 tcresult 0 0 0 tgresult 0 0 0 age 0 0 0 female 1 0 0 ageXfemale 0 1 0

Finally, if your model involves mean structure, an \(\mathbf{M}\) matrix must also be created to obtain the model-implied means. This is just a column vector of estimated means or intercepts for the observed and latent variables in the model. Unless explicitly estimated, you can assume all the latent variable means are 0.

. matrix M = (111.25\72.88\203.61\123.03\47.94\.52\24.84\0) . matrix rownames M = bpsystol bpdiast tcresult tgresult age female > ageXfemale Health . matrix list M M[8,1] c1 bpsystol 111.25 bpdiast 72.88 tcresult 203.61 tgresult 123.03 age 47.94 female .52 ageXfemale 24.84 Health 0

Now we can use matrix algebra to derive the model-implied means and covariance. The dimension of **I()** will need to be changed based on the number of observed and latent variables in the model. Here it’s 8 to include the seven observed variables and the one latent variable.

. *model-implied means . matrix mu = F*inv(I(8)-A)*M . matrix list mu mu[7,1] c1 bpsystol 129.9264 bpdiast 81.097616 tcresult 215.93642 tgresult 143.94757 age 47.94 female .52 ageXfemale 24.84 . *model-implied covariance . matrix Sigma = F*inv(I(8)-A)*S*inv(I(8)-A)'*F' . matrix list Sigma symmetric Sigma[7,7] bpsystol bpdiast tcresult tgresult age bpsystol 533.17918 bpdiast 204.84164 170.32032 tcresult 307.26246 135.19548 2286.3032 tgresult 521.41508 229.42264 344.13395 9308.8649 age 180.3901 79.371644 119.05747 202.03691 293.25 female -1.1083 -.487652 -.731478 -1.241296 .07 ageXfemale 65.3975 28.7749 43.16235 73.2452 155.36 female ageXfemale female .25 ageXfemale 12.01 729.2

We see the covariance matrix implied by the model with the supplied model parameters.