"""均值及秩检验，使用与研究设计匹配的比较和效应量。""" import numpy as np import pandas as pd from scipy import stats from statsmodels.stats.multitest import multipletests from flask_babel import gettext as _ from .contracts import Column, PaperTable, StatisticalResult, columns, numeric_frame, enough, number, choice, boolean def grouped_samples(data, variable, group, minimum_groups=2): columns(data, [group]) if variable == group: raise ValueError(_('分组变量不能与分析变量相同')) numeric = numeric_frame(data, [variable]) valid = numeric[variable].notna() & data[group].notna() values, names = [], [] for key in pd.unique(data.loc[valid, group]): mask = valid & (data[group] == key) names.append(str(key)) values.append(numeric.loc[mask, variable].to_numpy()) if len(values) < minimum_groups or len(values) > 50: raise ValueError(_('有效分组数不满足此方法的要求')) return names, values def mean_sd(values): return (float(np.mean(values)), float(np.std(values, ddof=1))) def median_iqr(values): return tuple(np.quantile(values, [.5, .25, .75])) def mean_interval(estimate, se, df, alpha, alternative): if alternative == 'less': return (-np.inf, estimate+stats.t.ppf(1-alpha, df)*se) if alternative == 'greater': return (estimate-stats.t.ppf(1-alpha, df)*se, np.inf) width = stats.t.ppf(1-alpha/2, df)*se return (estimate-width, estimate+width) def paired_columns(data, options): pairs = options.get('pairs') if not isinstance(pairs, list) or not pairs or len(pairs) > 50: raise ValueError(_('请至少设置一对配对变量')) return [columns(data, pair, 2, 2) for pair in pairs] def t_tests(data, options, method): alpha = number(options, 'alpha', .05, .0001, .25) alternative = choice(options, 'alternative', 'two-sided', ('two-sided', 'less', 'greater')) rows, details = [], {} headers = [Column('variable', _('变量'), 'text')] if method == 'stat_t_one_sample': selected = columns(data, options.get('variables')) numeric = numeric_frame(data, selected) test_value = number(options, 'test_value', 0) headers += [Column('n', 'N', 'integer'), Column('summary', _('均值±标准差'), 'mean_sd')] for name in selected: x = enough(numeric[name].dropna(), 2).to_numpy() n, sd = len(x), np.std(x, ddof=1) if sd <= 0: raise ValueError(_('变量“%(name)s”没有有效变异', name=name)) estimate, se, df = np.mean(x)-test_value, sd/np.sqrt(n), n-1 outcome = stats.ttest_1samp(x, test_value, alternative=alternative) rows.append({'variable': name, 'n': n, 'summary': mean_sd(x), 'difference': estimate, 't': outcome.statistic, 'df': df, 'p': outcome.pvalue, 'ci': mean_interval(estimate, se, df, alpha, alternative), 'effect': estimate/sd}) details['test_value'] = test_value effect_label = "Cohen's d" elif method == 'stat_t_independent': selected = columns(data, options.get('variables')) group = options.get('group') equal = boolean(options, 'equal_variance', False) headers += [Column('group_a', _('组1'), 'text'), Column('n_a', 'N', 'integer', _('组1')), Column('a', _('均值±标准差'), 'mean_sd', _('组1')), Column('group_b', _('组2'), 'text'), Column('n_b', 'N', 'integer', _('组2')), Column('b', _('均值±标准差'), 'mean_sd', _('组2'))] for name in selected: names, samples = grouped_samples(data, name, group) if len(samples) != 2: raise ValueError(_('独立样本t检验需要恰好两个有效组')) a, b = [enough(x, 2) for x in samples] na, nb = len(a), len(b) va, vb = np.var(a, ddof=1), np.var(b, ddof=1) pooled = ((na-1)*va+(nb-1)*vb)/(na+nb-2) if pooled <= 0: raise ValueError(_('两组都没有组内变异，无法进行t检验')) if equal: se, df = np.sqrt(pooled*(1/na+1/nb)), na+nb-2 else: se = np.sqrt(va/na+vb/nb) df = (va/na+vb/nb)**2/((va/na)**2/(na-1)+(vb/nb)**2/(nb-1)) estimate = np.mean(a)-np.mean(b) outcome = stats.ttest_ind(a, b, equal_var=equal, alternative=alternative) levene = stats.levene(a, b) details[name] = {'levene': {'statistic': levene.statistic, 'p': levene.pvalue}, 'groups': names} rows.append({'variable': name, 'group_a': names[0], 'group_b': names[1], 'n_a': na, 'n_b': nb, 'a': mean_sd(a), 'b': mean_sd(b), 'difference': estimate, 't': outcome.statistic, 'df': df, 'p': outcome.pvalue, 'ci': mean_interval(estimate, se, df, alpha, alternative), 'effect': estimate/np.sqrt(pooled)}) details['variance_assumption'] = 'equal' if equal else 'Welch' effect_label = "Cohen's d" else: pairs = paired_columns(data, options) headers += [Column('n', _('配对数'), 'integer'), Column('a', _('前项均值±标准差'), 'mean_sd'), Column('b', _('后项均值±标准差'), 'mean_sd')] for first, second in pairs: pair = enough(numeric_frame(data, [first, second]).dropna(), 2) a, b = pair[first].values, pair[second].values difference = a-b sd = np.std(difference, ddof=1) if sd <= 0: raise ValueError(_('配对差值没有变异，无法进行配对t检验')) estimate, se, df = np.mean(difference), sd/np.sqrt(len(pair)), len(pair)-1 outcome = stats.ttest_rel(a, b, alternative=alternative) rows.append({'variable': first+' − '+second, 'n': len(pair), 'a': mean_sd(a), 'b': mean_sd(b), 'difference': estimate, 't': outcome.statistic, 'df': df, 'p': outcome.pvalue, 'ci': mean_interval(estimate, se, df, alpha, alternative), 'effect': estimate/sd}) effect_label = "Cohen's dz" headers += [Column('difference', _('均值差')), Column('t', 't'), Column('df', 'df'), Column('p', 'p', 'p'), Column('ci', _('%(level)s%%置信区间', level='%g' % (100*(1-alpha))), 'interval'), Column('effect', effect_label)] details.update(alternative=alternative, alpha=alpha, results=rows) return StatisticalResult([PaperTable(_('t检验结果'), headers, rows)], details) def _dunn(groups): """Dunn比较使用合并样本秩与并列修正，而非重复两组秩检验。""" combined = np.concatenate(groups) ranks = stats.rankdata(combined) n = len(combined) unique_values, ties = np.unique(combined, return_counts=True) variance = n*(n+1)/12 - np.sum(ties**3-ties)/(12*(n-1)) if variance <= 0: raise ValueError(_('所有观测值相同，无法进行秩比较')) means, offset = [], 0 for group in groups: means.append(np.mean(ranks[offset:offset+len(group)])) offset += len(group) output = [] for i in range(len(groups)): for j in range(i+1, len(groups)): z = (means[i]-means[j])/np.sqrt(variance*(1/len(groups[i])+1/len(groups[j]))) output.append({'i': i, 'j': j, 'statistic': z, 'p': 2*stats.norm.sf(abs(z))}) return output def _wilcoxon(first, second, alternative): difference = first-second nonzero = difference[difference != 0] if not len(nonzero): return {'statistic': 0., 'p': 1., 'rank_biserial': 0., 'nonzero_n': 0} ranks = stats.rankdata(abs(nonzero)) positive, negative = ranks[nonzero > 0].sum(), ranks[nonzero < 0].sum() # 统一丢弃零差值后检验；并列或零差值存在时明确使用正态近似。 exact = len(nonzero) <= 50 and len(nonzero) == len(difference) and len(np.unique(abs(nonzero))) == len(nonzero) outcome = stats.wilcoxon(difference, zero_method='wilcox', correction=False, alternative=alternative, method='exact' if exact else 'approx') return {'statistic': outcome.statistic, 'p': outcome.pvalue, 'rank_biserial': (positive-negative)/(positive+negative), 'nonzero_n': len(nonzero), 'positive_rank_sum': positive, 'negative_rank_sum': negative, 'p_method': 'exact' if exact else 'normal_approximation'} def nonparametric(data, options, method): alternative = choice(options, 'alternative', 'two-sided', ('two-sided', 'less', 'greater')) followup = boolean(options, 'posthoc', False) correction = choice(options, 'adjustment', 'holm', ('holm', 'bonferroni')) rows, diagnostics, followups = [], {}, [] if method == 'stat_wilcoxon': for first, second in paired_columns(data, options): sample = enough(numeric_frame(data, [first, second]).dropna(), 3) outcome = _wilcoxon(sample[first].values, sample[second].values, alternative) rows.append({'variable': first+' − '+second, 'n': len(sample), 'group': _('配对样本'), 'summary': median_iqr(sample[first]-sample[second]), **outcome, 'effect': outcome['rank_biserial']}) diagnostics[first+' − '+second] = outcome headers = [Column('variable', _('变量对'), 'text'), Column('n', _('配对数'), 'integer'), Column('summary', _('差值中位数（四分位数）'), 'median_iqr'), Column('statistic', 'W'), Column('p', 'p', 'p'), Column('effect', _('秩二列相关'))] elif method == 'stat_friedman': selected = columns(data, options.get('variables'), 3, 50) sample = enough(numeric_frame(data, selected).dropna(), 3) outcome = stats.friedmanchisquare(*[sample[c] for c in selected]) if not np.isfinite(outcome.statistic): raise ValueError(_('数据缺乏有效的组内秩差异，无法进行Friedman检验')) kendall = outcome.statistic/(len(sample)*(len(selected)-1)) for i, name in enumerate(selected): rows.append({'variable': name, 'n': len(sample), 'summary': median_iqr(sample[name]), 'statistic': outcome.statistic if i == 0 else None, 'p': outcome.pvalue if i == 0 else None, 'effect': kendall if i == 0 else None}) if followup: for i, first in enumerate(selected): for second in selected[i+1:]: result = _wilcoxon(sample[first].values, sample[second].values, 'two-sided') followups.append({'variable': _('相关样本'), 'first': first, 'second': second, **result}) headers = [Column('variable', _('测量变量'), 'text'), Column('n', 'N', 'integer'), Column('summary', _('中位数（四分位数）'), 'median_iqr'), Column('statistic', 'χ²'), Column('p', 'p', 'p'), Column('effect', "Kendall's W")] diagnostics['df'] = len(selected)-1 else: selected = columns(data, options.get('variables')) for name in selected: names, groups = grouped_samples(data, name, options.get('group')) if method == 'stat_mann_whitney' and len(groups) != 2: raise ValueError(_('Mann–Whitney U检验需要两个有效组')) if method == 'stat_mann_whitney': a, b = groups outcome = stats.mannwhitneyu(a, b, alternative=alternative, method='auto') effect = 2*outcome.statistic/(len(a)*len(b))-1 else: if len(np.unique(np.concatenate(groups))) < 2: raise ValueError(_('所有观测值相同，无法进行Kruskal–Wallis检验')) outcome = stats.kruskal(*groups) effect = max(0., (outcome.statistic-len(groups)+1)/(sum(map(len, groups))-len(groups))) if sum(map(len, groups)) > len(groups) else None if followup: for comparison in _dunn(groups): followups.append({'variable': name, 'first': names[comparison['i']], 'second': names[comparison['j']], 'statistic': comparison['statistic'], 'p': comparison['p']}) for i, (label, group) in enumerate(zip(names, groups)): rows.append({'variable': name, 'group': label, 'n': len(group), 'summary': median_iqr(group), 'statistic': outcome.statistic if i == 0 else None, 'p': outcome.pvalue if i == 0 else None, 'effect': effect if i == 0 else None}) diagnostics[name] = {'statistic': outcome.statistic, 'p': outcome.pvalue, 'groups': names, 'n': list(map(len, groups)), 'effect': effect} headers = [Column('variable', _('变量'), 'text', merge=True), Column('group', _('组别'), 'text'), Column('n', 'N', 'integer'), Column('summary', _('中位数（四分位数）'), 'median_iqr'), Column('statistic', 'U' if method == 'stat_mann_whitney' else 'H'), Column('p', 'p', 'p'), Column('effect', _('秩二列相关') if method == 'stat_mann_whitney' else 'ε²')] tables = [PaperTable(_('非参数检验结果'), headers, rows)] if followups: # 每个结局的比较构成一个校正家族，保持不同分析变量相互独立。 for name in dict.fromkeys(row['variable'] for row in followups): family = [row for row in followups if row['variable'] == name] adjusted = multipletests([row['p'] for row in family], method=correction)[1] for row, value in zip(family, adjusted): row['adjusted_p'] = value tables.append(PaperTable(_('事后比较'), [Column('variable', _('变量'), 'text'), Column('first', _('比较组1'), 'text'), Column('second', _('比较组2'), 'text'), Column('statistic', _('统计量')), Column('adjusted_p', _('校正p值'), 'p')], followups, _('多重比较采用%(method)s校正。', method=correction))) diagnostics.update(alternative=alternative, posthoc=followups) return StatisticalResult(tables, diagnostics)