R描述性统计与假设检验：t检验、Wilcoxon与卡方检验
分析数据的第一步不是建模，而是把分布摸清楚。两个变量均值都是 20，可能一个是「所有人都挤在 20 附近」，另一个是「一半人 5、一半人 35」——后者用均值描述就丢掉了全部信息，后续该选哪种检验也完全不同。
描述性统计：先看分位数，再看均值
Section titled “描述性统计：先看分位数，再看均值”summary() 是成本最低的一步，一次给出最小值、四分位数、中位数、均值和最大值。
summary(mtcars[, c("mpg", "wt", "hp")]) mpg wt hp Min. :10.40 Min. :1.513 Min. : 52.0 1st Qu.:15.43 1st Qu.:2.581 1st Qu.: 96.5 Median :19.20 Median :3.325 Median :123.0 Mean :20.09 Mean :3.217 Mean :146.7 3rd Qu.:22.80 3rd Qu.:3.610 3rd Qu.:180.0 Max. :33.90 Max. :5.424 Max. :335.0看 mpg 这一列：均值 20.09，中位数 19.20，两者差不到 1。hp 就不一样了，均值 146.7 明显大于中位数 123.0，是右偏（right-skewed）分布——少数大马力车型把均值拉了上去。中间一半的数据落在 96.5 到 180 之间，而最大值 335 离第三四分位数很远，这个变量放进线性模型时值得留意。
均值与中位数的差是判断偏斜的粗略办法，base R 没有现成的偏度函数，e1071 包的 skewness() 可以直接给数：
library(e1071)
c(mpg = skewness(mtcars$mpg), hp = skewness(mtcars$hp)) mpg hp0.6106550 0.7260237偏度（skewness）为 0 表示对称，正值为右偏，负值为左偏。经验上绝对值小于 0.5 算轻微偏斜，0.5 到 1 之间是中度。两个变量都在中度右偏的区间，hp 更重一些。
这里有个容易漏掉的细节：mpg 的均值和中位数只差 0.9，看起来接近对称，但偏度系数是 0.61，并不小。均值减中位数只用到两个统计量，对尾部信息的利用很粗糙；上下两端的极端值互相抵消时，这个差会偏小。判断偏斜程度以偏度系数为准，均值减中位数只能当快速一瞥。
需要单独取值时：
mean(mtcars$mpg)sd(mtcars$mpg)quantile(mtcars$mpg, c(0, 0.25, 0.5, 0.75, 1))[1] 20.09062[1] 6.026948 0% 25% 50% 75% 100%10.400 15.425 19.200 22.800 33.900sd() 默认算样本标准差（分母 n−1），不是总体标准差。quantile() 的第二个参数是概率向量，返回值的名字就是百分位。按分组算，base R 里最顺手的是 aggregate() 和 tapply()：
aggregate(mpg ~ cyl, data = mtcars, FUN = mean)tapply(mtcars$mpg, mtcars$cyl, sd)tapply(mtcars$mpg, mtcars$cyl, length) cyl mpg1 4 26.663642 6 19.742863 8 15.10000
4 6 84.509828 1.453567 2.560048
4 6 811 7 14气缸数越多油耗越高，方向符合常识；但 6 缸组的标准差只有 1.45，8 缸组是 2.56，离散程度差了一倍。后面做方差分析时，这个差距会不会影响方差齐性假设，得先记在心里。
标准差与标准误不是一回事
Section titled “标准差与标准误不是一回事”标准差（standard deviation）描述样本内部的离散程度：一个人和另一个人差多少。标准误（standard error of the mean）描述「样本均值」这个估计量的不确定度：换一批样本，均值会晃多少。两个量都叫「标准」，用途完全不同。
m <- mtcars$mpgsd(m)sd(m) / sqrt(length(m))[1] 6.026948[1] 1.065424base R 没有 sem() 这样的函数，标准误要自己按定义写 sd(x) / sqrt(length(x))。按分组算同样用 tapply()：
tapply(mtcars$mpg, mtcars$cyl, function(x) sd(x) / sqrt(length(x))) 4 6 81.3597642 0.5493967 0.6842016标准误总比标准差小，比例是 1/√n：n 越大，均值这个估计越稳。比较三组时不能只看标准差——6 缸组的样本量最少（7 辆），标准差也最小（1.45），两个因素相乘之后它的标准误 0.55 反而是三组里最低的。
论文里写「均值 ± 数值」时，多数期刊要求的是标准误。这两种写法给读者的印象差别很大：mtcars 的 mpg 写成 20.09 ± 6.03 是标准差，写成 20.09 ± 1.07 是标准误，前者说的是「车与车之间差很多」，后者说的是「均值估得很准」。两者都对，但必须写清楚是哪一个。图表上的误差棒同理，画之前先确认图注的定义。
缺失值：na.rm 不是可选项
Section titled “缺失值：na.rm 不是可选项”airquality 里有真实缺失值，是练手缺失值处理的好材料。
summary(airquality$Ozone)mean(airquality$Ozone)mean(airquality$Ozone, na.rm = TRUE)colSums(is.na(airquality)) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 1.00 18.00 31.50 42.13 63.25 168.00 37
[1] NA[1] 42.12931 Ozone Solar.R Wind Temp Month Day 37 7 0 0 0 0summary() 会在末尾报出缺失值个数，所以它比 mean() 更安全。mean() 碰到 NA 直接返回 NA，这不是 bug，是 R 在提醒你「还没决定怎么处理缺失值」。加 na.rm = TRUE 表示跳过，但要想清楚：跳过意味着默认这些数据是随机缺失的（missing at random）。如果缺失和观测值本身相关——比如污染最严重的那几天仪器正好坏了——算出来的均值就是有偏的。
正态性检验：shapiro.test()
Section titled “正态性检验：shapiro.test()”t 检验和方差分析都假设数据来自正态分布，样本量小的时候值得查一下。
shapiro.test(mtcars$mpg)shapiro.test(airquality$Ozone) Shapiro-Wilk normality test
data: mtcars$mpgW = 0.94756, p-value = 0.1229
Shapiro-Wilk normality test
data: airquality$OzoneW = 0.87867, p-value = 2.79e-08原假设是「数据来自正态分布」，所以 p 值大才说明没有证据反对正态性。mpg 的 p = 0.12，不能拒绝；Ozone 的 p 值极小且 W 明显偏离 1，是显著右偏。这个函数要求样本量在 3 到 5000 之间，超出会报错；而大样本下它过于敏感，几乎总能拒绝，这时候该看 QQ 图而不是 p 值。
单样本检验问的是「总体均值是否等于某个给定值」。mtcars 的平均油耗 20.09 和公认的 20 有区别吗：
t.test(mtcars$mpg, mu = 20) One Sample t-test
data: mtcars$mpgt = 0.08506, df = 31, p-value = 0.9328alternative hypothesis: true mean is not equal to 2095 percent confidence interval: 17.91768 22.26357sample estimates:mean of x 20.09062p = 0.93，没有证据说明总体均值不是 20，置信区间 [17.92, 22.26] 也包含 20。t 统计量只有 0.085，因为分子（样本均值减假设值）才 0.09，而数据本身波动很大：sd ≈ 6、n = 32，标准误约 1.06。
双样本比较用公式写法最省事。mtcars 的 am 是变速箱类型（0 自动，1 手动）：
t.test(mpg ~ am, data = mtcars) Welch Two Sample t-test
data: mpg by amt = -3.7671, df = 18.332, p-value = 0.001374alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 095 percent confidence interval: -11.280194 -3.209684sample estimates:mean in group 0 mean in group 1 17.14737 24.39231R 默认跑的是 Welch t 检验，不假设两组方差相等，自由度因此不是整数（18.332 而不是 30）。手动挡平均比自动挡多跑约 7.2 英里/加仑，p = 0.0014。
配对设计（同一个体前后测量）要用 paired = TRUE，但公式写法不接受这个参数：
t.test(extra ~ group, data = sleep, paired = TRUE)Error in t.test.formula(extra ~ group, data = sleep, paired = TRUE) : cannot use 'paired' in formula method这是很容易踩的一个坑。sleep 是配对数据，每位受试者在两组里各有一条记录，必须先把两组拆成两个向量：
t.test(sleep$extra[sleep$group == 1], sleep$extra[sleep$group == 2], paired = TRUE) Paired t-test
data: sleep$extra[sleep$group == 1] and sleep$extra[sleep$group == 2]t = -4.0621, df = 9, p-value = 0.002833alternative hypothesis: true mean difference is not equal to 095 percent confidence interval: -2.4598858 -0.7001142sample estimates:mean difference -1.58同一批数据按独立样本处理，结论完全不同：t = -1.8608，p = 0.0794，卡在 0.05 边缘。配对检验先把「同一个人」内部的差异减掉，消除了个体间差异带来的噪声，功效高得多。方法选错，真实效应会被判成不显著。
把差值的均值与标准误单独算出来，能看清配对检验到底做了什么：
d <- sleep$extra[sleep$group == 1] - sleep$extra[sleep$group == 2]mean(d)sd(d)sd(d) / sqrt(length(d))mean(d) / (sd(d) / sqrt(length(d)))[1] -1.58[1] 1.229995[1] 0.3889587[1] -4.062128配对的 t 统计量就是「差值的均值除以差值的标准误」：−1.58 / 0.389 = −4.06，与 t.test() 的输出一致。自由度是 9（10 对减 1），不是 18。20 个观测在这一步被压成 10 个差值，个体间差异同时被减掉，这就是自由度变小而功效反而变高的原因。
差值向量 d 本身也值得看一眼：sd(d) 只有 1.23，而两组原始数据的标准差分别是 1.79 和 2.00。差值比原始观测更集中，说明「同一个人的两次测量」正相关，两组的相关系数是 0.795。配对结构是真实存在的，用独立样本检验就浪费掉了这部分信息。
sd(sleep$extra[sleep$group == 1])sd(sleep$extra[sleep$group == 2])cor(sleep$extra[sleep$group == 1], sleep$extra[sleep$group == 2])[1] 1.78901[1] 2.002249[1] 0.7951702方差齐性：var.test()
Section titled “方差齐性：var.test()”经典（非 Welch）t 检验和方差分析都要求各组方差相等。
var.test(mpg ~ am, data = mtcars) F test to compare two variances
data: mpg by amF = 0.38656, num df = 18, denom df = 12, p-value = 0.06691alternative hypothesis: true ratio of variances is not equal to 195 percent confidence interval: 0.1243721 1.0703429sample estimates:ratio of variances 0.3865615p = 0.067，没有证据说明两组方差不齐。但这不等于「两组方差相等」——不拒绝原假设只是证据不足，别把 p > 0.05 读成「假设成立」。
F 检验假定两组都服从正态分布，对偏离正态敏感。换一个对分布形状要求更宽松的 Levene 检验，同一批数据给出的是另一个答案：
library(car)
leveneTest(mpg ~ factor(am), data = mtcars)Levene's Test for Homogeneity of Variance (center = median) Df F value Pr(>F)group 1 4.1876 0.04957 * 30---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1p = 0.0496，刚好落在 0.05 以下。同一份数据，var.test() 说不显著，Levene 说显著。原因在于两者的构造不同：F 检验比较两个样本方差之比，前提是数据正态；Levene 先算每个观测到组内中心的绝对离差，再对这些离差做方差分析，把「方差是否相等」转化成了「离差的均值是否相等」，对分布形状的要求低得多。数据有偏态时优先信 Levene。
leveneTest() 的分组变量必须是因子。直接写 leveneTest(mpg ~ am, data = mtcars) 会报错：
leveneTest(mpg ~ am, data = mtcars)Error in leveneTest.formula(mpg ~ am, data = mtcars) : Levene's test is not appropriate with quantitative explanatory variables.mtcars$am 是数值型的 0/1，函数会把它当成连续的解释变量而不是分组变量，所以要先 factor() 转一下。
两个检验给出不同结论，正好说明把「先检验方差齐性、再决定用哪种 t 检验」当成固定流程是不可靠的：换一个方差齐性检验就能翻转下一步的选择，而整个流程的实际错误率并不等于 0.05。样本量不大时，直接上 Welch t 检验通常更稳妥。
非参数替代：wilcox.test()
Section titled “非参数替代：wilcox.test()”数据明显非正态，或者只有序数尺度时，用基于秩的 Wilcoxon 检验：
wilcox.test(mpg ~ am, data = mtcars) Wilcoxon rank sum test with continuity correction
data: mpg by amW = 42, p-value = 0.001871alternative hypothesis: true location shift is not equal to 0
Warning message:In wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...) : cannot compute exact p-value with ties结论和 t 检验一致（都显著），但 p 值不同，因为它检验的是位置移动而非均值差。警告说存在并列值（ties），无法精确计算 p 值，改用正态近似。对连续变量来说这几乎必然发生——mpg 只保留一位小数，撞值很正常，可以接受；但如果数据全是整数评分、并列值极多，正态近似就不太可靠了。
wilcox.test() 只返回 W 和 p 值，没有效应量，而 W 本身没有量纲、不可解读。转成秩双列相关（rank-biserial correlation）就有界了：
w <- wilcox.test(mpg ~ am, data = mtcars)n1 <- sum(mtcars$am == 0)n2 <- sum(mtcars$am == 1)1 - 2 * unname(w$statistic) / (n1 * n2)[1] 0.659919秩双列相关取值在 −1 到 1 之间，绝对值越接近 1，两组分离得越彻底；0.66 表示两组的秩分布已经分离得比较明显，和 t 检验给出的显著结论方向一致。公式里的 W 是 wilcox.test() 返回的统计量，n₁、n₂ 是两组样本量（这里 19 和 13）。符号取决于把哪一组当作第一组，所以报告时要说清分组顺序，或者只报绝对值。
p 值和效应量回答的是两个问题：p = 0.0019 说的是「有证据说明两组不同」，r = 0.66 说的是「不同到什么程度」。样本量足够大时，微小的位置移动也能得到很小的 p 值，只报 p 值会让读者高估实际差异。
配对版本写成 wilcox.test(x, y, paired = TRUE)，同样不能用公式加 paired。
卡方检验：两个分类变量是否独立
Section titled “卡方检验：两个分类变量是否独立”table(mtcars$am, mtcars$vs)chisq.test(table(mtcars$am, mtcars$vs)) 0 1 0 12 7 1 6 7
Pearson's Chi-squared test with Yates' continuity correction
data: table(mtcars$am, mtcars$vs)X-squared = 0.34754, df = 1, p-value = 0.5555vs 是发动机形状（0 为 V 型，1 为直列）。p = 0.56，看不出变速箱类型和发动机形状有关联。2×2 表默认加 Yates 连续性校正，想关掉要显式传 correct = FALSE。
卡方检验有个容易被忽略的前提：每个格子的期望频数最好都不小于 5。换个变量试试：
chisq.test(table(mtcars$cyl, mtcars$am)) Pearson's Chi-squared test
data: table(mtcars$cyl, mtcars$am)X-squared = 8.7407, df = 2, p-value = 0.01265
Warning message:In chisq.test(table(mtcars$cyl, mtcars$am)) : Chi-squared approximation may be incorrectp = 0.013 看着显著，但警告说近似可能不准。R 不会告诉你原因，得自己把期望频数取出来看：
tab <- table(mtcars$cyl, mtcars$am)
chisq.test(tab)$expectedround(chisq.test(tab)$expected, 2)sum(chisq.test(tab)$expected < 5) 0 1 4 6.53125 4.46875 6 4.15625 2.84375 8 8.31250 5.68750
0 1 4 6.53 4.47 6 4.16 2.84 8 8.31 5.69[1] 3这张 3×2 的表有 3 个格子的期望频数低于 5，最小的是 6 缸手动挡的 2.84，远低于门槛值。卡方统计量依赖「观测频数近似正态」这个前提，期望频数太小的时候这个近似不成立，p 值也就不可信。这种情况应该改用 Fisher 精确检验，或者合并相邻类别，而不是照抄 p 值。
chisq.test() 的结果对象里，$expected 是期望频数，$observed 是原始列联表，$residuals 是每个格子的 Pearson 残差——残差最大的那个格子就是「最出乎意料」的格子，定位差异来源时比看整张 p 值有用。
卡方没有现成的效应量，报告时需要自己补一个 Cramér’s V：
ct <- chisq.test(tab)sqrt(unname(ct$statistic) / (sum(tab) * (min(dim(tab)) - 1)))[1] 0.5226355Cramér’s V 的算法是卡方统计量除以「样本量 ×（行列数较小者减 1）」，再开平方，取值在 0 到 1 之间。0.52 说明气缸数和变速箱类型之间有明显的关联。它和 p 值回答的不是同一个问题：p = 0.013 说的是「有证据说明两者有关联」，V = 0.52 说的是「关联到什么程度」。样本量越大，p 值越小；V 不受样本量影响。
2×2 表还要注意连续性校正的影响。同一张表跑三遍：
t2 <- table(mtcars$am, mtcars$vs)chisq.test(t2)chisq.test(t2, correct = FALSE)fisher.test(t2) Pearson's Chi-squared test with Yates' continuity correction
data: t2X-squared = 0.34754, df = 1, p-value = 0.5555
Pearson's Chi-squared test
data: t2X-squared = 0.90688, df = 1, p-value = 0.3409
Fisher's Exact Test for Count Data
data: t2p-value = 0.4727alternative hypothesis: true odds ratio is not equal to 195 percent confidence interval: 0.3825342 10.5916087sample estimates:odds ratio 1.956055关掉校正后统计量从 0.348 涨到 0.907，p 值从 0.56 降到 0.34，差别不小。默认开启校正的原因是 Yates 校正让检验更保守，在 2×2 表里更接近精确检验的结果。自由度大于 1 的表不受影响，这个参数只对 2×2 生效。
三个 p 值方向一致，都不显著。样本量小、或者每格期望频数接近 5 时优先用 fisher.test()，它算的是精确概率，不依赖任何近似。fisher.test() 还会附上优势比（odds ratio）的估计和置信区间，这本身就是个效应量，2×2 场合下比 Cramér’s V 更常被报告。
p 值到底在说什么
Section titled “p 值到底在说什么”几种最常见的误读，值得单独说清楚：
- p 值不是「原假设为真的概率」，也不是「结果由随机造成的概率」。
- p < 0.05 不表示效应有 95% 的概率存在。它的含义是：如果原假设成立，出现当前或更极端结果的概率小于 5%。
- p 值大小和效应大小无关。样本量足够大时，0.01 的均值差异也能得到 p < 0.001。报告结果时请连置信区间和效应量一起给出。
- p = 0.06 和 p = 0.04 之间没有质的差别，别把一个说成「无效应」、另一个说成「有效应」。
配合 t.test() 输出的置信区间一起读，比盯着 p 值本身可靠得多。
基础检验覆盖的是「两组比较」，但研究里更常见的是三组以上、或者要同时控制多个变量。前者用方差分析，后者用线性回归，两者背后其实是同一个线性模型框架。如果你也用 Python，可以对照 /python/visualization/seaborn 画分布图与箱线图，先看图再决定用哪种检验；三组以上的均值比较见 /r/modeling/anova，连续变量之间的线性关系见 /r/modeling/linear-regression。