Difficulty testing measurement invariance across gender and using CFA factor scores for an interaction in Stata SEM

[Rebecca Oteng]

Hello everyone,

I am working on a PhD-level term paper using the Ghana Afrobarometer Round 10 dataset (N = 2,400). I am using Stata 16 and SEM/CFA to examine a multidimensional measure of empowerment.

My current measurement model uses six 5-point Likert-type items (Q86A–Q86F), specified as three correlated factors, with two indicators per factor:

• Reproductive agency: Q86A, Q86B
• Normative support/gender norms: Q86C, Q86D
• Decision-making autonomy: Q86E, Q86F

My hypotheses are:

• H1: Empowerment items load on ≥ 2 distinct factors rather than one.
• H2: The measurement model is invariant across gender (configural → metric → scalar), permitting valid score comparison.
• H3: Religion and region interact in predicting decision-autonomy, not merely add.

The three-factor CFA fits reasonably well when estimated for the full sample and when estimated separately for men and women. For example, the overall three-factor model gives approximately:

• χ²(6) = 58.93
• RMSEA = .061
• CFI = .981
• TLI = .953
• SRMR = .026

The separate gender models also converge and have similar fit.

However, I am having difficulty with H2. I tried estimating the three-factor model as a multiple-group SEM using gender (Q101) as the grouping variable. My understanding is that I should proceed sequentially from:

1. Configural invariance — same factor structure, parameters free across groups
2. Metric invariance — factor loadings constrained equal
3. Scalar invariance — factor loadings and intercepts constrained equal

I have tried commands such as:
sem (REPRODUCTIVEAGENCY -> Q86A_c Q86B_c) /// (NORMATIVESUPPORT -> Q86C_c Q86D_c) /// (DECISIONMAKING -> Q86E_c Q86F_c), /// group(Q101) ginvariant(none)
The problem is that the configural model itself fails to converge. Stata goes through many iterations and repeatedly reports (not concave). The same basic model converges very quickly when I estimate men and women separately, so I am unsure whether this is a problem with model identification, the fact that each latent factor has only two indicators, the ordinal nature of the five-point items, or the way I have specified the multiple-group model in Stata.

I have also tried more constrained specifications, but I now realize that a model with all parameters constrained equal is not a configural model.

I also conducted an exploratory factor analysis on the six items. A two-factor solution was suggested by the eigenvalues, with the rotated pattern showing two broad clusters. However, my theoretically specified three-factor CFA has good fit. A three-factor ML EFA also produced a Heywood-type problem, so I am not sure how much weight to give the EFA result.

My second problem concerns H3.

I want to use the observed factor score for the decision-making autonomy construct as the dependent variable and test whether religion and region interact:
regress decisionmaking_score i.religion_group##i.REGION
I would then test the interaction using:
testparm i.religion_group#i.REGION
and examine predicted values using:
margins religion_group#REGION marginsplot
My concern is that because I have not been able to establish measurement invariance across gender, I am unsure whether using factor scores from the CFA is appropriate, particularly if I want to make comparisons across groups.

I would therefore appreciate advice on two issues:

1. How should I correctly approach the measurement invariance problem in Stata when the three-factor model converges separately by gender but the configural multiple-group model does not converge? Is the two-indicator-per-factor specification likely to be causing the problem, and would you recommend changing the hypothesis/model?
2. Given this measurement problem, what is the most defensible way to generate and use factor scores for H3 and test the religion × region interaction? Should I use factor scores from the full-sample CFA, estimate the measurement model separately, or use another approach?

I am particularly interested in advice that is appropriate for Stata 16 and SEM/CFA, rather than general advice on SEM software.

[Erik Ruzek]

Almost certainly, two indicators per factor is the major limitation here. In order to run any model such as this (multiple 2-item factors), it's important to recognize that the individual factors themselves are not identified on their own. They borrow identification from their correlations with the other factors. When you throw in that you are estimating separate models by group, you make the problem even worse because you are decreasing the sample size and weird things can happen because you have so little information per parameter. That said, it is possible to estimate these models, but you have to do some constraining for any hope of it working. Below is the sequence I would use, and note that I am generating the data to be pretty clean (moderately correlated factors, relatively large loadings, and true mean differences). To the extent that your data has, for example, extremely low or high correlations, small loadings, or cross-loadings, you are likely to experience problems with this factor structure.

I also simulated a cross-loading pattern so you could see how to handle it. Please note that you can probably only get away with one cross-loading item and the cross-loading item cannot be the scaling indicator for either factor it loads on. You must explicitly do so in your code, as shown below, as Stata will not automatically do this for you. I do not show invariance testing in the cross-loading case, but it only adds to the loss of information problem you have with two-item factors in a multigroup setting.

The code below was generated in collaboration with Claude Sonnet 5. Claude is a bit wordy in its code but I think it is helpful for those trying to learn what the code is doing.

*-----------------------------------------------------------------
* Toy simulation: 3-factor CFA, 2 indicators per factor
* Illustrates identification via the "two-indicator rule"
*-----------------------------------------------------------------
version 16
clear all
set seed 20260915
set obs 500

*-----------------------------------------------------------------
* 1. Simulate three correlated latent factors
*    (variances = 1, moderate positive correlations among them --
*     these correlations are what identify the 2-indicator factors)
*-----------------------------------------------------------------
matrix C = (1.00, 0.40, 0.30 \ ///
            0.40, 1.00, 0.35 \ ///
            0.30, 0.35, 1.00)

drawnorm f1 f2 f3, corr(C) means(0 0 0) sds(1 1 1)

*-----------------------------------------------------------------
* 1b. Generate a grouping variable and impose population mean
*     differences on the factors: cat==1 is higher on F1 and F2,
*     lower on F3. Loadings/residual variances are left identical
*     across groups, so the population model is configurally,
*     metrically, AND scalarly invariant.
*-----------------------------------------------------------------
gen byte cat = (runiform() < 0.5)

replace f1 = f1 + 0.5 if cat == 1
replace f2 = f2 + 0.4 if cat == 1
replace f3 = f3 - 0.3 if cat == 1

*-----------------------------------------------------------------
* 2. Generate two indicators per factor
*    item = loading*factor + unique residual
*    (loadings chosen to give roughly 0.6-0.8 standardized loadings)
*-----------------------------------------------------------------
gen y1 = 0.8*f1 + rnormal(0, 0.6)   // factor 1 indicators
gen y2 = 0.7*f1 + rnormal(0, 0.7)

gen y3 = 0.75*f2 + rnormal(0, 0.65)  // factor 2 indicators
gen y4 = 0.65*f2 + rnormal(0, 0.75)

gen y5 = 0.8*f3 + rnormal(0, 0.6)   // factor 3 indicators
gen y6 = 0.7*f3 + rnormal(0, 0.7)

*-----------------------------------------------------------------
* 2b. A cross-loading version of y2: loads on F1 (primary, 0.7,
*     same as the clean y2 above) AND on F2 (secondary, 0.3).
*     Used later (section 7) to illustrate fitting a CFA with one
*     cross-loading; y2 itself is left untouched so the clean
*     baseline/invariance models above are unaffected.
*-----------------------------------------------------------------
gen y2_cl = 0.7*f1 + 0.3*f2 + rnormal(0, 0.65)

drop f1 f2 f3   // pretend these are unobserved, as in real data

*-----------------------------------------------------------------
* 3. Fit the CFA
*    - marker-indicator scaling: first loading per factor fixed to 1
*    - must freely estimate three factor covariances for
*      identification
*    - no cross-loadings, no correlated residuals
*-----------------------------------------------------------------
sem (F1 -> y1 y2) (F2 -> y3 y4) (F3 -> y5 y6), ///
    latent(F1 F2 F3) cov(F1*F2 F1*F3 F2*F3) standardized

estat gof, stats(all)     // overall fit; df should be 6 for this setup
estat mindices            // check for obvious misspecification

*-----------------------------------------------------------------
* 6. Measurement invariance across cat: configural -> metric -> scalar
*    Marker-variable scaling (y1, y3, y5 fixed to 1) throughout, so
*    ginvariant() controls only loadings (mcoefficients) and item
*    intercepts (mcuts). Latent means are fixed at 0 in the
*    reference group (cat==0) and freely estimated in cat==1 once
*    intercepts are constrained equal (scalar step) -- that's what
*    lets you recover the population mean differences you built in.
*-----------------------------------------------------------------

* --- Configural invariance: same pattern of free/fixed parameters
*     in both groups, nothing constrained equal across groups
*     NOTE: means(F1@0 F2@0 F3@0) fixes the latent means at 0 in
*     BOTH groups. This is not optional at this stage -- with item
*     intercepts still free across groups, a free latent mean and
*     the group-specific intercepts are perfectly aliased (an
*     observed mean difference can be explained by either one).
*     Latent means only become identifiable once
*     intercepts are constrained equal, at the scalar step below.
sem (F1 -> y1 y2) (F2 -> y3 y4) (F3 -> y5 y6), ///
    latent(F1 F2 F3) cov(F1*F2 F1*F3 F2*F3) group(cat) ginvariant(none) ///
    means(F1@0 F2@0 F3@0)
estat gof, stats(all)
estimates store m_configural

* --- Metric invariance: factor loadings held equal across groups
*     (latent means still fixed at 0 in both groups, same reason)
sem (F1 -> y1 y2) (F2 -> y3 y4) (F3 -> y5 y6), ///
    latent(F1 F2 F3) cov(F1*F2 F1*F3 F2*F3) group(cat) ///
    ginvariant(mcoef) means(F1@0 F2@0 F3@0)
estat gof, stats(all)
estimates store m_metric

* --- Scalar invariance: loadings AND item intercepts held equal;
*     this is what identifies the latent mean differences below
sem (F1 -> y1 y2) (F2 -> y3 y4) (F3 -> y5 y6), ///
    latent(F1 F2 F3) cov(F1*F2 F1*F3 F2*F3) group(cat) ///
    ginvariant(mcoef mcons)
estat gof, stats(all)
estimates store m_scalar

* --- Nested-model chi-square difference tests
lrtest m_configural m_metric
lrtest m_metric m_scalar

*-----------------------------------------------------------------
* 7. Single-group CFA with one cross-loading
*    y2_cl loads on both F1 (primary, population value 0.7) and F2
*    (secondary, population value 0.3). y1 stays the marker for F1
*    and y3 stays the marker for F2 -- neither marker is the
*    cross-loading item, so scale-setting and the cross-loading
*    estimate stay cleanly separated. F1 keeps y1 as a pure
*    indicator and F2 keeps y3/y4 as pure indicators, so both
*    factors retain at least one indicator that loads on them
*    exclusively.
*-----------------------------------------------------------------

sem (F1 -> y1 y2_cl) (F2 -> y3@1 y4 y2_cl) (F3 -> y5 y6), ///
    latent(F1 F2 F3) cov(F1*F2 F1*F3 F2*F3) means(F1@0 F2@0 F3@0) ///
    standardized

estat gof, stats(all)     // df should be 5 here vs. 6 in the clean model
estimates store m_crossload