"""横截面统计回归：显式编码参照组、共同样本与论文式结果。""" import warnings import numpy as np import pandas as pd import statsmodels.api as sm from scipy import stats from statsmodels.stats.outliers_influence import variance_inflation_factor from statsmodels.stats.stattools import durbin_watson, jarque_bera from statsmodels.tools.sm_exceptions import PerfectSeparationError from flask_babel import gettext as _ from .contracts import Column, PaperTable, StatisticalResult, columns, numeric_frame, enough, number, choice, boolean, varying from .descriptive import _same_code def unique_names(names): return list(dict.fromkeys(names)) def prepare_sample(data, outcome, predictors, options, binary=False): columns(data, [outcome]) predictors = columns(data, predictors) if outcome in predictors: raise ValueError(_('因变量不能同时作为自变量')) categorical = options.get('categorical') or [] if not isinstance(categorical, list) or any(name not in predictors for name in categorical): raise ValueError(_('分类变量必须来自已选择的解释变量')) numeric_names = [name for name in predictors if name not in categorical] if not binary: numeric_names += [outcome] sample = data[[outcome]+predictors].copy() if numeric_names: sample[numeric_names] = numeric_frame(data, numeric_names) sample = enough(sample.dropna(), 5) if binary: levels = list(pd.unique(sample[outcome])) if len(levels) != 2: raise ValueError(_('二元Logistic回归的因变量必须恰好有两个有效类别')) event = options.get('event', 1) matches = [value for value in levels if _same_code(value, event)] if len(matches) != 1: raise ValueError(_('事件编码未唯一匹配因变量的类别，请重新选择')) sample[outcome] = (sample[outcome] == matches[0]).astype(float) else: varying(sample[outcome].values, outcome) return sample def design_matrix(sample, predictors, options): """矩阵列使用内部名称，用户列名只作为展示元数据。""" categorical = options.get('categorical') or [] references = options.get('references') or {} if not isinstance(references, dict): raise ValueError(_('参照组配置必须是对象')) matrix = [np.ones(len(sample))] terms = [{'key': 'intercept', 'label': _('常数项'), 'variable': None, 'kind': 'constant'}] used_references = {} for name in predictors: if name in categorical: levels = list(pd.unique(sample[name])) if not 2 <= len(levels) <= 30: raise ValueError(_('分类变量“%(name)s”需要2至30个有效类别', name=name)) reference = levels[0] if name in references: matching = [value for value in levels if _same_code(value, references[name])] if len(matching) != 1: raise ValueError(_('变量“%(name)s”的参照组不在有效样本中', name=name)) reference = matching[0] used_references[name] = reference for value in levels: if value == reference: continue matrix.append((sample[name] == value).to_numpy(dtype=float)) terms.append({'key': 'v'+str(len(terms)), 'label': name+'='+str(value), 'variable': name, 'kind': 'category', 'level': value, 'reference': reference}) else: values = sample[name].to_numpy(dtype=float) varying(values, name) matrix.append(values) terms.append({'key': 'v'+str(len(terms)), 'label': name, 'variable': name, 'kind': 'continuous'}) return np.column_stack(matrix), terms, used_references def checked_design(matrix): if matrix.shape[0] <= matrix.shape[1]+1: raise ValueError(_('有效样本不足以估计当前模型，请减少解释变量或增加样本')) if np.linalg.matrix_rank(matrix) < matrix.shape[1]: raise ValueError(_('解释变量存在完全共线性，请检查重复变量、常数列或分类编码')) def fit_ols(y, matrix, covariance): checked_design(matrix) model = sm.OLS(y, matrix) result = model.fit() if covariance == 'standard' else model.fit(cov_type=covariance, use_t=True) if not np.all(np.isfinite(result.params)) or not np.all(np.isfinite(result.bse)): raise ValueError(_('模型估计不稳定，请检查数据尺度及共线性')) return result def coefficient_rows(fitted, matrix, terms, alpha, standardize=True): rows = [] intervals = fitted.conf_int(alpha=alpha) y_sd = np.std(fitted.model.endog, ddof=1) for i, term in enumerate(terms): beta = fitted.params[i]*np.std(matrix[:, i], ddof=1)/y_sd if i and standardize else None vif = variance_inflation_factor(matrix, i) if i and matrix.shape[1] > 2 else 1. if i else None rows.append({'term': term['label'], 'term_key': term['key'], 'b': fitted.params[i], 'se': fitted.bse[i], 'beta': beta, 'statistic': fitted.tvalues[i], 'p': fitted.pvalues[i], 'ci': tuple(intervals[i]), 'vif': vif}) return rows def model_details(result): jb, jb_p, skew, kurtosis = jarque_bera(result.resid) return {'n': int(result.nobs), 'df_model': result.df_model, 'df_resid': result.df_resid, 'r_squared': result.rsquared, 'adjusted_r_squared': result.rsquared_adj, 'f': result.fvalue, 'f_p': result.f_pvalue, 'aic': result.aic, 'bic': result.bic, 'residual_sd': np.sqrt(result.mse_resid), 'durbin_watson': durbin_watson(result.resid), 'residual_normality': {'jarque_bera': jb, 'p': jb_p, 'skewness': skew, 'kurtosis': kurtosis}} def regression_headers(alpha, vif=True): head = [Column('term', _('变量'), 'text'), Column('b', 'B'), Column('se', _('标准误')), Column('beta', 'β'), Column('statistic', 't'), Column('p', 'p', 'p'), Column('ci', _('%(level)s%%置信区间', level='%g' % (100*(1-alpha))), 'interval')] if vif: head.append(Column('vif', 'VIF')) return head def linear_regression(data, options, hierarchical=False, simple=False): outcome = options.get('outcome') controls = options.get('controls') or [] if controls: columns(data, controls) alpha = number(options, 'alpha', .05, .0001, .25) covariance = choice(options, 'covariance', 'standard', ('standard', 'HC1', 'HC3')) if hierarchical: blocks = options.get('blocks') if not isinstance(blocks, list) or not 2 <= len(blocks) <= 10: raise ValueError(_('层级回归请设置2至10个依次加入的变量组')) blocks = [columns(data, block) for block in blocks] all_predictors = controls + [name for block in blocks for name in block] else: predictors = columns(data, options.get('predictors'), 1, 1 if simple else 100) if simple and controls: raise ValueError(_('简单线性回归只包含一个解释因素；加入控制变量请使用多元回归')) blocks = [predictors] all_predictors = controls+predictors if len(unique_names(all_predictors)) != len(all_predictors): raise ValueError(_('变量组之间或自变量与控制变量之间不能重复选择')) sample = prepare_sample(data, outcome, all_predictors, options) full_matrix, all_terms, references = design_matrix(sample, all_predictors, options) y = sample[outcome].values included, fitted_models, coefficient_sets, diagnostics = list(controls), [], [], [] for block in blocks: included.extend(block) positions = [i for i, term in enumerate(all_terms) if i == 0 or term['variable'] in included] matrix = full_matrix[:, positions] terms = [all_terms[i] for i in positions] fit = fit_ols(y, matrix, covariance) coeffs = coefficient_rows(fit, matrix, terms, alpha) detail = model_details(fit) detail.update(variables=list(included), coefficients=coeffs) if fitted_models: previous = fitted_models[-1] detail['r_squared_change'] = fit.rsquared-previous.rsquared if covariance == 'standard': f_change, p_change, df_change = fit.compare_f_test(previous) else: new_positions = [j for j, term in enumerate(terms) if term['variable'] in block] restriction = np.eye(matrix.shape[1])[new_positions] test = fit.f_test(restriction) f_change, p_change, df_change = float(test.fvalue), float(test.pvalue), len(new_positions) detail.update(change_f=f_change, change_p=p_change, change_df=df_change, change_test='partial_F' if covariance == 'standard' else 'robust_Wald_F') else: detail.update(r_squared_change=fit.rsquared, change_f=fit.fvalue, change_p=fit.f_pvalue, change_df=fit.df_model) fitted_models.append(fit) coefficient_sets.append(coeffs) diagnostics.append(detail) if not hierarchical: note = _('N=%(n)s；R²=%(r2)s；调整R²=%(adjusted)s。', n=len(sample), r2='%.3f' % fit.rsquared, adjusted='%.3f' % fit.rsquared_adj) table = PaperTable(_('线性回归结果'), regression_headers(alpha), coefficient_sets[0], note) else: head = [Column('term', _('变量'), 'text')] rows = [{'term': term['label'], 'term_key': term['key']} for term in all_terms] for i, coeffs in enumerate(coefficient_sets): label = _('模型%(n)s', n=i+1) head += [Column('b'+str(i), 'B (SE)', 'coefficient', label), Column('beta'+str(i), 'β', 'number', label), Column('p'+str(i), 'p', 'p', label)] lookup = {row['term_key']: row for row in coeffs} for row in rows: coefficient = lookup.get(row['term_key']) if coefficient: row.update({'b'+str(i): (coefficient['b'], coefficient['se']), 'beta'+str(i): coefficient['beta'], 'p'+str(i): coefficient['p']}) for key, label in [('n', 'N'), ('r_squared', 'R²'), ('adjusted_r_squared', _('调整R²')), ('r_squared_change', 'ΔR²'), ('f', 'F'), ('change_f', _('增量F'))]: row = {'term': label} for i, detail in enumerate(diagnostics): row['b'+str(i)] = detail[key] if key in ('f', 'change_f'): row['p'+str(i)] = detail['f_p' if key == 'f' else 'change_p'] rows.append(row) table = PaperTable(_('层级回归结果'), head, rows) return StatisticalResult([table], {'outcome': outcome, 'models': diagnostics, 'references': references, 'covariance': covariance, 'missing_policy': 'common_complete_sample'}) def fit_logistic(y, matrix): checked_design(matrix) try: with warnings.catch_warnings(record=True) as recorded: warnings.simplefilter('always') fitted = sm.Logit(y, matrix).fit(disp=False, maxiter=200) if not fitted.mle_retvals.get('converged') or not np.all(np.isfinite(fitted.bse)): raise ValueError(_('Logistic模型未稳定收敛，请检查完全分离、稀疏类别或共线性')) information = -fitted.model.hessian(fitted.params) if np.min(np.linalg.eigvalsh(information)) <= 1e-10 or np.max(abs(fitted.params)) > 500: raise ValueError(_('Logistic估计接近分离或数值奇异，不能报告可靠结果')) if any('separation' in str(w.message).lower() for w in recorded): raise ValueError(_('数据存在完全或近完全分离，普通Logistic估计不适用')) return fitted except (PerfectSeparationError, np.linalg.LinAlgError): raise ValueError(_('Logistic模型存在完全分离或奇异信息矩阵，请调整变量')) def logistic_regression(data, options, univariate=False): predictors = columns(data, options.get('predictors')) controls = options.get('controls') or [] if univariate and controls: raise ValueError(_('单因素Logistic分析不同时加入控制变量')) if controls: columns(data, controls) all_names = predictors+controls if len(all_names) != len(set(all_names)): raise ValueError(_('自变量与控制变量不能重复')) categorical = options.get('categorical') or [] if not isinstance(categorical, list) or any(name not in all_names for name in categorical): raise ValueError(_('分类变量必须来自已选择的解释变量')) alpha = number(options, 'alpha', .05, .0001, .25) models = [[name] for name in predictors] if univariate else [all_names] rows, diagnostics = [], [] for names in models: model_options = {**options, 'categorical': [name for name in options.get('categorical', []) if name in names]} sample = prepare_sample(data, options.get('outcome'), names, model_options, binary=True) matrix, terms, references = design_matrix(sample, names, model_options) y = sample[options['outcome']].values fitted = fit_logistic(y, matrix) intervals = fitted.conf_int(alpha=alpha) if np.max(abs(intervals)) > 700: raise ValueError(_('Logistic置信区间过大，估计不稳定，请检查稀疏类别或分离')) for i, term in enumerate(terms): if univariate and i == 0: continue rows.append({'term': term['label'], 'n': len(sample), 'b': fitted.params[i], 'se': fitted.bse[i], 'wald': fitted.tvalues[i]**2, 'p': fitted.pvalues[i], 'or': np.exp(fitted.params[i]), 'ci': tuple(np.exp(intervals[i]))}) joint_tests = {} for name in names: indexes = [i for i, term in enumerate(terms) if term['variable'] == name] wald = fitted.wald_test(np.eye(matrix.shape[1])[indexes], scalar=True) joint_tests[name] = {'wald_chi2': float(wald.statistic), 'df': len(indexes), 'p': float(wald.pvalue)} predicted = fitted.predict(matrix) diagnostics.append({'variables': names, 'n': len(sample), 'events': int(y.sum()), 'references': references, 'event': options.get('event', 1), 'coefficients': fitted.params, 'covariance': fitted.cov_params(), 'log_likelihood': fitted.llf, 'lr_chi2': fitted.llr, 'lr_p': fitted.llr_pvalue, 'pseudo_r2_mcfadden': fitted.prsquared, 'aic': fitted.aic, 'bic': fitted.bic, 'factor_tests': joint_tests, 'classification_threshold': .5, 'accuracy': np.mean((predicted >= .5) == y)}) head = [Column('term', _('变量'), 'text'), Column('n', 'N', 'integer'), Column('b', 'B'), Column('se', _('标准误')), Column('wald', 'Wald χ²'), Column('p', 'p', 'p'), Column('or', 'OR'), Column('ci', _('OR的%(level)s%%置信区间', level='%g' % (100*(1-alpha))), 'interval')] return StatisticalResult([PaperTable(_('单因素Logistic回归') if univariate else _('多因素Logistic回归'), head, rows)], {'models': diagnostics, 'missing_policy': 'per_factor_complete_sample' if univariate else 'common_complete_sample'}) def analyze_regression(data, method, options): if method in ('stat_linear_simple', 'stat_linear_multiple', 'stat_linear_hierarchical'): return linear_regression(data, options, hierarchical=method == 'stat_linear_hierarchical', simple=method == 'stat_linear_simple') if method in ('stat_logistic_univariate', 'stat_logistic_multivariate'): return logistic_regression(data, options, univariate=method == 'stat_logistic_univariate') if method.startswith('stat_mediation_'): from .mechanisms import mediation return mediation(data, method, options) from .mechanisms import moderation return moderation(data, method, options)