"""量表信度、题项鉴别和因子探索；保留量表维度及样本口径。""" import numpy as np from scipy import stats from factor_analyzer import FactorAnalyzer from factor_analyzer.rotator import Rotator from flask_babel import gettext as _ from .contracts import Column, PaperTable, StatisticalResult, columns, numeric_frame, enough, number, choice, boolean, varying def dimension_sets(data, options, minimum_items=2): dimensions = options.get('dimensions') or [] if not dimensions: return [(_('总量表'), columns(data, options.get('variables'), minimum_items))] if not isinstance(dimensions, list) or len(dimensions) > 30: raise ValueError(_('维度设置格式错误或维度过多')) groups, assigned, labels = [], set(), set() for dimension in dimensions: if not isinstance(dimension, dict): raise ValueError(_('每个维度都需要名称和题项')) name = str(dimension.get('name') or '').strip() selected = columns(data, dimension.get('variables'), minimum_items) if not name or name in labels or len(name) > 80: raise ValueError(_('维度名称不能为空、重复或超过80个字符')) if assigned.intersection(selected): raise ValueError(_('同一个题项不能同时分配给多个维度')) assigned.update(selected) labels.add(name) groups.append((name, selected)) if boolean(options, 'include_total', False) and len(groups) > 1: groups.append((_('总量表'), [name for label, selected in groups for name in selected])) return groups def alpha_coefficient(matrix): p = matrix.shape[1] if p < 2: return None covariance = np.cov(matrix, rowvar=False, ddof=1) denominator = np.sum(covariance) if denominator <= np.finfo(float).eps: return None return p/(p-1)*(1-np.trace(covariance)/denominator) def item_statistics(matrix, selected): output = [] for i, name in enumerate(selected): values = matrix[:, i] remaining = np.delete(matrix, i, axis=1) rest_sum = remaining.sum(axis=1) citc = stats.pearsonr(values, rest_sum)[0] if np.std(values) > 0 and np.std(rest_sum) > 0 else None output.append({'item': name, 'mean': np.mean(values), 'sd': np.std(values, ddof=1), 'citc': citc, 'alpha_deleted': alpha_coefficient(remaining)}) return output def reliability(data, options, split=False): rows, item_rows, details = [], [], {} for label, selected in dimension_sets(data, options): frame = enough(numeric_frame(data, selected).dropna(), 3) values = frame.values for name in selected: varying(frame[name].values, name) coefficient = alpha_coefficient(values) if coefficient is None: raise ValueError(_('维度“%(name)s”的总分没有有效变异', name=label)) correlation = np.corrcoef(values, rowvar=False) count = len(selected) standardized = count/(count-1)*(1-count/np.sum(correlation)) if np.sum(correlation) > 0 else None row = {'dimension': label, 'n': len(frame), 'items': count, 'alpha': coefficient, 'standard_alpha': standardized} detail = {'n': len(frame), 'items': selected, 'alpha': coefficient, 'standardized_alpha': standardized, 'covariance': np.cov(values, rowvar=False), 'correlation': correlation} if split: split_method = choice(options, 'split_method', 'odd_even', ('odd_even', 'first_last')) indexes_a = np.arange(0, count, 2) if split_method == 'odd_even' else np.arange((count+1)//2) indexes_b = np.array([i for i in range(count) if i not in indexes_a]) first, second = values[:, indexes_a].sum(axis=1), values[:, indexes_b].sum(axis=1) varying(first, label) varying(second, label) r = float(stats.pearsonr(first, second)[0]) if r <= -1+1e-12: raise ValueError(_('两个半量表完全负相关，请检查反向计分')) proportion = len(indexes_a)/count equal_sb = 2*r/(1+r) # 不等长校正使用题项数比例，等长时退化为常规Spearman–Brown公式。 unequal_sb = 2*r/(r+np.sqrt(r*r+4*proportion*(1-proportion)*(1-r*r))) total_var = np.var(first+second, ddof=1) guttman = 2*(1-(np.var(first, ddof=1)+np.var(second, ddof=1))/total_var) row.update(first_items=len(indexes_a), second_items=len(indexes_b), half_r=r, equal_sb=equal_sb, unequal_sb=unequal_sb, guttman=guttman) detail.update(split_method=split_method, first_items=[selected[i] for i in indexes_a], second_items=[selected[i] for i in indexes_b], half_r=r, equal_sb=equal_sb, unequal_sb=unequal_sb, guttman=guttman) else: items = item_statistics(values, selected) item_rows.extend({'dimension': label, **item} for item in items) detail['item_statistics'] = items rows.append(row) details[label] = detail headers = [Column('dimension', _('量表或维度'), 'text'), Column('n', 'N', 'integer'), Column('items', _('题项数'), 'integer')] if split: headers += [Column('first_items', _('前半题项数'), 'integer'), Column('second_items', _('后半题项数'), 'integer'), Column('half_r', _('两半相关')), Column('equal_sb', _('等长Spearman–Brown')), Column('unequal_sb', _('不等长Spearman–Brown')), Column('guttman', _('Guttman分半系数'))] else: headers += [Column('alpha', "Cronbach's α"), Column('standard_alpha', _('标准化α'))] tables = [PaperTable(_('分半信度分析') if split else _('内部一致性信度分析'), headers, rows)] if not split and boolean(options, 'item_details', False): tables.append(PaperTable(_('题项信度明细'), [Column('dimension', _('维度'), 'text'), Column('item', _('题项'), 'text'), Column('mean', _('均值')), Column('sd', _('标准差')), Column('citc', 'CITC'), Column('alpha_deleted', _('删除该题项后的α'))], item_rows)) return StatisticalResult(tables, details) def item_analysis(data, options): fraction = number(options, 'extreme_fraction', .27, .1, .49) rows, details = [], {} for label, selected in dimension_sets(data, options): frame = enough(numeric_frame(data, selected).dropna(), 8) values = frame.values scores = values.sum(axis=1) low_cut = np.quantile(scores, fraction) high_cut = np.quantile(scores, 1-fraction) low, high = scores <= low_cut, scores >= high_cut if low_cut >= high_cut or np.any(low & high): raise ValueError(_('维度“%(name)s”的高低分组因并列分数重叠，请调整比例或题项', name=label)) if min(low.sum(), high.sum()) < 2: raise ValueError(_('高低分组各需要至少两条有效记录')) for i, item in enumerate(item_statistics(values, selected)): a, b = values[high, i], values[low, i] if np.var(a)+np.var(b) == 0: # 常数题项保留在项目分析表中，避免自动删除掩盖问题。 t, p = (0., 1.) if np.mean(a) == np.mean(b) else (None, None) else: outcome = stats.ttest_ind(a, b, equal_var=False) t, p = outcome.statistic, outcome.pvalue rows.append({'dimension': label, **item, 'low': (np.mean(b), np.std(b, ddof=1)), 'high': (np.mean(a), np.std(a, ddof=1)), 't': t, 'p': p}) details[label] = {'n': len(values), 'high_n': int(high.sum()), 'low_n': int(low.sum()), 'high_cutoff': high_cut, 'low_cutoff': low_cut, 'tie_policy': 'include_boundary_ties', 'test': 'Welch', 'item_order': selected} table = PaperTable(_('项目分析'), [Column('dimension', _('维度'), 'text', merge=True), Column('item', _('题项'), 'text'), Column('mean', _('均值')), Column('sd', _('标准差')), Column('low', _('低分组'), 'mean_sd'), Column('high', _('高分组'), 'mean_sd'), Column('t', 't'), Column('p', 'p', 'p'), Column('citc', 'CITC'), Column('alpha_deleted', _('删除后α'))], rows) return StatisticalResult([table], details) def sampling_adequacy(correlation, n): p = len(correlation) sign, logdet = np.linalg.slogdet(correlation) if sign <= 0 or np.min(np.linalg.eigvalsh(correlation)) <= 1e-10: raise ValueError(_('相关矩阵奇异或近似奇异，请检查重复题项、共线性及样本量')) inverse = np.linalg.inv(correlation) partial = -inverse/np.sqrt(np.outer(np.diag(inverse), np.diag(inverse))) r = correlation.copy() np.fill_diagonal(r, 0) np.fill_diagonal(partial, 0) square, partial_square = r*r, partial*partial denominator = np.sum(square)+np.sum(partial_square) kmo = np.sum(square)/denominator if denominator > 0 else None per_item_denominator = np.sum(square+partial_square, axis=0) per_item = np.divide(np.sum(square, axis=0), per_item_denominator, out=np.full(p, np.nan), where=per_item_denominator > 0) statistic = -(n-1-(2*p+5)/6)*logdet df = p*(p-1)/2 return {'kmo': kmo, 'kmo_by_item': per_item, 'bartlett_chi2': statistic, 'bartlett_df': df, 'bartlett_p': stats.chi2.sf(statistic, df)} def principal_axis(correlation, factors): """主轴因子法迭代共同度；不将PCA或其他提取法伪装为PAF。""" communalities = 1-1/np.diag(np.linalg.inv(correlation)) for iteration in range(1000): reduced = correlation.copy() np.fill_diagonal(reduced, communalities) eigenvalues, eigenvectors = np.linalg.eigh(reduced) indices = np.argsort(eigenvalues)[::-1][:factors] if np.min(eigenvalues[indices]) <= 0: raise ValueError(_('指定的公共因子数超过有效正特征根数量')) loadings = eigenvectors[:, indices]*np.sqrt(eigenvalues[indices]) updated = np.sum(loadings*loadings, axis=1) if np.max(updated) > 1.000001: raise ValueError(_('因子解出现大于1的共同度，请调整因子数或检查题项')) if np.max(abs(updated-communalities)) < 1e-7: return loadings, iteration+1 communalities = updated raise ValueError(_('主轴因子迭代未收敛，请调整因子数或数据')) def rotate_loadings(loadings, method): """多起点正交旋转避免对称载荷在鞍点停住，再求斜交目标。""" factors = loadings.shape[1] random = np.random.default_rng(60926) starts = [np.eye(factors)] for index in range(8): orthogonal, upper = np.linalg.qr(random.normal(size=(factors, factors))) starts.append(orthogonal) best, best_score = None, -np.inf for initial in starts: candidate = Rotator(method='varimax', normalize=True, max_iter=1000, tol=1e-8).fit_transform(loadings @ initial) norm = np.sqrt(np.sum(candidate**2, axis=1)) scaled = candidate/np.where(norm > 0, norm, 1)[:, None] score = np.sum(scaled**4)-np.sum(np.sum(scaled**2, axis=0)**2)/len(scaled) if score > best_score: best, best_score = candidate, score if method == 'varimax': return best, np.eye(factors) rotator = Rotator(method='promax', normalize=True, max_iter=1000, tol=1e-8) pattern = rotator.fit_transform(best) return pattern, rotator.phi_ def factor_analysis(data, options, pca=False): selected = columns(data, options.get('variables'), 3, 80) frame = enough(numeric_frame(data, selected).dropna(), len(selected)+2) for name in selected: varying(frame[name].values, name) correlation = frame.corr().to_numpy() adequacy = sampling_adequacy(correlation, len(frame)) eigenvalues, eigenvectors = np.linalg.eigh(correlation) order = np.argsort(eigenvalues)[::-1] eigenvalues, eigenvectors = eigenvalues[order], eigenvectors[:, order] factors = number(options, 'factors', 0, 0, len(selected), integer=True) selection = 'manual' if factors else 'eigenvalue_greater_than_one' factors = factors or max(1, int(np.sum(eigenvalues > 1))) extraction = 'pca' if pca else choice(options, 'extraction', 'paf', ('paf', 'ml')) if not pca and ((len(selected)-factors)**2-(len(selected)+factors))/2 < 0: raise ValueError(_('指定公共因子数导致模型欠识别，请减少因子数')) iterations = None if pca: loadings = eigenvectors[:, :factors]*np.sqrt(eigenvalues[:factors]) elif extraction == 'paf': loadings, iterations = principal_axis(correlation, factors) else: model = FactorAnalyzer(n_factors=factors, rotation=None, method='ml', is_corr_matrix=True, bounds=(.005, 1)) model.fit(correlation) loadings = model.loadings_ if np.any(model.get_uniquenesses() <= .00501): raise ValueError(_('最大似然因子解触及独特性边界，请检查模型与题项')) unrotated = loadings.copy() rotation = choice(options, 'rotation', 'varimax', ('none', 'varimax', 'promax')) phi = np.eye(factors) if rotation != 'none' and factors > 1: loadings, phi = rotate_loadings(loadings, rotation) # 统一因子方向，斜交因子相关矩阵随载荷同时变号。 signs = np.where(loadings[np.argmax(abs(loadings), axis=0), np.arange(factors)] < 0, -1, 1) loadings = loadings*signs phi = phi*np.outer(signs, signs) communalities = np.diag(loadings @ phi @ loadings.T) headers = [Column('item', _('题项'), 'text')] headers += [Column('factor'+str(i), _('成分%(n)s', n=i+1) if pca else _('因子%(n)s', n=i+1)) for i in range(factors)] headers.append(Column('communality', _('共同度'))) rows = [{'item': name, 'communality': communalities[i], **{'factor'+str(j): loadings[i, j] for j in range(factors)}} for i, name in enumerate(selected)] raw = {'n': len(frame), 'variables': selected, 'correlation': correlation, **adequacy, 'eigenvalues': eigenvalues, 'factor_count': factors, 'factor_selection': selection, 'extraction': extraction, 'rotation': rotation, 'unrotated_loadings': unrotated, 'pattern_matrix': loadings, 'structure_matrix': loadings @ phi, 'factor_correlation': phi, 'communalities': communalities, 'unrotated_explained_variance': np.sum(unrotated**2, axis=0)/len(selected), 'iterations': iterations} note = _('N=%(n)s；提取方法：%(method)s；旋转：%(rotation)s。', n=len(frame), method=extraction.upper(), rotation=rotation) table = PaperTable(_('主成分载荷矩阵') if pca else _('探索性因子载荷矩阵'), headers, rows, note) from .plots import scree_plot return StatisticalResult([table], raw, [scree_plot(eigenvalues, factors)]) def analyze_measurement(data, method, options): if method == 'stat_reliability_alpha': return reliability(data, options) if method == 'stat_reliability_split': return reliability(data, options, split=True) if method == 'stat_item_analysis': return item_analysis(data, options) return factor_analysis(data, options, pca=method == 'stat_pca')