Reviewer 2 Criticized Your Shortened or Modified Survey Scale? Here Is How to Defend It

R

library(lavaan)
library(psych)   # psych::omega also needs the GPArotation package installed

# 1. Compute McDonald's Omega Total & Cronbach's Alpha (one-factor solution)
omega_res <- suppressWarnings(psych::omega(short_scale_df, nfactors = 1, plot = FALSE))
cat(sprintf("Alpha: %.3f | McDonald's Omega: %.3f\n", omega_res$alpha, omega_res$omega.tot))

# 2. Fit Robust 1-Factor Confirmatory Factor Analysis (MLR)
cfa_model <- paste0("F1 =~ ", paste(colnames(short_scale_df), collapse = " + "))
fit <- cfa(cfa_model, data = short_scale_df, std.lv = TRUE, estimator = "MLR")
print(fitMeasures(fit, c("cfi.robust", "tli.robust", "rmsea.robust", "srmr")))

Python

import numpy as np
import pandas as pd

def posthoc_short_form_audit(df: pd.DataFrame, max_iter: int = 500, tol: float = 1e-8) -> dict:
    """
    Computes Cronbach's alpha, McDonald's omega total (one-factor model fitted by
    iterated principal-axis factoring), and the largest absolute residual correlation
    (a model-based stand-in for Yen's Q3) for an abbreviated scale.
    """
    X = df.dropna().astype(float).to_numpy()
    n_obs, k_items = X.shape

    # 1. Cronbach's alpha: (K / (K - 1)) * (1 - sum(item_vars) / total_var)
    item_vars = X.var(axis=0, ddof=1)
    total_var = X.sum(axis=1).var(ddof=1)
    alpha = (k_items / (k_items - 1.0)) * (1.0 - item_vars.sum() / total_var)

    # 2. One-factor loadings by iterated principal-axis factoring:
    #    put the communalities on the diagonal, take the top eigenvector, repeat until stable
    R = np.corrcoef(X, rowvar=False)
    h2 = 1.0 - 1.0 / np.diag(np.linalg.inv(R))          # start: squared multiple correlations
    for _ in range(max_iter):
        R_red = R.copy()
        np.fill_diagonal(R_red, h2)
        vals, vecs = np.linalg.eigh(R_red)
        lam = vecs[:, -1] * np.sqrt(max(vals[-1], 0.0))
        lam = lam if lam.sum() >= 0 else -lam            # fix the arbitrary sign
        h2_new = lam ** 2
        done = np.max(np.abs(h2_new - h2)) < tol
        h2 = h2_new
        if done:
            break

    # McDonald's omega total = (sum(lambda))^2 / [ (sum(lambda))^2 + sum(uniqueness) ]
    omega_total = lam.sum() ** 2 / (lam.sum() ** 2 + (1.0 - h2).sum())

    # 3. Local independence: residual correlations = observed - model-implied (lambda * lambda^T)
    resid = R - np.outer(lam, lam)
    np.fill_diagonal(resid, 0.0)
    max_abs_resid = float(np.abs(resid).max())

    return {
        "N": n_obs,
        "K_items": k_items,
        "cronbach_alpha": round(float(alpha), 3),
        "mcdonalds_omega": round(float(omega_total), 3),
        "max_abs_residual": round(max_abs_resid, 3),
        "local_independence_pass": max_abs_resid < 0.20,
    }