"""独立组单因素/双因素方差分析及与模型一致的事后比较。""" import itertools import numpy as np import pandas as pd from scipy import stats import statsmodels.api as sm from statsmodels.formula.api import ols from statsmodels.stats.multicomp import pairwise_tukeyhsd 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 from .comparisons import grouped_samples, mean_sd def posthoc_comparisons(names, samples, method, alpha): rows = [] k = len(samples) pairs = k*(k-1)//2 sizes = np.array([len(x) for x in samples], dtype=float) means = np.array([np.mean(x) for x in samples]) variances = np.array([np.var(x, ddof=1) for x in samples]) residual_df = sum(sizes)-k mse = np.sum((sizes-1)*variances)/residual_df tukey = None if method == 'tukey': tukey = pairwise_tukeyhsd(np.concatenate(samples), np.repeat(np.arange(k), sizes.astype(int)), alpha=alpha) for position, (i, j) in enumerate(itertools.combinations(range(k), 2)): difference = means[i]-means[j] if method == 'tukey': p = float(tukey.pvalues[position]) ci = (-tukey.confint[position][1], -tukey.confint[position][0]) se = np.sqrt(mse*(1/sizes[i]+1/sizes[j])) df = residual_df elif method == 'bonferroni': se = np.sqrt(mse*(1/sizes[i]+1/sizes[j])) df = residual_df p = min(1., 2*stats.t.sf(abs(difference)/se, df)*pairs) half = stats.t.ppf(1-alpha/(2*pairs), df)*se ci = (difference-half, difference+half) elif method == 'games_howell': a, b = variances[i]/sizes[i], variances[j]/sizes[j] se = np.sqrt(a+b) if se <= 0: raise ValueError(_('所比较的两组都没有组内变异，无法进行Games–Howell比较')) df = (a+b)**2/(a*a/(sizes[i]-1)+b*b/(sizes[j]-1)) p = stats.studentized_range.sf(abs(difference)*np.sqrt(2)/se, k, df) half = stats.studentized_range.ppf(1-alpha, k, df)*se/np.sqrt(2) ci = (difference-half, difference+half) else: raise ValueError(_('不支持的事后比较方法')) rows.append({'first': names[i], 'second': names[j], 'difference': difference, 'se': se, 'df': df, 'p': p, 'ci': ci}) return rows def oneway(data, options): selected = columns(data, options.get('variables')) group_name = columns(data, [options.get('group')])[0] alpha = number(options, 'alpha', .05, .0001, .25) method = choice(options, 'anova_method', 'classic', ('classic', 'welch')) posthoc = choice(options, 'posthoc', 'none', ('none', 'tukey', 'bonferroni', 'games_howell')) if method == 'welch' and posthoc in ('tukey', 'bonferroni'): raise ValueError(_('Welch方差分析请配合Games–Howell比较，或不进行事后比较')) rows, comparisons, details = [], [], {} all_groups = list(dict.fromkeys(str(v) for v in data[group_name].dropna())) head = [Column('variable', _('变量'), 'text')] for i, label in enumerate(all_groups): head.extend([Column('n'+str(i), 'N', 'integer', label), Column('s'+str(i), _('均值±标准差'), 'mean_sd', label)]) head += [Column('f', 'F'), Column('df1', _('组间df')), Column('df2', _('误差df')), Column('p', 'p', 'p'), Column('eta2', 'η²')] for variable in selected: names, groups = grouped_samples(data, variable, group_name) groups = [enough(x, 2) for x in groups] n = np.array([len(x) for x in groups], dtype=float) means = np.array([np.mean(x) for x in groups]) variances = np.array([np.var(x, ddof=1) for x in groups]) k = len(groups) grand_mean = np.sum(n*means)/sum(n) between = np.sum(n*(means-grand_mean)**2) within = np.sum((n-1)*variances) if within <= 0: raise ValueError(_('各组均缺乏组内变异，无法进行方差分析')) if method == 'classic': df1, df2 = k-1, sum(n)-k f = between/df1/(within/df2) else: if np.any(variances <= 0): raise ValueError(_('Welch方差分析要求各组存在有效组内变异')) weights = n/variances weight_sum = sum(weights) weighted_mean = sum(weights*means)/weight_sum correction = sum((1-weights/weight_sum)**2/(n-1)) df1, df2 = k-1, (k*k-1)/(3*correction) f = sum(weights*(means-weighted_mean)**2)/(k-1)/(1+2*(k-2)*correction/(k*k-1)) p = stats.f.sf(f, df1, df2) eta2 = between/(between+within) row = {'variable': variable, 'f': f, 'df1': df1, 'df2': df2, 'p': p, 'eta2': eta2} for label, values in zip(names, groups): position = all_groups.index(label) row['n'+str(position)], row['s'+str(position)] = len(values), mean_sd(values) rows.append(row) levene = stats.levene(*groups, center='median') details[variable] = {'ss_between': between, 'ss_within': within, 'method': method, 'levene': {'statistic': levene.statistic, 'p': levene.pvalue}, 'n': n} if posthoc != 'none': comparisons.extend({'variable': variable, **row} for row in posthoc_comparisons(names, groups, posthoc, alpha)) tables = [PaperTable(_('单因素方差分析'), head, rows)] if comparisons: tables.append(PaperTable(_('事后多重比较'), [Column('variable', _('变量'), 'text'), Column('first', _('组1'), 'text'), Column('second', _('组2'), 'text'), Column('difference', _('均值差')), Column('se', _('标准误')), Column('p', _('校正p值'), 'p'), Column('ci', _('%(level)s%%置信区间', level='%g' % (100*(1-alpha))), 'interval')], comparisons, _('比较方法：%(method)s。', method=posthoc))) details['posthoc'] = comparisons return StatisticalResult(tables, details) def twoway(data, options): selected = columns(data, options.get('variables')) factors = columns(data, [options.get('group'), options.get('second_group')], 2, 2) if set(selected).intersection(factors): raise ValueError(_('因变量不能同时作为分组因素')) ss_type = number(options, 'ss_type', 3, 2, 3, integer=True) simple = boolean(options, 'simple_effects', False) alpha = number(options, 'alpha', .05, .0001, .25) rows, details, simple_rows = [], {}, [] for name in selected: # 使用内部列别名，用户的中文列名及特殊字符不会成为公式代码。 numeric = numeric_frame(data, [name]) frame = pd.DataFrame({'response': numeric[name], 'factor_a': data[factors[0]], 'factor_b': data[factors[1]]}).dropna() enough(frame, 5) levels = [list(pd.unique(frame[key])) for key in ('factor_a', 'factor_b')] if min(map(len, levels)) < 2: raise ValueError(_('双因素方差分析的每个因素都需要至少两个有效水平')) cell_counts = pd.crosstab(frame.factor_a, frame.factor_b) if (cell_counts == 0).any().any(): raise ValueError(_('因素组合存在空单元格，请调整分组或分析范围')) fitted = ols('response ~ C(factor_a, Sum) * C(factor_b, Sum)', frame).fit() if fitted.df_resid <= 0 or fitted.mse_resid <= 0 or np.linalg.matrix_rank(fitted.model.exog) < fitted.model.exog.shape[1]: raise ValueError(_('模型无法估计独立的误差变异，请增加每个组合的有效观测')) result = sm.stats.anova_lm(fitted, typ=ss_type) names = {'C(factor_a, Sum)': factors[0], 'C(factor_b, Sum)': factors[1], 'C(factor_a, Sum):C(factor_b, Sum)': factors[0]+' × '+factors[1], 'Residual': _('误差')} for term, label in names.items(): item = result.loc[term] rows.append({'variable': name, 'source': label, 'ss': item['sum_sq'], 'df': item['df'], 'ms': item['sum_sq']/item['df'], 'f': item.get('F'), 'p': item.get('PR(>F)'), 'partial_eta2': item['sum_sq']/(item['sum_sq']+fitted.ssr) if term != 'Residual' else None}) groups = [part.response.values for key, part in frame.groupby(['factor_a', 'factor_b'], observed=True)] levene = stats.levene(*groups, center='median') details[name] = {'n': len(frame), 'ss_type': ss_type, 'anova': result, 'cell_descriptives': frame.groupby(['factor_a', 'factor_b'], observed=True).response.agg(['count', 'mean', 'std']).reset_index(), 'levene': {'statistic': levene.statistic, 'p': levene.pvalue}} if simple: # 简单效应沿用完整双因素模型的误差均方与自由度。 family = [] for varying_factor, fixed_factor in [('factor_a', 'factor_b'), ('factor_b', 'factor_a')]: for fixed_level, part in frame.groupby(fixed_factor, sort=False, observed=True): split = [g.response.values for key, g in part.groupby(varying_factor, sort=False, observed=True)] sizes, means = np.array([len(v) for v in split]), np.array([np.mean(v) for v in split]) average = np.average(means, weights=sizes) ss = np.sum(sizes*(means-average)**2) df = len(split)-1 f = (ss/df)/fitted.mse_resid family.append({'variable': name, 'factor': factors[0 if varying_factor == 'factor_a' else 1], 'condition': str(factors[1 if fixed_factor == 'factor_b' else 0])+'='+str(fixed_level), 'f': f, 'df1': df, 'df2': fitted.df_resid, 'p': stats.f.sf(f, df, fitted.df_resid)}) adjusted = multipletests([row['p'] for row in family], alpha=alpha, method='holm')[1] for row, p in zip(family, adjusted): row['adjusted_p'] = p simple_rows.extend(family) tables = [PaperTable(_('双因素方差分析'), [Column('variable', _('因变量'), 'text', merge=True), Column('source', _('变异来源'), 'text'), Column('ss', _('平方和')), Column('df', 'df'), Column('ms', _('均方')), Column('f', 'F'), Column('p', 'p', 'p'), Column('partial_eta2', _('偏η²'))], rows, _('采用第%(type)s类平方和。', type=ss_type))] if simple_rows: tables.append(PaperTable(_('简单效应检验'), [Column('variable', _('因变量'), 'text'), Column('factor', _('因素'), 'text'), Column('condition', _('条件'), 'text'), Column('f', 'F'), Column('df1', _('效应df')), Column('df2', _('误差df')), Column('adjusted_p', _('校正p值'), 'p')], simple_rows, _('采用完整模型的误差项及Holm校正。'))) details['simple_effects'] = simple_rows return StatisticalResult(tables, details)