Metric Invariance in Exploratory Graph Analysis via Permutation Testing

Measurement Invariance in Traditional Psychometrics

Exploratory Graph Analysis

Network psychometric methods are an alternative to latent variable modeling. Networks represent variables as nodes (circles) and their relationships (e.g., partial correlations) as edges (lines). Because the relationships between variables are not known a priori, they must be estimated. There are many methods to estimate a network with the graphical least absolute shrinkage and selection operator (GLASSO; Epskamp & Fried, 2018; Friedman et al., 2008) being one of the most common.

A key feature of network models is that each node is (usually) not connected to all other nodes (known as sparsity). Often, some nodes are more densely connected to each other relative to other nodes in the network. These sets of connected nodes are often referred to as communities, which are consistent with latent factors when data are generated from a factor model (Golino & Epskamp, 2017). Community detection algorithms are a common, data-driven way to identify communities in networks (Fortunato, 2010). The combination of the GLASSO with the Walktrap community detection algorithm (Pons & Latapy, 2006) has been labeled as Exploratory Graph Analysis (Golino et al., 2020; EGA; Golino & Epskamp, 2017) in the network psychometrics literature.

Across the broader field of network psychometrics, several methods have been developed to identify differences in network structure (Van Borkulo et al., 2022; Williams et al., 2020) and sub-groups (Danaher et al., 2014; Haslbeck & Bork, 2022; Jones et al., 2020). Although these methods aim to identify differences between networks, they all tend to treat the networks as unidimensional—that is, these methods do not account for the community structure of the network. Therefore, unless the construct is assumed to be unidimensional, the detected differences are unlikely to parallel traditional measurement invariance procedures. Establishing community structure and a metric consistent with factor loadings is key for developing such a comparison method.

Recent work has demonstrated that a node’s strength or absolute sum of a node’s connections to other nodes is related to confirmatory factor analysis (CFA) loadings (Hallquist et al., 2021). Hallquist et al. (2021) found, however, that the strength of a node is comprised of both dominant and cross loadings. To circumvent this issue, Christensen and Golino (2021b) proposed a measure called $networkloadings$ that splits a node’s strength based on the dimensions identified by EGA. In their simulation, they found that this measure was consistent with factor loadings when data were generated by a factor model. The development of network loadings opened the door for broader measurement evaluation within network psychometrics such as item selection, weighted between-person scores, and hierarchical dimensionality assessment (Christensen & Golino, 2021b; Jiménez et al., 2023).

The goal of this study is to leverage these network loadings to establish a method within the network psychometric framework to test for measurement invariance. The extent of measurement invariance within the network psychometrics framework only includes configural and metric invariance because latent variables are not estimated using networks. Consequently, intercepts and residual variances are not feasible because there are no latent factors created. Therefore, measurement invariance using network psychometrics, like network loadings, is a heuristic for configural and metric invariance in latent factors rather than a direct equivalent.

Present Research

SEM Procedure

To test metric and partial metric invariance using SEM, we estimated two models: a configural, unconstrained model, see Equation (4), and a model with loadings constrained to equality across k populations, constrained metric model, see Equation (6). In order to directly compare to our proposed method, we tested for partial metric invariance. Testing for partial metric invariance was conducted using three methods: Free, Fixed, or Wald. The Free method follows the method proposed by Raykov et al. (2013). Using the {semTools} package (Version 0.5-6; Jorgensen et al., 2022), the Fixed and Wald methods are run simultaneously with Free. Since it is a common software method for researchers to use in practice, we evaluated the results of all three approaches. In all methods, an original model was chosen to be either the constrained or unconstrained model. After, the loadings were either fixed or freed iteratively to create a new model which was then compared to the original model. Using these methods circumvented a common problem in many approaches to invariance testing: we did not exclude the testing of any variables by fixing the loading of one variable per factor to 1. In this way, we could make direct comparisons between the SEM and proposed methods.

The Free method uses the constrained model as the original model. Iteratively, each variable j is freed in the matrix Λ to create J models. Each model is then compared to the original model using a likelihood ratio test and an assessment of CFI for a total of J tests. The Fixed method uses the unconstrained model as the original model. Iteratively, each variable j is constrained to be equal across populations k to create J models. Each model is then compared to the original model using a likelihood ratio test and an assessment of CFI for a total of J tests. The Wald method is similar to Free. It uses the constrained model as the original model, but rather than iteratively freeing each variable j and conducting likelihood ratio tests, it uses a multivariate Wald test. Nonetheless, multiple hypotheses are being tested. These methods do not adjust for Type I error rate and so it is often necessary to apply a multiple comparison test (Raykov et al., 2013).

Multiple Comparison Problem

Within both partial invariance frameworks (EGA and SEM), multiple hypotheses are tested which can artificially inflate the Type I error rate. To adjust for this inflation, a multiple comparison procedure (MCP) can be applied (Raykov et al., 2013; Steinberg, 2001). To select which MCP to apply, it’s important to consider consequences in the trade off of identifying (non)invariant items. Most MCPs focus on controlling the Family Wise Error Rate (FWER). FWER attempts to avoid making any Type I error and is inclined toward the notion that a Type I error is a serious issue.

In the context of partial invariance, a Type I error would suggest that a variable is noninvariant when it is truly invariant. In most research contexts, the cost of falsely identifying an item as invariant is greater than falsely identifying an item as noninvariant, particularly if the construct will be used to compare across groups (Shi et al., 2019). Therefore, FWER may suggest that more variables are invariant than there truly are, potentially leading to more costly consequences than using an uncorrected p-value. An alternative and less conservative MCP is the Benjamini-Hochberg procedure (BH-procedure; Benjamini & Hochberg, 1995) which controls the False Discovery Rate (FDR). FDR takes a more balanced approach to the multiple comparison problem by using the expected number of falsely rejected null hypotheses if any null hypotheses are rejected to determine its correction. Formally, FDR works to control ϕ:

where V is the number of falsely rejected null hypotheses and S is the number of correctly rejected null hypotheses out of the set of all hypotheses tested. The BH-procedure provides adequate control over false positives while showing marked improvements in power above and beyond traditional MCP methods (e.g., Tukey, Bonferroni, Scheffe; Benjamini & Hochberg, 1995).

The BH-procedure works by sorting individual p-values in ascending order and assigning them a rank. The adjusted p-value is computed using:

where mr represents the total number of p-values and r represents to rank for the corresponding p-value j.

Raykov et al. (2013) first proposed the use of the BH-procedure to test partial invariance due to its more liberal approach of lowering the risk of false positive noninvariant variables relative to FWER’s focus on lowering the risk of false positive noninvariant variables entirely. Given that the consequences are usually more dire when not correctly detecting noninvariant variables, the BH-procedure was preferred over FWER as our MCP. All p-values calculated using the BH-procedure, hereafter, are referred to as corrected p-values.

Data Generation

Data was generated following a common factor model, as was done by Golino et al. (2020). We begin by computing a population correlation matrix for each group, ${\mathbf{R}}_{{\mathbf{R}}_{\mathbf{G}}}$, with communalities in the diagonal,

where ${\mathbf{R}}_{{\mathbf{R}}_{\mathbf{G}}}$ is the reproduced population correlation matrix for each group G, ${\mathrm{\Lambda }}_{\mathrm{G}}$ is a k × m factor loading matrix for k variables and m factors for each group G, and Phi (Φ) is the structure matrix of the latent variables (i.e., a m × m matrix of correlations among factors). The population does not contain any correlated residuals and therefore no minor factors.

Then, by inserting unities in the diagonal of ${\mathrm{R}}_{{\mathrm{R}}_{\mathrm{G}}}$ it becomes a full rank matrix and is now population correlation matrix ${\mathrm{R}}_{{\mathrm{P}}_{\mathrm{G}}}$. Each group in G is assigned a ${\mathrm{R}}_{{\mathrm{P}}_{\mathrm{G}}}$ matrix. A Cholesky decomposition is performed on each ${\mathrm{R}}_{{\mathrm{P}}_{\mathrm{G}}}$:

If any ${\mathbf{R}}_{{\mathbf{P}}_{\mathbf{G}}}$ is not semi-positive definite or an item’s communality is greater than 0.90, then a new ${\mathbf{R}}_{{\mathbf{P}}_{\mathbf{G}}}$ matrix is constructed. From this, the sample data matrix (continuous variables) can be computed as:

where ${\mathbf{Z}}_{\mathbf{G}}$ is a matrix of random standard normal deviates with rows equal to the sample size and columns equal to the number of variables.

Design

The overall design of the simulation study followed closely that of Kim and Yoon (2011) with a few modifications. A two-factor model was simulated, each factor containing six variables, similar to Kim and Yoon (2011) and Yoon and Millsap (2007). Typically, simulation studies investigating invariance methods use unidimensional models; however, we decided to simulate two factors. This approach allowed us to manipulate interfactor correlation and investigate whether it impacted the power of the proposed method. Only one variable in one factor was simulated to have unequal dominant loadings across group. Since our main goal was to assess each method’s ability to identify noninvariant items correctly, having only one noninvariant item allows for a direct estimate of the true positive rate. Additionally, it allows us to compare the ability of each method to detect invariant items within factors both with and without noninvariant items.

For simplicity, we only simulated two groups. Factor loadings were set to be the same across factors for each respective variable (0.80, 0.70, 0.60, 0.80, 0.70, 0.60). Keeping high, static factor loadings allowed us to make sure configural invariance was not negatively impacted, particularly for data conditions with a high difference in loadings and/or a high interfactor correlation. Similar to Golino and Epskamp (2017), the correlation between factors was set to be low (0.30), medium (0.50), or high (0.70).

The loading of Variable 5 in Factor 1 (0.70) was decreased in G₁ by either 0.20 (small difference) or 0.40 (large difference) as was done in Kim and Yoon (2011). Static factor loadings ensured that the magnitude of loading differences will have the same interpretation across data conditions (Yoon & Millsap, 2007). Per group, there was either the same sample size per group (500 or 1000 in both G₁ and G₂) or different sample sizes per group (500 in G₁ and 1000 in G₂). This design allowed us to compare the new method’s ability to detect noninvariant items in conditions that traditional methods currently usually struggle (i.e., disparate and/or small sample). This resulted in 18 separate conditions. For each condition, 500 datasets were simulated.

Measurement invariance was tested on each simulated dataset using both EGA in the {EGAnet} package (Version 1.1.1; Golino & Christensen, 2022) and SEM using the {lavaan} (Version 0.6.17; Rosseel, 2012) and {semTools} in R. All analyses were conducted in R and full code can be found in Jamison, Golino, and Christensen (2024).

Data Analysis

To assess the accuracy of each model’s (non)invariance detection, we used confusion matrix metrics. Because loadings were only changed for one variable (Variable 5 in Factor 1) and all other variables had equivalent loadings in the population, noninvariance should only be detected for Variable 5. Therefore, a true positive (TP) occurs when the model identifies noninvariance in Variable 5, and a false positive (FP) occurs when any other variable is identified as noninvariant. A false negative (FN) occurs when the model identifies Variable 5 as invariant and a true negative (TN) occurs when any other variable is identified as invariant. An item is considered noninvariant if its significance is $p<0.05$ and invariant if it’s significance is $p\ge 0.05$.

The following confusion matrix metrics were used to identify more specific measures of accuracy: $HitRate$, Sensitivity, Specificity, and F1. The {caret} package (version 6.0.94; Kuhn, 2022) in R was used to calculate Sensitivity, Specificity, and F1. All metrics were calculated separately using both uncorrected and corrected (using the BH-procedure) p-values.

$HitRate$, or $\frac{TP+TN}{TP+FP+TN+FN}$, provides a straightforward, overall assessment of method accuracy of correctly identifying invariance or noninvariance. Sensitivity or $\frac{TP}{TP+FN}$ represents the proportion of true positives correctly identified by the method out of all the truly noninvariant items. Specificity, or $\frac{TN}{TN+FP}$, represents the proportion of true negatives correctly identified by the method out of all true null hypotheses. F1, or $\frac{2TP}{2TP+FP+FN}$, provides a similar metric to Sensitivity but places greater emphasis on identifying Variable 5 as noninvariant relative to identifying the other variables as invariant.

It’s important to contextualize these measures in our current simulation. Because there is only one possible TP or FN (i.e., Variable 5), Sensitivity breaks down to TP. $HitRate$ either breaks down to $\frac{TN}{FP+TN+FN}$ if Variable 5 is identified as invariant or $\frac{TP+TN}{TP+FP+TN}$ when Variable 5 is identified as noninvariant. Similarly, F1 either breaks down to zero if Variable 5 is identified as invariant or $\frac{2TP}{2TP+FP}$ when Variable 5 is identified as noninvariant. Therefore, F1 is weighted toward identifying the noninvariant variable while lowering the (relative) cost of a FP. Finally, Specificity is a pure measure of the extent to which all invariant variables are correctly identified as invariant. Within the context of our study, greater weight should be given to detecting noninvariance over invariance. Therefore, Sensitivity and F1 should be given greater weight over Specificity and $HitRate$.

Effect of MCP on p-Values

Configural invariance was recovered in 99.73% of the simulated datasets using EGA and 100% using SEM. To provide a direct, full comparison between both methods (rather than removing items showing configural noninvariance), only those datasets where configural invariance was found for EGA were retained and the others were discarded. All datasets were retained for analysis using SEM methods. Within each method, the accuracy of metric invariance methods were assessed. Figures 1 and 2 show the mean and 95% confidence interval across all datasets of the p-values split by method, sample size, correlation between factors, and loading difference. A dashed line intercepts the y-axis at .05 representing the α level. The mean p-value is represented for each variable.

Click to enlarge

Figure 1

Mean Uncorrected p-Value by Variable

Click to enlarge

Figure 2

Mean Corrected p-Value by Variable

In both Figures 1 and 2, Variable 5 is the only variable which should be significant or should show a mean p-value consistently below .05. Across both uncorrected (Figure 1) and corrected (Figure 2) p-values for all four methods, regardless of condition, the lowest mean p-value across variables is indeed Variable 5. This is in line with the manipulation used, changing the loading between groups by either 0.2 or 0.4 for only Variable 5. When the p-value is not corrected, as in Figure 1, the Free method has a lower mean p-value for all variables in Factor 1 (where noninvariance was simulated to exist), but not in Factor 2 where no noninvariance was simulated. This pattern is not present for the other methods. When the p-value is corrected (see Figure 2), the average p-value is higher than when it is uncorrected regardless of whether an item is invariant or not. All 3 SEM methods have a more noticeable increase in average p-value for Variable 5 when the difference in loadings for Variable 5 is set to 0.2 and sample size is either different or 500. Under these same conditions, this same trend in EGA is only noticeable when the correlation between factors increases to 0.7.

Hit Rate

In almost all cases, corrected p-values produce a higher mean $HitRate$ than uncorrected p-Values (Figure 3). When the difference in loadings is 0.4, EGA, Fixed, and Wald all have almost perfect $HitRate$ across all variables. In this condition, the same trend arises in Free as was seen in Figures 1 and 2: mean $HitRate$ is lower in general for items in Factor 1, however its level of mean $HitRate$ for Factor 2 when noninvariance is not present, is more similar to that of the other three methods.

Click to enlarge

Figure 3

Hit Rate by Condition

When the difference in loading is set to 0.2, for Variable 5, all four methods experience lower mean $HitRate$ when the p-value is corrected as compared to the uncorrected p-value. This trend is most notable for Fixed, Free, and Wald when sample size is “Different” or 500, but does not appear when sample size is 1000. EGA only shows this trend when sample size is 500 and gradually becomes more disparate as the correlation between factors increases from 0.3 to 0.7. The magnitude of this effect is the same for Fixed, Free, and Wald regardless of the correlation between factors. This indicates that EGA’s ability to correctly identify noninvariant variables is not as heavily influenced by data structures as Fixed, Free, and Wald. The Free method is better able to accurately identify invariant variables when noninvariant items are not present in the same factor.

Overall Metrics

Looking at the accuracy methods overall (not split by simulation condition), when a correction is applied, an interesting pattern appears (Table 1). Both F1 and Specificity increase for all four methods, but Sensitivity decreases. Using uncorrected p-values, Sensitivity is nearly 1 for all four methods, EGA being the highest at 0.99 and Fixed the lowest at 0.96. Once the BH-procedure is applied, Sensitivity decreases for all four methods, most dramatically for Fixed and Wald, falling below 0.90. When using corrected p-values, EGA has the highest values for F1 (0.91) and is tied for the highest Specificity with Fixed and Wald at 0.99. Free has the lowest values (using corrected p-values) of all of four methods for both F1 (0.83) and Sensitivity (0.97). For all four methods, Sensitivity increased slightly by applying the BH-procedure, while F1 values dramatically increase by applying the BH-procedure, going up on average by 0.14.

Sensitivity

When the difference in loadings is 0.4, all methods in all conditions have perfect Sensitivity regardless of whether or not the p-value is corrected (Figure 4). When the difference in loadings is 0.2, uncorrected p-values lead to a higher level of Sensitivity. In this condition, almost perfect Sensitivity is achieved using corrected p-values when sample size size is 1000 for all methods. When the difference in loadings is set to 0.2, corrected p-values are used, and sample size is either “Different” or 500, EGA and Free are performing better than Fixed and Wald. However, EGA is more heavily influenced by the increase in correlation between factors; when the correlation between factors reaches 0.7, EGA’s performance falls below Free’s to the same level as Fixed and Wald. Though when the correlation between factors is 0.3 or 0.5, EGA outperforms Free. All in all, setting the difference in loadings to 0.4 does not affect the ability of any of the methods to identify TP’s (noninvariant variables). However, when the difference is lower, correcting p-values lowers the Sensitivity for all the methods. EGA is, again, less affected by this difference and sample size in its ability to detect TP’s except when the correlation between factors is high.

Click to enlarge

Figure 4

Sensitivity by Condition

F1

In all conditions and across all four methods, corrected p-values produce higher F1 values than uncorrected p-values (Figure 5). When the difference between loadings is set to 0.4, EGA, Fixed, and Wald have similar (and nearly perfect) F1 values. Free, however, has lower F1 values in this condition than the other three methods, particularly when the sample size is increased to 1000. When the difference between loadings is set to 0.2 and F1 is calculated using corrected p-values, a similar pattern arises that was seen in Sensitivity. EGA outperforms the other three methods when sample size is “Different” or 500. However, EGA is more heavily influenced by the increase in correlation between factors; when the correlation between factors reached 0.7, EGA’s performance falls below Free’s to the same level as Fixed and Wald. When the correlation between factors is 0.3 or 0.5, EGA outperforms Free.

Click to enlarge

Figure 5

F1 by Condition

Specificity

Across all these conditions, Specificity calculated using corrected p-values is higher than uncorrected p-values (Figure 6). All methods have consistently high and comparable levels of Specificity, except for the same trend that has been appearing for Free. When the difference in loadings increases from 0.2 to 0.4, the Specificity for the Free method decreases. Altogether, this indicates that each method is able to comparably recover TN’s or invariant items (except for the Free method in one condition).

Click to enlarge

Figure 6

Specificity by Condition

Applying the Test for Metric Invariance to the BAPQ

To demonstrate a substantive application of this approach, we apply our proposed test for metric invariance to the Broad Autism Phenotype Questionnaire (BAPQ; Hurley et al., 2007). Appendix B contains the results from the application of the traditional partial invariance SEM method as implemented by the partialInvariance() function in {semTools}. Data was obtained from the Simons Foundation Powering Autism Research for Knowledge (SPARK) of the Simons Foundation Autism Research Initiative (SFARI), a large research initiative which has collected data from over 50,000 individuals with autism and their families (Feliciano et al., 2018). The BAPQ is a 36-item questionnaire designed to assess autism-related traits in adults. Participants are asked to rate the how often a statement applies to them on a 6-point Likert scale ranging from (1) Very Rarely to (6) Very Often. Items were intended to relate to one of three domains: aloofness, rigid personality, or pragmatic language.

This questionnaire was given to the parents (either mother or father) of an autistic child to assess their phenotypic level of autistic traits. We begin assessing measurement invariance between mothers and fathers by establishing configural invariance. To do so we apply EGA separately to the data on mothers and the data on fathers and compare their community structures. For this example, we are using the {EGAnet} package (Version 2.0.6; Golino & Christensen, 2024).

# Load EGAnet Package
library(EGAnet)

# Load in the Data
load("../2. Data/bapq.all.RData")

# Set mother indices
mother <- bapq.all$Parent == "Mother"

# Extract items only
items <- bapq.all[,4:39]

## Mother
ega.mother <- EGA (data = items[mother,])

## Father
ega.father <- EGA(data = items[!mother,])

Visually we can see that the two graphs contain nonequivalent community structures (Figure 7). We can apply the invariance() function to the data, which will first identify a common structure that exists using bootEGA(), removing item stabilities less than 0.70, to establish configural invariance. After, establishing configural invariance, the procedure will continue to test metric invariance.¹

# Perform invariance
bapq_invariance <- invariance(
  data = items, group = bapq.all$Parent, 
  ncores = 8, seed = 1, loading.method = "experimental"
)

Click to enlarge

Figure 7

Comparison of the EGA Networks for the Mothers (Left) and Fathers (Right) With Color Denoting the Community for Each Node

The function will print out how many items were identified for configural invariance, for example: Configural invariance was found with 32 variables. To view which items were removed from the the original 36, the following object can be accessed:

	[1] "q12" "q23" "q25" "q28"

In Figure 8, the before and after item stability are plotted with the latter being accessed in the results using plot(bapq_invariance$configural.results$item_stability).

Click to enlarge

Figure 8

Comparison of the BAPQ Item Stability Before (Left) and After (Right) Configural Invariance

Evaluating each group separately, we can see that the both groups have equivalent structures (Figure 9):

## Father
ega.father <- EGA(data = items[!mother, stable_names])

## Mother
ega.mother <- EGA(data = items[mother, stable_names])

Click to enlarge

Figure 9

Comparison of the EGA Networks for the Fathers (Left) and Mothers (Right) With Color Denoting the Community for Each Node

Finally, we can print the results of invariance to see a table that breaks down the metric invariance for each item:

# Print summary
summary(bapq_invariance)

      Membership   Difference     p       p_BH   sig  Direction
q01       1          -0.006     0.662    0.850
q05       1           0.028     0.038    0.174    *   Father > Mother
q09       1           0.022     0.208    0.428
q16       1           0.027     0.054    0.192    .
q18       1          -0.001     0.948    0.979
q27       1           0.030     0.076    0.243    .
q31       1           0.033     0.030    0.174    *   Father > Mother
q36       1          -0.018     0.354    0.539
q02       2           0.017     0.296    0.515
q04       2           0.003     0.870    0.960
q07       2          -0.021     0.214    0.428
q10       2          -0.017     0.332    0.531
q11       2           0.038     0.038    0.174    *   Father > Mother
q14       2          -0.004     0.810    0.926
q17       2          -0.016     0.306    0.515
q20       2          -0.059     0.004    0.064   **   Father < Mother
q21       2          -0.043     0.004    0.064   **   Father < Mother
q29       2           0.044     0.006    0.064   **   Father > Mother
q32       2           0.019     0.214    0.428
q34       2           0.024     0.152    0.428
q03       3          -0.032     0.044    0.176    *   Father < Mother
q06       3           0.001     0.922    0.979
q08       3          -0.004     0.798    0.926
q13       3           0.010     0.564    0.785
q15       3          -0.017     0.264    0.497
q19       3          -0.007     0.664    0.850
q22       3           0.039     0.030    0.174    *   Father > Mother 
q24       3           0.022     0.182    0.428
q26       3          -0.019     0.186    0.428
q30       3           0.000     0.984    0.984
q33       3           0.005     0.708    0.871
q35       3          -0.011     0.422    0.614
----
Signif. code: 0 ’***’  0.001 ’**’  0.01 ’*’  0.05 ’.’  0.1 ’n.s.’ 1

The item text for the 6 items showing metric noninvariance using uncorrected p-values are displayed in Table 2. Using the corrected p-values, there were no items that were detected as noninvariant. The noninvariant items detected with the uncorrected p-values spanned each dimension of the BAPQ. All differences between mothers and fathers corresponded with deficits in fathers relative to mothers—that is, there was larger loadings for items related to deficits or behaviors contrary to norms in the general population or smaller loadings for items related to norms in the general population.

These results can further be visualized by using the plot() function on the output of invariance (Figure 10).

Click to enlarge

Figure 10

Comparison of the EGA Networks for the Fathers (Left) and Mothers (Right) With Color Denoting the Community for Each Node and Transparency Denoting (Non)invariance

Conclusion

Ensuring the equivalence of a measure across assessment groups is vital to the efficacy of group comparison. Although many methods have been proposed to improve measurement invariance in the SEM framework, many unresolved problems still remain. The EGA approach proposed in this study performed comparably to existing SEM methods and, in several conditions, outperformed them with the aim of detecting noninvariant items. The method was then applied to a substantive dataset to demonstrate its assessment of metric invariance in a real-world dataset, finding important differences in the BAPQ inventory that exist between mothers and fathers of children with ASD.

Appendix A

Hit Rate

Click to enlarge

Figure 11

Mean Hit Rate by Condition

Sensitivity

Click to enlarge

Figure 12

Sensitivity by Condition

F1

Click to enlarge

Figure 13

F1 by Condition

Specificity

Click to enlarge

Figure 14

Specificity by Condition

Appendix B

#Loading libraries
library (semTools)
library (lavaan)

# Confirmatory Three Factor Model
cfa.model <- "
aloof =~ q01 + q05 + q09 + q12 + q16 + q18 + q23 + q25 + q27 + q28 + q31 + q36
pragmatic =~ q02 + q04 + q07 + q10 + q11 + q14 + q17 + q20 + q21 + q29 + q32 + q34
rigid =~ q03 + q06 + q08 + q13 + q15 + q19 + q22 + q24 + q26 + q30 + q33 + q35"

cfa_fit <- cfa(cfa.model, data = bapq.all)

# Obtaining item names
items <- bapq.all[,4:39]
var.names <- names(items)

# Assessing EFA from 1 to 4 factors
fit <- efa(data = bapq.all[,var.names], nfactors = 1:4)

# Four-Factor configural invariance model
conf <- "
f1 =~ q02 + q04 + q14 + q17 + q20 + q29 + q32
f2 =~ q03 + q06 + q08 + q13 + q15 + q19 + q22 + q24 + q26 + q30 + q33 + q35
f3 =~ q01 + q05 + q09 + q10 + q11 + q12 + q16 + q18 + q23 + q25 + q27 + q28 + q31 + q36
f4 =~ q07 + q21 + q34"

configural <- cfa(conf, data = bapq.all, std.lv = TRUE, group = "Parent")

# Two-Factor configural invariance model
conf <- "
f2 =~ q03 + q08 + q13 + q19 + q22 + q24
f3 =~ q01 + q09 + q16 + q23 + q25 + q36"

configural  <- cfa (conf,  data =  bapq.all, std.lv = TRUE, group = "Parent")

# Metric invariance model
weak <- "
f2 =~ q03 + q08 + q13 + q19 + q22 + q24
f3 =~ q01 + q09 + q16 + q23 + q25 + q36
f2 ~~ c(1, NA)*f2
f3 ~~ c(1, NA)*f3"
weak <- cfa(weak, data = bapq.all, group="Parent", group.equal="loadings"
models <- list(fit.configural = configural, fit.loadings = weak)

# Partial invariance models
pi_model <- partialInvariance(models, "metric")
pi_model.h <- partialInvariance(models, "metric", p.adjust = "hochberg")