In jsakaluk/dySEM: Dyadic Structural Equation Modeling

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Fair Use of this Tutorial

This tutorial is a supplemental material from the following article:

Sakaluk, J. K., & Camanto, O. J. (2025). Dyadic Data Analysis via Structural Equation Modeling with Latent Variables: A Tutorial with the dySEM package for R.

This article is intended to serve as the primary citation of record for the dySEM package's functionality for the techniques described in this tutorial. If this tutorial has informed your modeling strategy (including but not limited to your use of dySEM), please cite this article.

The citation of research software aiding in analyses--like the use of dySEM--is a required practice, according to the Journal Article Reporting Standards (JARS) for Quantitative Research in Psychology (Appelbaum et al., 2018).

Furthermore, citations remain essential for our development team to demonstrate the impact of our work to our local institutions and our funding sources, and with your support, we will be able to continue justifying our time and efforts spent improving dySEM and expanding its reach.

Overview

The CDFM is a "uni-construct" dyadic SEM (i.e., used to represent dyadic data about one--and only one--construct, like "relationship satisfaction").

It contains:

• two latent variables/factors, onto which
• each partner's observed variables discriminantly load (i.e., one partner's observed variables onto their factor, and the other partner's observed variables onto their factor)

It also features two varieties of covariances (or correlations, depending on scale-setting/output standardization):

• one between the two latent variables (effectively, the latent intraclass correlation coefficient), and
• several between the residual variances of the same observed variables across each partner (e.g., between Item 1 for Partner A and Partner B; another between Item 2 for Partner A and Partner B, etc.,).

Owing to it serving as one half of the measurement model of a latent-APIM--and the APIM being the clear default structural model in relationship science (Kenny, 2018)--the CDFM could (taciltly) be considered relationship science's "default"--for better and for worse--uni-construct dyadic SEM.

Packages, Data, and dySEM Overview

This exemplar makes use of the tidverse meta-package, and the gt, dySEM, and lavaan (Rosseel, 2012) packages.

library(dplyr) #for data management
library(gt) #for reproducible tabling
library(dySEM) #for dyadic SEM scripting and outputting
library(lavaan) #for fitting dyadic SEMs

For the exemplar, we will use the built-in commitmentM dataset from dySEM, which includes ratings of satisfaction and committment from 282 mixed-sex dyads. These data were collected using the "global" version of items for relationship satisfaction from the Rusbult et al. (1998) Investment Model Scale. More information about this data set can be found in Sakaluk et al. (2021).

This dataset already possesses many of the desirable features of a dataset for use with dySEM, as it is in wide format, and the variable names have a repetitive predictable structure (to learn more, see our tutorial on variable name structure). Specificially:

• The satisfaction items follow a "sip" pattern (Stem, Item, Partner) with a stem of "sat.g", followed immediately by a number for the item (i.e., no delimeter between these elements), followed by a delimiting "_", and then either the character "f" or "m to indicate to which partner the item refers. For example, the first satisfaction item for partner 1 is "sat.g1_f", and the first satisfaction item for partner 2 is "sat.g1_m".

• Likewise, the commitment items follow a similar "sip" pattern, albeit with a different stem, "com".

For this exemplar of dyadic invariance testing, we will focus on the satisfaction items, and assign them into a stand-alone demonstration data frame called my_dat:

my_dat <- commitmentM

#show satisfaction item previews
my_dat |> 
  select(starts_with("sat.g")) |> 
  as_tibble() 

As described in the paper, and the Get Started tutorial for dySEM , the steps we must take to carry out dyadic invariance testing with dySEM are:

1. Scrape the relevant items from the data set, to be included in each CDFM
2. Script the lavaan syntax for the dyadic configural, loading, interecept, and residual invariance models
3. Fit the invariance models with lavaan
4. generate reproducible output of our analyses, and
5. consider the use of optional functionality for supplementary indexes and tests that dySEM can compute (we will demonstrate a couple)

1. Scrape the Variables

The scrapeVarCross from dySEM is the appropriate function to scrape the variable information needed for automating the scripting of latent dyadic models of cross-sectional data. In this instance, we are only scraping information about indicators about one (pair of) dyadic latent variable(s): men and women's latent relationship satisfaction.

As described before, the satisfaction items follow a "sip" pattern (Stem, Item, Partner) with a stem of "sat.g", followed immediately by a number for the item (i.e., no delimeter between these elements), followed by a delimiting "_", and then either the character "f" or "m to indicate to which partner the item refers.

We therefore refer to our data frame ("my_dat") for the dat argument, and the pattern, stem, delimeters (if any), and distinguishing characters for the x_order, x_stem, x_delim1, x_delim2, distinguish_1, and distinguish_2 arguments, respectively. In this instance, our stem is "sat.g"; there is no first delimeter, and an "_" for the second delimeter; and the distinguishing characters are "f" and "m".

We assign the output of this function (a list of variable names an information) to an object we are (arbitrarily) calling "sat_dvn"; we often refer to this list as a "dyad variable names" list, or "dvn" for short.

sat_dvn <- scrapeVarCross(dat = my_dat, x_order = "sip", x_stem = "sat.g", x_delim1 = "", x_delim2="_", distinguish_1="f", distinguish_2="m")

sat_dvn

As can be seen, the information contained in the dvn is rather mundane: the indicator names for the "f" partners (e.g., "sat.g1_f"), the indicator names for the "m" partners (e.g., "sat.g5_m"), the number of indicators per partner (in this case, 5), the distinguishing characters ("f" and "m"), and the total number of indicators (10) in the entire set. Yet this information is all that is needed for dySEM to dramatically simplify the process of scripting dyadic SEM models, including dyadic invariance models.

2. Script the Model(s)

Next, we will script the dyadic configural, loading, intercept, and residual invariance models. We will use the scriptCor function from dySEM to do this. This function requires the dvn object we created in the prior step, an arbitrary name for the latent variable for lavaan, and the constr_dy_meas and constr_dy_struct arguments to indicate which measurement (i.e., "meas") and structural (i.e., "struct") parameters we want to constraint (i.e., "constr") across dyad members (i.e., "dy"). The scaleset argument automatically defaults to the fixed-factor scale-setting and identification method, prefered for invariance testing . However, we encourage users who are just beginning to include this argument explicitly, so they (and others) can readily determine their model specification approach.

Below, we script the sequence of four models to be fit and compared; we forego the use of the writeTo and fileName arguments, which are used to save the .txt of scripts to a file. Finally, we use an intuitive naming structure to keep track of which model script resides in which R object.

sat.config.script <- scriptCor(sat_dvn, lvname = "Sat", scaleset = "FF", constr_dy_meas = "none", constr_dy_struct = "none")

sat.load.script <- scriptCor(sat_dvn, lvname = "Sat", scaleset = "FF", constr_dy_meas = c("loadings"), constr_dy_struct = "none")

sat.int.script <- scriptCor(sat_dvn, lvname = "Sat", scaleset = "FF", constr_dy_meas = c("loadings", "intercepts"), constr_dy_struct = "none")

sat.resid.script <- scriptCor(sat_dvn, lvname = "Sat", scaleset = "FF", constr_dy_meas = c("loadings", "intercepts", "residuals"), constr_dy_struct = "none")

In their stored form, these scripts are not particularly readable, but as we shall soon see, they are immediately passable to lavaan:

sat.resid.script

However, when exported (via the writeTo and fileName arguments) or concatenated, they can be presented in a much more readable format that mirrors what users expect from a manually written lavaan script:

cat(sat.resid.script)

3. Fit the Model(s)

All of these models can now be fit with whichever lavaan wrapper the user prefers, and with their chosen analytic options (e.g., estimator, missing data treatment, etc)

sat.config.mod <- cfa(sat.config.script, data = my_dat)

sat.load.mod <- cfa(sat.load.script, data = my_dat)

sat.int.mod <- cfa(sat.int.script, data = my_dat)

sat.resid.mod <- cfa(sat.resid.script, data = my_dat)

4. Output and Interpret the Model(s)

Users could then inspect these models and compare these models using their preferred strategy(ies), such as likelihood ratio test, and/or changes in certain fit measures and/or information criteria. Given the analytic researcher degrees of freedom and traditions of preference herein, we leave it to researchers to apply (and preferably also preregister) their preferred approach.

For example, researchers could use the built in anova() function as below:

anova(sat.config.mod, sat.load.mod, sat.int.mod, sat.resid.mod)

Conversely, they could instead consider the use of outputInvarCompTab(), a function we have written to wrap some of the same comparative output in a more reproducible-reporting-friendly format:

mods <- list(sat.config.mod, sat.load.mod, sat.int.mod, sat.resid.mod)

outputInvarCompTab(mods)

Without additional arguemntation, outputInvarCompTab() (and other outputters in dySEM) return data frames; these are the most flexible output types to supply to other tabling and/or visualization functionality. However, users looking for a more expedient way to "clean up" reporting from outputters can either make use of the gtTab argument; TRUE will return a basic gt table which can then be further modified using additional gt layers. This is a good option for reproducible reporting document formats (like .Rmd and .Qmd):

outputInvarCompTab(mods, gtTab = TRUE)

Conversely, if users prefer to export the output to a file, they can set the additional writeTo argument to a directory path (a path of "." will write to the user's current working directory, like that associated with an .Rproj file), and the fileName argument to a desired base name for the output file. This will save the output as an .rtf file in that directory, which can then be opened in Microsoft Word or other word processing software:

outputInvarCompTab(mods, gtTab = TRUE, writeTo = ".", fileName = "sat_invar_comp")

comp <- anova(sat.config.mod, sat.load.mod, sat.int.mod, sat.resid.mod)
comp <- janitor::clean_names(comp)

chisq_diff <- round(comp$chisq_diff, 2)
p <- round(comp$pr_chisq, 2)

Were one strictly to follow only the likelihood ratio test statistic as a guide, they might feel comfortable claiming they have reasonable support for most forms of dyadic invariance, with the exception of significant residual noninvariances, $\chi^2$ (r comp$df_diff[4]) = r chisq_diff[4], p = r p[4].

Users can then request that dySEM provide them with reproducible output via either the outputParamTab() (for tabular output) or the outputParamFig() (for a path diagram, back-ended by semPlot). For example, if one wanted a table of measurement model output from the dyadic intercept invariance model they could use:

outputParamTab(sat_dvn, model = "cfa", fit = sat.int.mod, 
               tabletype = "measurement") |> 
  gt()

The model and tabletype arguments are the main providers of optionality to be aware of, as outputParamTab is also used to provide tabular output from other model types (e.g., APIMs), and users may be most interested in the output from one portion of the model (e.g., the measurement model) versus another (e.g., the structural model).

Path diagrams, meanwhile, are created with the outputParamFig() function, which is used in a similar manner to outputParamTab() (and includes the same writeTo and fileName arguments if you wish to export a .png). Otherwise, the primary argument to be aware of is figtype, which allows users to specify what kind of parameter estimates they wish to be depicted in the path diagram.

outputParamFig(fit = sat.int.mod, figtype = "standardized")

params <- parameterEstimates(sat.int.mod, standardized = TRUE)

params <- params |> 
  filter(op == "~~" & lhs == "Satf" & rhs == "Satm") |> 
  select(std.all, z, pvalue, ci.lower, ci.upper) |> 
  mutate_if(is.numeric, round, 2)

Thus, in little more than a dozen lines of code, researchers can script, fit, compare, and output a series of dyadic invariance comparisons using the Correlated Dyadic Factors Model. From this, researchers would learn that the satisfaction measure demonstrates reasonable configural, loading, and intercept invariance between partners, and that the latent ICC between partner's satisfaction levels is r params$std.all, z = r params$z, p r if(params$p < .001){paste("<.001")} (95% CI: r params$ci.lower, r params$ci.upper).

5. Consider Optional Indexes and Output

A remaining question would be: if our chosen model fails to demonstrate residual invariance, which indicator(s) are responsible for introducing the noninvariance? The outputConstraintTab() function from dySEM, while not strictly necessary, offers supplemental functionality to help answer this question.

Users simply enter the lavaan model responsible for the statistically significant degradation of model fit, and are provided with a table of each dyadic invariance constraint, by indicator, and its Langrange multiplier test statistic and corresponding test:

outputConstraintTab(sat.resid.mod) |> 
  gt()

In this particular instance, significant noninvariance was only introduced by forcing the constraint on the residual variances for the 1st and 5th items

jsakaluk/dySEM documentation built on June 11, 2025, 2:16 p.m.