Skip to main content Skills Marketplace Descubre y explora habilidades de IA creadas por la comunidad.
Instalar con Codex o Claude Copia este prompt, pégalo en Codex, Claude u otro asistente, y deja que revise la página de la skill y la instale por ti.
Mostrar detalles del prompt Un comando directo omite el prompt de revisión. Revisa el origen antes de ejecutarlo.
npx skills add https://github.com/xjtulyc/awesome-rosetta-skills --skill psychometricsEl comando permanece en una sola línea. Desplázate horizontalmente para revisarlo antes de copiarlo.
¿Prefieres una copia local? Descarga los archivos que SkillsMP tiene disponibles ahora.
Descargar Zip Descargando... Ocupaciones relacionadas SOC
Basado en la clasificación ocupacional SOC
name psychometrics description Psychometric analysis in Python: CTT item stats, Cronbach's alpha, EFA, CFA with semopy, 2PL IRT, measurement invariance, and Mantel-Haenszel DIF detection.
tags ["psychology","psychometrics","factor-analysis","item-response-theory","python","scale-development"] version 1.0.0 authors [{"name":"awesome-rosetta-skills contributors","github":"@awesome-rosetta-skills"}] license MIT platforms ["claude-code","codex","gemini-cli","cursor"] dependencies ["pandas>=2.0.0","numpy>=1.24.0","scipy>=1.10.0","factor_analyzer>=0.5.0","semopy>=2.3.0","matplotlib>=3.7.0","scikit-learn>=1.3.0"] last_updated 2026-03-17
Psychometric Analysis in Python
This skill covers the full psychometric validation pipeline:
Classical Test Theory (CTT) item statistics, reliability (Cronbach's alpha),
Exploratory Factor Analysis (EFA) with parallel analysis, Confirmatory Factor
Analysis (CFA) with semopy, Item Response Theory (IRT) 2PL models,
measurement invariance testing, scale development workflow, and differential
item functioning (DIF) via the Mantel-Haenszel procedure.
1. Setup
pip install pandas numpy scipy factor_analyzer semopy matplotlib scikit-learn
import numpy as np
import pandas as pd
from scipy import stats
from scipy.optimize import minimize
from scipy.special import expit
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore" )
from factor_analyzer import FactorAnalyzer
from factor_analyzer.factor_analyzer import calculate_kmo, calculate_bartlett_sphericity
import semopy
np.random.seed(42 )
N_RESPONDENTS = 500
N_ITEMS = 20
F1 = np.random.normal(0 , 1 , N_RESPONDENTS)
F2 = np.random.normal(0 , 1 , N_RESPONDENTS)
loadings_F1 = np.array([0.75
0.72
0.68
0.80
0.65
0.70
0.73
0.60
0.78
0.55
0.70
0.74
0.65
0.80
0.72
0.68
0.76
0.62
0.55
0.78
f"item_{i+1 :02d} "
for
in
range
for
in
range
10
1
2
0
1
for
in
range
10
20
10
1
10
2
0
1
for
in
True
5
1
2
3
4
5
float
"group"
1
0.5
print
f"Synthetic data: {responses.shape[0 ]} respondents × {N_ITEMS} items"
2. Classical Test Theory Item Statistics def compute_item_stats (response_matrix: pd.DataFrame ) -> pd.DataFrame:
"""
Compute Classical Test Theory item-level statistics.
Parameters
----------
response_matrix : pd.DataFrame
Rows = respondents, columns = items. Values are ordinal scores.
Returns
-------
pd.DataFrame with columns:
mean, std, p_value (difficulty), skewness, kurtosis,
item_total_r, corrected_item_total_r, alpha_if_deleted.
"""
data = response_matrix.copy()
n_items = data.shape[1 ]
total = data.sum (axis=1 )
rows = []
for col in data.columns:
item_scores = data[col]
rest_scores = total - item_scores
itc, _ = stats.pearsonr(item_scores, total)
citc, _ = stats.pearsonr(item_scores, rest_scores)
max_score = item_scores.max ()
min_score = item_scores.min ()
p_val = (item_scores.mean() - min_score) / max (max_score - min_score, 1e-9 )
remaining = data.drop(columns=[col])
alpha_del = _cronbach_alpha_raw(remaining)
rows.append({
"item" : col,
"mean" : round (item_scores.mean(), 3 ),
"std" : round (item_scores.std(), 3 ),
"p_value" : round (p_val, 3 ),
"skewness" : round (float (stats.skew(item_scores)), 3 ),
"kurtosis" : round (float (stats.kurtosis(item_scores)), 3 ),
"item_total_r" : round (itc, 3 ),
"corrected_item_total_r" : round (citc, 3 ),
"alpha_if_deleted" : round (alpha_del, 3 ),
})
return pd.DataFrame(rows).set_index("item" )
def _cronbach_alpha_raw (data: pd.DataFrame ) -> float :
"""Internal helper: compute alpha from a DataFrame of item scores."""
k = data.shape[1 ]
if k < 2 :
return np.nan
item_var = data.var(ddof=1 , axis=0 ).sum ()
total_var = data.sum (axis=1 ).var(ddof=1 )
if total_var == 0 :
return np.nan
return (k / (k - 1 )) * (1 - item_var / total_var)
3. Cronbach's Alpha and Omega def cronbach_alpha (data: pd.DataFrame, ci: float = 0.95 ) -> dict :
"""
Compute Cronbach's alpha with bootstrapped confidence interval.
Parameters
----------
data : pd.DataFrame
Item score matrix (respondents × items).
ci : float
Confidence level for bootstrap CI.
Returns
-------
dict with keys: alpha, ci_lower, ci_upper, n_items, n_respondents.
"""
alpha = _cronbach_alpha_raw(data)
n = data.shape[0 ]
bs_alphas = []
for _ in range (1000 ):
idx = np.random.choice(n, size=n, replace=True )
bs_alphas.append(_cronbach_alpha_raw(data.iloc[idx]))
lo = np.percentile(bs_alphas, (1 - ci) / 2 * 100 )
hi = np.percentile(bs_alphas, (1 + ci) / 2 * 100 )
result = {
"alpha" : round (alpha, 4 ),
"ci_lower" : round (lo, 4 ),
"ci_upper" : round (hi, 4 ),
"n_items" : data.shape[1 ],
"n_respondents" : data.shape[0 ],
"interpretation" : (
"Excellent (≥0.90)" if alpha >= 0.90 else
"Good (0.80-0.89)" if alpha >= 0.80 else
"Acceptable (0.70-0.79)" if alpha >= 0.70 else
"Questionable (0.60-0.69)" if alpha >= 0.60 else
"Poor (<0.60)"
),
}
print (f"Cronbach's alpha = {alpha:.4 f} [{lo:.4 f} , {hi:.4 f} ] — {result['interpretation' ]} " )
return result
4. Exploratory Factor Analysis def run_efa (
data: pd.DataFrame,
n_factors: int | None = None ,
rotation: str = "oblimin" ,
method: str = "minres" ,
parallel_analysis: bool = True ,
) -> dict :
"""
Run Exploratory Factor Analysis with optional parallel analysis.
Parameters
----------
data : pd.DataFrame
Item response matrix.
n_factors : int or None
Number of factors. If None, determined by parallel analysis or Kaiser.
rotation : str
'oblimin' (oblique) or 'varimax' (orthogonal).
method : str
Extraction method: 'minres', 'ml', 'principal'.
parallel_analysis : bool
Use parallel analysis to select n_factors if n_factors is None.
Returns
-------
dict with keys: loadings, communalities, eigenvalues, n_factors, variance.
"""
chi2, p_bartlett = calculate_bartlett_sphericity(data)
kmo_all, kmo_model = calculate_kmo(data)
print (f"Bartlett's test: chi2={chi2:.2 f} , p={p_bartlett:.4 f} " )
print (f"KMO measure: {kmo_model:.3 f} (>0.60 acceptable, >0.80 meritorious)" )
fa_initial = FactorAnalyzer(n_factors=data.shape[1 ], rotation=None )
fa_initial.fit(data)
ev, _ = fa_initial.get_eigenvalues()
if n_factors is None and parallel_analysis:
random_evs = []
for _ in range (100 ):
rnd = pd.DataFrame(
np.random.normal(size=data.shape),
columns=data.columns,
)
fa_rnd = FactorAnalyzer(n_factors=data.shape[1 ], rotation=None )
fa_rnd.fit(rnd)
ev_rnd, _ = fa_rnd.get_eigenvalues()
random_evs.append(ev_rnd)
mean_rnd_ev = np.mean(random_evs, axis=0 )
n_factors = int ((ev > mean_rnd_ev).sum ())
n_factors = max (1 , n_factors)
print (f"Parallel analysis suggests {n_factors} factor(s)." )
elif n_factors is None :
n_factors = int ((ev > 1 ).sum ())
print (f"Kaiser criterion suggests {n_factors} factor(s)." )
fa = FactorAnalyzer(n_factors=n_factors, rotation=rotation, method=method)
fa.fit(data)
loadings = pd.DataFrame(
fa.loadings_,
index=data.columns,
columns=[f"F{i+1 } " for i in range (n_factors)],
)
communalities = pd.Series(fa.get_communalities(), index=data.columns, name="h2" )
variance_df = pd.DataFrame(
fa.get_factor_variance(),
index=["SS Loadings" , "Proportion Var" , "Cumulative Var" ],
columns=[f"F{i+1 } " for i in range (n_factors)],
)
print ("\n=== Factor Loadings (|loading| > 0.30 shown) ===" )
print (loadings.where(loadings.abs () > 0.30 , "" ).to_string())
print ("\n=== Variance Explained ===" )
print (variance_df.round (4 ).to_string())
fig, ax = plt.subplots(figsize=(8 , 4 ))
ax.plot(range (1 , len (ev) + 1 ), ev, "bo-" , linewidth=1.5 , label="Actual eigenvalues" )
if parallel_analysis:
ax.plot(range (1 , len (mean_rnd_ev) + 1 ), mean_rnd_ev, "r--" ,
linewidth=1.2 , label="Parallel analysis (random)" )
ax.axhline(1 , color="gray" , linestyle=":" , linewidth=0.8 , label="Kaiser criterion" )
ax.set_xlabel("Factor Number" )
ax.set_ylabel("Eigenvalue" )
ax.set_title("Scree Plot with Parallel Analysis" , fontweight="bold" )
ax.legend()
ax.set_xlim(0.5 , min (20 , len (ev)) + 0.5 )
plt.tight_layout()
plt.savefig("scree_plot.png" , dpi=150 )
plt.show()
return {
"loadings" : loadings,
"communalities" : communalities,
"eigenvalues" : ev,
"n_factors" : n_factors,
"variance" : variance_df,
"fa_object" : fa,
}
5. Confirmatory Factor Analysis with semopy def run_cfa (
data: pd.DataFrame,
model_str: str ,
estimator: str = "MLW" ,
) -> dict :
"""
Run a CFA model using semopy and return key fit indices.
Parameters
----------
data : pd.DataFrame
Item response matrix.
model_str : str
semopy model specification. Example:
'F1 =~ item_01 + item_02 + item_03
F2 =~ item_11 + item_12 + item_13'
estimator : str
'MLW' (default, robust ML), 'ULS', 'DWLS'.
Returns
-------
dict with keys: cfi, rmsea, srmr, tli, chi2, df, p_chi2, model, results.
"""
mod = semopy.Model(model_str)
res = mod.fit(data, solver=estimator)
stats_obj = semopy.calc_stats(mod)
cfi = float (stats_obj.loc["CFI" , "Value" ]) if "CFI" in stats_obj.index else np.nan
tli = float (stats_obj.loc["TLI" , "Value" ]) if "TLI" in stats_obj.index else np.nan
rmsea= float (stats_obj.loc["RMSEA" ,"Value" ]) if "RMSEA" in stats_obj.index else np.nan
srmr = float (stats_obj.loc["SRMR" , "Value" ]) if "SRMR" in stats_obj.index else np.nan
chi2 = float (stats_obj.loc["chi2" , "Value" ]) if "chi2" in stats_obj.index else np.nan
df_ = float (stats_obj.loc["DoF" , "Value" ]) if "DoF" in stats_obj.index else np.nan
pval = float (stats_obj.loc["chi2 p-value" , "Value" ]) if "chi2 p-value" in stats_obj.index else np.nan
print ("\n=== CFA Fit Indices ===" )
print (f" chi2({df_:.0 f} ) = {chi2:.3 f} , p = {pval:.4 f} " )
print (f" CFI = {cfi:.3 f} (≥0.95 acceptable)" )
print (f" TLI = {tli:.3 f} (≥0.95 acceptable)" )
print (f" RMSEA= {rmsea:.3 f} (≤0.06 close fit, ≤0.08 acceptable)" )
print (f" SRMR = {srmr:.3 f} (≤0.08 acceptable)" )
return {
"cfi" : cfi,
"tli" : tli,
"rmsea" : rmsea,
"srmr" : srmr,
"chi2" : chi2,
"df" : df_,
"p_chi2" : pval,
"model" : mod,
"results" : res,
"stats" : stats_obj,
}
6. IRT 2-Parameter Logistic Model def run_irt_2pl (
responses: np.ndarray,
n_iter: int = 300 ,
tol: float = 1e-5 ,
) -> dict :
"""
Estimate the 2PL IRT model via marginal maximum likelihood.
Model: P(X_ij = 1 | theta_i) = sigmoid(a_j * (theta_i - b_j))
where a_j = discrimination, b_j = difficulty.
Parameters
----------
responses : np.ndarray
Binary response matrix (n_respondents × n_items). 1 = correct/endorsed.
n_iter : int
Maximum EM iterations.
tol : float
Convergence tolerance.
Returns
-------
dict with keys: a (discrimination), b (difficulty), theta (ability),
ll_history.
"""
R = responses.astype(float )
N, J = R.shape
a = np.ones(J)
b = np.zeros(J)
row_prop = R.mean(axis=1 )
row_prop = np.clip(row_prop, 0.01 , 0.99 )
theta = stats.norm.ppf(row_prop)
ll_history = []
for iteration in range (n_iter):
eta = a[None , :] * (theta[:, None ] - b[None , :])
P = expit(eta)
P = np.clip(P, 1e-9 , 1 - 1e-9 )
ll = (R * np.log(P) + (1 - R) * np.log(1 - P)).sum ()
ll_history.append(ll)
W = P * (1 - P)
for j in range (J):
r_j = R[:, j]
p_j = P[:, j]
w_j = W[:, j]
res_j = r_j - p_j
grad_b = -a[j] * res_j.sum ()
hess_b = -a[j]**2 * w_j.sum ()
if abs (hess_b) > 1e-9 :
b[j] -= grad_b / hess_b
theta_b = theta - b[j]
grad_a = (res_j * theta_b).sum ()
hess_a = -(w_j * theta_b**2 ).sum ()
if abs (hess_a) > 1e-9 :
a[j] -= grad_a / hess_a
a[j] = max (a[j], 0.01 )
for i in range (N):
p_i = P[i]
w_i = W[i]
res_i = R[i] - p_i
grad_th = (a * res_i).sum ()
hess_th = -(a**2 * w_i).sum ()
if abs (hess_th) > 1e-9 :
theta[i] -= grad_th / hess_th
theta = (theta - theta.mean()) / max (theta.std(), 1e-9 )
if iteration > 0 and abs (ll_history[-1 ] - ll_history[-2 ]) < tol:
print (f"IRT 2PL converged at iteration {iteration} " )
break
item_df = pd.DataFrame({"a_discrimination" : a, "b_difficulty" : b},
index=[f"item_{j+1 :02d} " for j in range (J)])
return {
"a" : a,
"b" : b,
"theta" : theta,
"item_params" : item_df,
"ll_history" : ll_history,
}
7. Differential Item Functioning (Mantel-Haenszel) def detect_dif (
responses: pd.DataFrame,
group: pd.Series,
n_score_strata: int = 5 ,
alpha: float = 0.05 ,
) -> pd.DataFrame:
"""
Detect DIF using the Mantel-Haenszel (1959) procedure.
Items are flagged if the MH chi-squared test is significant after
Bonferroni correction.
Parameters
----------
responses : pd.DataFrame
Binary or polytomous item responses (respondents × items).
group : pd.Series
Binary group indicator (0 = reference, 1 = focal).
n_score_strata : int
Number of total-score strata for conditioning.
alpha : float
Significance level (before Bonferroni correction).
Returns
-------
pd.DataFrame with columns: item, mh_chi2, p_value, mh_odds_ratio, dif_flag.
"""
items = [c for c in responses.columns if c != group.name]
total_score = responses[items].sum (axis=1 )
strata = pd.qcut(total_score, q=n_score_strata, labels=False , duplicates="drop" )
n_items = len (items)
bonferroni = alpha / n_items
rows = []
for item in items:
A_list, B_list, C_list, D_list = [], [], [], []
for s in strata.unique():
mask = strata == s
ref = group[mask] == 0
foc = group[mask] == 1
threshold = responses[item][mask].median()
endorsed = responses[item][mask] >= threshold
A = int ((ref & endorsed)[mask].sum ())
B = int ((ref & ~endorsed)[mask].sum ())
C = int ((foc & endorsed)[mask].sum ())
D = int ((foc & ~endorsed)[mask].sum ())
sub_ref_end = ((group[mask] == 0 ) & (responses[item][mask] >= threshold)).sum ()
sub_ref_not = ((group[mask] == 0 ) & (responses[item][mask] < threshold)).sum ()
sub_foc_end = ((group[mask] == 1 ) & (responses[item][mask] >= threshold)).sum ()
sub_foc_not = ((group[mask] == 1 ) & (responses[item][mask] < threshold)).sum ()
A_list.append(sub_ref_end)
B_list.append(sub_ref_not)
C_list.append(sub_foc_end)
D_list.append(sub_foc_not)
A, B, C, D = map (np.array, [A_list, B_list, C_list, D_list])
N = A + B + C + D + 1e-9
numerator = (A * D / N).sum ()
denominator = (B * C / N).sum ()
mh_or = numerator / (denominator + 1e-12 )
E_A = (A + B) * (A + C) / N
var = (A + B) * (C + D) * (A + C) * (B + D) / (N**2 * (N - 1 ) + 1e-9 )
chi2 = (abs ((A - E_A).sum ()) ** 2 ) / max (var.sum (), 1e-9 )
p_val = 1 - stats.chi2.cdf(chi2, df=1 )
rows.append({
"item" : item,
"mh_chi2" : round (chi2, 4 ),
"p_value" : round (p_val, 6 ),
"mh_odds_ratio" : round (mh_or, 4 ),
"dif_flag" : "DIF" if p_val < bonferroni else "No DIF" ,
"dif_magnitude" : (
"C (large)" if abs (np.log(mh_or)) > np.log(2.35 ) else
"B (moderate)" if abs (np.log(mh_or)) > np.log(1.5 ) else
"A (negligible)"
),
})
df = pd.DataFrame(rows).sort_values("mh_chi2" , ascending=False ).reset_index(drop=True )
n_dif = (df["dif_flag" ] == "DIF" ).sum ()
print (f"\nDIF Analysis: {n_dif} /{n_items} items flagged (Bonferroni α={bonferroni:.4 f} )" )
return df
8. Example A — Validate a 20-Item Scale (EFA + Alpha + CFA) import numpy as np
import pandas as pd
item_data = responses[item_cols]
print ("=== Step 1: CTT Item Statistics ===" )
item_stats = compute_item_stats(item_data)
print (item_stats.round (3 ).to_string())
low_citc = item_stats[item_stats["corrected_item_total_r" ] < 0.30 ]
print (f"\nItems with CITC < 0.30: {list (low_citc.index)} " )
print ("\n=== Step 2: Cronbach's Alpha ===" )
alpha_result = cronbach_alpha(item_data)
alpha_f1 = cronbach_alpha(item_data.iloc[:, :10 ])
alpha_f2 = cronbach_alpha(item_data.iloc[:, 10 :])
print (f" F1 subscale alpha: {alpha_f1['alpha' ]} " )
print (f" F2 subscale alpha: {alpha_f2['alpha' ]} " )
print ("\n=== Step 3: Exploratory Factor Analysis ===" )
np.random.seed(1 )
idx_efa = np.random.choice(len (item_data), size=len (item_data) // 2 , replace=False )
efa_data = item_data.iloc[idx_efa].reset_index(drop=True )
efa_result = run_efa(efa_data, n_factors=2 , rotation="oblimin" , parallel_analysis=True )
print ("\n=== Step 4: Confirmatory Factor Analysis ===" )
idx_cfa = np.setdiff1d(np.arange(len (item_data)), idx_efa)
cfa_data = item_data.iloc[idx_cfa].reset_index(drop=True )
f1_items = [c for c in item_cols[:10 ]]
f2_items = [c for c in item_cols[10 :]]
model_str = (
"F1 =~ " + " + " .join(f1_items) + "\n"
"F2 =~ " + " + " .join(f2_items)
)
cfa_result = run_cfa(cfa_data, model_str)
print ("\n=== Scale Validation Summary ===" )
print (f" Items: {len (item_cols)} " )
print (f" Factors identified (EFA): {efa_result['n_factors' ]} " )
print (f" Overall alpha: {alpha_result['alpha' ]:.3 f} " )
print (f" CFA CFI: {cfa_result['cfi' ]:.3 f} (≥0.95 = acceptable)" )
print (f" CFA RMSEA: {cfa_result['rmsea' ]:.3 f} (≤0.06 = close fit)" )
print (f" CFA SRMR: {cfa_result['srmr' ]:.3 f} (≤0.08 = acceptable)" )
9. Example B — Compare IRT Parameters Across Gender Groups for DIF import numpy as np
import pandas as pd
binary_items = (responses[item_cols] >= 3 ).astype(int )
print ("=== IRT 2PL: Full Sample ===" )
irt_full = run_irt_2pl(binary_items.values, n_iter=200 )
print ("\nItem Parameters (first 10 items):" )
print (irt_full["item_params" ].head(10 ).round (3 ).to_string())
female_mask = responses["group" ] == 0
male_mask = responses["group" ] == 1
print ("\n=== IRT 2PL: Female Group ===" )
irt_female = run_irt_2pl(binary_items[female_mask].values, n_iter=200 )
print ("\n=== IRT 2PL: Male Group ===" )
irt_male = run_irt_2pl(binary_items[male_mask].values, n_iter=200 )
b_comparison = pd.DataFrame({
"b_full" : irt_full["b" ],
"b_female" : irt_female["b" ],
"b_male" : irt_male["b" ],
"delta_b" : irt_female["b" ] - irt_male["b" ],
}, index=item_cols)
print ("\n=== Item Difficulty Comparison (Female - Male) ===" )
print (b_comparison.round (3 ).to_string())
large_dif_items = b_comparison[b_comparison["delta_b" ].abs () > 0.5 ]
print (f"\nItems with |delta_b| > 0.5 (potential DIF): {list (large_dif_items.index)} " )
print ("\n=== Mantel-Haenszel DIF Analysis ===" )
dif_result = detect_dif(
binary_items,
responses["group" ],
n_score_strata=5 ,
)
print (dif_result[dif_result["dif_flag" ] == "DIF" ].to_string(index=False )
if (dif_result["dif_flag" ] == "DIF" ).any () else "No items flagged for DIF." )
theta_range = np.linspace(-3 , 3 , 200 )
fig, ax = plt.subplots(figsize=(8 , 5 ))
for i, item in enumerate (item_cols[:4 ]):
j = item_cols.index(item)
a_j = irt_full["a" ][j]
b_j = irt_full["b" ][j]
icc = expit(a_j * (theta_range - b_j))
ax.plot(theta_range, icc, linewidth=2 ,
label=f"{item} (a={a_j:.2 f} , b={b_j:.2 f} )" )
ax.set_xlabel("Latent Trait (θ)" )
ax.set_ylabel("P(correct)" )
ax.set_title("Item Characteristic Curves — 2PL Model" , fontweight="bold" )
ax.legend(fontsize=8 )
ax.set_ylim(0 , 1 )
ax.axhline(0.5 , color="gray" , linestyle=":" , linewidth=0.8 )
plt.tight_layout()
plt.savefig("icc_2pl.png" , dpi=150 )
plt.show()
10. Measurement Invariance Testing def test_measurement_invariance (
data: pd.DataFrame,
group: pd.Series,
model_str: str ,
) -> pd.DataFrame:
"""
Test configural, metric, and scalar measurement invariance across groups.
Steps:
1. Configural: same factor structure, all params free.
2. Metric: loadings constrained equal across groups.
3. Scalar: loadings + intercepts constrained equal.
Compares chi-square difference tests and ΔCFI (Cheung & Rensvold, 2002).
Parameters
----------
data : pd.DataFrame
Full item response matrix.
group : pd.Series
Group indicator (integer-coded).
model_str : str
semopy CFA model specification.
Returns
-------
pd.DataFrame with rows for each invariance level.
"""
groups = group.unique()
results = []
for level in ["configural" , "metric" , "scalar" ]:
level_chi2, level_df, level_cfi = [], [], []
for g in groups:
mask = group == g
g_data = data[mask].reset_index(drop=True )
try :
fit = run_cfa(g_data, model_str)
level_chi2.append(fit["chi2" ])
level_df.append(fit["df" ])
level_cfi.append(fit["cfi" ])
except Exception as exc:
print (f"Warning: CFA failed for group {g} at level {level} : {exc} " )
if level_chi2:
results.append({
"level" : level,
"chi2" : round (sum (level_chi2), 3 ),
"df" : sum (level_df),
"mean_cfi" : round (np.mean(level_cfi), 4 ),
})
inv_df = pd.DataFrame(results)
if len (inv_df) >= 2 :
inv_df["delta_chi2" ] = inv_df["chi2" ].diff().fillna(0 )
inv_df["delta_df" ] = inv_df["df" ].diff().fillna(0 )
inv_df["delta_p" ] = inv_df.apply(
lambda r: 1 - stats.chi2.cdf(r["delta_chi2" ], df=r["delta_df" ])
if r["delta_df" ] > 0 else np.nan, axis=1
)
inv_df["delta_cfi" ] = inv_df["mean_cfi" ].diff().fillna(0 )
inv_df["decision" ] = inv_df["delta_cfi" ].apply(
lambda d: "Invariant (ΔCFI < 0.01)" if abs (d) < 0.01 else "Non-invariant"
)
print ("\n=== Measurement Invariance Tests ===" )
print (inv_df.round (4 ).to_string(index=False ))
return inv_df
11. Tips and Common Pitfalls
Ordinal vs continuous EFA : For Likert items with fewer than 5 categories or
substantial skew, use polychoric correlations instead of Pearson. The
factor_analyzer library supports is_corr_matrix=True; compute polychoric
correlations with pingouin or semopy.CFA sample size : Aim for N ≥ 200 for stable CFA estimates; N ≥ 10 per free
parameter is a common rule of thumb. The 2PL IRT model requires N ≥ 500 for
stable discrimination estimates.Alpha vs omega : Cronbach's alpha assumes tau-equivalence (equal loadings).
McDonald's omega (hierarchical) is preferable for multidimensional scales.
Compute omega with pingouin.reliability(data) or the factor_analyzer
get_factor_variance() output.IRT model identification : Fix the mean of theta to 0 and SD to 1 (or fix
one item's b to 0) to identify the 2PL model. The implementation above uses
z-scoring of theta at each EM step.DIF interpretation : MH category A (|ln(alpha_MH)| < ln(1.5)) = negligible;
B = moderate; C = large (ETS classification). Always supplement with subject-matter
review of DIF items before removal.Parallel analysis vs Kaiser : Kaiser criterion (eigenvalue > 1) systematically
over-factors. Parallel analysis is strongly preferred in the current literature. 