"""单组一阶CFA及递归SEM，直接估计协方差结构而非因子分数回归。""" import itertools import warnings import numpy as np import pandas as pd from scipy import stats from semopy import Model, calc_stats from flask_babel import gettext as _ from .contracts import Column, PaperTable, StatisticalResult, columns, numeric_frame, enough, number, boolean, varying from .measurement import dimension_sets from .plots import structural_plot def validate_paths(paths, nodes): if not isinstance(paths, list) or not paths or len(paths) > 100: raise ValueError(_('请设置至少一条有效的结构路径')) edges = [] for edge in paths: if not isinstance(edge, dict): raise ValueError(_('每条路径都需要起点和终点')) source, target = edge.get('from'), edge.get('to') if source not in nodes or target not in nodes or source == target: raise ValueError(_('路径节点不存在或包含自循环，请检查模型设置')) if (source, target) in edges: raise ValueError(_('模型中不能重复设置同一条路径')) edges.append((source, target)) # 拓扑排序同时检查有向循环，首版仅支持递归模型。 pending = set(nodes) ordered = [] while pending: roots = [name for name in nodes if name in pending and not any(target == name and source in pending for source, target in edges)] if not roots: raise ValueError(_('首版结构方程支持递归模型，路径中不能形成有向循环')) ordered.extend(roots) pending.difference_update(roots) participating = {name for edge in edges for name in edge} if participating != set(nodes): raise ValueError(_('部分节点未连接到任何结构路径，请补充路径或移除这些节点')) return edges, ordered def model_description(data, method, options): is_path = method == 'stat_path_analysis' if is_path: observed = columns(data, options.get('variables'), 2, 30) dimensions = [] nodes = list(observed) else: if method == 'stat_cfa_single': dimensions = [(_('潜变量1'), columns(data, options.get('variables'), 3, 50))] else: dimensions = dimension_sets(data, {**options, 'include_total': False}, minimum_items=3) if len(dimensions) < 2: raise ValueError(_('此方法需要至少两个维度，每个维度至少三个题项')) observed = [item for name, items in dimensions for item in items] nodes = [name for name, items in dimensions] if len(observed) > 100 or len(nodes) > 15: raise ValueError(_('首版模型最多支持15个潜变量和100个观测题项')) observed_aliases = {name: 'observed'+str(i) for i, name in enumerate(observed)} node_aliases = observed_aliases if is_path else {name: 'latent'+str(i) for i, name in enumerate(nodes)} inverse = {value: key for key, value in observed_aliases.items()} inverse.update({value: key for key, value in node_aliases.items()}) lines = [] for name, items in dimensions: lines.append(node_aliases[name]+' =~ '+' + '.join(observed_aliases[item] for item in items)) edges, order = ([], list(nodes)) if method in ('stat_path_analysis', 'stat_sem'): edges, order = validate_paths(options.get('paths'), nodes) for i, (source, target) in enumerate(edges): lines.append(node_aliases[target]+' ~ path'+str(i)+'*'+node_aliases[source]) roots = [name for name in nodes if not any(target == name for source, target in edges)] correlated = boolean(options, 'correlate_exogenous', True) for name in nodes: # 显式估计外生节点方差，防止固定样本协方差改变自由参数计数。 alias = node_aliases[name] lines.append(alias+' ~~ '+alias) for first, second in itertools.combinations(nodes, 2): coefficient = '' if first in roots and second in roots and correlated else '0*' lines.append(node_aliases[first]+' ~~ '+coefficient+node_aliases[second]) return {'description': '\n'.join(lines), 'observed': observed, 'dimensions': dimensions, 'nodes': nodes, 'edges': edges, 'order': order, 'observed_aliases': observed_aliases, 'node_aliases': node_aliases, 'inverse': inverse, 'roots': roots} def finite_estimate(value): try: number = float(value) return number if np.isfinite(number) else None except (TypeError, ValueError): return None def fit_description(fit): pieces = [] for name, key in [('χ²', 'chi2'), ('df', 'DoF'), ('CFI', 'CFI'), ('TLI', 'TLI'), ('RMSEA', 'RMSEA'), ('SRMR', 'SRMR')]: value = fit.get(key) if value is not None and np.isfinite(value): pieces.append(name+'='+('%d' % value if key == 'DoF' else '%.3f' % value)) return '；'.join(pieces) def effect_decomposition(model, spec, estimates, covariance, alpha): nodes = spec['nodes'] positions = {name: i for i, name in enumerate(nodes)} coefficients = np.zeros((len(nodes), len(nodes))) active_names = [name for name, parameter in model.parameters.items() if parameter.active] active_positions = {name: i for i, name in enumerate(active_names)} path_parameters = [] for i, (source, target) in enumerate(spec['edges']): subset = estimates[(estimates.op == '~') & (estimates.lval == spec['node_aliases'][target]) & (estimates.rval == spec['node_aliases'][source])] value = float(subset.iloc[0]['Estimate']) row, column = positions[target], positions[source] coefficients[row, column] = value path_parameters.append((row, column, active_positions['path'+str(i)])) inverse = np.linalg.inv(np.eye(len(nodes))-coefficients) total = inverse-np.eye(len(nodes)) rows = [] for source in nodes: reachable = {target for start, target in spec['edges'] if start == source} for step in range(len(nodes)): reachable.update(target for start, target in spec['edges'] if start in reachable) for target in nodes: i, j = positions[target], positions[source] if i == j or target not in reachable: continue gradient_total, gradient_direct = np.zeros(len(active_names)), np.zeros(len(active_names)) for row, col, parameter in path_parameters: gradient_total[parameter] = inverse[i, row]*inverse[col, j] gradient_direct[parameter] = 1 if row == i and col == j else 0 values = [('direct', _('直接效应'), coefficients[i, j], gradient_direct), ('indirect', _('间接效应'), total[i, j]-coefficients[i, j], gradient_total-gradient_direct), ('total', _('总效应'), total[i, j], gradient_total)] for key, label, estimate, gradient in values: if not np.any(gradient): continue variance = float(gradient @ covariance @ gradient) se = np.sqrt(max(variance, 0)) z = estimate/se if se > 0 else None half = stats.norm.ppf(1-alpha/2)*se rows.append({'path': source+' → '+target, 'effect': label, 'effect_code': key, 'estimate': estimate, 'se': se, 'z': z, 'p': 2*stats.norm.sf(abs(z)) if z is not None else None, 'ci': (estimate-half, estimate+half)}) return rows def analyze_structural(data, method, options): alpha = number(options, 'alpha', .05, .0001, .25) spec = model_description(data, method, options) sample = enough(numeric_frame(data, spec['observed']).dropna(), max(20, len(spec['observed'])+2)) for name in spec['observed']: varying(sample[name].values, name) correlation = sample.corr().to_numpy() if np.min(np.linalg.eigvalsh(correlation)) < 1e-10: raise ValueError(_('观测变量相关矩阵奇异，请检查重复指标及完全共线性')) renamed = sample.rename(columns=spec['observed_aliases']) model = Model(spec['description']) free_count = sum(parameter.active for parameter in model.parameters.values()) moments = len(sample.columns)*(len(sample.columns)+1)//2 if free_count > moments: raise ValueError(_('自由参数超过可用协方差矩，模型欠识别')) with warnings.catch_warnings(record=True) as warning_records: warnings.simplefilter('always') result = model.fit(renamed, obj='MLW', solver='SLSQP', options={'maxiter': 2000, 'ftol': 1e-10}) if not result.success or not np.all(np.isfinite(model.param_vals)): raise ValueError(_('结构方程模型未收敛，请检查量表结构、路径或数据尺度')) information = model.calc_fim() diagonal = np.diag(information) if np.any(diagonal <= 0): raise ValueError(_('模型信息矩阵无效，无法计算可靠的标准误')) normalized = information/np.sqrt(np.outer(diagonal, diagonal)) if np.min(np.linalg.eigvalsh(normalized)) <= 1e-9: raise ValueError(_('模型局部欠识别或估计不稳定，请减少冗余参数')) covariance = np.linalg.inv(information) estimates = model.inspect(std_est=True) fitted_stats = calc_stats(model).iloc[0].to_dict() variances = estimates[(estimates.op == '~~') & (estimates.lval == estimates.rval)] if (variances['Estimate'] <= 1e-9).any(): raise ValueError(_('模型出现接近零或负的误差方差，请检查不当解与题项结构')) sigma, ignored_auxiliary = model.calc_sigma() # SRMR以样本方差标准化协方差残差，包含唯一协方差元素。 scale = np.sqrt(np.outer(np.diag(model.mx_cov), np.diag(model.mx_cov))) standardized_residuals = (model.mx_cov-sigma)/scale lower = np.tril_indices_from(standardized_residuals) fitted_stats['SRMR'] = np.sqrt(np.mean(standardized_residuals[lower]**2)) fit = {key: finite_estimate(value) for key, value in fitted_stats.items()} if fit['DoF'] == 0: for key in ('RMSEA', 'TLI', 'AGFI', 'chi2 p-value'): fit[key] = None raw = {'n': len(sample), 'estimator': 'MLW', 'missing_policy': 'listwise', 'model_specification': spec['description'], 'alias_mapping': spec['inverse'], 'fit': fit, 'estimates': estimates, 'parameter_covariance': covariance, 'sample_covariance': model.mx_cov, 'implied_covariance': sigma, 'observed_order': model.vars['observed'], 'warnings': [str(w.message) for w in warning_records], 'free_parameters': free_count} measurements, validity, graph_edges = [], [], [] for factor, items in spec['dimensions']: factor_alias = spec['node_aliases'][factor] factor_rows = [] for item in items: item_alias = spec['observed_aliases'][item] estimate = estimates[(estimates.op == '~') & (estimates.lval == item_alias) & (estimates.rval == factor_alias)].iloc[0] loading = float(estimate['Est. Std']) if loading*loading >= 1-1e-8: raise ValueError(_('标准化因子载荷触及不当解边界，请检查模型')) row = {'dimension': factor, 'item': item, 'loading': loading, 'estimate': float(estimate['Estimate']), 'se': finite_estimate(estimate['Std. Err']), 'z': finite_estimate(estimate['z-value']), 'p': finite_estimate(estimate['p-value'])} factor_rows.append(row) graph_edges.append({'from': factor_alias, 'to': item_alias, 'label': '%.2f' % loading}) loadings = np.array([row['loading'] for row in factor_rows]) errors = 1-loadings**2 cr = np.sum(loadings)**2/(np.sum(loadings)**2+np.sum(errors)) ave = np.sum(loadings**2)/(np.sum(loadings**2)+np.sum(errors)) factor_rows[0].update(cr=cr, ave=ave) measurements.extend(factor_rows) validity.append({'factor': factor, 'cr': cr, 'ave': ave}) raw['measurement'] = measurements raw['convergent_validity'] = validity if spec['dimensions']: inner_inverse = np.linalg.inv(np.eye(len(model.mx_beta))-model.mx_beta) inner_covariance = inner_inverse @ model.mx_psi @ inner_inverse.T aliases = list(model.names_beta[0]) indexes = [aliases.index(spec['node_aliases'][name]) for name, items in spec['dimensions']] latent_covariance = inner_covariance[np.ix_(indexes, indexes)] latent_correlation = latent_covariance/np.sqrt(np.outer(np.diag(latent_covariance), np.diag(latent_covariance))) discriminant = latent_correlation.copy() np.fill_diagonal(discriminant, np.sqrt([row['ave'] for row in validity])) raw['factor_correlation'] = latent_correlation raw['fornell_larcker'] = {'factors': [name for name, items in spec['dimensions']], 'matrix': discriminant} note = _('N=%(n)s；%(fit)s。', n=len(sample), fit=fit_description(fit)) if method.startswith('stat_cfa_'): tables = [PaperTable(_('验证性因子分析'), [Column('dimension', _('维度'), 'text', merge=True), Column('item', _('题项'), 'text'), Column('loading', _('标准化载荷')), Column('estimate', _('非标准化载荷')), Column('se', _('标准误')), Column('z', 'z'), Column('p', 'p', 'p'), Column('cr', 'CR'), Column('ave', 'AVE')], measurements, note)] if boolean(options, 'show_discriminant_validity', False): factors = [name for name, items in spec['dimensions']] rows = [{'factor': name, **{'c'+str(j): discriminant[i, j] for j in range(len(factors))}} for i, name in enumerate(factors)] tables.append(PaperTable(_('区分效度'), [Column('factor', _('维度'), 'text')]+[Column('c'+str(j), name) for j, name in enumerate(factors)], rows, _('对角线为AVE平方根，其余为潜变量相关系数。'))) else: rows = [] for source, target in spec['edges']: estimate = estimates[(estimates.op == '~') & (estimates.lval == spec['node_aliases'][target]) & (estimates.rval == spec['node_aliases'][source])].iloc[0] b, se = float(estimate['Estimate']), float(estimate['Std. Err']) half = stats.norm.ppf(1-alpha/2)*se rows.append({'path': source+' → '+target, 'estimate': b, 'standardized': float(estimate['Est. Std']), 'se': se, 'z': float(estimate['z-value']), 'p': float(estimate['p-value']), 'ci': (b-half, b+half)}) graph_edges.append({'from': spec['node_aliases'][source], 'to': spec['node_aliases'][target], 'label': '%.2f' % float(estimate['Est. Std'])}) tables = [PaperTable(_('结构路径估计'), [Column('path', _('路径'), 'text'), Column('estimate', 'B'), Column('se', _('标准误')), Column('standardized', 'β'), Column('z', 'z'), Column('p', 'p', 'p'), Column('ci', _('%(level)s%%置信区间', level='%g' % (100*(1-alpha))), 'interval')], rows, note)] effects = effect_decomposition(model, spec, estimates, covariance, alpha) raw['effects'] = effects raw['effect_se_method'] = 'delta_method_full_parameter_covariance' if boolean(options, 'show_effects', True) and any(row['effect_code'] == 'indirect' for row in effects): tables.append(PaperTable(_('效应分解'), [Column('path', _('路径'), 'text'), Column('effect', _('效应'), 'text'), Column('estimate', _('效应值')), Column('se', _('标准误')), Column('p', 'p', 'p'), Column('ci', _('%(level)s%%置信区间', level='%g' % (100*(1-alpha))), 'interval')], effects, _('效应标准误与区间采用Delta法。'))) if measurements and boolean(options, 'show_measurement', False): tables.append(PaperTable(_('测量模型'), [Column('dimension', _('维度'), 'text'), Column('item', _('题项'), 'text'), Column('loading', _('标准化载荷')), Column('cr', 'CR'), Column('ave', 'AVE')], measurements)) latent_aliases = set(spec['node_aliases'].values()) if spec['dimensions'] else set() image = structural_plot(spec['inverse'], graph_edges, latent_aliases) return StatisticalResult(tables, raw, [image])