最近整理数据时又跟列联表打了几场交道发现很多人用 R语言做独立性检验翻来覆去只会chisq.test()一行代码。可实际数据一旦遇到小样本、稀疏单元格、分层混杂、因子有空白水平这行代码就会给出误导性结果。这篇笔记是 R语言笔记系列的第十篇专门把独立性检验函数从原理到实战捋一遍内容包括chisq.test()、fisher.test()、mantelhaen.test()、loglm()以及配套的效应量计算和批量检验小工具。无论你是在 RStudio、Positron 还是纯命令行环境里跑代码都能直接用。1. 先想清楚你手里的是原始记录还是汇总好的频数表1.1 chisq.test() 的两种等价写法独立性检验最常见的场景是判断两个分类变量之间是否存在关联。比如你想知道“用户选择哪个方案”和“是否最终购买”有没有关系数据结构通常有两种一种是原始观测记录一行为一个样本每列是一个变量。这种数据在 R语言里做检验十分直接set.seed(2024) df - data.frame( 方案 factor(sample(c(方案A, 方案B, 方案C), 200, replace TRUE)), 是否购买 factor(sample(c(购买, 未购买), 200, replace TRUE)) ) tb - with(df, table(方案, 是否购买)) tb chisq.test(tb)另一种是已经汇总好的频数矩阵。比如从数据库里直接group by两个字段后拿到了一个 2×3 或 3×3 的频数表。这时不需要还原原始数据直接把矩阵或table对象丢进去chisq.test(tb)很多人不知道chisq.test()还支持直接用两个因子向量做参数chisq.test(df$方案, df$是否购买)两种写法最终的统计量和 p 值一致区别在于内部实现传两个向量的写法会先调用table()生成列联表再走同样的检验过程。如果数据本来就以data.frame形式存在我会直接传两个因子向量如果是从论文表格里复制出来的频数就直接传矩阵。这里有一个很容易踩的坑chisq.test()的第一个参数如果是矩阵第二个参数就不能再传向量否则函数会把它当成“两个变量分别传”的接口结果跟你预期完全不一样。新手最容易在“已经拿到频数表、又额外传了一个向量”的时候翻车。1.2 预期频数不够卡方结果就不可信chisq.test()默认计算的是 Pearson 卡方统计量理论上统计量服从卡方分布的前提是“样本量足够大、单元格预期频数不能太小”。经验法则一般有两个门槛表中不能有预期频数小于 1 的单元格预期频数小于 5 的单元格占比不能超过 20%。这个“预期频数”不是表格里实际看到的频数而是在“两个变量独立”的假设下推算出来的理论频数。chisq.test()并不会阻止你计算但会在结果最后打印一行警告Warning message: In chisq.test(tb) : Chi-squared approximation may be incorrect看到这行警告第一反应不是去网上搜“怎么关掉warning”而是要判断自己的表是不是太稀疏。检查方法很简单expected_tb - chisq.test(tb)$expected expected_tb min(expected_tb) # 最小预期频数 mean(expected_tb 5) # 小于5的单元格占比如果违反条件优先考虑合并分类、增大样本量或者改用第 2 节要讲的fisher.test()而不是硬着头皮用 p 值汇报结论。2. 卡方统计量背后的计算以及什么情况必须换 Fisher 精确检验2.1 期望频数、自由度、连续性校正到底在算什么理解卡方独立性检验不用去背公式推导但至少要知道三个数字从哪来。第一个是预期频数。对于一个 I 行 J 列的列联表每个单元格的预期频数是E_ij (第 i 行合计 × 第 j 列合计) / 总样本数独立假设下行变量和列变量互不影响所以每个格子里的“理论人数”应该按行、列边际比例去分配。实际频数与预期频数差得越多说明“独立”这个假设越站不住脚。第二个是自由度。自由度不是(行数-1) × (列数-1)那么简单的一句话它反映的是在外围行、列合计都固定时表里可以自由变动的单元格数量。比如 2×2 表四个格子只要知道了一个剩下三个都会被行合计和列合计锁死所以自由度是 1。这个大原则决定了后续你在chisq.test()输出里看到的df参数。第三个是统计量本身X² Σ (观测频数 - 预期频数)² / 预期频数这个 X² 越大p 值越小。R语言里chisq.test()默认还有个correct TRUE参数这个参数只对 2×2 表生效意思是把每个|O - E|先减去 0.5 再做平方让统计量偏保守一点减少小样本下的第一类错误。很多老教程会建议把它关掉但我个人习惯是**事先不确定的时候不手动改尤其不要为了“让 p 值变小”而设成correct FALSE。**如果你面对的是大于 2×2 的表格这个参数设不设都一样不必纠结。2.2 fisher.test() 的适用场景与代价当预期频数条件不满足时R语言里最稳妥的替代方案就是fisher.test()。它在行列合计固定的条件下把所有可能出现的列联表全部枚举一遍逐张计算概率然后把“当前表以及比当前表更极端”的表概率累加得到精确 p 值。操作非常直接small_tb - matrix( c(8, 1, 2, 9), nrow 2, dimnames list( 处理 c(处理组, 对照组), 改善 c(改善, 无改善) ) ) fisher.test(small_tb)这组数据如果硬用卡方预期频数里会出现小于 5 的格子卡方近似可能失效Fisher 精确检验则没有这个问题。但代价也很明显枚举所有可能表意味着计算量会随表尺寸指数级增长。如果是大一点的表直接跑fisher.test()可能卡很久甚至内存不足。R语言里对应的参数是workspace默认 20 万实在跑不动可以调大也可以在大型表上使用蒙特卡洛模拟 p 值fisher.test(large_tb, simulate.p.value TRUE, B 10000)注意simulate.p.value TRUE得到的不再是精确 p 值而是基于随机抽样计算的近似 p 值。代码里B是模拟次数B10000是最低可接受的量级报告结果时要说明自己用的是模拟 p 值。另外fisher.test()的输出里除了 p 值还有odds ratio和置信区间。当某个单元格为 0 时优势比可能直接是 0 或无穷大这时候置信区间会异常宽解读时要格外小心。3. 当表从二维变成三维分层表与条件独立性3.1 mantelhaen.test() 处理 2×2×k 表很多数据分析场景不止是两个变量的关系第三个变量往往是“需要控制”的混杂因素。比如看新药是否有效你收集了两个研究中心的数据如果直接把两个中心的数据合并成一张 2×2 表很可能会出问题中心之间患者基线不同合并后可能把真实关联掩盖掉甚至出现方向反转这就是常说的辛普森悖论。R语言里处理这种情况的经典函数是mantelhaen.test()它的官方定位是“2×2×k 分层表的条件独立性检验”。也就是说它回答的问题是在控制了分层变量 Z 以后X 和 Y 是否仍然有关联。举个例子mytable - array( c(30, 20, 5, 15, 18, 12, 7, 31), dim c(2, 2, 2), dimnames list( 治疗 c(新药, 对照), 有效 c(是, 否), 研究中心 c(中心1, 中心2) ) ) mantelhaen.test(mytable)如果 p 值不显著说明在控制研究中心后治疗与效果之间没有统计学上的独立关系被推翻如果显著则说明关联在各层内依然存在。**这个函数只认 2×2×k 的array所以新建数据时dim千万不能写错。**我见过不少人在用array()构造三维表时把维数顺序写反导致行、列、层彻底错乱跑出来的 p 值自然没有任何意义。3.2 更广义的独立性模型用 loglm() 做对数线性模型比较mantelhaen.test()的局限性很明显它只处理二维分层表。如果变量不止三个或者不是标准的 2×2×k 结构可以换对数线性模型的思路。推荐用的是MASS包里的loglm()这个函数的基础用法和线性模型很像但拟合的对象是列联表频数的对数library(MASS) # 两变量独立性检验等价于卡方检验 loglm(~方案 是否购买, data tb)拟合“治疗方案、是否有效、研究中心”三分组数据时如果中心是混杂变量我们就想检验“治疗和有效在给定研究中心后是否独立”。对应模型可以写成loglm(~治疗 有效 研究中心 治疗:研究中心 有效:研究中心, data mytable)这个模型里没有治疗:有效这个交互项意思就是“治疗和有效不在同一个条件下直接关联”只允许它们各自跟研究中心有关。loglm()输出里的似然比卡方和 Pearson 卡方都能用来判断模型拟合是否够好p 值小说明模型与数据不符也就是“条件独立”的假设不成立。用对数线性模型的好处是它把独立性检验放到了一个统一的框架里想加入更多控制变量就继续在公式里加交互项灵活性比单独的函数高很多。缺点是理解门槛高新手很容易把公式里的.和:弄混。我平时只要场景还够用会先用mantelhaen.test()快速看结果如果是要写论文再用loglm()做完整的模型比较顺便把效应量和置信区间一起报告。4. 几个容易翻车的执行细节因子水平、缺失值与模拟 p 值4.1 因子里藏着的空白水平会悄悄改变自由度R语言里最隐蔽的坑之一就是因子变量明明只出现了三种取值水平定义里却有五个。table()默认会把这个因子所有水平都展示出来没出现的水平对应频数为 0。df$方案 - factor(df$方案, levels c(方案A, 方案B, 方案C, 方案D)) table(df$方案, df$是否购买)这时表里会多出一行全 0 的“方案D”自由度从(3-1)×(2-1)2变成了(4-1)×(2-1)3。如果你没看table()直接跑检验卡方统计量和 p 值都会是另一个东西。处理办法是在建模前把空白水平清掉df$方案 - droplevels(df$方案)或者在使用table()后删掉全零行和全零列tb - tb[rowSums(tb) 0, colSums(tb) 0, drop FALSE]另外数据里有缺失值时table()默认不把NA算进去但如果某个变量大量缺失而你又不检查样本量其实已经在“看不见的地方”缩水了。我建议先summary(df)看一眼各列缺失情况再用na.omit()处理如果缺失比例很高单靠删除样本会影响检验功效要先想清楚缺失机制再说。4.2 simulate.p.value TRUE 不是万能药当预期频数条件不满足、Fisher 精确检验又跑不动时chisq.test()提供了一个折中方案set.seed(1) chisq.test(tb, simulate.p.value TRUE, B 10000)它会根据当前表的结构随机生成大量独立同分布的列联表统计这些模拟表中出现更极端统计量的频率把它当作 p 值。注意三点第一这个 p 值不是确定性的。同样一份数据跑两次结果会略有波动所以必须用set.seed()保证可复现。第二B不能设得太小。B2000的模拟 p 值分辨率太粗我一般至少用 5000 到 10000。第三模拟 p 值替代不了精确检验。如果你的表只有 2×2 且样本量不大fisher.test()十秒内就能出结果没必要用模拟只有在表很大、Fisher 已经不可行时才考虑simulate.p.valueTRUE。4.3 配对设计千万别用独立性检验这是独立性检验函数最常见的误用场景。比如你研究同一批患者治疗前后的改善情况或者同一个产品在两种包装下的用户偏好数据是配对的两个人或两次测量之间并不独立。这时候再用chisq.test()去检验关联等于把“同一只羊”当成两只羊来数p 值会严重失真。配对二分类数据应该用mcnemar.test()paired_tb - matrix(c(35, 10, 15, 40), nrow 2, dimnames list(前后 c(改善, 未改善), 实际 c(改善, 未改善))) mcnemar.test(paired_tb)mcnemar.test()检验的是两个相关比例是否一致跟独立性检验是两回事。往数据框里塞一条“患者ID”列并试图做配对校正也不是chisq.test()能解决的。这一点如果写论文时搞错评审一句话就能把结论推翻。4.4 不要只盯 p 值效应量同样重要独立性检验得到显著 p 值只说明“关联不太可能是随机造成的”并不说明关联有多大。大样本下哪怕实际差异极小p 值也可能小于 0.05。R语言里看二分类变量的关联强度最常用的是vcd包的assocstats()install.packages(vcd) library(vcd) assocstats(tb)它会同时输出 phi 系数、Cramers V 和似然比卡方等。对于 2×2 表phi 系数就是相关系数范围在 -1 到 1 之间对更大的表Cramers V 通常按 0.1、0.3、0.5 分别看作小、中、大效应。虽然这个分界标准算不上公认铁律但在论文里报告比只写 p 值要有说服力得多。5. 把独立性检验做成可复用的“小工具”5.1 一个兼顾统计量和预期频数的包装函数项目里如果反复要做几十张列联表每次都手写chisq.test()再手动看$expected太低效。我习惯写一个简单的包装函数把最关心的信息一次性收集起来ind_summary - function(x, y NULL, simulate FALSE, B 10000) { tb - as.table(if (is.null(y)) x else table(x, y)) tb - tb[rowSums(tb) 0, colSums(tb) 0, drop FALSE] chi - suppressWarnings( chisq.test(tb, simulate.p.value simulate, B B) ) df_value - if (length(chi$parameter) 0) NA else unname(chi$parameter) data.frame( X2 unname(chi$statistic), df df_value, p.value chi$p.value, min.expected min(chi$expected), nrow nrow(tb), ncol ncol(tb) ) } ind_summary(df$方案, df$是否购买)函数里的两个关键设计是先把全零行和全零列删掉避免空白因子水平干扰自由度然后通过参数simulate手动决定是否启用蒙特卡洛模拟。返回值里我特意保留了min.expected这样每次跑完都能快速发现哪些表违反了预期频数条件。写这个函数不是为了炫技而是减少重复劳动和手滑。实际项目里数据列名可能很乱变量类型也可能不是因子这时候统一走同一个入口至少能保证每一步都用了同样一套判断逻辑。5.2 批量跑多变量组合并校正 p 值如果你手里有一个包含多个分类变量的数据框想快速筛查“哪些变量之间有潜在关联”可以用combn()把变量两两组合然后批量跑独立性检验var_names - c(方案, 是否购买, 年龄段, 性别) combos - combn(var_names, 2, simplify FALSE) batch_p - lapply(combos, function(v) { chi - chisq.test(df[[v[1]]], df[[v[2]]], simulate.p.value TRUE, B 5000) data.frame(x v[1], y v[2], p chi$p.value) }) p_vals - do.call(rbind, batch_p) p_vals$p_adjust - p.adjust(p_vals$p, method BH)最后一步p.adjust()很重要。做几十组两两检验时按 0.05 的阈值看纯靠随机也可能出现三五次“显著”所以必须做多重比较校正。method BH是控制假发现率的常用选择在筛查场景里比Bonferroni更不容易漏掉真实信号。批量检验适合做“初筛”它不能替代你针对核心假设做的正式分析。筛出有信号的变量对之后我会回到原始数据画一张分组条形图或者马赛克图看一眼方向再决定下一步用什么模型。最后再多说一句独立性检验函数在 R语言里只是工具箱的一小部分。真正写论文或者做业务复盘时千万不要把一个 p 值当作全部证据结合预期频数、效应量、分层控制才能让结论经得起问。遇到小样本和稀疏表优先考虑 Fisher 精确检验遇到复杂分层结构多想想mantelhaen.test()和loglm()。这些函数单个看都不难难的是在正确的地方调用它们。