执行Two-Way ANOVA的Tukey HSD检验时遇TypeError求助

问题原因 #

报错核心是value列包含字符串类型数据，导致计算均值时无法转换为数值。从报错信息看，value列中出现了AAAAAAAAA这类字符串，说明数据读取或转换过程中响应变量被错误识别为字符串类型；此外，bioinfokit的tukey_hsd函数在计算均值时，默认对整个DataFrame（包含字符串类型的Genotype和years列）求均值，会将字符串拼接，进一步引发类型错误。

解决方法 #

方法1：修复响应变量的数据类型 #

在数据melt后，强制将value列转换为数值类型，去除无效值：

import pandas as pd
import scipy.stats as stats
import statsmodels.api as sm
from statsmodels.formula.api import ols
from bioinfokit.analys import stat

d = pd.read_csv("https://reneshbedre.github.io/assets/posts/anova/twowayanova.txt", sep="\t")

d_melt = pd.melt(d, id_vars=['Genotype'], value_vars=['1_year', '2_year', '3_year'])
d_melt.columns = ['Genotype', 'years', 'value']

# 关键修复：转换value列为数值类型，剔除无效值
d_melt['value'] = pd.to_numeric(d_melt['value'], errors='coerce')
d_melt = d_melt.dropna(subset=['value'])

# 重新执行Tukey HSD
res = stat()
res.tukey_hsd(df=d_melt, res_var='value', xfac_var='Genotype', anova_model='value~C(Genotype)+C(years)+C(Genotype):C(years)')
print(res.tukey_summary)

方法2：改用statsmodels的Tukey HSD实现（更稳定） #

避免bioinfokit的潜在bug，直接用statsmodels的官方实现：

import pandas as pd
from statsmodels.stats.multicomp import pairwise_tukeyhsd

d = pd.read_csv("https://reneshbedre.github.io/assets/posts/anova/twowayanova.txt", sep="\t")
d_melt = pd.melt(d, id_vars=['Genotype'], value_vars=['1_year', '2_year', '3_year'])
d_melt.columns = ['Genotype', 'years', 'value']

# 转换value列为数值类型
d_melt['value'] = pd.to_numeric(d_melt['value'], errors='coerce')
d_melt = d_melt.dropna(subset=['value'])

# 对Genotype主效应做Tukey HSD
tukey_genotype = pairwise_tukeyhsd(endog=d_melt['value'], groups=d_melt['Genotype'], alpha=0.05)
print("Genotype主效应Tukey HSD结果：")
print(tukey_genotype.summary())

# 对years主效应做Tukey HSD
tukey_years = pairwise_tukeyhsd(endog=d_melt['value'], groups=d_melt['years'], alpha=0.05)
print("\nyears主效应Tukey HSD结果：")
print(tukey_years.summary())

# 对交互项做Tukey HSD（组合两个因子为分组）
d_melt['interaction_group'] = d_melt['Genotype'] + "_" + d_melt['years']
tukey_interaction = pairwise_tukeyhsd(endog=d_melt['value'], groups=d_melt['interaction_group'], alpha=0.05)
print("\n交互项Tukey HSD结果：")
print(tukey_interaction.summary())

说明 #

• 方法1针对bioinfokit的使用场景修复，确保响应变量为纯数值；
• 方法2使用更成熟的statsmodels工具，避免第三方库的潜在问题，同时支持主效应和交互项的多重比较。

内容的提问来源于stack exchange，提问作者Fulgence