"""观测变量中介及调节：共同样本、个案重抽样与条件效应协方差。""" import numpy as np from scipy import stats from flask_babel import gettext as _ from .contracts import Column, PaperTable, StatisticalResult, columns, number, choice, boolean from .regression import prepare_sample, design_matrix, fit_ols, coefficient_rows, regression_headers, model_details, unique_names from .plots import simple_slopes_plot def mediation(data, method, options): predictor = columns(data, options.get('predictors'), 1, 1)[0] outcome = options.get('outcome') mediators = columns(data, options.get('mediators'), 1, 10) if method == 'stat_mediation_simple' and len(mediators) != 1: raise ValueError(_('简单中介需要一个中介变量')) if method == 'stat_mediation_parallel' and len(mediators) < 2: raise ValueError(_('平行中介至少需要两个中介变量')) if method == 'stat_mediation_serial' and len(mediators) != 2: raise ValueError(_('链式中介首版需要按先后顺序指定两个中介变量')) controls = options.get('controls') or [] if controls: columns(data, controls) names = [predictor]+mediators+controls if len(names) != len(set(names)): raise ValueError(_('自变量、中介变量和控制变量不能重复')) categorical = options.get('categorical') or [] if predictor in categorical or any(name in categorical for name in mediators): raise ValueError(_('首版中介分析的自变量和中介变量须按连续变量处理')) sample = prepare_sample(data, outcome, names, options) alpha = number(options, 'alpha', .05, .0001, .25) repetitions = number(options, 'bootstrap_samples', 1000, 200, 5000, integer=True) seed = number(options, 'seed', 20260926, 0, 4294967295, integer=True) # 同一个编码矩阵用于全部方程和重抽样，分类控制变量的参照组不会漂移。 matrix, terms, references = design_matrix(sample, names, options) position = {term['variable']: i for i, term in enumerate(terms) if term['kind'] == 'continuous'} control_positions = [i for i, term in enumerate(terms) if i == 0 or term['variable'] in controls] x_position = position[predictor] mediator_positions = [position[name] for name in mediators] equations = [('total', outcome, control_positions+[x_position])] for i, mediator in enumerate(mediators): previous = mediator_positions[:i] if method == 'stat_mediation_serial' else [] equations.append(('mediator'+str(i), mediator, control_positions+[x_position]+previous)) equations.append(('direct', outcome, control_positions+[x_position]+mediator_positions)) fits = {} diagnostics = {} path_rows = [] for key, response, indexes in equations: fit = fit_ols(sample[response].values, matrix[:, indexes], 'standard') fits[key] = fit coeffs = coefficient_rows(fit, matrix[:, indexes], [terms[i] for i in indexes], alpha) diagnostics[key] = {'outcome': response, **model_details(fit), 'coefficients': coeffs} path_rows.extend({'equation': response, **row} for row in coeffs) effect_names = [_('总效应'), _('直接效应'), _('总间接效应')] labels = [predictor+' → '+mediator+' → '+outcome for mediator in mediators] if method == 'stat_mediation_serial': labels.append(predictor+' → '+mediators[0]+' → '+mediators[1]+' → '+outcome) effect_names += labels def estimates(index=None): coefficients = {} for key, response, positions in equations: if index is None: beta = fits[key].params else: design = matrix[index][:, positions] if np.linalg.matrix_rank(design) < len(positions): raise np.linalg.LinAlgError('重抽样设计矩阵秩不足') beta = np.linalg.lstsq(design, sample[response].values[index], rcond=None)[0] coefficients[key] = dict(zip(positions, beta)) total = coefficients['total'][x_position] direct = coefficients['direct'][x_position] indirect = [coefficients['mediator'+str(i)][x_position]*coefficients['direct'][mediator_positions[i]] for i in range(len(mediators))] if method == 'stat_mediation_serial': indirect.append(coefficients['mediator0'][x_position]*coefficients['mediator1'][mediator_positions[0]]*coefficients['direct'][mediator_positions[1]]) return np.array([total, direct, sum(indirect)]+indirect) original = estimates() random = np.random.default_rng(seed) draws, failures = [], 0 progress = options.get('_progress') for iteration in range(repetitions): index = random.integers(0, len(sample), len(sample)) try: draw = estimates(index) if not np.all(np.isfinite(draw)): raise np.linalg.LinAlgError('非有限估计') draws.append(draw) except np.linalg.LinAlgError: failures += 1 if callable(progress) and (iteration % 20 == 0 or iteration+1 == repetitions): progress(iteration+1, repetitions) if len(draws) < max(180, repetitions*.9): raise ValueError(_('有效Bootstrap抽样不足90%，请检查样本量、稀疏类别或共线性')) bootstrap = np.asarray(draws) intervals = np.quantile(bootstrap, [alpha/2, 1-alpha/2], axis=0) standard_errors = np.std(bootstrap, axis=0, ddof=1) rows = [{'effect': name, 'estimate': original[i], 'se': standard_errors[i], 'ci': (intervals[0, i], intervals[1, i])} for i, name in enumerate(effect_names)] table = PaperTable(_('中介效应检验'), [Column('effect', _('效应或路径'), 'text'), Column('estimate', _('效应值')), Column('se', 'Boot SE'), Column('ci', _('Bootstrap %(level)s%%置信区间', level='%g' % (100*(1-alpha))), 'interval')], rows, _('N=%(n)s；个案Bootstrap有效抽样%(b)s次，百分位置信区间。', n=len(sample), b=len(draws))) tables = [table] if boolean(options, 'show_path_regressions', False): tables.append(PaperTable(_('路径回归'), [Column('equation', _('因变量'), 'text')]+regression_headers(alpha), path_rows)) return StatisticalResult(tables, {'sample_n': len(sample), 'outcome': outcome, 'predictor': predictor, 'mediators': mediators, 'controls': controls, 'references': references, 'equations': diagnostics, 'effects': rows, 'bootstrap_requested': repetitions, 'bootstrap_valid': len(draws), 'bootstrap_failed': failures, 'seed': seed, 'ci_method': 'percentile_case_bootstrap', 'decomposition_error': original[0]-original[1]-original[2]}) def moderation(data, method, options): predictor = columns(data, options.get('predictors'), 1, 1)[0] moderator = columns(data, [options.get('moderator')])[0] outcome = options.get('outcome') controls = options.get('controls') or [] if controls: columns(data, controls) names = [predictor, moderator]+controls if len(names) != len(set(names)): raise ValueError(_('自变量、调节变量和控制变量不能重复')) categorical = options.get('categorical') or [] if predictor in categorical: raise ValueError(_('首版调节分析的焦点自变量须按连续变量处理')) is_category = method == 'stat_moderation_categorical' if not is_category and moderator in categorical: raise ValueError(_('连续调节分析不能将调节变量按分类变量编码')) options = {**options, 'categorical': unique_names(categorical+([moderator] if is_category else []))} sample = prepare_sample(data, outcome, names, options) alpha = number(options, 'alpha', .05, .0001, .25) covariance = choice(options, 'covariance', 'standard', ('standard', 'HC1', 'HC3')) centered = boolean(options, 'center', True) predictor_center = float(sample[predictor].mean()) if centered else 0. moderator_center = float(sample[moderator].mean()) if centered and not is_category else 0. work = sample.copy() work[predictor] -= predictor_center if not is_category: work[moderator] -= moderator_center matrix, terms, references = design_matrix(work, names, options) x_index = next(i for i, term in enumerate(terms) if term['variable'] == predictor) w_indexes = [i for i, term in enumerate(terms) if term['variable'] == moderator] interaction_indexes = [] for i in w_indexes: interaction_indexes.append(matrix.shape[1]) matrix = np.column_stack([matrix, matrix[:, x_index]*matrix[:, i]]) terms.append({'key': 'interaction'+str(i), 'label': predictor+' × '+terms[i]['label'], 'variable': predictor+' × '+moderator, 'kind': 'interaction'}) fit = fit_ols(work[outcome].values, matrix, covariance) rows = coefficient_rows(fit, matrix, terms, alpha) restrictions = np.eye(matrix.shape[1])[interaction_indexes] joint = fit.f_test(restrictions) if is_category: values = [references[moderator]]+[terms[i]['level'] for i in w_indexes] else: values = options.get('moderator_values') if values is None or values == []: mean, sd = sample[moderator].mean(), sample[moderator].std(ddof=1) values = [mean-sd, mean, mean+sd] if not isinstance(values, (list, tuple)) or not 1 <= len(values) <= 8: raise ValueError(_('请设置1至8个调节变量取值')) try: values = [float(v) for v in values] except (TypeError, ValueError): raise ValueError(_('调节变量取值必须是数字')) if not np.all(np.isfinite(values)): raise ValueError(_('调节变量取值必须是有限数字')) conditional, lines = [], [] x_grid = np.linspace(sample[predictor].min(), sample[predictor].max(), 45) covariance_matrix = fit.cov_params() for value in values: contrast = np.zeros(matrix.shape[1]) contrast[x_index] = 1. if is_category: for index, interaction in zip(w_indexes, interaction_indexes): contrast[interaction] = 1. if value == terms[index]['level'] else 0. else: contrast[interaction_indexes[0]] = value-moderator_center estimate = float(contrast @ fit.params) variance = float(contrast @ covariance_matrix @ contrast) if variance <= 0: raise ValueError(_('条件效应的方差无法估计，请检查数据和模型')) se = np.sqrt(variance) t = estimate/se half = stats.t.ppf(1-alpha/2, fit.df_resid)*se label = moderator+'='+str(value) conditional.append({'term': label, 'b': estimate, 'se': se, 'statistic': t, 'p': 2*stats.t.sf(abs(t), fit.df_resid), 'ci': (estimate-half, estimate+half)}) plot_matrix = np.tile(matrix.mean(axis=0), (len(x_grid), 1)) plot_matrix[:, 0] = 1 plot_matrix[:, x_index] = x_grid-predictor_center for index, interaction in zip(w_indexes, interaction_indexes): plot_matrix[:, index] = (1. if value == terms[index]['level'] else 0.) if is_category else value-moderator_center plot_matrix[:, interaction] = plot_matrix[:, x_index]*plot_matrix[:, index] lines.append({'label': label, 'x': x_grid, 'y': plot_matrix @ fit.params}) # 同一张表将回归系数与条件效应并列呈现，复杂诊断留在原始输出。 for row in conditional: rows.append({**row, 'term': _('条件效应：%(value)s', value=row['term'])}) table = PaperTable(_('调节效应分析'), regression_headers(alpha), rows, _('N=%(n)s；R²=%(r2)s。', n=len(sample), r2='%.3f' % fit.rsquared)) return StatisticalResult([table], {'model': model_details(fit), 'conditional_effects': conditional, 'interaction_joint_f': float(joint.fvalue), 'interaction_joint_p': float(joint.pvalue), 'interaction_df': len(interaction_indexes), 'references': references, 'centering': {predictor: predictor_center, moderator: moderator_center}, 'covariance': covariance, 'parameter_covariance': covariance_matrix}, [simple_slopes_plot(lines, predictor, outcome)])