Purpose of simstandard

In the figure below, you can see a hypothetical structural model with its standardized loadings and path coefficients.

Suppose you need to simulate multivariate normal data based on this model, but you do not know the error variances and the latent disturbance variances needed to make your model produce standardized data. It is often difficult to find such values algebraically and instead they must be found iteratively.

The simstandard package finds the standardized variances and creates standardized multivariate normal data using lavaan syntax. It can also create latent variable scores, error terms, disturbance terms, estimated factor scores, and equally weighted composite scores for each latent variable.

Generate Model-Based Multivariate Data

A1 A2 A3 B1 B2 B3 A B e_A1 e_A2 e_A3 e_B1 e_B2 e_B3 d_B
1.51 0.88 1.03 2.38 3.58 2.87 1.31 2.79 0.60 -0.17 -0.15 0.04 1.35 0.36 2.00
1.57 0.59 0.51 0.33 0.51 1.12 1.05 0.68 0.83 -0.26 -0.43 -0.46 -0.03 0.51 0.05
-0.50 0.59 0.76 0.36 -0.31 0.70 0.04 0.17 -0.52 0.56 0.73 0.23 -0.45 0.55 0.15
-0.27 -0.94 -0.76 -0.29 0.22 -0.68 -0.86 -0.46 0.33 -0.26 0.01 0.29 0.58 -0.27 0.06
-0.02 -0.81 0.45 -0.23 0.72 0.20 1.04 -0.27 -0.75 -1.65 -0.49 -0.35 0.94 0.44 -0.89
-1.44 -0.19 -0.79 -1.13 -0.44 -0.82 -0.97 -0.67 -0.77 0.58 0.07 -0.37 0.09 -0.22 -0.09

Let’s make a function to display correlations and covariance matrices:

Because the data are standardized, the covariance matrix of the observed and latent variables should be nearly identical to a correlation matrix. The error and disturbance terms are not standardized.

cov(d) %>% 
  ggcor

To return only the observed variables

A1 A2 A3 B1 B2 B3
1.61 1.15 2.06 0.90 -0.05 -0.08
1.46 1.50 0.68 0.12 -0.48 0.34
-0.89 -1.19 -1.22 -0.52 0.40 -0.73
0.13 -1.55 -0.79 -1.00 -0.17 -1.18
-0.93 -1.64 -0.96 -1.23 -1.51 -1.01
-0.43 -0.24 -0.10 -0.87 -0.78 -1.21

Comparison with lavaan::simulateData

I love the lavaan package. However, one aspect of one function in lavaan is not quite right yet. lavaan’s simulateData function is known to generate non-standardized data, even when the standardized parameter is set to TRUE. See how it creates variables B1, B2, and B3 with variances much higher than 1. Furthermore, it only produces observed variables.

Inspecting model matrices

You can inspect the matrices that simstandard uses to create the data by calling simstandardized_matrices.

The A matrix contains all the asymmetric path coefficients (i.e., the loadings and the structural coefficients). These coefficients are specified in the lavaan model syntax.

The S matrix contains all the symmetric path coefficients (i.e., the variances and correlations of the observed and latent variables). For endogenous variables, the variances and correlations refer to the variance and correlations of the variable’s associated error or disturbance term. In this case, A is the only endogenous variable, and thus its variance on the diagonal of the S matrix is 1.

Thus, we can use these results to insert the missing values from the path diagram at the beginning of this tutorial

Estimated Factor Scores

If you want to estimate factor scores using the regression method (i.e., Thurstone’s method), set factor_scores to TRUE. All scores ending in *_FS* are factor score estimates.

A1 A2 A3 A e_A1 e_A2 e_A3 A_FS e_A1_FS e_A2_FS e_A3_FS
0.54 -0.76 0.54 -0.81 1.26 -0.11 1.10 0.20 0.82 -1.53 0.56
-3.67 -3.10 -2.72 -3.67 -0.37 -0.17 -0.15 -3.50 -1.19 -0.50 -0.38
0.13 -1.01 -1.10 -0.26 0.37 -0.79 -0.92 -0.39 1.11 -1.16 -1.16
1.87 1.10 1.35 0.79 1.16 0.46 0.79 1.65 0.90 -0.36 0.27
-0.72 -0.21 -0.39 -0.33 -0.42 0.06 -0.16 -0.55 -0.52 0.39 0.00
0.44 0.97 0.91 0.82 -0.31 0.31 0.33 0.68 -0.41 0.71 0.60

Adding factor scores to new data

Suppose you have some new data and which to add estimated factor scores to it. The add_factor_scores function will take your data and return your data with the estimated factors added to it.

A1 A2 A3 A
2 2.5 1.3 2.1
-1 -1.5 -2.1 -1.4

Composite Scores

If you want to calculate equally-weighted composite scores based on the indicators of each latent variable, set `composites = TRUE’.

A1 A2 A3 A e_A1 e_A2 e_A3 A_Composite
0.03 -0.30 -0.94 -0.55 0.53 0.13 -0.55 -0.46
-0.11 0.04 1.05 0.01 -0.12 0.03 1.04 0.38
-2.43 -1.75 -1.29 -2.01 -0.62 -0.15 0.12 -2.09
-0.88 -1.26 -1.76 -0.70 -0.25 -0.70 -1.28 -1.49
-0.15 -0.13 -0.04 -0.21 0.04 0.04 0.10 -0.12
0.93 0.94 1.71 1.41 -0.34 -0.18 0.72 1.37

Return lavaan syntax with all parameters set free

Suppose that we want to verify that the data generated by the sim_standardized function is correct. We will need an analogous model, but with all the fixed parameters set free. We could manually remove the fixed parameter values, but with large models the process is tedious and introduces a risk of error. The fixed2free function painlessly removes the fixed parameters values from the model.

Now let’s use lavaan to see if the observed data in d conform to the model in m_free.

# Set the random number generator for reproducible results
set.seed(12)
# Generate data based on model m
d <- sim_standardized(
  m,
  n = 100000,
  latent = FALSE,
  errors = FALSE)

# Evaluate the fit of model m_free on data d
library(lavaan)
lav_results <- sem(
  model = m_free, 
  data = d)

# Display summary of model
summary(
  lav_results, 
  standardized = TRUE, 
  fit.measures = TRUE)
#> lavaan 0.6-3 ended normally after 27 iterations
#> 
#>   Optimization method                           NLMINB
#>   Number of free parameters                         14
#> 
#>   Number of observations                        100000
#> 
#>   Estimator                                         ML
#>   Model Fit Test Statistic                       7.493
#>   Degrees of freedom                                 7
#>   P-value (Chi-square)                           0.379
#> 
#> Model test baseline model:
#> 
#>   Minimum Function Test Statistic           371352.125
#>   Degrees of freedom                                15
#>   P-value                                        0.000
#> 
#> User model versus baseline model:
#> 
#>   Comparative Fit Index (CFI)                    1.000
#>   Tucker-Lewis Index (TLI)                       1.000
#> 
#> Loglikelihood and Information Criteria:
#> 
#>   Loglikelihood user model (H0)             -666610.982
#>   Loglikelihood unrestricted model (H1)     -666607.236
#> 
#>   Number of free parameters                         14
#>   Akaike (AIC)                              1333249.965
#>   Bayesian (BIC)                            1333383.146
#>   Sample-size adjusted Bayesian (BIC)       1333338.653
#> 
#> Root Mean Square Error of Approximation:
#> 
#>   RMSEA                                          0.001
#>   90 Percent Confidence Interval          0.000  0.004
#>   P-value RMSEA <= 0.05                          1.000
#> 
#> Standardized Root Mean Square Residual:
#> 
#>   SRMR                                           0.001
#> 
#> Parameter Estimates:
#> 
#>   Information                                 Expected
#>   Information saturated (h1) model          Structured
#>   Standard Errors                             Standard
#> 
#> Latent Variables:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>   A =~                                                                  
#>     A1                1.000                               0.703    0.703
#>     A2                1.142    0.005  231.116    0.000    0.803    0.800
#>     A3                1.284    0.005  247.676    0.000    0.903    0.901
#>     B1                0.427    0.004  114.547    0.000    0.300    0.300
#>   B =~                                                                  
#>     B1                1.000                               0.701    0.700
#>     B2                1.139    0.005  238.139    0.000    0.798    0.798
#>     B3                1.288    0.005  247.366    0.000    0.902    0.901
#> 
#> Regressions:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>   B ~                                                                   
#>     A                 0.597    0.004  141.196    0.000    0.599    0.599
#> 
#> Variances:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>    .A1                0.507    0.003  194.579    0.000    0.507    0.506
#>    .A2                0.362    0.002  165.904    0.000    0.362    0.359
#>    .A3                0.188    0.002   97.276    0.000    0.188    0.188
#>    .B1                0.170    0.001  130.982    0.000    0.170    0.169
#>    .B2                0.362    0.002  178.775    0.000    0.362    0.362
#>    .B3                0.189    0.002  106.863    0.000    0.189    0.189
#>     A                 0.495    0.004  121.519    0.000    1.000    1.000
#>    .B                 0.315    0.003  124.248    0.000    0.641    0.641

# Extract RAM paths
RAM <- lav2ram(lav_results)

# Display asymmetric paths (i.e., single-headed arrows for 
# loadings and structure coefficients)
RAM$A %>% ggcor()

As can be seen, all the fit measures indicate a near-perfect fit, and the parameter estimates are within rounding error of the fixed parameters in model m.

Return lavaan syntax for a model with standardized variances specified

Although the simstandardized function will generate data for you, you might want to use a function from a different package instead, such as lavaan::simulateData or simsem::sim. In this case, you can use the model_complete function to output the lavaan syntax for a standardized model with all standardized variances specified.

Return lavaan syntax from matrices

Suppose that a research article provides model coefficients in a table. We could spend time creating lavaan syntax by hand, but such work can be tedious. The matrix2lavaan function can help save time when the models are already specified in matrix form.

The measurement model

The measurement model can be specified with a matrix in which the column names are latent variables and the row names are indicator variables.

Here we have three latent variables, Vocabulary, Working Memory Capacity, and Reading, each defined by three indicator variables.

The structural model

The structural model can be specified with a matrix in which the predictors are the column names and the criterion variables are the row names.

Here we have Vocabulary and Working Memory Capacity predicting Reading Scores.

m_struct <- matrix(
  c(0.4,0.3), 
    ncol = 2,
    dimnames = list(
      "Reading",
      c("Vocabulary", "WorkingMemory"))) 

This could have been a 3 by 3 matrix with zeroes (which are ignored).

m_struct <- matrix(c(
    0,   0,   0,  # Vocabulary
    0,   0,   0,  # WorkingMemory
    0.4, 0.3, 0), # Reading
    nrow = 3, 
    byrow = TRUE) 
rownames(m_struct) <- c("Vocabulary", "WorkingMemory", "Reading")
colnames(m_struct) <- c("Vocabulary", "WorkingMemory", "Reading")

Covariances

The variances and covariances must be specified as a symmetric matrix, though variables can be omitted.

Here we specify that the latent variables Vocabulary and Working Memory Capacity are correlated.

m_cov <- matrix(c(
    1,   0.5, 
    0.5, 1), 
    nrow = 2,
    dimnames = list(
      c("Vocabulary", "WorkingMemory"),
      c("Vocabulary", "WorkingMemory"))) 

Using the matrix2lavaan function

The matrix2lavaan function takes arguments for the measurement model, structural model, and covariances. Any of the three matrices can be omitted.