双因素方差分析

比较多个独立组的均值，并在适用时分析因素交互。

计算口径

本入口为独立组方差分析，不是重复测量设计。需关注组内残差、方差齐性、各组合样本量与空单元格。

同时指定两个分组变量，模型包含两个主效应及交互。第Ⅱ/Ⅲ类平方和在不平衡设计中的问题不同。简单效应使用模型合并误差，并报告Holm校正。

总体F检验显著不代表每两组都不同。需要定位差异时选择相匹配的事后比较，并解释校正后的p值。

同源实现

以下片段来自 core/statistical/anova.py 的 twoway，由镜像脚本按语法树提取。共享辅助函数和分发逻辑包含在完整下载包中。

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)

查看完整模块 · 下载算法包 · 复现指南

复现本方法

在解压目录安装 requirements.txt 后执行。样例为固定种子的模拟数据，仅供验证；输出不得冒充真实研究结果。

import examples._bootstrap
from examples.statistics_cases import build_data, method_cases
from core.statistical.runner import analyze

method = "stat_anova_twoway"
options = dict(method_cases())[method]
result = analyze(build_data(), method, options)
print(result.tables[0].html())

页面操作与解读教程