R语言统计分析
小天使格礼
在Rstudio中安装本节课所有需要的包，运行以下指令：
install.packages(c(
"dplyr", "tidyr", "ggplot2",
"psych", "pastecs", "e1071",
"boot", "pwr", "reshape2",
"car", "DescTools", "MASS",
"vcd", "exact2x2", "coin",
"multcomp", "emmeans", "afex",
"ggpubr", "rstatix",
"corrplot", "ppcor",
"lm.beta", "broom"
))目录
| 章节 | 主题 | 核心内容 |
|---|---|---|
| 第一章 | 描述性统计基础 | 中心趋势、离散程度、分布形状度量 |
| 第二章 | 概率分布与随机变量 | 离散分布、连续分布、正态性检验 |
| 第三章 | 参数估计 | 点估计、置信区间、自举法 |
| 第四章 | 假设检验基础 | 原假设、p值、两类错误 |
| 第五章 | 常见参数检验 | t检验、比例检验、方差检验 |
| 第六章 | 非参数检验 | 秩检验、卡方检验、Fisher检验 |
| 第七章 | 方差分析 | 单因素、多因素、重复测量ANOVA |
| 第八章 | 相关分析与基础回归 | 相关系数、线性回归 |
| 第九章 | 样本量与功效分析 | 功效分析、样本量计算 |
| 第十章 | 统计模拟与重采样方法 | 蒙特卡洛模拟、自举法、置换检验 |
第一章：描述性统计基础
描述性统计是统计分析的第一步，它帮助我们了解数据的基本特征。在深入进行推断统计之前，我们需要先”认识”我们的数据——数据的中心在哪里？数据有多分散？数据的分布形状如何？这些问题都需要通过描述性统计来回答。
为什么描述性统计如此重要？
- 数据质量检查：通过描述性统计可以发现异常值、缺失值等问题
- 选择合适的统计方法：不同的数据分布需要不同的分析方法
- 结果解读的基础：了解数据的基本特征有助于正确解读分析结果
- 数据报告的必要组成部分：任何研究都需要报告样本的基本特征
1.1 中心趋势度量
中心趋势度量描述数据的”典型值”或”代表值”，告诉我们数据集中在哪个位置。不同的中心趋势度量适用于不同的数据类型和分布特征。
数学基础
设我们有一组数据 \(x_1, x_2, ..., x_n\)，各中心趋势度量的定义如下：
均值（Mean）： \[\bar{x} = \frac{1}{n}\sum_{i=1}^{n}x_i\]
中位数（Median）： 将数据从小到大排列后，位于中间位置的值。若n为奇数，中位数为第 \(\frac{n+1}{2}\) 个值；若n为偶数，中位数为第 \(\frac{n}{2}\) 和第 \(\frac{n}{2}+1\) 个值的平均。
众数（Mode）： 出现次数最多的值。
截尾均值（Trimmed Mean）： 去掉最大和最小的p%数据后计算均值： \[\bar{x}_{trimmed} = \frac{1}{n(1-2p)}\sum_{i=[np]+1}^{n-[np]}x_{(i)}\] 其中 \(x_{(i)}\) 表示排序后的第i个值。
使用场景与意义
| 度量指标 | 适用场景 | 优点 | 缺点 |
|---|---|---|---|
| 均值 | 对称分布、无极端值 | 利用所有信息、数学性质好 | 受极端值影响大 |
| 中位数 | 偏态分布、有极端值 | 不受极端值影响 | 信息利用不充分 |
| 众数 | 分类数据、多峰分布 | 适用于分类数据 | 可能不存在或不唯一 |
| 截尾均值 | 有少量极端值 | 兼顾稳健性和信息利用 | 需要选择截尾比例 |
# 创建示例数据：某班级学生考试成绩
# 包含一些极端值（非常低或非常高的分数）
scores <- c(85, 78, 92, 88, 76, 95, 82, 79, 91, 87,
84, 89, 77, 93, 86, 81, 90, 83, 88, 85,
35, 98) # 35和98是极端值
# 计算均值：所有数据的算术平均
# 均值 = 总和 / 个数
mean_score <- mean(scores)
mean_score## [1] 83.72727# 计算中位数：将数据排序后取中间值
# 中位数不受极端值影响，更能代表"典型"水平
median_score <- median(scores)
median_score## [1] 85.5# R基础包没有直接计算众数的函数
# 我们可以自定义一个函数来计算众数
get_mode <- function(x) {
# 计算每个值出现的频数
freq_table <- table(x)
# 找出出现次数最多的值
mode_value <- as.numeric(names(freq_table)[which.max(freq_table)])
return(mode_value)
}
# 对于连续数据，众数可能不太有意义
# 让我们创建一个有明确众数的数据
grades <- c("A", "B", "B", "C", "B", "A", "B", "C", "B", "A")
# 计算众数（出现次数最多的等级）
mode_grade <- get_mode(grades)## Warning in get_mode(grades): NAs introduced by coercionmode_grade## [1] NA# 截尾均值：去掉两端各10%的数据后计算均值
# 这可以减少极端值的影响
# trim参数指定截尾比例（两端各去掉的比例）
trimmed_mean <- mean(scores, trim = 0.1)
trimmed_mean## [1] 85.44444# 比较不同的中心趋势度量
cat("原始均值:", mean_score, "\n")## 原始均值: 83.72727cat("截尾均值:", trimmed_mean, "\n")## 截尾均值: 85.44444cat("中位数:", median_score, "\n")## 中位数: 85.5cat("均值与中位数的差异:", mean_score - median_score, "\n")## 均值与中位数的差异: -1.772727# 差异较大说明数据可能存在偏态或极端值解读要点： - 当均值 ≈ 中位数时，数据分布近似对称 - 当均值 > 中位数时，数据右偏（正偏），存在较大的极端值 - 当均值 < 中位数时，数据左偏（负偏），存在较小的极端值
1.2 离散程度度量
离散程度度量描述数据的”分散程度”或”变异性”，告诉我们数据围绕中心趋势的分布范围。仅知道中心趋势是不够的，我们还需要了解数据的波动情况。
数学基础
方差（Variance）： \[s^2 = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2\]
标准差（Standard Deviation）： \[s = \sqrt{s^2} = \sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2}\]
极差（Range）： \[R = x_{max} - x_{min}\]
四分位距（Interquartile Range, IQR）： \[IQR = Q_3 - Q_1\] 其中 \(Q_1\) 是第25百分位数，\(Q_3\) 是第75百分位数。
变异系数（Coefficient of Variation, CV）： \[CV = \frac{s}{\bar{x}} \times 100\%\]
使用场景与意义
| 度量指标 | 适用场景 | 特点 |
|---|---|---|
| 方差 | 统计推断的基础 | 单位是原单位的平方，不直观 |
| 标准差 | 描述数据分散程度 | 与原数据同单位，最常用 |
| 极差 | 快速了解数据范围 | 只利用两端信息，不稳定 |
| IQR | 有极端值时的稳健度量 | 不受极端值影响 |
| CV | 比较不同量纲数据的变异 | 无量纲，便于比较 |
# 使用之前的考试成绩数据
scores <- c(85, 78, 92, 88, 76, 95, 82, 79, 91, 87,
84, 89, 77, 93, 86, 81, 90, 83, 88, 85,
35, 98)
# 计算方差：衡量数据偏离均值的程度
# 方差越大，数据越分散
var_scores <- var(scores)
var_scores## [1] 152.684# 计算标准差：方差的平方根
# 标准差与原数据同单位，更易解释
sd_scores <- sd(scores)
sd_scores## [1] 12.35654# 计算极差：最大值减最小值
# 极差只反映数据的范围，不反映中间分布
range_scores <- diff(range(scores))
range_scores## [1] 63# 或者使用
range_result <- range(scores)
range_result## [1] 35 98# 计算四分位距：第75百分位数 - 第25百分位数
# IQR是稳健的离散度量，不受极端值影响
iqr_scores <- IQR(scores)
iqr_scores## [1] 8.5# 也可以用quantile函数计算
q1 <- quantile(scores, 0.25)
q3 <- quantile(scores, 0.75)
iqr_manual <- q3 - q1
iqr_manual## 75%
## 8.5# 计算变异系数：标准差/均值
# CV是无量纲的，可用于比较不同量纲数据的变异程度
cv_scores <- sd(scores) / mean(scores) * 100
cv_scores## [1] 14.75808# 示例：比较两组不同量纲数据的变异程度
# 身高数据（厘米）
height <- c(170, 175, 168, 172, 180, 165, 178, 173)
# 体重数据（千克）
weight <- c(65, 70, 62, 68, 75, 60, 72, 67)
# 直接比较标准差没有意义（单位不同）
sd(height)## [1] 5.012484sd(weight)## [1] 5.012484# 使用变异系数比较
cv_height <- sd(height) / mean(height) * 100
cv_weight <- sd(weight) / mean(weight) * 100
cat("身高CV:", cv_height, "%\n")## 身高CV: 2.903684 %cat("体重CV:", cv_weight, "%\n")## 体重CV: 7.43968 %# CV较大说明相对变异程度更大解读要点： - 标准差与均值结合使用：对于正态分布，约68%的数据在均值±1个标准差范围内 - IQR与中位数结合使用：对于任何分布，约50%的数据在Q1到Q3范围内 - CV > 30%通常表示数据变异较大
1.3 分布形状度量
分布形状度量描述数据分布的对称性和尖锐程度，帮助我们判断数据是否符合正态分布等理论分布。
数学基础
偏度（Skewness）： \[g_1 = \frac{n}{(n-1)(n-2)}\sum_{i=1}^{n}\left(\frac{x_i - \bar{x}}{s}\right)^3\]
偏度衡量分布的对称性： - \(g_1 = 0\)：对称分布 - \(g_1 > 0\)：右偏（正偏），右侧有长尾 - \(g_1 < 0\)：左偏（负偏），左侧有长尾
峰度（Kurtosis）： \[g_2 = \frac{n(n+1)}{(n-1)(n-2)(n-3)}\sum_{i=1}^{n}\left(\frac{x_i - \bar{x}}{s}\right)^4 - \frac{3(n-1)^2}{(n-2)(n-3)}\]
峰度衡量分布的”尖峭”程度（相对于正态分布）： - \(g_2 = 0\)：与正态分布相同（正态峰） - \(g_2 > 0\)：尖峰分布（尾部更厚） - \(g_2 < 0\)：平峰分布（尾部更薄）
使用场景与意义
# 创建不同分布特征的数据进行演示
# 对称分布数据
set.seed(123)
symmetric_data <- rnorm(1000, mean = 50, sd = 10)
# 右偏分布数据（如收入分布）
right_skewed <- rexp(1000, rate = 0.1)
# 左偏分布数据
left_skewed <- 100 - rexp(1000, rate = 0.1)
# 使用e1071包计算偏度和峰度
# type = 2 是默认类型，适用于小样本
# 偏度
skew_symmetric <- e1071::skewness(symmetric_data)
skew_right <- e1071::skewness(right_skewed)
skew_left <- e1071::skewness(left_skewed)
cat("对称分布偏度:", skew_symmetric, "\n")## 对称分布偏度: 0.065196cat("右偏分布偏度:", skew_right, "\n")## 右偏分布偏度: 1.68554cat("左偏分布偏度:", skew_left, "\n")## 左偏分布偏度: -2.170876# 峰度
kurt_symmetric <- e1071::kurtosis(symmetric_data)
kurt_right <- e1071::kurtosis(right_skewed)
kurt_left <- e1071::kurtosis(left_skewed)
cat("对称分布峰度:", kurt_symmetric, "\n")## 对称分布峰度: -0.08010201cat("右偏分布峰度:", kurt_right, "\n")## 右偏分布峰度: 3.074591cat("左偏分布峰度:", kurt_left, "\n")## 左偏分布峰度: 6.699674# 可视化不同分布的形状
par(mfrow = c(1, 3))
# 对称分布
hist(symmetric_data, main = "对称分布",
xlab = "值", col = "lightblue", breaks = 30, freq = FALSE)
curve(dnorm(x, mean = mean(symmetric_data), sd = sd(symmetric_data)),
add = TRUE, col = "red", lwd = 2)
# 右偏分布
hist(right_skewed, main = paste("右偏分布\n偏度 =", round(skew_right, 2)),
xlab = "值", col = "lightgreen", breaks = 30, freq = FALSE)
# 左偏分布
hist(left_skewed, main = paste("左偏分布\n偏度 =", round(skew_left, 2)),
xlab = "值", col = "lightyellow", breaks = 30, freq = FALSE)
par(mfrow = c(1, 1))解读要点： - 偏度绝对值 > 1 表示严重偏态，可能需要数据转换 - 峰度绝对值 > 1 表示与正态分布差异较大 - 偏态会影响均值作为中心趋势度量的代表性
1.4 分位数与百分位数
分位数是将数据等分的分割点，百分位数是分位数的特例（将数据分为100等份）。它们在描述数据分布和确定参考范围时非常有用。
数学基础
第p百分位数 \(P_p\) 是满足以下条件的值： \[P(X \leq P_p) \geq \frac{p}{100} \text{ 且 } P(X \geq P_p) \geq \frac{100-p}{100}\]
常用的分位数包括： - 四分位数：将数据分为4等份（\(Q_1=25\%, Q_2=50\%, Q_3=75\%\)） - 十分位数：将数据分为10等份 - 百分位数：将数据分为100等份
使用场景与意义
分位数在医学参考值范围的确定中尤为重要。例如，儿童生长发育标准常用百分位数表示。
# 使用考试成绩数据
scores <- c(85, 78, 92, 88, 76, 95, 82, 79, 91, 87,
84, 89, 77, 93, 86, 81, 90, 83, 88, 85,
35, 98)
# 计算四分位数
quartiles <- quantile(scores, probs = c(0.25, 0.5, 0.75))
quartiles## 25% 50% 75%
## 81.25 85.50 89.75# 计算常用的百分位数
percentiles <- quantile(scores, probs = c(0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95))
percentiles## 5% 10% 25% 50% 75% 90% 95%
## 76.05 77.10 81.25 85.50 89.75 92.90 94.90# 计算任意百分位数
# 例如，第10百分位数和第90百分位数
p10 <- quantile(scores, 0.10)
p90 <- quantile(scores, 0.90)
cat("第10百分位数:", p10, "\n")## 第10百分位数: 77.1cat("第90百分位数:", p90, "\n")## 第90百分位数: 92.9# 五数概括：最小值、Q1、中位数、Q3、最大值
fivenum(scores)## [1] 35.0 81.0 85.5 90.0 98.0# 使用summary函数获取基本统计量
summary(scores)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 35.00 81.25 85.50 83.73 89.75 98.00# 医学应用示例：儿童身高参考值
# 假设这是某年龄组男孩的身高数据（厘米）
set.seed(456)
boys_height <- rnorm(500, mean = 120, sd = 8)
# 计算常用百分位数（用于生长曲线）
height_percentiles <- quantile(boys_height, probs = c(0.03, 0.10, 0.25, 0.50, 0.75, 0.90, 0.97))
height_percentiles## 3% 10% 25% 50% 75% 90% 97%
## 107.4552 110.9130 115.1824 120.6781 126.5114 130.9721 134.9521# 解释：第3百分位数意味着只有3%的儿童身高低于此值
# 常用于判断是否低于正常范围解读要点： - 第50百分位数 = 中位数 - IQR = Q3 - Q1，反映中间50%数据的分布范围 - 医学上常用第5和第95百分位数（或第3和第97百分位数）作为参考范围
1.5 汇总统计函数
R提供了多种汇总统计函数，可以一次性获取多个统计量。常用的有summary()、describe()（psych包）和stat.desc()（pastecs包）。
使用场景与意义
不同的汇总函数提供不同详细程度和类型的统计量：
| 函数 | 来源 | 特点 |
|---|---|---|
| summary() | 基础R | 提供基本六数概括 |
| describe() | psych包 | 提供偏度、峰度等更多统计量 |
| stat.desc() | pastecs包 | 提供最全面的描述统计 |
# 创建示例数据框
# 模拟一项医学研究的基线数据
set.seed(789)
medical_data <- data.frame(
patient_id = 1:50,
age = round(rnorm(50, mean = 55, sd = 12)),
bmi = round(rnorm(50, mean = 25, sd = 3), 1),
sbp = round(rnorm(50, mean = 130, sd = 15)), # 收缩压
dbp = round(rnorm(50, mean = 85, sd = 10)), # 舒张压
glucose = round(rnorm(50, mean = 5.5, sd = 1), 1), # 血糖
group = sample(c("Treatment", "Control"), 50, replace = TRUE)
)
# 查看数据结构
str(medical_data)## 'data.frame': 50 obs. of 7 variables:
## $ patient_id: int 1 2 3 4 5 6 7 8 9 10 ...
## $ age : num 61 28 55 57 51 49 47 53 43 64 ...
## $ bmi : num 24.9 20.4 27.4 24.6 22.9 26.8 28.3 22.9 28.5 28.7 ...
## $ sbp : num 122 124 139 132 132 139 141 114 140 127 ...
## $ dbp : num 79 80 87 99 77 86 81 78 84 92 ...
## $ glucose : num 2.3 3.8 4.3 7.5 5.1 4.7 5.4 6.4 3.8 5.7 ...
## $ group : chr "Control" "Control" "Treatment" "Treatment" ...# 使用summary()：基础R函数
# 提供最小值、Q1、中位数、均值、Q3、最大值
summary(medical_data)## patient_id age bmi sbp
## Min. : 1.00 Min. :18.00 Min. :17.10 Min. : 85.0
## 1st Qu.:13.25 1st Qu.:47.00 1st Qu.:22.95 1st Qu.:120.5
## Median :25.50 Median :53.00 Median :25.45 Median :132.0
## Mean :25.50 Mean :53.30 Mean :25.20 Mean :129.6
## 3rd Qu.:37.75 3rd Qu.:61.75 3rd Qu.:27.40 3rd Qu.:138.8
## Max. :50.00 Max. :76.00 Max. :32.00 Max. :176.0
## dbp glucose group
## Min. : 63.00 Min. :2.300 Length:50
## 1st Qu.: 77.00 1st Qu.:4.825 Class :character
## Median : 84.00 Median :5.350 Mode :character
## Mean : 83.76 Mean :5.378
## 3rd Qu.: 91.50 3rd Qu.:6.250
## Max. :105.00 Max. :7.500# 使用psych包的describe()函数
# 提供更多统计量，包括标准差、偏度、峰度等
psych::describe(medical_data[, c("age", "bmi", "sbp", "dbp", "glucose")])## vars n mean sd median trimmed mad min max range skew
## age 1 50 53.30 11.02 53.00 53.70 11.12 18.0 76.0 58.0 -0.52
## bmi 2 50 25.20 3.09 25.45 25.26 3.63 17.1 32.0 14.9 -0.17
## sbp 3 50 129.58 15.55 132.00 130.15 10.38 85.0 176.0 91.0 -0.19
## dbp 4 50 83.76 10.59 84.00 83.72 10.38 63.0 105.0 42.0 0.07
## glucose 5 50 5.38 1.02 5.35 5.42 0.96 2.3 7.5 5.2 -0.46
## kurtosis se
## age 0.84 1.56
## bmi -0.43 0.44
## sbp 1.20 2.20
## dbp -0.74 1.50
## glucose 0.16 0.14# 使用pastecs包的stat.desc()函数
# 提供最全面的描述统计，包括标准误、置信区间等
pastecs::stat.desc(medical_data[, c("age", "bmi", "sbp", "dbp", "glucose")])## age bmi sbp dbp glucose
## nbr.val 50.0000000 50.0000000 50.0000000 50.0000000 50.0000000
## nbr.null 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## nbr.na 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## min 18.0000000 17.1000000 85.0000000 63.0000000 2.3000000
## max 76.0000000 32.0000000 176.0000000 105.0000000 7.5000000
## range 58.0000000 14.9000000 91.0000000 42.0000000 5.2000000
## sum 2665.0000000 1260.1000000 6479.0000000 4188.0000000 268.9000000
## median 53.0000000 25.4500000 132.0000000 84.0000000 5.3500000
## mean 53.3000000 25.2020000 129.5800000 83.7600000 5.3780000
## SE.mean 1.5587148 0.4366336 2.1992745 1.4979060 0.1443717
## CI.mean.0.95 3.1323547 0.8774481 4.4196075 3.0101548 0.2901258
## var 121.4795918 9.5324449 241.8404082 112.1861224 1.0421592
## std.dev 11.0217781 3.0874658 15.5512189 10.5917951 1.0208620
## coef.var 0.2067876 0.1225088 0.1200125 0.1264541 0.1898219# stat.desc()的norm = TRUE参数可进行正态性检验
pastecs::stat.desc(medical_data[, c("age", "bmi", "sbp")], norm = TRUE)## age bmi sbp
## nbr.val 50.0000000 50.0000000 50.0000000
## nbr.null 0.0000000 0.0000000 0.0000000
## nbr.na 0.0000000 0.0000000 0.0000000
## min 18.0000000 17.1000000 85.0000000
## max 76.0000000 32.0000000 176.0000000
## range 58.0000000 14.9000000 91.0000000
## sum 2665.0000000 1260.1000000 6479.0000000
## median 53.0000000 25.4500000 132.0000000
## mean 53.3000000 25.2020000 129.5800000
## SE.mean 1.5587148 0.4366336 2.1992745
## CI.mean.0.95 3.1323547 0.8774481 4.4196075
## var 121.4795918 9.5324449 241.8404082
## std.dev 11.0217781 3.0874658 15.5512189
## coef.var 0.2067876 0.1225088 0.1200125
## skewness -0.5182621 -0.1734366 -0.1870188
## skew.2SE -0.7698469 -0.2576296 -0.2778051
## kurtosis 0.8369006 -0.4282681 1.1970114
## kurt.2SE 0.6321876 -0.3235101 0.9042123
## normtest.W 0.9679493 0.9891814 0.9610137
## normtest.p 0.1907601 0.9252156 0.0980297# skew.2SE和kurt.2SE列：如果绝对值>1，说明偏度/峰度显著解读要点： -
describe()中的se是标准误，反映均值估计的精度 -
stat.desc()中的skew.2SE和kurt.2SE绝对值>1表示显著偏离正态
1.6 分组汇总
在实际研究中，我们经常需要按组计算统计量，例如比较不同治疗组的基线特征。R提供了多种分组汇总的方法。
使用场景与意义
分组汇总是数据探索和结果报告的基础，常用于： - 比较不同组别的基线特征 - 计算各组的描述性统计 - 生成汇总表格用于论文报告
# 使用之前创建的医学数据
# 按组别计算各变量的描述性统计
# 方法1：使用tapply函数（基础R）
# tapply(数据, 分组变量, 统计函数)
# 计算各组的平均年龄
tapply(medical_data$age, medical_data$group, mean)## Control Treatment
## 51.08333 55.34615# 计算各组的年龄标准差
tapply(medical_data$age, medical_data$group, sd)## Control Treatment
## 10.02569 11.68569# 计算多个统计量
tapply(medical_data$age, medical_data$group, function(x) {
c(mean = mean(x), sd = sd(x), min = min(x), max = max(x))
})## $Control
## mean sd min max
## 51.08333 10.02569 28.00000 74.00000
##
## $Treatment
## mean sd min max
## 55.34615 11.68569 18.00000 76.00000# 方法2：使用aggregate函数（基础R）
# aggregate(公式, 数据, 统计函数)
# 计算各组的平均年龄
aggregate(age ~ group, data = medical_data, mean)## group age
## 1 Control 51.08333
## 2 Treatment 55.34615# 计算多个变量的均值
aggregate(cbind(age, bmi, sbp) ~ group, data = medical_data, mean)## group age bmi sbp
## 1 Control 51.08333 24.52500 124.5000
## 2 Treatment 55.34615 25.82692 134.2692# 方法3：使用dplyr包（推荐）
# dplyr提供更灵活和可读性更好的语法
# group_by() + summarise() 组合
medical_data %>%
dplyr::group_by(group) %>%
dplyr::summarise(
n = dplyr::n(), # 样本量
mean_age = mean(age), # 平均年龄
sd_age = sd(age), # 年龄标准差
mean_bmi = mean(bmi), # 平均BMI
sd_bmi = sd(bmi), # BMI标准差
mean_sbp = mean(sbp), # 平均收缩压
median_glucose = median(glucose) # 血糖中位数
)## # A tibble: 2 × 8
## group n mean_age sd_age mean_bmi sd_bmi mean_sbp median_glucose
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Control 24 51.1 10.0 24.5 2.61 124. 5.3
## 2 Treatment 26 55.3 11.7 25.8 3.40 134. 5.6# 更复杂的分组汇总
# 计算各组的详细描述统计
group_summary <- medical_data %>%
dplyr::group_by(group) %>%
dplyr::summarise(
across(c(age, bmi, sbp, dbp, glucose),
list(mean = ~mean(.), sd = ~sd(.), median = ~median(.)),
.names = "{.col}_{.fn}")
)
# 查看结果
group_summary## # A tibble: 2 × 16
## group age_mean age_sd age_median bmi_mean bmi_sd bmi_median sbp_mean sbp_sd
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Control 51.1 10.0 50.5 24.5 2.61 24.0 124. 14.7
## 2 Treatme… 55.3 11.7 58 25.8 3.40 26.6 134. 15.1
## # ℹ 7 more variables: sbp_median <dbl>, dbp_mean <dbl>, dbp_sd <dbl>,
## # dbp_median <dbl>, glucose_mean <dbl>, glucose_sd <dbl>,
## # glucose_median <dbl># 使用psych包的describeBy函数
# 可以一次性获取各组的完整描述统计
psych::describeBy(medical_data[, c("age", "bmi", "sbp", "glucose")],
group = medical_data$group)##
## Descriptive statistics by group
## group: Control
## vars n mean sd median trimmed mad min max range skew
## age 1 24 51.08 10.03 50.50 51.00 6.67 28.0 74.0 46.0 0.11
## bmi 2 24 24.52 2.61 24.05 24.37 2.67 20.4 30.5 10.1 0.53
## sbp 3 24 124.50 14.70 127.50 125.95 14.08 85.0 143.0 58.0 -0.92
## glucose 4 24 5.25 1.01 5.30 5.34 0.82 2.3 6.6 4.3 -1.00
## kurtosis se
## age 0.04 2.05
## bmi -0.46 0.53
## sbp 0.15 3.00
## glucose 0.84 0.21
## ------------------------------------------------------------
## group: Treatment
## vars n mean sd median trimmed mad min max range skew
## age 1 26 55.35 11.69 58.0 56.09 9.64 18.0 76.0 58.0 -1.03
## bmi 2 26 25.83 3.40 26.6 26.04 2.97 17.1 32.0 14.9 -0.69
## sbp 3 26 134.27 15.09 135.5 133.95 10.38 103.0 176.0 73.0 0.35
## glucose 4 26 5.50 1.03 5.6 5.50 1.26 3.8 7.5 3.7 0.01
## kurtosis se
## age 1.77 2.29
## bmi -0.09 0.67
## sbp 0.84 2.96
## glucose -1.11 0.20解读要点： -
dplyr方法语法清晰，适合复杂数据处理管道 -
describeBy()适合快速生成报告用的描述统计表 -
分组汇总结果是论文”表1”的主要内容
1.7 频数表与列联表
频数表和列联表用于描述分类变量的分布情况，是分类数据分析的基础。
数学基础
频数（Frequency）：各类别的出现次数 相对频数（Relative Frequency）：各类别的比例 \[f_{rel} = \frac{n_i}{n}\]
列联表（Contingency Table）：两个或多个分类变量的交叉频数表
使用场景与意义
频数表和列联表在以下场景中常用： - 描述样本的基本特征（如性别、疾病分期分布） - 分析两个分类变量的关联性 - 为卡方检验等统计检验做准备
# 创建分类变量数据
set.seed(111)
survey_data <- data.frame(
gender = sample(c("Male", "Female"), 200, replace = TRUE),
age_group = sample(c("Young", "Middle", "Senior"), 200, replace = TRUE),
treatment = sample(c("Drug A", "Drug B", "Placebo"), 200, replace = TRUE),
outcome = sample(c("Improved", "No Change", "Worsened"), 200, replace = TRUE)
)
# 单变量频数表
# 使用table()函数
gender_table <- table(survey_data$gender)
gender_table##
## Female Male
## 100 100# 添加边际合计
addmargins(gender_table)##
## Female Male Sum
## 100 100 200# 计算比例（相对频数）
prop.table(gender_table)##
## Female Male
## 0.5 0.5# 转换为百分比
prop.table(gender_table) * 100##
## Female Male
## 50 50# 二维列联表（两个分类变量）
# 治疗方案与疗效的关系
treatment_outcome <- table(survey_data$treatment, survey_data$outcome)
treatment_outcome##
## Improved No Change Worsened
## Drug A 20 19 18
## Drug B 21 26 29
## Placebo 20 26 21# 添加边际合计
addmargins(treatment_outcome)##
## Improved No Change Worsened Sum
## Drug A 20 19 18 57
## Drug B 21 26 29 76
## Placebo 20 26 21 67
## Sum 61 71 68 200# 计算行百分比（每组各类结果的比例）
# margin = 1 表示按行计算
prop.table(treatment_outcome, margin = 1) * 100##
## Improved No Change Worsened
## Drug A 35.08772 33.33333 31.57895
## Drug B 27.63158 34.21053 38.15789
## Placebo 29.85075 38.80597 31.34328# 计算列百分比（各结果中不同组的比例）
# margin = 2 表示按列计算
prop.table(treatment_outcome, margin = 2) * 100##
## Improved No Change Worsened
## Drug A 32.78689 26.76056 26.47059
## Drug B 34.42623 36.61972 42.64706
## Placebo 32.78689 36.61972 30.88235# 计算总百分比
prop.table(treatment_outcome) * 100##
## Improved No Change Worsened
## Drug A 10.0 9.5 9.0
## Drug B 10.5 13.0 14.5
## Placebo 10.0 13.0 10.5# 使用xtabs()函数创建列联表（公式语法）
# 语法：xtabs(频数 ~ 行变量 + 列变量, data)
xtabs(~ treatment + outcome, data = survey_data)## outcome
## treatment Improved No Change Worsened
## Drug A 20 19 18
## Drug B 21 26 29
## Placebo 20 26 21# 三维列联表
three_way <- xtabs(~ gender + treatment + outcome, data = survey_data)
three_way## , , outcome = Improved
##
## treatment
## gender Drug A Drug B Placebo
## Female 9 7 8
## Male 11 14 12
##
## , , outcome = No Change
##
## treatment
## gender Drug A Drug B Placebo
## Female 9 17 10
## Male 10 9 16
##
## , , outcome = Worsened
##
## treatment
## gender Drug A Drug B Placebo
## Female 11 19 10
## Male 7 10 11# 查看特定维度的边际表
# 按性别分层查看治疗与疗效的关系
ftable(three_way) # 扁平化显示## outcome Improved No Change Worsened
## gender treatment
## Female Drug A 9 9 11
## Drug B 7 17 19
## Placebo 8 10 10
## Male Drug A 11 10 7
## Drug B 14 9 10
## Placebo 12 16 11# 使用dplyr创建频数表
freq_table <- survey_data %>%
dplyr::count(treatment, outcome) %>%
dplyr::group_by(treatment) %>%
dplyr::mutate(
percent = n / sum(n) * 100,
cum_percent = cumsum(percent)
)
freq_table## # A tibble: 9 × 5
## # Groups: treatment [3]
## treatment outcome n percent cum_percent
## <chr> <chr> <int> <dbl> <dbl>
## 1 Drug A Improved 20 35.1 35.1
## 2 Drug A No Change 19 33.3 68.4
## 3 Drug A Worsened 18 31.6 100
## 4 Drug B Improved 21 27.6 27.6
## 5 Drug B No Change 26 34.2 61.8
## 6 Drug B Worsened 29 38.2 100
## 7 Placebo Improved 20 29.9 29.9
## 8 Placebo No Change 26 38.8 68.7
## 9 Placebo Worsened 21 31.3 100# 创建美观的汇总表
# 使用table()和prop.table()组合
create_crosstab <- function(data, row_var, col_var) {
# 创建频数表
freq <- table(data[[row_var]], data[[col_var]])
# 创建行百分比表
row_pct <- prop.table(freq, 1) * 100
# 合并为一个表格
result <- matrix(paste0(freq, " (", round(row_pct, 1), "%)"),
nrow = nrow(freq))
rownames(result) <- rownames(freq)
colnames(result) <- colnames(freq)
return(result)
}
create_crosstab(survey_data, "treatment", "outcome")## Improved No Change Worsened
## Drug A "20 (35.1%)" "19 (33.3%)" "18 (31.6%)"
## Drug B "21 (27.6%)" "26 (34.2%)" "29 (38.2%)"
## Placebo "20 (29.9%)" "26 (38.8%)" "21 (31.3%)"解读要点： - 行百分比用于比较不同组的分布差异 - 列百分比用于了解某结果来自哪个组 - 列联表是卡方检验的基础
1.8 多变量描述：协方差与相关系数
当研究多个连续变量之间的关系时，协方差和相关系数是重要的描述性工具。
数学基础
协方差（Covariance）： \[Cov(X,Y) = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})\]
协方差的符号表示关系的方向： - \(Cov > 0\)：正相关（同向变化） - \(Cov < 0\)：负相关（反向变化） - \(Cov = 0\)：无线性相关
Pearson相关系数： \[r = \frac{Cov(X,Y)}{s_X \cdot s_Y} = \frac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum(x_i-\bar{x})^2 \cdot \sum(y_i-\bar{y})^2}}\]
相关系数的范围：\(-1 \leq r \leq 1\) - \(|r| \approx 1\)：强相关 - \(|r| \approx 0\)：弱相关或无相关
使用场景与意义
| 指标 | 范围 | 特点 |
|---|---|---|
| 协方差 | \((-\infty, +\infty)\) | 受量纲影响，难以直接解释 |
| 相关系数 | \([-1, 1]\) | 无量纲，便于比较和解释 |
# 使用医学数据演示协方差和相关系数
# 计算两个变量的协方差
# 收缩压和舒张压的协方差
cov_sb_db <- cov(medical_data$sbp, medical_data$dbp)
cov_sb_db## [1] -11.38857# 计算相关系数
cor_sb_db <- cor(medical_data$sbp, medical_data$dbp)
cor_sb_db## [1] -0.06914093# 相关系数的解释
if (abs(cor_sb_db) >= 0.8) {
cat("强相关\n")
} else if (abs(cor_sb_db) >= 0.5) {
cat("中等相关\n")
} else if (abs(cor_sb_db) >= 0.3) {
cat("弱相关\n")
} else {
cat("极弱或无相关\n")
}## 极弱或无相关# 计算相关系数矩阵
# 多个变量之间的两两相关
cor_matrix <- cor(medical_data[, c("age", "bmi", "sbp", "dbp", "glucose")])
cor_matrix## age bmi sbp dbp glucose
## age 1.000000000 0.003160536 0.01098979 -0.007587041 0.04993343
## bmi 0.003160536 1.000000000 0.10542960 -0.104017226 0.06573473
## sbp 0.010989793 0.105429600 1.00000000 -0.069140935 -0.11950275
## dbp -0.007587041 -0.104017226 -0.06914093 1.000000000 -0.11091203
## glucose 0.049933435 0.065734733 -0.11950275 -0.110912034 1.00000000# 计算协方差矩阵
cov_matrix <- cov(medical_data[, c("age", "bmi", "sbp", "dbp", "glucose")])
cov_matrix## age bmi sbp dbp glucose
## age 121.4795918 0.1075510 1.883673 -0.8857143 0.5618367
## bmi 0.1075510 9.5324449 5.062082 -3.4015510 0.2071878
## sbp 1.8836735 5.0620816 241.840408 -11.3885714 -1.8971837
## dbp -0.8857143 -3.4015510 -11.388571 112.1861224 -1.1992653
## glucose 0.5618367 0.2071878 -1.897184 -1.1992653 1.0421592# 可视化相关系数矩阵
# 使用corrplot包（如果已安装）
# 这里用基础方法演示
# 将相关矩阵转换为长格式
cor_df <- as.data.frame(as.table(cor_matrix))
names(cor_df) <- c("Var1", "Var2", "Correlation")
# 绘制热图
ggplot2::ggplot(cor_df, ggplot2::aes(x = Var1, y = Var2, fill = Correlation)) +
ggplot2::geom_tile() +
ggplot2::scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0, limits = c(-1, 1)) +
ggplot2::geom_text(ggplot2::aes(label = round(Correlation, 2)), size = 3) +
ggplot2::theme_minimal() +
ggplot2::labs(title = "变量相关系数矩阵热图")
# 散点图展示两个变量的关系
ggplot2::ggplot(medical_data, ggplot2::aes(x = sbp, y = dbp)) +
ggplot2::geom_point(alpha = 0.6) +
ggplot2::geom_smooth(method = "lm", se = TRUE, color = "red") +
ggplot2::labs(title = paste("收缩压与舒张压的关系 (r =", round(cor_sb_db, 3), ")"),
x = "收缩压 (mmHg)",
y = "舒张压 (mmHg)") +
ggplot2::theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
# 不同类型的相关系数
# Pearson相关（默认）：适用于连续变量、线性关系
cor_pearson <- cor(medical_data$sbp, medical_data$dbp, method = "pearson")
# Spearman相关：适用于等级数据或非线性单调关系
cor_spearman <- cor(medical_data$sbp, medical_data$dbp, method = "spearman")
# Kendall相关：适用于小样本或有序数据
cor_kendall <- cor(medical_data$sbp, medical_data$dbp, method = "kendall")
cat("Pearson相关:", cor_pearson, "\n")## Pearson相关: -0.06914093cat("Spearman相关:", cor_spearman, "\n")## Spearman相关: -0.009209166cat("Kendall相关:", cor_kendall, "\n")## Kendall相关: -0.01499386解读要点： - 相关系数只反映线性关系，不反映因果关系 - 相关系数受极端值影响，分析前应检查异常值 - Pearson相关假设数据近似正态分布，否则考虑Spearman相关
1.9 数据标准化与中心化
标准化和中心化是数据预处理的重要步骤，使不同量纲的变量具有可比性。
数学基础
中心化（Centering）： \[x_{centered} = x - \bar{x}\] 中心化后的数据均值为0。
标准化（Standardization/Z-score）： \[z = \frac{x - \bar{x}}{s}\] 标准化后的数据均值为0，标准差为1。
使用场景与意义
| 应用场景 | 说明 |
|---|---|
| 变量比较 | 消除量纲影响，使不同变量可比 |
| 回归分析 | 减少多重共线性，便于比较系数 |
| 聚类分析 | 防止大量纲变量主导距离计算 |
| 主成分分析 | 使各变量贡献相等 |
# 使用医学数据演示标准化
# 原始数据的均值和标准差
cat("原始数据:\n")## 原始数据:cat("年龄 - 均值:", mean(medical_data$age), "标准差:", sd(medical_data$age), "\n")## 年龄 - 均值: 53.3 标准差: 11.02178cat("BMI - 均值:", mean(medical_data$bmi), "标准差:", sd(medical_data$bmi), "\n")## BMI - 均值: 25.202 标准差: 3.087466cat("收缩压 - 均值:", mean(medical_data$sbp), "标准差:", sd(medical_data$sbp), "\n")## 收缩压 - 均值: 129.58 标准差: 15.55122# 使用scale()函数进行标准化
# scale()默认进行中心化和标准化
age_scaled <- scale(medical_data$age)
# 标准化后的均值和标准差
cat("\n标准化后:\n")##
## 标准化后:cat("年龄 - 均值:", mean(age_scaled), "标准差:", sd(age_scaled), "\n")## 年龄 - 均值: 2.339101e-16 标准差: 1# 对多个变量进行标准化
scaled_vars <- scale(medical_data[, c("age", "bmi", "sbp", "dbp", "glucose")])
head(scaled_vars)## age bmi sbp dbp glucose
## [1,] 0.6986169 -0.09781485 -0.4874216 -0.4494045 -3.0150991
## [2,] -2.2954554 -1.55532089 -0.3588143 -0.3549918 -1.5457525
## [3,] 0.1542401 0.71191073 0.6057403 0.3058972 -1.0559704
## [4,] 0.3356990 -0.19498192 0.1556148 1.4388496 2.0786355
## [5,] -0.2086778 -0.74559531 0.1556148 -0.6382299 -0.2723189
## [6,] -0.3901367 0.51757659 0.6057403 0.2114845 -0.6641446# 检查标准化结果
cat("\n标准化后各变量的均值:\n")##
## 标准化后各变量的均值:colMeans(scaled_vars)## age bmi sbp dbp glucose
## 2.339795e-16 -5.087597e-16 -7.920054e-16 -4.864165e-16 -1.386391e-16cat("\n标准化后各变量的标准差:\n")##
## 标准化后各变量的标准差:apply(scaled_vars, 2, sd)## age bmi sbp dbp glucose
## 1 1 1 1 1# 仅中心化（不减标准差）
# scale()的center和scale参数可单独控制
age_centered <- scale(medical_data$age, center = TRUE, scale = FALSE)
cat("\n仅中心化后的年龄均值:", mean(age_centered), "\n")##
## 仅中心化后的年龄均值: 2.842344e-15cat("中心化后的年龄标准差:", sd(age_centered), "\n")## 中心化后的年龄标准差: 11.02178# 仅标准化（不中心化）
# 通常不推荐，但某些情况下有用
age_scaled_only <- scale(medical_data$age, center = FALSE, scale = TRUE)
cat("\n仅标准化（不中心化）的均值:", mean(age_scaled_only), "\n")##
## 仅标准化（不中心化）的均值: 0.9698371# 将标准化结果添加到数据框
medical_data_scaled <- medical_data
medical_data_scaled$age_z <- scale(medical_data$age)
medical_data_scaled$bmi_z <- scale(medical_data$bmi)
medical_data_scaled$sbp_z <- scale(medical_data$sbp)
head(medical_data_scaled)## patient_id age bmi sbp dbp glucose group age_z bmi_z
## 1 1 61 24.9 122 79 2.3 Control 0.6986169 -0.09781485
## 2 2 28 20.4 124 80 3.8 Control -2.2954554 -1.55532089
## 3 3 55 27.4 139 87 4.3 Treatment 0.1542401 0.71191073
## 4 4 57 24.6 132 99 7.5 Treatment 0.3356990 -0.19498192
## 5 5 51 22.9 132 77 5.1 Control -0.2086778 -0.74559531
## 6 6 49 26.8 139 86 4.7 Treatment -0.3901367 0.51757659
## sbp_z
## 1 -0.4874216
## 2 -0.3588143
## 3 0.6057403
## 4 0.1556148
## 5 0.1556148
## 6 0.6057403# 标准化的实际应用示例
# 比较不同量纲变量的"异常程度"
# Z分数的绝对值越大，表示离均值越远
# 找出年龄最极端的5个样本
medical_data$age_z <- scale(medical_data$age)
medical_data %>%
dplyr::arrange(dplyr::desc(abs(age_z))) %>%
dplyr::select(patient_id, age, age_z) %>%
head(5)## patient_id age age_z
## 1 47 18 -3.202750
## 2 2 28 -2.295455
## 3 42 76 2.059559
## 4 43 74 1.878100
## 5 45 70 1.515182# 标准化在回归分析中的应用
# 当变量量纲差异很大时，标准化系数便于比较重要性
# 例如：年龄每增加1个标准差 vs BMI每增加1个标准差对血压的影响
# 创建标准化数据框用于回归
reg_data <- as.data.frame(scale(medical_data[, c("sbp", "age", "bmi")]))
head(reg_data)## sbp age bmi
## 1 -0.4874216 0.6986169 -0.09781485
## 2 -0.3588143 -2.2954554 -1.55532089
## 3 0.6057403 0.1542401 0.71191073
## 4 0.1556148 0.3356990 -0.19498192
## 5 0.1556148 -0.2086778 -0.74559531
## 6 0.6057403 -0.3901367 0.51757659解读要点： - Z分数表示数据点距离均值有多少个标准差 - |Z| > 2 或 |Z| > 3 通常被视为异常值 - 标准化不改变变量之间的相关关系
第一章小结
| 类别 | 主要指标 | R函数 | 适用场景 |
|---|---|---|---|
| 中心趋势 | 均值、中位数、众数、截尾均值 | mean(), median() | 描述数据典型值 |
| 离散程度 | 方差、标准差、IQR、CV | var(), sd(), IQR() | 描述数据分散程度 |
| 分布形状 | 偏度、峰度 | e1071::skewness(), kurtosis() | 判断正态性 |
| 分位数 | 百分位数、四分位数 | quantile() | 确定参考范围 |
| 汇总统计 | 多指标汇总 | summary(), describe(), stat.desc() | 快速了解数据 |
| 分组汇总 | 按组计算统计量 | tapply(), aggregate(), dplyr | 比较组间差异 |
| 频数表 | 频数、比例 | table(), prop.table() | 分类变量描述 |
| 相关分析 | 协方差、相关系数 | cov(), cor() | 变量关系描述 |
| 数据变换 | 标准化、中心化 | scale() | 数据预处理 |
第二章：概率分布与随机变量
概率分布是统计推断的理论基础。理解各种概率分布的特征和应用场景，对于正确选择统计方法至关重要。本章将介绍常见的离散分布和连续分布，以及如何在R中进行相关计算。
为什么学习概率分布？
- 统计推断的基础：假设检验和置信区间的推导依赖于概率分布
- 数据建模：了解数据可能服从的分布有助于选择合适的分析方法
- 模拟研究：生成符合特定分布的随机数据用于模拟和验证
- 风险评估：计算特定事件发生的概率
2.1 离散分布
离散分布描述离散随机变量的概率分布，即随机变量只能取有限或可数无限个值。
数学基础
概率质量函数（PMF）：离散随机变量X取值为x的概率 \[P(X = x) = f(x)\]
累积分布函数（CDF）：随机变量小于等于某值的概率 \[F(x) = P(X \leq x) = \sum_{t \leq x}f(t)\]
2.1.1 二项分布（Binomial Distribution）
定义：n次独立伯努利试验中成功次数X的分布
\[P(X = k) = \binom{n}{k}p^k(1-p)^{n-k}\]
其中： - n：试验次数 - p：每次试验成功的概率 - k：成功次数（k = 0, 1, 2, …, n）
期望与方差： \[E(X) = np, \quad Var(X) = np(1-p)\]
使用场景： - n次治疗中治愈的患者数 - n次试验中成功的次数 - 样本中某种疾病的患者数
# 二项分布示例：药物临床试验
# 假设某药物的有效率为60%（p = 0.6）
# 在20名患者中，有多少人会有效？
# 参数设置
n <- 20 # 试验次数（患者数）
p <- 0.6 # 成功概率（有效率）
# 计算恰好有k人有效的概率
# dbinom(x, size, prob) 计算概率质量函数
# 恰好有12人有效的概率
prob_12 <- dbinom(12, size = n, prob = p)
prob_12## [1] 0.1797058# 计算不同成功次数的概率分布
k_values <- 0:n
probs <- dbinom(k_values, size = n, prob = p)
# 可视化二项分布
plot(k_values, probs, type = "h", lwd = 3, col = "steelblue",
main = "二项分布概率质量函数\n(n=20, p=0.6)",
xlab = "有效人数", ylab = "概率")
points(k_values, probs, pch = 16, col = "steelblue")
# 计算累积概率
# pbinom(q, size, prob) 计算累积分布函数
# 最多有10人有效的概率
pbinom(10, size = n, prob = p)## [1] 0.2446628# 至少有15人有效的概率
1 - pbinom(14, size = n, prob = p)## [1] 0.125599# 计算分位数
# qbinom(p, size, prob) 计算分位数
# 95%的患者有效人数不超过多少？
qbinom(0.95, size = n, prob = p)## [1] 16# 生成随机数
# rbinom(n, size, prob) 生成随机数
# 模拟100次试验，每次20名患者
set.seed(123)
sim_results <- rbinom(100, size = n, prob = p)
head(sim_results)## [1] 13 10 13 9 9 16# 计算模拟结果的均值和理论期望
cat("模拟均值:", mean(sim_results), "\n")## 模拟均值: 11.99cat("理论期望:", n * p, "\n")## 理论期望: 122.1.2 泊松分布（Poisson Distribution）
定义：单位时间/空间内稀有事件发生次数的分布
\[P(X = k) = \frac{\lambda^k e^{-\lambda}}{k!}\]
其中λ是单位时间/空间内事件的平均发生次数。
期望与方差： \[E(X) = Var(X) = \lambda\]
使用场景： - 单位时间内到达医院的急诊患者数 - 单位面积内某种疾病的发病数 - DNA序列中突变的次数
# 泊松分布示例：急诊科患者到达
# 假设平均每小时有5名急诊患者到达
lambda <- 5 # 平均到达率
# 计算恰好有k名患者到达的概率
# dpois(x, lambda) 计算概率质量函数
# 恰好有3名患者到达的概率
dpois(3, lambda = lambda)## [1] 0.1403739# 计算不同患者数的概率分布
k_values <- 0:15
probs <- dpois(k_values, lambda = lambda)
# 可视化泊松分布
plot(k_values, probs, type = "h", lwd = 3, col = "darkgreen",
main = "泊松分布概率质量函数\n(λ=5)",
xlab = "患者数", ylab = "概率")
points(k_values, probs, pch = 16, col = "darkgreen")
# 计算累积概率
# ppois(q, lambda) 计算累积分布函数
# 最多有3名患者到达的概率
ppois(3, lambda = lambda)## [1] 0.2650259# 超过8名患者到达的概率
1 - ppois(8, lambda = lambda)## [1] 0.06809363# 生成随机数
set.seed(456)
poisson_sim <- rpois(100, lambda = lambda)
head(poisson_sim)## [1] 2 3 6 7 7 4# 比较模拟均值与理论期望
cat("模拟均值:", mean(poisson_sim), "\n")## 模拟均值: 5.3cat("理论期望:", lambda, "\n")## 理论期望: 52.1.3 几何分布与超几何分布
几何分布：首次成功所需的试验次数
\[P(X = k) = (1-p)^{k-1}p\]
超几何分布：不放回抽样中成功次数的分布
\[P(X = k) = \frac{\binom{K}{k}\binom{N-K}{n-k}}{\binom{N}{n}}\]
# 几何分布示例：首次成功所需的试验次数
# 假设治疗成功率为30%，首次成功需要多少次治疗？
p_success <- 0.3
# 计算第3次才成功的概率
# dgeom(x, prob) 中x是失败次数
# 第3次成功意味着前2次失败
dgeom(2, prob = p_success)## [1] 0.147# 可视化几何分布
k_values <- 0:15
probs <- dgeom(k_values, prob = p_success)
plot(k_values, probs, type = "h", lwd = 3, col = "purple",
main = "几何分布概率质量函数\n(p=0.3)",
xlab = "失败次数（首次成功前）", ylab = "概率")
# 超几何分布示例：不放回抽样
# 假设一批药品中有10件次品，随机抽取5件，其中有多少次品？
# 参数
m <- 10 # 总体中成功的数量（次品数）
n <- 40 # 总体中失败的数量（合格品数）
k <- 5 # 抽取数量
# 计算恰好有2件次品的概率
# dhyper(x, m, n, k)
dhyper(2, m = m, n = n, k = k)## [1] 0.2098397# 计算概率分布
x_values <- 0:min(5, m)
probs <- dhyper(x_values, m = m, n = n, k = k)
# 可视化
plot(x_values, probs, type = "h", lwd = 3, col = "orange",
main = "超几何分布概率质量函数\n(次品10件，合格40件，抽5件)",
xlab = "抽到的次品数", ylab = "概率")
2.2 连续分布
连续分布描述连续随机变量的概率分布，随机变量可以取某一区间内的任意值。
数学基础
概率密度函数（PDF）：描述连续随机变量在各点的密度 \[P(a \leq X \leq b) = \int_a^b f(x)dx\]
性质： - \(f(x) \geq 0\) - \(\int_{-\infty}^{+\infty}f(x)dx = 1\)
2.2.1 正态分布（Normal Distribution）
正态分布（又称高斯分布）是统计学中最重要的分布，许多自然现象都近似服从正态分布。
概率密度函数： \[f(x) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}\]
记号：\(X \sim N(\mu, \sigma^2)\)
标准正态分布：\(\mu = 0, \sigma = 1\)，记为\(Z \sim N(0,1)\)
使用场景： - 身高、体重等生理指标 - 测量误差 - 大样本均值的抽样分布
# 正态分布示例：成人血压分布
# 假设收缩压服从均值为120，标准差为15的正态分布
mu <- 120 # 均值
sigma <- 15 # 标准差
# 计算概率密度
# dnorm(x, mean, sd) 计算概率密度函数
x_values <- seq(60, 180, by = 1)
density <- dnorm(x_values, mean = mu, sd = sigma)
# 可视化正态分布
plot(x_values, density, type = "l", lwd = 2, col = "blue",
main = "正态分布概率密度函数\n(μ=120, σ=15)",
xlab = "收缩压 (mmHg)", ylab = "概率密度")
# 添加均值线
abline(v = mu, col = "red", lty = 2, lwd = 2)
legend("topright", legend = c("概率密度", "均值"),
col = c("blue", "red"), lty = c(1, 2), lwd = 2)
# 计算累积概率
# pnorm(q, mean, sd) 计算累积分布函数
# 收缩压低于100的概率
pnorm(100, mean = mu, sd = sigma)## [1] 0.09121122# 收缩压在100到140之间的概率
pnorm(140, mean = mu, sd = sigma) - pnorm(100, mean = mu, sd = sigma)## [1] 0.8175776# 收缩压高于150的概率
1 - pnorm(150, mean = mu, sd = sigma)## [1] 0.02275013# 计算分位数
# qnorm(p, mean, sd) 计算分位数
# 收缩压的第95百分位数
qnorm(0.95, mean = mu, sd = sigma)## [1] 144.6728# 计算参考范围（95%参考范围）
lower <- qnorm(0.025, mean = mu, sd = sigma)
upper <- qnorm(0.975, mean = mu, sd = sigma)
cat("95%参考范围:", lower, "-", upper, "\n")## 95%参考范围: 90.60054 - 149.3995# 生成随机数
set.seed(789)
bp_sim <- rnorm(1000, mean = mu, sd = sigma)
# 比较模拟结果与理论分布
hist(bp_sim, breaks = 30, probability = TRUE,
main = "模拟数据与理论分布比较",
xlab = "收缩压", col = "lightblue")
curve(dnorm(x, mean = mu, sd = sigma), add = TRUE, col = "red", lwd = 2)
legend("topright", legend = c("模拟数据", "理论分布"),
col = c("lightblue", "red"), lty = c(NA, 1), lwd = c(NA, 2),
pch = c(15, NA))
2.2.2 其他重要连续分布
# t分布：小样本推断
# t分布比正态分布有更厚的尾部，适用于小样本
# 自由度越大，t分布越接近正态分布
par(mfrow = c(2, 2))
# 比较不同自由度的t分布与标准正态分布
x <- seq(-4, 4, by = 0.01)
plot(x, dnorm(x), type = "l", lwd = 2, col = "black",
main = "t分布与正态分布比较", ylab = "密度", xlab = "x")
lines(x, dt(x, df = 1), col = "red", lwd = 2)
lines(x, dt(x, df = 5), col = "blue", lwd = 2)
lines(x, dt(x, df = 30), col = "green", lwd = 2)
legend("topright", legend = c("N(0,1)", "t(1)", "t(5)", "t(30)"),
col = c("black", "red", "blue", "green"), lwd = 2)
# 卡方分布：方差检验、拟合优度检验
# 自由度决定分布形状
x_chi <- seq(0, 20, by = 0.01)
plot(x_chi, dchisq(x_chi, df = 3), type = "l", lwd = 2, col = "red",
main = "卡方分布", ylab = "密度", xlab = "x",
xlim = c(0, 20), ylim = c(0, 0.3))
lines(x_chi, dchisq(x_chi, df = 5), col = "blue", lwd = 2)
lines(x_chi, dchisq(x_chi, df = 10), col = "green", lwd = 2)
legend("topright", legend = c("df=3", "df=5", "df=10"),
col = c("red", "blue", "green"), lwd = 2)
# F分布：方差分析、方差比较
# 两个自由度参数
x_f <- seq(0, 5, by = 0.01)
plot(x_f, df(x_f, df1 = 3, df2 = 10), type = "l", lwd = 2, col = "red",
main = "F分布", ylab = "密度", xlab = "x")
lines(x_f, df(x_f, df1 = 5, df2 = 10), col = "blue", lwd = 2)
lines(x_f, df(x_f, df1 = 10, df2 = 20), col = "green", lwd = 2)
legend("topright", legend = c("F(3,10)", "F(5,10)", "F(10,20)"),
col = c("red", "blue", "green"), lwd = 2)
# 均匀分布：随机数生成的基础
x_unif <- seq(0, 10, by = 0.01)
plot(x_unif, dunif(x_unif, min = 2, max = 8), type = "l", lwd = 2, col = "purple",
main = "均匀分布 U(2,8)", ylab = "密度", xlab = "x")
par(mfrow = c(1, 1))2.2.3 指数分布
概率密度函数： \[f(x) = \lambda e^{-\lambda x}, \quad x \geq 0\]
使用场景： - 事件发生的时间间隔 - 设备的寿命 - 患者等待时间
# 指数分布示例：患者等待时间
# 假设平均等待时间为10分钟（λ = 0.1）
lambda_exp <- 0.1 # 速率参数（1/平均时间）
# 计算概率密度
x_exp <- seq(0, 50, by = 0.1)
density_exp <- dexp(x_exp, rate = lambda_exp)
# 可视化
plot(x_exp, density_exp, type = "l", lwd = 2, col = "darkorange",
main = "指数分布概率密度函数\n(λ=0.1, 平均等待10分钟)",
xlab = "等待时间（分钟）", ylab = "概率密度")
# 计算等待时间少于5分钟的概率
pexp(5, rate = lambda_exp)## [1] 0.3934693# 等待时间超过20分钟的概率
1 - pexp(20, rate = lambda_exp)## [1] 0.1353353# 等待时间在5到15分钟之间的概率
pexp(15, rate = lambda_exp) - pexp(5, rate = lambda_exp)## [1] 0.3834005# 中位数等待时间
qexp(0.5, rate = lambda_exp)## [1] 6.9314722.3 概率分布函数总结
R中概率分布函数的命名规则：
| 前缀 | 功能 | 示例 |
|---|---|---|
| d | 概率密度/质量函数 (PDF/PMF) | dnorm(), dbinom() |
| p | 累积分布函数 (CDF) | pnorm(), pbinom() |
| q | 分位数函数 | qnorm(), qbinom() |
| r | 随机数生成 | rnorm(), rbinom() |
# 综合示例：演示四种函数类型
# 以标准正态分布为例
# d - 密度函数：计算某点的密度值
dnorm(0, mean = 0, sd = 1) # 均值处的密度最大## [1] 0.3989423# p - 分布函数：计算累积概率
pnorm(1.96, mean = 0, sd = 1) # Z=1.96以下的概率## [1] 0.9750021# q - 分位数函数：计算给定概率对应的值
qnorm(0.975, mean = 0, sd = 1) # 97.5%分位数## [1] 1.959964# r - 随机数：生成随机样本
set.seed(111)
random_normal <- rnorm(5, mean = 0, sd = 1)
random_normal## [1] 0.2352207 -0.3307359 -0.3116238 -2.3023457 -0.1708760# 常用分布函数汇总表
dist_summary <- data.frame(
分布 = c("正态分布", "t分布", "卡方分布", "F分布",
"二项分布", "泊松分布", "指数分布", "均匀分布"),
R函数 = c("norm", "t", "chisq", "f", "binom", "pois", "exp", "unif"),
主要参数 = c("mean, sd", "df", "df", "df1, df2",
"size, prob", "lambda", "rate", "min, max"),
应用场景 = c("连续变量、抽样分布", "小样本推断",
"方差检验、拟合优度", "方差分析、方差比较",
"成功次数", "计数数据", "等待时间", "随机数生成")
)
dist_summary## 分布 R函数 主要参数 应用场景
## 1 正态分布 norm mean, sd 连续变量、抽样分布
## 2 t分布 t df 小样本推断
## 3 卡方分布 chisq df 方差检验、拟合优度
## 4 F分布 f df1, df2 方差分析、方差比较
## 5 二项分布 binom size, prob 成功次数
## 6 泊松分布 pois lambda 计数数据
## 7 指数分布 exp rate 等待时间
## 8 均匀分布 unif min, max 随机数生成2.4 正态性检验
许多统计方法假设数据服从正态分布，因此正态性检验是统计分析的重要步骤。
数学基础
Shapiro-Wilk检验： - 原假设H₀：数据来自正态分布 - 备择假设H₁：数据不来自正态分布 - p值 < 0.05 拒绝原假设，认为数据非正态
# 创建正态和非正态数据
set.seed(222)
normal_data <- rnorm(100, mean = 50, sd = 10)
skewed_data <- rexp(100, rate = 0.1) # 指数分布，右偏
# Shapiro-Wilk检验
# 正态数据的检验
shapiro.test(normal_data)##
## Shapiro-Wilk normality test
##
## data: normal_data
## W = 0.97674, p-value = 0.07387# 偏态数据的检验
shapiro.test(skewed_data)##
## Shapiro-Wilk normality test
##
## data: skewed_data
## W = 0.86574, p-value = 4.738e-08# Kolmogorov-Smirnov检验
# 将数据标准化后与标准正态分布比较
ks.test(scale(normal_data), "pnorm", mean = 0, sd = 1)##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: scale(normal_data)
## D = 0.077204, p-value = 0.5902
## alternative hypothesis: two-sidedks.test(scale(skewed_data), "pnorm", mean = 0, sd = 1)##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: scale(skewed_data)
## D = 0.12809, p-value = 0.07513
## alternative hypothesis: two-sided# Q-Q图（分位数-分位数图）
# 正态数据的Q-Q图
par(mfrow = c(1, 2))
qqnorm(normal_data, main = "正态数据Q-Q图", col = "blue")
qqline(normal_data, col = "red", lwd = 2)
qqnorm(skewed_data, main = "偏态数据Q-Q图", col = "orange")
qqline(skewed_data, col = "red", lwd = 2)
par(mfrow = c(1, 1))
# 使用car包绘制更详细的Q-Q图
# 包含置信区间
car::qqPlot(normal_data, main = "正态数据Q-Q图（带置信区间）")
## [1] 70 74解读要点： - Q-Q图中点近似落在直线上表示正态 - 两端偏离直线表示有偏态或重尾 - 检验的p值 > 0.05 不能拒绝正态假设 - 样本量很大时，检验可能过于敏感，应结合图形判断
2.5 中心极限定理与抽样分布
数学基础
中心极限定理（CLT）： 设 \(X_1, X_2, ..., X_n\) 是独立同分布的随机变量，具有有限的均值μ和方差σ²，则当n足够大时： \[\bar{X} \sim N\left(\mu, \frac{\sigma^2}{n}\right)\]
意义：无论总体分布如何，样本均值的分布在大样本时近似正态分布。
抽样分布：统计量的概率分布，如样本均值的分布。
# 演示中心极限定理
# 从非正态分布（指数分布）抽样，观察样本均值的分布
set.seed(333)
# 总体：指数分布（明显右偏）
population <- rexp(100000, rate = 0.5)
# 绘制总体分布
par(mfrow = c(2, 2))
hist(population, breaks = 50, main = "总体分布（指数分布）",
xlab = "值", col = "lightgreen", probability = TRUE)
# 不同样本量下样本均值的分布
sample_sizes <- c(5, 30, 100)
n_simulations <- 1000 # 模拟次数
for (i in seq_along(sample_sizes)) {
n <- sample_sizes[i]
# 进行多次抽样，计算样本均值
sample_means <- replicate(n_simulations, {
sample_data <- sample(population, size = n)
mean(sample_data)
})
# 绘制样本均值的分布
hist(sample_means, breaks = 30,
main = paste("样本均值分布 (n =", n, ")"),
xlab = "样本均值", col = "lightblue", probability = TRUE)
# 叠加正态分布曲线
curve(dnorm(x, mean = mean(population),
sd = sd(population)/sqrt(n)),
add = TRUE, col = "red", lwd = 2)
}
par(mfrow = c(1, 1))
# 理论值与模拟值比较
cat("总体均值:", mean(population), "\n")## 总体均值: 1.992382cat("总体标准差:", sd(population), "\n")## 总体标准差: 1.998025cat("样本均值的理论标准误:", sd(population)/sqrt(30), "\n")## 样本均值的理论标准误: 0.3647878cat("模拟的样本均值标准差:", sd(replicate(1000, mean(sample(population, 30)))), "\n")## 模拟的样本均值标准差: 0.36811162.6 大数定律
数学基础
大数定律：当样本量n趋向无穷时，样本均值依概率收敛于总体均值。
\[\bar{X}_n \xrightarrow{P} \mu \quad \text{当} \quad n \to \infty\]
意义：样本量越大，样本统计量越接近总体参数。
# 演示大数定律
set.seed(444)
# 模拟掷骰子（期望值为3.5）
n_rolls <- 10000
rolls <- sample(1:6, n_rolls, replace = TRUE)
# 计算累积均值
cumulative_mean <- cumsum(rolls) / (1:n_rolls)
# 绘制累积均值随样本量的变化
plot(1:n_rolls, cumulative_mean, type = "l",
main = "大数定律演示：掷骰子",
xlab = "投掷次数", ylab = "累积平均",
col = "blue", lwd = 1)
abline(h = 3.5, col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("累积平均", "理论期望 (3.5)"),
col = c("blue", "red"), lty = c(1, 2), lwd = c(1, 2))
# 放大前500次
plot(1:500, cumulative_mean[1:500], type = "l",
main = "大数定律演示（前500次）",
xlab = "投掷次数", ylab = "累积平均",
col = "blue", lwd = 2, ylim = c(2, 5))
abline(h = 3.5, col = "red", lwd = 2, lty = 2)
第三章：参数估计
参数估计是统计推断的核心内容之一，它利用样本数据来推断总体参数。参数估计分为点估计和区间估计两种形式，本章将详细介绍这两种方法及其在R中的实现。
为什么需要参数估计？
- 总体参数未知：我们通常无法获得整个总体的数据，只能通过样本推断
- 科学决策：基于数据做出合理的决策需要了解参数的可能范围
- 研究结论的可靠性：区间估计提供了估计的不确定性信息
3.1 点估计
点估计是用样本统计量直接估计总体参数的方法。
3.1.1 矩估计法
数学基础：
矩估计的基本思想是用样本矩估计总体矩。设总体有k个未知参数，用前k阶样本矩等于相应的总体矩，解方程组得到参数估计。
一阶矩（均值）： \[\bar{X} = \frac{1}{n}\sum_{i=1}^{n}X_i \rightarrow E(X) = \mu\]
二阶矩（方差）： \[S^2 = \frac{1}{n}\sum_{i=1}^{n}(X_i - \bar{X})^2 \rightarrow Var(X) = \sigma^2\]
# 矩估计示例：正态分布参数估计
set.seed(555)
# 从已知参数的正态分布生成样本
true_mu <- 100
true_sigma <- 15
n <- 50
sample_data <- rnorm(n, mean = true_mu, sd = true_sigma)
# 矩估计：用样本均值估计总体均值
mu_moment <- mean(sample_data)
mu_moment## [1] 100.3589# 矩估计：用样本方差估计总体方差
# 注意：样本方差有两种计算方式
# 有偏估计（总体方差）：除以n
sigma2_biased <- mean((sample_data - mu_moment)^2)
# 无偏估计：除以n-1
sigma2_unbiased <- var(sample_data)
cat("真实均值:", true_mu, "\n")## 真实均值: 100cat("矩估计均值:", mu_moment, "\n")## 矩估计均值: 100.3589cat("真实方差:", true_sigma^2, "\n")## 真实方差: 225cat("矩估计方差（有偏）:", sigma2_biased, "\n")## 矩估计方差（有偏）: 222.5963cat("矩估计方差（无偏）:", sigma2_unbiased, "\n")## 矩估计方差（无偏）: 227.13913.1.2 极大似然估计（MLE）
数学基础：
极大似然估计寻找使观测数据出现概率最大的参数值。
似然函数： \[L(\theta) = \prod_{i=1}^{n}f(x_i|\theta)\]
对数似然函数： \[\ell(\theta) = \sum_{i=1}^{n}\ln f(x_i|\theta)\]
MLE估计量： \[\hat{\theta}_{MLE} = \arg\max_{\theta} L(\theta)\]
对于正态分布 \(X \sim N(\mu, \sigma^2)\)： - \(\hat{\mu}_{MLE} = \bar{X}\) - \(\hat{\sigma}^2_{MLE} = \frac{1}{n}\sum_{i=1}^{n}(X_i - \bar{X})^2\)
# 极大似然估计示例
# 方法1：手动计算正态分布的MLE
# 对于正态分布，均值的MLE就是样本均值
mu_mle <- mean(sample_data)
# 标准差的MLE（有偏版本）
sigma_mle <- sqrt(mean((sample_data - mu_mle)^2))
cat("MLE均值估计:", mu_mle, "\n")## MLE均值估计: 100.3589cat("MLE标准差估计:", sigma_mle, "\n")## MLE标准差估计: 14.91966# 方法2：使用MASS包的fitdistr函数进行MLE
# 自动计算多种分布的MLE
mle_result <- MASS::fitdistr(sample_data, "normal")
mle_result## mean sd
## 100.358878 14.919662
## ( 2.109959) ( 1.491966)# 方法3：自定义似然函数进行优化
# 定义正态分布的对数似然函数
log_likelihood <- function(params, data) {
mu <- params[1]
sigma <- params[2]
# 约束sigma > 0
if (sigma <= 0) return(-Inf)
# 计算对数似然
ll <- sum(dnorm(data, mean = mu, sd = sigma, log = TRUE))
return(ll)
}
# 使用optim函数最大化对数似然
# 初始值
init_params <- c(mean(sample_data), sd(sample_data))
# 优化（最大化似然等价于最小化负似然）
mle_optim <- optim(init_params, function(p) -log_likelihood(p, sample_data),
method = "L-BFGS-B", lower = c(-Inf, 0.001))
cat("\n自定义优化结果:\n")##
## 自定义优化结果:cat("均值:", mle_optim$par[1], "\n")## 均值: 100.3589cat("标准差:", mle_optim$par[2], "\n")## 标准差: 14.91966# 比较不同方法的结果
comparison <- data.frame(
方法 = c("真实值", "矩估计", "MLE(fitdistr)", "MLE(optim)"),
均值 = c(true_mu, mu_moment, mle_result$estimate[1], mle_optim$par[1]),
标准差 = c(true_sigma, sqrt(sigma2_unbiased), mle_result$estimate[2], mle_optim$par[2])
)
comparison## 方法 均值 标准差
## 1 真实值 100.0000 15.00000
## 2 矩估计 100.3589 15.07113
## 3 MLE(fitdistr) 100.3589 14.91966
## 4 MLE(optim) 100.3589 14.919663.2 估计量的评价标准
不同的估计方法可能得到不同的估计量，我们需要评价估计量的优劣。
数学基础
无偏性（Unbiasedness）： \[E(\hat{\theta}) = \theta\] 估计量的期望等于真实参数值。
有效性（Efficiency）： 在所有无偏估计量中，方差最小的估计量最有效。
一致性（Consistency）： \[\hat{\theta}_n \xrightarrow{P} \theta \quad \text{当} \quad n \to \infty\] 样本量增大时，估计量收敛于真实参数。
# 演示无偏性和一致性
set.seed(666)
# 比较样本方差的有偏和无偏估计
n_simulations <- 10000
true_var <- 225 # 真实方差 (15^2)
# 存储不同样本量下的估计结果
sample_sizes <- c(10, 30, 100, 500)
results <- data.frame()
for (n in sample_sizes) {
biased_estimates <- numeric(n_simulations)
unbiased_estimates <- numeric(n_simulations)
for (i in 1:n_simulations) {
sample_data <- rnorm(n, mean = 100, sd = 15)
# 有偏估计（除以n）
biased_estimates[i] <- mean((sample_data - mean(sample_data))^2)
# 无偏估计（除以n-1）
unbiased_estimates[i] <- var(sample_data)
}
results <- rbind(results, data.frame(
n = n,
方法 = "有偏估计",
均值 = mean(biased_estimates),
偏差 = mean(biased_estimates) - true_var,
方差 = var(biased_estimates)
))
results <- rbind(results, data.frame(
n = n,
方法 = "无偏估计",
均值 = mean(unbiased_estimates),
偏差 = mean(unbiased_estimates) - true_var,
方差 = var(unbiased_estimates)
))
}
# 查看结果
results## n 方法 均值 偏差 方差
## 1 10 有偏估计 202.0247 -22.97534650 8993.0528
## 2 10 无偏估计 224.4718 -0.52816278 11102.5343
## 3 30 有偏估计 217.2684 -7.73163421 3267.1957
## 4 30 无偏估计 224.7604 -0.23962160 3496.4044
## 5 100 有偏估计 222.6066 -2.39335395 1004.0881
## 6 100 无偏估计 224.8552 -0.14480197 1024.4752
## 7 500 有偏估计 224.4865 -0.51350947 203.2422
## 8 500 无偏估计 224.9364 -0.06363675 204.0576# 可视化：随着样本量增加，两种估计都趋于真值
library(ggplot2)
ggplot2::ggplot(results, ggplot2::aes(x = n, y = 均值, color = 方法)) +
ggplot2::geom_point(size = 3) +
ggplot2::geom_line() +
ggplot2::geom_hline(yintercept = true_var, linetype = "dashed", color = "red") +
ggplot2::labs(title = "方差估计的无偏性与一致性",
x = "样本量", y = "估计均值") +
ggplot2::theme_minimal() +
ggplot2::annotate("text", x = 300, y = true_var + 10, label = "真实方差", color = "red")
3.3 置信区间基本原理
数学基础
置信区间定义： 对于参数θ，如果统计量 \(\hat{\theta}_L\) 和 \(\hat{\theta}_U\) 满足： \[P(\hat{\theta}_L \leq \theta \leq \hat{\theta}_U) = 1 - \alpha\] 则称 \([\hat{\theta}_L, \hat{\theta}_U]\) 为θ的 \((1-\alpha)\) 置信区间。
正确理解置信区间： - 不是”参数有95%的概率落在区间内” - 而是”如果重复抽样多次，95%的区间会包含真实参数”
# 演示置信区间的含义
set.seed(777)
# 参数设置
true_mu <- 50
true_sigma <- 10
n <- 30
n_intervals <- 100 # 生成100个置信区间
# 存储结果
intervals <- matrix(NA, nrow = n_intervals, ncol = 3)
colnames(intervals) <- c("lower", "mean", "upper")
# 生成多个置信区间
for (i in 1:n_intervals) {
sample_data <- rnorm(n, mean = true_mu, sd = true_sigma)
sample_mean <- mean(sample_data)
sample_se <- sd(sample_data) / sqrt(n)
# 95%置信区间
intervals[i, ] <- c(
sample_mean - qt(0.975, df = n-1) * sample_se,
sample_mean,
sample_mean + qt(0.975, df = n-1) * sample_se
)
}
# 计算有多少区间包含真实均值
contains_true <- (intervals[, "lower"] <= true_mu) & (intervals[, "upper"] >= true_mu)
coverage_rate <- mean(contains_true)
cat("覆盖率:", coverage_rate * 100, "%\n")## 覆盖率: 94 %# 可视化前20个区间
par(mar = c(5, 4, 4, 2))
plot(1, type = "n", xlim = c(35, 65), ylim = c(0, 21),
main = "95%置信区间演示（前20个）",
xlab = "值", ylab = "区间编号")
abline(v = true_mu, col = "red", lwd = 2, lty = 2)
for (i in 1:20) {
color <- ifelse(intervals[i, "lower"] <= true_mu & intervals[i, "upper"] >= true_mu,
"blue", "red")
segments(intervals[i, "lower"], i, intervals[i, "upper"], i, col = color, lwd = 2)
points(intervals[i, "mean"], i, pch = 16, col = color)
}
legend("topright", legend = c("包含真值", "不包含真值", "真实均值"),
col = c("blue", "red", "red"), lty = c(1, 1, 2), lwd = c(2, 2, 2))
3.4 单个总体均值的置信区间
3.4.1 Z区间（总体方差已知或大样本）
数学公式： \[\bar{X} \pm z_{\alpha/2} \cdot \frac{\sigma}{\sqrt{n}}\]
适用条件： - 总体方差已知，或 - 大样本（n ≥ 30）可用样本标准差代替
# Z区间示例
# 假设已知总体标准差为10，估计均值的置信区间
set.seed(888)
sample_data <- rnorm(50, mean = 100, sd = 10)
n <- length(sample_data)
sample_mean <- mean(sample_data)
known_sigma <- 10 # 已知总体标准差
# 手动计算95%置信区间
z_critical <- qnorm(0.975)
se <- known_sigma / sqrt(n)
ci_lower <- sample_mean - z_critical * se
ci_upper <- sample_mean + z_critical * se
cat("样本均值:", sample_mean, "\n")## 样本均值: 99.40236cat("95%置信区间: [", ci_lower, ",", ci_upper, "]\n")## 95%置信区间: [ 96.63055 , 102.1742 ]# 大样本情况（用样本标准差）
sample_sd <- sd(sample_data)
se_sample <- sample_sd / sqrt(n)
ci_lower_large <- sample_mean - z_critical * se_sample
ci_upper_large <- sample_mean + z_critical * se_sample
cat("\n大样本近似（用样本标准差）:\n")##
## 大样本近似（用样本标准差）:cat("95%置信区间: [", ci_lower_large, ",", ci_upper_large, "]\n")## 95%置信区间: [ 96.3753 , 102.4294 ]3.4.2 t区间（总体方差未知，小样本）
数学公式： \[\bar{X} \pm t_{\alpha/2, n-1} \cdot \frac{s}{\sqrt{n}}\]
适用条件： - 总体方差未知 - 小样本（n < 30） - 总体近似正态分布
# t区间示例
set.seed(999)
small_sample <- rnorm(15, mean = 50, sd = 8)
n <- length(small_sample)
# 手动计算
sample_mean <- mean(small_sample)
sample_sd <- sd(small_sample)
t_critical <- qt(0.975, df = n - 1)
se <- sample_sd / sqrt(n)
ci_lower <- sample_mean - t_critical * se
ci_upper <- sample_mean + t_critical * se
cat("样本均值:", sample_mean, "\n")## 样本均值: 48.35829cat("样本标准差:", sample_sd, "\n")## 样本标准差: 7.76903cat("t临界值:", t_critical, "\n")## t临界值: 2.144787cat("95%置信区间: [", ci_lower, ",", ci_upper, "]\n")## 95%置信区间: [ 44.05595 , 52.66064 ]# 使用R内置函数t.test计算
t_result <- t.test(small_sample)
t_result##
## One Sample t-test
##
## data: small_sample
## t = 24.107, df = 14, p-value = 8.435e-13
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 44.05595 52.66064
## sample estimates:
## mean of x
## 48.35829# 提取置信区间
t_result$conf.int## [1] 44.05595 52.66064
## attr(,"conf.level")
## [1] 0.953.5 单个总体比例的置信区间
数学基础
正态近似法： \[\hat{p} \pm z_{\alpha/2}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}\]
适用条件：\(n\hat{p} \geq 5\) 且 \(n(1-\hat{p}) \geq 5\)
Wilson区间（更精确）： \[\frac{\hat{p} + \frac{z^2}{2n} \pm z\sqrt{\frac{\hat{p}(1-\hat{p})}{n} + \frac{z^2}{4n^2}}}{1 + \frac{z^2}{n}}\]
# 比例的置信区间示例
# 某药物临床试验：100名患者中75人有效
n <- 100
successes <- 75
p_hat <- successes / n
# 方法1：正态近似
z_critical <- qnorm(0.975)
se <- sqrt(p_hat * (1 - p_hat) / n)
ci_normal <- c(p_hat - z_critical * se, p_hat + z_critical * se)
cat("正态近似法:\n")## 正态近似法:cat("有效率:", p_hat, "\n")## 有效率: 0.75cat("95%置信区间: [", ci_normal[1], ",", ci_normal[2], "]\n")## 95%置信区间: [ 0.6651311 , 0.8348689 ]# 方法2：Wilson区间（推荐）
# 使用prop.test函数
prop_result <- prop.test(successes, n, correct = FALSE)
prop_result$conf.int## [1] 0.6569554 0.8245479
## attr(,"conf.level")
## [1] 0.95# 方法3：精确二项分布区间（Clopper-Pearson）
binom_result <- binom.test(successes, n)
binom_result$conf.int## [1] 0.6534475 0.8312203
## attr(,"conf.level")
## [1] 0.95# 比较三种方法
comparison_prop <- data.frame(
方法 = c("正态近似", "Wilson", "Clopper-Pearson"),
下限 = c(ci_normal[1], prop_result$conf.int[1], binom_result$conf.int[1]),
上限 = c(ci_normal[2], prop_result$conf.int[2], binom_result$conf.int[2])
)
comparison_prop## 方法 下限 上限
## 1 正态近似 0.6651311 0.8348689
## 2 Wilson 0.6569554 0.8245479
## 3 Clopper-Pearson 0.6534475 0.83122033.6 单个总体方差的置信区间
数学基础
基于卡方分布： \[\left[\frac{(n-1)s^2}{\chi^2_{\alpha/2, n-1}}, \frac{(n-1)s^2}{\chi^2_{1-\alpha/2, n-1}}\right]\]
# 方差的置信区间示例
set.seed(1111)
sample_data <- rnorm(30, mean = 100, sd = 15)
n <- length(sample_data)
sample_var <- var(sample_data)
# 手动计算
chi_lower <- qchisq(0.975, df = n - 1)
chi_upper <- qchisq(0.025, df = n - 1)
ci_var_lower <- (n - 1) * sample_var / chi_lower
ci_var_upper <- (n - 1) * sample_var / chi_upper
cat("样本方差:", sample_var, "\n")## 样本方差: 306.5221cat("方差95%置信区间: [", ci_var_lower, ",", ci_var_upper, "]\n")## 方差95%置信区间: [ 194.4159 , 553.9416 ]cat("标准差95%置信区间: [", sqrt(ci_var_lower), ",", sqrt(ci_var_upper), "]\n")## 标准差95%置信区间: [ 13.94331 , 23.53596 ]# 使用DescTools包的VarCI函数
DescTools::VarCI(sample_data, method = "classic")## var lwr.ci upr.ci
## 306.5221 194.4159 553.94163.7 两个总体均值差的置信区间
3.7.1 独立样本
方差相等时（ pooled variance ）： \[\bar{X}_1 - \bar{X}_2 \pm t_{\alpha/2, n_1+n_2-2} \cdot s_p\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}\]
其中 \(s_p^2 = \frac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1+n_2-2}\)
方差不等时（Welch校正）： \[\bar{X}_1 - \bar{X}_2 \pm t_{\alpha/2, \nu} \cdot \sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}}\]
# 两个独立样本均值差的置信区间
set.seed(2222)
group1 <- rnorm(25, mean = 100, sd = 15)
group2 <- rnorm(30, mean = 110, sd = 15)
# 使用t.test计算
# 方差相等假设
t_equal <- t.test(group1, group2, var.equal = TRUE)
t_equal##
## Two Sample t-test
##
## data: group1 and group2
## t = -0.92681, df = 53, p-value = 0.3582
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -11.612622 4.272472
## sample estimates:
## mean of x mean of y
## 104.0225 107.6926# 方差不等假设（Welch校正，默认）
t_welch <- t.test(group1, group2, var.equal = FALSE)
t_welch##
## Welch Two Sample t-test
##
## data: group1 and group2
## t = -0.89579, df = 41.033, p-value = 0.3756
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -11.943996 4.603847
## sample estimates:
## mean of x mean of y
## 104.0225 107.6926# 提取置信区间
cat("\n方差相等假设下的均值差置信区间:\n")##
## 方差相等假设下的均值差置信区间:t_equal$conf.int## [1] -11.612622 4.272472
## attr(,"conf.level")
## [1] 0.95cat("\n方差不等假设下的均值差置信区间:\n")##
## 方差不等假设下的均值差置信区间:t_welch$conf.int## [1] -11.943996 4.603847
## attr(,"conf.level")
## [1] 0.953.7.2 配对样本
数学公式： \[\bar{d} \pm t_{\alpha/2, n-1} \cdot \frac{s_d}{\sqrt{n}}\]
其中d为配对差值。
# 配对样本均值差的置信区间
set.seed(3333)
# 治疗前后的血压数据
before <- rnorm(20, mean = 140, sd = 10)
after <- before - rnorm(20, mean = 8, sd = 5) # 治疗后血压下降
# 使用t.test（配对）
t_paired <- t.test(before, after, paired = TRUE)
t_paired##
## Paired t-test
##
## data: before and after
## t = 10.167, df = 19, p-value = 4.026e-09
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## 7.06303 10.72480
## sample estimates:
## mean difference
## 8.893917# 手动计算
differences <- before - after
mean_diff <- mean(differences)
sd_diff <- sd(differences)
n <- length(differences)
t_crit <- qt(0.975, df = n - 1)
ci_diff <- mean_diff + c(-1, 1) * t_crit * sd_diff / sqrt(n)
cat("\n手动计算结果:\n")##
## 手动计算结果:cat("平均差值:", mean_diff, "\n")## 平均差值: 8.893917cat("差值95%置信区间: [", ci_diff[1], ",", ci_diff[2], "]\n")## 差值95%置信区间: [ 7.06303 , 10.7248 ]3.8 两个总体比例差的置信区间
数学基础
正态近似： \[(\hat{p}_1 - \hat{p}_2) \pm z_{\alpha/2}\sqrt{\frac{\hat{p}_1(1-\hat{p}_1)}{n_1} + \frac{\hat{p}_2(1-\hat{p}_2)}{n_2}}\]
# 两个比例差的置信区间
# 两组治疗效果比较
# 治疗组：80人有效，共100人
# 对照组：60人有效，共100人
# 使用prop.test
prop_comparison <- prop.test(c(80, 60), c(100, 100))
prop_comparison##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(80, 60) out of c(100, 100)
## X-squared = 8.5952, df = 1, p-value = 0.00337
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## 0.06604099 0.33395901
## sample estimates:
## prop 1 prop 2
## 0.8 0.6# 提取比例差置信区间
prop_comparison$conf.int## [1] 0.06604099 0.33395901
## attr(,"conf.level")
## [1] 0.95# 手动计算
p1 <- 80/100
p2 <- 60/100
n1 <- n2 <- 100
diff <- p1 - p2
se_diff <- sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)
z_crit <- qnorm(0.975)
ci_diff_prop <- diff + c(-1, 1) * z_crit * se_diff
cat("\n手动计算:\n")##
## 手动计算:cat("比例差:", diff, "\n")## 比例差: 0.2cat("95%置信区间: [", ci_diff_prop[1], ",", ci_diff_prop[2], "]\n")## 95%置信区间: [ 0.07604099 , 0.323959 ]3.9 自举法（Bootstrap）置信区间
数学基础
自举法是一种非参数方法，通过有放回重抽样来估计统计量的分布。
基本步骤： 1. 从原始样本中有放回抽取n个观测值（重抽样） 2. 计算感兴趣的统计量 3. 重复步骤1-2 B次（如1000次） 4. 用统计量的经验分布构建置信区间
适用场景： - 统计量的理论分布难以推导 - 小样本 - 复杂统计量（如中位数、相关系数）
# Bootstrap置信区间示例
set.seed(4444)
original_sample <- rnorm(30, mean = 50, sd = 10)
# 使用boot包进行bootstrap
# 定义统计量函数
mean_fun <- function(data, indices) {
return(mean(data[indices]))
}
# 进行bootstrap
boot_result <- boot::boot(data = original_sample,
statistic = mean_fun,
R = 1000)
boot_result##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot::boot(data = original_sample, statistic = mean_fun, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 51.37265 0.05731496 1.982868# 计算不同类型的置信区间
boot_ci <- boot::boot.ci(boot_result, type = c("norm", "perc", "bca"))
boot_ci## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot::boot.ci(boot.out = boot_result, type = c("norm", "perc",
## "bca"))
##
## Intervals :
## Level Normal Percentile BCa
## 95% (47.43, 55.20 ) (47.70, 55.29 ) (47.71, 55.30 )
## Calculations and Intervals on Original Scale# 可视化bootstrap分布
hist(boot_result$t, breaks = 30, main = "Bootstrap分布（均值）",
xlab = "均值", col = "lightblue", probability = TRUE)
abline(v = mean(original_sample), col = "red", lwd = 2, lty = 2)
abline(v = boot_ci$percent[4:5], col = "blue", lwd = 2)
legend("topright", legend = c("样本均值", "95% CI"),
col = c("red", "blue"), lty = c(2, 1), lwd = 2)
# Bootstrap用于中位数（理论分布复杂）
median_fun <- function(data, indices) {
return(median(data[indices]))
}
boot_median <- boot::boot(data = original_sample,
statistic = median_fun,
R = 1000)
boot_median_ci <- boot::boot.ci(boot_median, type = "perc")
cat("\n中位数的Bootstrap置信区间:\n")##
## 中位数的Bootstrap置信区间:boot_median_ci## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot::boot.ci(boot.out = boot_median, type = "perc")
##
## Intervals :
## Level Percentile
## 95% (46.07, 57.11 )
## Calculations and Intervals on Original Scale3.10 贝叶斯可信区间简介
数学基础
贝叶斯公式： \[P(\theta|data) = \frac{P(data|\theta)P(\theta)}{P(data)}\]
可信区间： 参数θ的95%可信区间表示：参数有95%的概率落在此区间内。
与置信区间的区别： - 置信区间：频率学派，针对抽样过程 - 可信区间：贝叶斯学派，针对参数本身
# 贝叶斯可信区间简单示例
# 使用共轭先验估计正态分布均值
# 假设：已知总体方差σ² = 100
# 先验：μ ~ N(μ₀, τ₀²)，其中μ₀ = 50, τ₀² = 25
# 数据
set.seed(5555)
sample_data <- rnorm(20, mean = 60, sd = 10)
n <- length(sample_data)
sample_mean <- mean(sample_data)
# 已知参数
sigma2 <- 100 # 总体方差
mu0 <- 50 # 先验均值
tau0_2 <- 25 # 先验方差
# 后验分布参数（共轭更新）
# 后验均值是先验均值和样本均值的加权平均
posterior_var <- 1 / (1/tau0_2 + n/sigma2)
posterior_mean <- posterior_var * (mu0/tau0_2 + n*sample_mean/sigma2)
cat("先验均值:", mu0, "\n")## 先验均值: 50cat("样本均值:", sample_mean, "\n")## 样本均值: 60.30519cat("后验均值:", posterior_mean, "\n")## 后验均值: 58.58766cat("后验标准差:", sqrt(posterior_var), "\n")## 后验标准差: 2.041241# 95%可信区间
bayesian_ci <- posterior_mean + c(-1, 1) * qnorm(0.975) * sqrt(posterior_var)
cat("95%贝叶斯可信区间: [", bayesian_ci[1], ",", bayesian_ci[2], "]\n")## 95%贝叶斯可信区间: [ 54.5869 , 62.58842 ]# 与频率学派置信区间比较
freq_ci <- sample_mean + c(-1, 1) * qt(0.975, df = n-1) * sd(sample_data)/sqrt(n)
cat("95%置信区间: [", freq_ci[1], ",", freq_ci[2], "]\n")## 95%置信区间: [ 56.26757 , 64.34282 ]# 可视化先验、似然和后验
x_range <- seq(40, 80, by = 0.1)
prior <- dnorm(x_range, mean = mu0, sd = sqrt(tau0_2))
likelihood <- dnorm(x_range, mean = sample_mean, sd = sqrt(sigma2/n))
posterior <- dnorm(x_range, mean = posterior_mean, sd = sqrt(posterior_var))
# 标准化便于比较
prior <- prior / max(prior)
likelihood <- likelihood / max(likelihood)
posterior <- posterior / max(posterior)
plot(x_range, prior, type = "l", col = "blue", lwd = 2,
main = "贝叶斯推断：先验、似然与后验",
xlab = "μ", ylab = "密度（标准化）")
lines(x_range, likelihood, col = "green", lwd = 2)
lines(x_range, posterior, col = "red", lwd = 2)
legend("topright", legend = c("先验", "似然", "后验"),
col = c("blue", "green", "red"), lwd = 2)
第四章：假设检验基础
假设检验是统计推断的核心方法之一，它帮助我们基于样本数据对总体参数做出判断。本章将介绍假设检验的基本概念、原理和实施步骤。
为什么需要假设检验？
- 科学决策：判断某种治疗是否有效
- 验证理论：检验理论预测是否符合实际数据
- 质量控制：判断产品是否符合标准
- 研究结论：为研究发现提供统计学依据
4.1 原假设与备择假设
数学基础
原假设（Null Hypothesis, H₀）： - 通常表示”无效应”或”无差异” - 是我们试图拒绝的假设 - 例如：\(H_0: \mu = \mu_0\)（总体均值等于某值）
备择假设（Alternative Hypothesis, H₁或Hₐ）： - 表示存在效应或差异 - 是我们试图支持的假设 - 例如：\(H_1: \mu \neq \mu_0\)（总体均值不等于某值）
假设检验的逻辑
假设检验采用”反证法”思想： 1. 假设H₀为真 2. 计算在H₀成立条件下，观测到当前数据（或更极端情况）的概率 3. 如果这个概率很小，则拒绝H₀
# 假设检验基本概念示例
# 场景：某药物声称能降低血压
# 原假设 H₀: 药物无效（血压变化均值 = 0）
# 备择假设 H₁: 药物有效（血压变化均值 ≠ 0）
set.seed(6666)
# 假设我们观测到以下血压变化数据
bp_change <- c(-8, -5, -12, -3, -7, -10, -6, -4, -9, -11,
-2, -8, -7, -5, -13, -6, -9, -4, -8, -10)
# 样本统计量
n <- length(bp_change)
sample_mean <- mean(bp_change)
sample_sd <- sd(bp_change)
cat("样本量:", n, "\n")## 样本量: 20cat("样本均值:", sample_mean, "\n")## 样本均值: -7.35cat("样本标准差:", sample_sd, "\n")## 样本标准差: 3.013566# 进行单样本t检验
# H₀: μ = 0 vs H₁: μ ≠ 0
t_test_result <- t.test(bp_change, mu = 0)
t_test_result##
## One Sample t-test
##
## data: bp_change
## t = -10.907, df = 19, p-value = 1.276e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -8.760392 -5.939608
## sample estimates:
## mean of x
## -7.35# 提取关键信息
cat("\n假设检验结果解读:\n")##
## 假设检验结果解读:cat("t统计量:", t_test_result$statistic, "\n")## t统计量: -10.90741cat("自由度:", t_test_result$parameter, "\n")## 自由度: 19cat("p值:", t_test_result$p.value, "\n")## p值: 1.276243e-09cat("结论:", ifelse(t_test_result$p.value < 0.05,
"拒绝H₀，药物有效",
"不能拒绝H₀，证据不足"), "\n")## 结论: 拒绝H₀，药物有效4.2 检验统计量、拒绝域与显著性水平
数学基础
检验统计量： 将样本数据转换为标准化的统计量，用于判断是否拒绝H₀。
对于单样本t检验： \[t = \frac{\bar{X} - \mu_0}{s/\sqrt{n}}\]
拒绝域： 当检验统计量落入的区域，我们拒绝H₀。
显著性水平（α）： 犯第一类错误的最大允许概率，通常取0.05或0.01。
临界值： 拒绝域的边界值，由显著性水平决定。
# 演示检验统计量、拒绝域和显著性水平
# 参数设置
alpha <- 0.05
df <- 19 # n-1 = 20-1
# 计算临界值（双侧检验）
t_critical <- qt(1 - alpha/2, df = df)
cat("显著性水平 α =", alpha, "\n")## 显著性水平 α = 0.05cat("双侧检验临界值: ±", t_critical, "\n")## 双侧检验临界值: ± 2.093024# 可视化拒绝域
x <- seq(-4, 4, by = 0.01)
y <- dt(x, df = df)
plot(x, y, type = "l", lwd = 2, col = "black",
main = "t分布与拒绝域（α = 0.05，双侧）",
xlab = "t值", ylab = "密度")
# 标记拒绝域
x_left <- seq(-4, -t_critical, by = 0.01)
x_right <- seq(t_critical, 4, by = 0.01)
polygon(c(-4, x_left, -t_critical),
c(0, dt(x_left, df = df), 0),
col = "red", border = NA)
polygon(c(t_critical, x_right, 4),
c(0, dt(x_right, df = df), 0),
col = "red", border = NA)
# 添加临界值线
abline(v = c(-t_critical, t_critical), col = "red", lty = 2, lwd = 2)
abline(v = 0, col = "gray", lty = 3)
# 标记观测的t值
observed_t <- t_test_result$statistic
abline(v = observed_t, col = "blue", lwd = 2)
points(observed_t, 0, pch = 18, col = "blue", cex = 2)
legend("topright",
legend = c("拒绝域", "观测t值", "临界值"),
col = c("red", "blue", "red"),
lty = c(NA, 1, 2), pch = c(15, 18, NA), lwd = c(NA, 2, 2))
# 判断是否拒绝H₀
cat("\n观测到的t统计量:", observed_t, "\n")##
## 观测到的t统计量: -10.90741cat("是否在拒绝域内:", ifelse(abs(observed_t) > t_critical, "是，拒绝H₀", "否，不能拒绝H₀"), "\n")## 是否在拒绝域内: 是，拒绝H₀4.3 p值的概念与解释
数学基础
p值定义： 在H₀为真的条件下，观测到当前统计量或更极端值的概率。
\[p = P(T \geq |t_{obs}| | H_0)\]
p值的解释： - p值越小，反对H₀的证据越强 - p < α：拒绝H₀ - p ≥ α：不能拒绝H₀
重要提醒： - p值不是H₀为真的概率 - p值不是效应大小的度量 - 统计显著 ≠ 实际重要
# p值计算与可视化
# 使用之前的t检验结果
observed_t <- as.numeric(t_test_result$statistic)
df <- 19
# 手动计算p值（双侧）
p_value <- 2 * pt(-abs(observed_t), df = df)
cat("手动计算的p值:", p_value, "\n")## 手动计算的p值: 1.276243e-09cat("t.test输出的p值:", t_test_result$p.value, "\n")## t.test输出的p值: 1.276243e-09# 可视化p值
# 确保x轴范围足够大以包含观测t值
x_range <- max(4, abs(observed_t) + 1)
x <- seq(-x_range, x_range, by = 0.01)
y <- dt(x, df = df)
plot(x, y, type = "l", lwd = 2, col = "black",
main = paste("p值可视化 (p =", round(p_value, 6), ")"),
xlab = "t值", ylab = "密度")
# 标记p值对应的区域
# 使用正确的seq方向
t_crit <- abs(observed_t)
if (t_crit < x_range) {
x_p_left <- seq(-x_range, -t_crit, length.out = 100)
x_p_right <- seq(t_crit, x_range, length.out = 100)
polygon(c(-x_range, x_p_left, -t_crit),
c(0, dt(x_p_left, df = df), 0),
col = "orange", border = NA)
polygon(c(t_crit, x_p_right, x_range),
c(0, dt(x_p_right, df = df), 0),
col = "orange", border = NA)
}
# 添加观测值线
abline(v = c(-t_crit, t_crit), col = "blue", lwd = 2, lty = 2)
abline(v = 0, col = "gray", lty = 3)
legend("topright",
legend = c("p值区域", "观测t值"),
col = c("orange", "blue"),
lty = c(NA, 2), pch = c(15, NA), lwd = c(NA, 2))
# p值与显著性水平比较
cat("\np值与显著性水平的比较:\n")##
## p值与显著性水平的比较:cat("p值 =", p_value, "\n")## p值 = 1.276243e-09cat("α = 0.05\n")## α = 0.05cat("结论:", ifelse(p_value < 0.05, "p < α，拒绝H₀", "p ≥ α，不能拒绝H₀"), "\n")## 结论: p < α，拒绝H₀4.4 第I类错误与第II类错误
数学基础
| H₀为真 | H₀为假 | |
|---|---|---|
| 不拒绝H₀ | 正确决策 (1-α) | 第II类错误 (β) |
| 拒绝H₀ | 第I类错误 (α) | 正确决策 (1-β) |
第I类错误（Type I Error）： - 当H₀为真时错误地拒绝H₀ - 概率 = 显著性水平α - 也称为”假阳性”
第II类错误（Type II Error）： - 当H₀为假时未能拒绝H₀ - 概率 = β - 也称为”假阴性”
检验功效（Power）： - 正确拒绝错误H₀的概率 - Power = 1 - β
# 演示两类错误
set.seed(7777)
# 设置参数
n <- 30
true_mu_null <- 100 # H₀下的均值
true_mu_alt <- 105 # 实际均值（H₁为真）
sigma <- 15
alpha <- 0.05
n_sim <- 1000
# 存储结果
reject_when_H0_true <- numeric(n_sim) # H₀为真时的决策
reject_when_H1_true <- numeric(n_sim) # H₁为真时的决策
for (i in 1:n_sim) {
# H₀为真时的样本
sample_H0 <- rnorm(n, mean = true_mu_null, sd = sigma)
t_result_H0 <- t.test(sample_H0, mu = true_mu_null)
reject_when_H0_true[i] <- ifelse(t_result_H0$p.value < alpha, 1, 0)
# H₁为真时的样本
sample_H1 <- rnorm(n, mean = true_mu_alt, sd = sigma)
t_result_H1 <- t.test(sample_H1, mu = true_mu_null)
reject_when_H1_true[i] <- ifelse(t_result_H1$p.value < alpha, 1, 0)
}
# 计算各类错误率
type1_error_rate <- mean(reject_when_H0_true)
type2_error_rate <- 1 - mean(reject_when_H1_true)
power <- mean(reject_when_H1_true)
cat("模拟结果（1000次）:\n")## 模拟结果（1000次）:cat("第I类错误率（假阳性）:", type1_error_rate, "\n")## 第I类错误率（假阳性）: 0.044cat("理论显著性水平 α:", alpha, "\n")## 理论显著性水平 α: 0.05cat("第II类错误率（假阴性）:", type2_error_rate, "\n")## 第II类错误率（假阴性）: 0.583cat("检验功效:", power, "\n")## 检验功效: 0.417# 可视化两类错误
par(mfrow = c(1, 2))
# H₀分布
x <- seq(70, 130, by = 0.1)
y_null <- dnorm(x, mean = true_mu_null, sd = sigma/sqrt(n))
plot(x, y_null, type = "l", lwd = 2, col = "blue",
main = "H₀为真时的分布",
xlab = "样本均值", ylab = "密度")
critical_value <- true_mu_null + qt(0.975, df = n-1) * sigma/sqrt(n)
x_reject <- seq(critical_value, 130, by = 0.1)
polygon(c(critical_value, x_reject, 130),
c(0, dnorm(x_reject, mean = true_mu_null, sd = sigma/sqrt(n)), 0),
col = "red", border = NA)
abline(v = critical_value, col = "red", lty = 2)
legend("topright", legend = c("H₀分布", "拒绝域"),
col = c("blue", "red"), lty = c(1, NA), pch = c(NA, 15))
# H₁分布
y_alt <- dnorm(x, mean = true_mu_alt, sd = sigma/sqrt(n))
plot(x, y_alt, type = "l", lwd = 2, col = "green",
main = "H₁为真时的分布",
xlab = "样本均值", ylab = "密度")
x_accept <- seq(70, critical_value, by = 0.1)
polygon(c(70, x_accept, critical_value),
c(0, dnorm(x_accept, mean = true_mu_alt, sd = sigma/sqrt(n)), 0),
col = "orange", border = NA)
abline(v = critical_value, col = "red", lty = 2)
legend("topright", legend = c("H₁分布", "接受域"),
col = c("green", "orange"), lty = c(1, NA), pch = c(NA, 15))
par(mfrow = c(1, 1))4.5 检验功效与样本量关系
数学基础
功效与样本量的关系： - 样本量越大，检验功效越高 - 效应量越大，检验功效越高 - 显著性水平越宽松（α越大），功效越高
功效分析公式（单样本t检验）： \[n = \left(\frac{(z_{1-\alpha/2} + z_{1-\beta})\sigma}{\delta}\right)^2\]
其中δ是期望检测的差异。
# 功效与样本量关系演示
# 使用pwr包进行功效分析
# 单样本t检验
# 参数设置
effect_size <- 0.5 # 中等效应量 (Cohen's d)
alpha <- 0.05
power <- 0.8
# 计算所需样本量
sample_size_needed <- pwr::pwr.t.test(d = effect_size,
sig.level = alpha,
power = power,
type = "one.sample")
sample_size_needed##
## One-sample t test power calculation
##
## n = 33.36713
## d = 0.5
## sig.level = 0.05
## power = 0.8
## alternative = two.sided# 不同样本量下的功效
sample_sizes <- seq(10, 100, by = 5)
power_values <- sapply(sample_sizes, function(n) {
result <- pwr::pwr.t.test(d = effect_size,
sig.level = alpha,
n = n,
type = "one.sample")
return(result$power)
})
# 可视化功效曲线
plot(sample_sizes, power_values, type = "b", pch = 16, col = "blue",
main = "功效曲线（效应量d = 0.5）",
xlab = "样本量", ylab = "检验功效")
abline(h = 0.8, col = "red", lty = 2, lwd = 2)
abline(v = sample_size_needed$n, col = "green", lty = 2, lwd = 2)
legend("bottomright",
legend = c("功效曲线", "目标功效(0.8)", paste("所需样本量:", round(sample_size_needed$n))),
col = c("blue", "red", "green"), lty = c(1, 2, 2), pch = c(16, NA, NA))
# 不同效应量下的功效曲线
effect_sizes <- c(0.2, 0.5, 0.8) # 小、中、大效应量
colors <- c("red", "blue", "green")
labels <- c("小效应(d=0.2)", "中效应(d=0.5)", "大效应(d=0.8)")
plot(NULL, xlim = c(10, 200), ylim = c(0, 1),
main = "不同效应量下的功效曲线",
xlab = "样本量", ylab = "检验功效")
for (i in seq_along(effect_sizes)) {
power_vals <- sapply(sample_sizes, function(n) {
result <- pwr::pwr.t.test(d = effect_sizes[i],
sig.level = alpha,
n = n,
type = "one.sample")
return(result$power)
})
lines(seq(10, 100, by = 5), power_vals, col = colors[i], lwd = 2)
}
abline(h = 0.8, col = "gray", lty = 2)
legend("bottomright", legend = labels, col = colors, lwd = 2)
4.6 单侧检验与双侧检验
数学基础
双侧检验： - H₁: μ ≠ μ₀ - 拒绝域在分布两端 - 适用于不确定效应方向的情况
单侧检验： - H₁: μ > μ₀ 或 H₁: μ < μ₀ - 拒绝域在分布一端 - 适用于有明确方向假设的情况
# 单侧与双侧检验比较
set.seed(8888)
sample_data <- rnorm(25, mean = 105, sd = 15)
mu_null <- 100
# 双侧检验
t_two_sided <- t.test(sample_data, mu = mu_null, alternative = "two.sided")
cat("双侧检验:\n")## 双侧检验:cat("H₀: μ =", mu_null, "vs H₁: μ ≠", mu_null, "\n")## H₀: μ = 100 vs H₁: μ ≠ 100cat("t统计量:", t_two_sided$statistic, "\n")## t统计量: 2.787575cat("p值:", t_two_sided$p.value, "\n")## p值: 0.01021967# 单侧检验（大于）
t_greater <- t.test(sample_data, mu = mu_null, alternative = "greater")
cat("\n单侧检验（大于）:\n")##
## 单侧检验（大于）:cat("H₀: μ =", mu_null, "vs H₁: μ >", mu_null, "\n")## H₀: μ = 100 vs H₁: μ > 100cat("t统计量:", t_greater$statistic, "\n")## t统计量: 2.787575cat("p值:", t_greater$p.value, "\n")## p值: 0.005109834# 单侧检验（小于）
t_less <- t.test(sample_data, mu = mu_null, alternative = "less")
cat("\n单侧检验（小于）:\n")##
## 单侧检验（小于）:cat("H₀: μ =", mu_null, "vs H₁: μ <", mu_null, "\n")## H₀: μ = 100 vs H₁: μ < 100cat("t统计量:", t_less$statistic, "\n")## t统计量: 2.787575cat("p值:", t_less$p.value, "\n")## p值: 0.9948902# 可视化比较
par(mfrow = c(1, 2))
# 双侧检验
x <- seq(-4, 4, by = 0.01)
y <- dt(x, df = 24)
plot(x, y, type = "l", lwd = 2, main = "双侧检验",
xlab = "t值", ylab = "密度")
t_crit_two <- qt(0.975, df = 24)
x_left <- seq(-4, -t_crit_two, by = 0.01)
x_right <- seq(t_crit_two, 4, by = 0.01)
polygon(c(-4, x_left, -t_crit_two), c(0, dt(x_left, 24), 0), col = "red")
polygon(c(t_crit_two, x_right, 4), c(0, dt(x_right, 24), 0), col = "red")
abline(v = t_two_sided$statistic, col = "blue", lwd = 2)
# 单侧检验
plot(x, y, type = "l", lwd = 2, main = "单侧检验（大于）",
xlab = "t值", ylab = "密度")
t_crit_one <- qt(0.95, df = 24)
x_right <- seq(t_crit_one, 4, by = 0.01)
polygon(c(t_crit_one, x_right, 4), c(0, dt(x_right, 24), 0), col = "red")
abline(v = t_greater$statistic, col = "blue", lwd = 2)
par(mfrow = c(1, 1))4.7 效应量
数学基础
效应量衡量效应的实际大小，独立于样本量。常用的效应量包括：
Cohen’s d（均值差异）： \[d = \frac{\bar{X}_1 - \bar{X}_2}{s_{pooled}}\]
Cramér’s V（分类变量关联）： \[V = \sqrt{\frac{\chi^2}{n(k-1)}}\]
η²（Eta-squared）（方差分析）： \[\eta^2 = \frac{SS_{between}}{SS_{total}}\]
效应量解释标准
| Cohen’s d | 效应大小 |
|---|---|
| 0.2 | 小 |
| 0.5 | 中 |
| 0.8 | 大 |
# 效应量计算示例
# Cohen's d 计算
# 单样本
sample_mean <- mean(sample_data)
sample_sd <- sd(sample_data)
mu_null <- 100
cohens_d <- (sample_mean - mu_null) / sample_sd
cat("Cohen's d (单样本):", cohens_d, "\n")## Cohen's d (单样本): 0.557515# 使用effectsize包（如果可用）或手动计算
# 双样本Cohen's d
set.seed(9999)
group1 <- rnorm(30, mean = 100, sd = 15)
group2 <- rnorm(30, mean = 110, sd = 15)
# 计算合并标准差
n1 <- length(group1)
n2 <- length(group2)
s1 <- sd(group1)
s2 <- sd(group2)
s_pooled <- sqrt(((n1-1)*s1^2 + (n2-1)*s2^2) / (n1 + n2 - 2))
# Cohen's d
d <- (mean(group2) - mean(group1)) / s_pooled
cat("Cohen's d (双样本):", d, "\n")## Cohen's d (双样本): 0.127303# 效应量解释
if (abs(d) < 0.2) {
effect_interpretation <- "极小效应"
} else if (abs(d) < 0.5) {
effect_interpretation <- "小效应"
} else if (abs(d) < 0.8) {
effect_interpretation <- "中等效应"
} else {
effect_interpretation <- "大效应"
}
cat("效应量解释:", effect_interpretation, "\n")## 效应量解释: 极小效应# 使用DescTools包计算效应量
# Cohen's d
DescTools::CohenD(group1, group2)## [1] -0.127303
## attr(,"magnitude")
## [1] "negligible"# 相关性效应量（r）
# r = d / sqrt(d^2 + 4)
r_effect <- d / sqrt(d^2 + 4)
cat("\n相关性效应量 r:", r_effect, "\n")##
## 相关性效应量 r: 0.06352294# 分类变量的效应量（Cramér's V）
# 创建列联表
contingency_table <- matrix(c(30, 10, 15, 45), nrow = 2)
rownames(contingency_table) <- c("治疗组", "对照组")
colnames(contingency_table) <- c("有效", "无效")
contingency_table## 有效 无效
## 治疗组 30 15
## 对照组 10 45# 卡方检验
chi_test <- chisq.test(contingency_table)
chi_test##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: contingency_table
## X-squared = 22.264, df = 1, p-value = 2.376e-06# 计算Cramér's V
n_total <- sum(contingency_table)
min_dim <- min(dim(contingency_table)) - 1
cramers_v <- sqrt(chi_test$statistic / (n_total * min_dim))
cat("Cramér's V:", cramers_v, "\n")## Cramér's V: 0.4718507# 使用vcd包计算
vcd::assocstats(contingency_table)## X^2 df P(> X^2)
## Likelihood Ratio 25.161 1 5.2745e-07
## Pearson 24.242 1 8.4941e-07
##
## Phi-Coefficient : 0.492
## Contingency Coeff.: 0.442
## Cramer's V : 0.4924.8 多重比较问题与校正
数学基础
问题：进行多次检验时，至少犯一次第I类错误的概率增加。
族错误率（Family-wise Error Rate, FWER）： \[FWER = 1 - (1-\alpha)^m\] 其中m是检验次数。
校正方法：
Bonferroni校正： \[\alpha_{adjusted} = \frac{\alpha}{m}\]
FDR（False Discovery Rate）： 控制错误发现的比例，更宽松但更常用。
# 多重比较问题演示
set.seed(11111)
# 模拟场景：检验5个变量是否与结果相关
# 假设所有变量实际上都与结果无关（H₀都为真）
n_tests <- 5
n_sim <- 1000
alpha <- 0.05
# 存储每次模拟中至少拒绝一次的比例
any_rejection <- numeric(n_sim)
for (i in 1:n_sim) {
# 生成5个独立的随机变量
p_values <- numeric(n_tests)
for (j in 1:n_tests) {
# 生成随机数据（H₀为真）
x <- rnorm(50)
y <- rnorm(50)
# 计算相关性检验的p值
p_values[j] <- cor.test(x, y)$p.value
}
# 是否至少有一个显著
any_rejection[i] <- ifelse(any(p_values < alpha), 1, 0)
}
# 计算实际犯第I类错误的概率
fwer_simulated <- mean(any_rejection)
fwer_theoretical <- 1 - (1 - alpha)^n_tests
cat("模拟的族错误率:", fwer_simulated, "\n")## 模拟的族错误率: 0.225cat("理论族错误率:", fwer_theoretical, "\n")## 理论族错误率: 0.2262191cat("单次检验显著性水平:", alpha, "\n")## 单次检验显著性水平: 0.05# 多重比较校正示例
# 假设我们进行了5次检验，得到以下p值
observed_pvalues <- c(0.01, 0.03, 0.04, 0.12, 0.25)
# Bonferroni校正
bonferroni_adjusted <- p.adjust(observed_pvalues, method = "bonferroni")
cat("\nBonferroni校正:\n")##
## Bonferroni校正:cat("原始p值:", observed_pvalues, "\n")## 原始p值: 0.01 0.03 0.04 0.12 0.25cat("校正后p值:", bonferroni_adjusted, "\n")## 校正后p值: 0.05 0.15 0.2 0.6 1# FDR校正（Benjamini-Hochberg）
fdr_adjusted <- p.adjust(observed_pvalues, method = "BH")
cat("\nFDR校正:\n")##
## FDR校正:cat("原始p值:", observed_pvalues, "\n")## 原始p值: 0.01 0.03 0.04 0.12 0.25cat("校正后p值:", fdr_adjusted, "\n")## 校正后p值: 0.05 0.06666667 0.06666667 0.15 0.25# 比较表格
comparison_table <- data.frame(
检验 = paste("检验", 1:5),
原始p值 = observed_pvalues,
Bonferroni = bonferroni_adjusted,
FDR = fdr_adjusted,
Bonferroni显著 = ifelse(bonferroni_adjusted < 0.05, "是", "否"),
FDR显著 = ifelse(fdr_adjusted < 0.05, "是", "否")
)
comparison_table## 检验 原始p值 Bonferroni FDR Bonferroni显著 FDR显著
## 1 检验 1 0.01 0.05 0.05000000 否 否
## 2 检验 2 0.03 0.15 0.06666667 否 否
## 3 检验 3 0.04 0.20 0.06666667 否 否
## 4 检验 4 0.12 0.60 0.15000000 否 否
## 5 检验 5 0.25 1.00 0.25000000 否 否第五章：常见参数检验
参数检验是在总体分布形式已知（通常假设正态分布）的条件下进行的假设检验。本章将介绍最常用的参数检验方法及其在R中的实现。
参数检验的前提条件： 1. 数据来自正态分布（或样本量足够大） 2. 样本是独立随机样本 3. 方差齐性（对于多组比较）
5.1 单样本t检验
数学基础
检验样本均值是否与已知的总体均值有显著差异。
检验统计量： \[t = \frac{\bar{X} - \mu_0}{s/\sqrt{n}}\]
假设： - H₀: μ = μ₀ - H₁: μ ≠ μ₀（双侧）或 μ > μ₀ / μ < μ₀（单侧）
适用条件： - 样本来自正态分布总体 - 或样本量足够大（n ≥ 30）
# 单样本t检验示例
# 场景：检验某地区成人平均收缩压是否为120 mmHg
set.seed(12121)
# 生成样本数据
sbp_data <- rnorm(35, mean = 125, sd = 12)
# 检查正态性
par(mfrow = c(1, 2))
hist(sbp_data, breaks = 10, main = "收缩压分布", xlab = "收缩压", col = "lightblue")
qqnorm(sbp_data, main = "Q-Q图")
qqline(sbp_data)
par(mfrow = c(1, 1))
# Shapiro-Wilk正态性检验
shapiro.test(sbp_data)##
## Shapiro-Wilk normality test
##
## data: sbp_data
## W = 0.96795, p-value = 0.3894# 单样本t检验
# H₀: μ = 120 vs H₁: μ ≠ 120
t_test_one <- t.test(sbp_data, mu = 120)
t_test_one##
## One Sample t-test
##
## data: sbp_data
## t = 0.33033, df = 34, p-value = 0.7432
## alternative hypothesis: true mean is not equal to 120
## 95 percent confidence interval:
## 116.5646 124.7689
## sample estimates:
## mean of x
## 120.6668# 提取关键信息
cat("\n检验结果解读:\n")##
## 检验结果解读:cat("样本均值:", mean(sbp_data), "\n")## 样本均值: 120.6668cat("t统计量:", t_test_one$statistic, "\n")## t统计量: 0.3303258cat("自由度:", t_test_one$parameter, "\n")## 自由度: 34cat("p值:", t_test_one$p.value, "\n")## p值: 0.7431807cat("95%置信区间:", t_test_one$conf.int, "\n")## 95%置信区间: 116.5646 124.7689cat("结论:", ifelse(t_test_one$p.value < 0.05,
"拒绝H₀，收缩压均值显著不同于120",
"不能拒绝H₀"), "\n")## 结论: 不能拒绝H₀# 效应量（Cohen's d）
cohens_d <- (mean(sbp_data) - 120) / sd(sbp_data)
cat("Cohen's d:", cohens_d, "\n")## Cohen's d: 0.055835255.2 独立样本t检验
数学基础
检验两个独立样本的均值是否有显著差异。
方差相等时（Student’s t-test）： \[t = \frac{\bar{X}_1 - \bar{X}_2}{s_p\sqrt{\frac{1}{n_1} + \frac{1}{n_2}}}\]
其中 \(s_p = \sqrt{\frac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1 + n_2 - 2}}\)
方差不等时（Welch’s t-test）： \[t = \frac{\bar{X}_1 - \bar{X}_2}{\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}}\]
# 独立样本t检验示例
# 场景：比较新药组和安慰剂组的血压降低值
set.seed(22222)
# 生成数据
drug_group <- rnorm(25, mean = 15, sd = 6) # 新药组血压降低值
placebo_group <- rnorm(25, mean = 8, sd = 6) # 安慰剂组血压降低值
# 创建数据框
bp_data <- data.frame(
reduction = c(drug_group, placebo_group),
group = factor(rep(c("Drug", "Placebo"), each = 25))
)
# 描述性统计
bp_data %>%
dplyr::group_by(group) %>%
dplyr::summarise(
n = dplyr::n(),
mean = mean(reduction),
sd = sd(reduction),
min = min(reduction),
max = max(reduction)
)## # A tibble: 2 × 6
## group n mean sd min max
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Drug 25 14.2 6.28 5.06 26.3
## 2 Placebo 25 8.98 5.60 1.77 25.0# 可视化
ggplot2::ggplot(bp_data, ggplot2::aes(x = group, y = reduction, fill = group)) +
ggplot2::geom_boxplot(alpha = 0.7) +
ggplot2::geom_jitter(width = 0.2, alpha = 0.5) +
ggplot2::labs(title = "两组血压降低值比较",
x = "组别", y = "血压降低值 (mmHg)") +
ggplot2::theme_minimal()
# 检查正态性（每组）
by(bp_data$reduction, bp_data$group, shapiro.test)## bp_data$group: Drug
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.92617, p-value = 0.07094
##
## ------------------------------------------------------------
## bp_data$group: Placebo
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.8888, p-value = 0.01055# 检查方差齐性
var_test <- var.test(reduction ~ group, data = bp_data)
var_test##
## F test to compare two variances
##
## data: reduction by group
## F = 1.2592, num df = 24, denom df = 24, p-value = 0.5768
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5548954 2.8574999
## sample estimates:
## ratio of variances
## 1.259211# 独立样本t检验
# 方差相等假设
t_equal <- t.test(reduction ~ group, data = bp_data, var.equal = TRUE)
t_equal##
## Two Sample t-test
##
## data: reduction by group
## t = 3.0893, df = 48, p-value = 0.003333
## alternative hypothesis: true difference in means between group Drug and group Placebo is not equal to 0
## 95 percent confidence interval:
## 1.815993 8.585922
## sample estimates:
## mean in group Drug mean in group Placebo
## 14.181663 8.980705# Welch校正（默认，更稳健）
t_welch <- t.test(reduction ~ group, data = bp_data, var.equal = FALSE)
t_welch##
## Welch Two Sample t-test
##
## data: reduction by group
## t = 3.0893, df = 47.376, p-value = 0.003352
## alternative hypothesis: true difference in means between group Drug and group Placebo is not equal to 0
## 95 percent confidence interval:
## 1.814842 8.587074
## sample estimates:
## mean in group Drug mean in group Placebo
## 14.181663 8.980705# 效应量（Cohen's d）
# 使用DescTools包
DescTools::CohenD(drug_group, placebo_group)## [1] 0.8737906
## attr(,"magnitude")
## [1] "large"# 手动计算
s_pooled <- sqrt(((25-1)*sd(drug_group)^2 + (25-1)*sd(placebo_group)^2) / (25 + 25 - 2))
d_manual <- (mean(drug_group) - mean(placebo_group)) / s_pooled
cat("手动计算的Cohen's d:", d_manual, "\n")## 手动计算的Cohen's d: 0.87379065.3 配对样本t检验
数学基础
检验配对数据（如治疗前后）的差异是否显著。
检验统计量： \[t = \frac{\bar{d}}{s_d/\sqrt{n}}\]
其中d是配对差值，\(\bar{d}\)是差值均值，\(s_d\)是差值标准差。
适用场景： - 治疗前后测量 - 同一受试者的两种测量方法 - 配对设计的实验
# 配对样本t检验示例
# 场景：患者治疗前后的血压变化
set.seed(33333)
n <- 20
# 治疗前血压
before <- rnorm(n, mean = 145, sd = 12)
# 治疗后血压（假设平均下降10 mmHg）
after <- before - rnorm(n, mean = 10, sd = 5)
# 创建数据框
paired_data <- data.frame(
patient = 1:n,
before = before,
after = after,
diff = before - after
)
# 查看数据
head(paired_data)## patient before after diff
## 1 1 158.6004 148.9063 9.694084
## 2 2 139.4294 128.3045 11.124910
## 3 3 127.3783 114.0198 13.358465
## 4 4 149.1485 130.5847 18.563822
## 5 5 141.2193 128.8387 12.380648
## 6 6 148.9844 134.7064 14.277963# 描述性统计
cat("治疗前均值:", mean(before), "\n")## 治疗前均值: 141.9287cat("治疗后均值:", mean(after), "\n")## 治疗后均值: 131.4311cat("平均差值:", mean(before - after), "\n")## 平均差值: 10.49755# 可视化配对数据
# 差值的分布
par(mfrow = c(1, 2))
hist(paired_data$diff, breaks = 10, main = "治疗前后差值分布",
xlab = "差值 (mmHg)", col = "lightgreen")
abline(v = 0, col = "red", lwd = 2, lty = 2)
# 配对连线图
plot(rep(1, n), before, xlim = c(0.5, 2.5), ylim = range(c(before, after)),
main = "治疗前后血压变化", xlab = "", ylab = "血压 (mmHg)",
xaxt = "n", pch = 16, col = "blue")
points(rep(2, n), after, pch = 16, col = "green")
segments(rep(1, n), before, rep(2, n), after, col = "gray")
axis(1, at = c(1, 2), labels = c("治疗前", "治疗后"))
legend("topright", legend = c("治疗前", "治疗后"),
col = c("blue", "green"), pch = 16)
par(mfrow = c(1, 1))
# 检查差值的正态性
shapiro.test(paired_data$diff)##
## Shapiro-Wilk normality test
##
## data: paired_data$diff
## W = 0.93545, p-value = 0.1965# 配对t检验
# 方法1：使用两个向量
t_paired <- t.test(before, after, paired = TRUE)
t_paired##
## Paired t-test
##
## data: before and after
## t = 9.2511, df = 19, p-value = 1.816e-08
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## 8.122531 12.872569
## sample estimates:
## mean difference
## 10.49755# 方法2：使用差值
t_diff <- t.test(paired_data$diff, mu = 0)
t_diff##
## One Sample t-test
##
## data: paired_data$diff
## t = 9.2511, df = 19, p-value = 1.816e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 8.122531 12.872569
## sample estimates:
## mean of x
## 10.49755# 效应量
# 对于配对数据，使用差值的标准差
cohens_d_paired <- mean(paired_data$diff) / sd(paired_data$diff)
cat("配对数据的Cohen's d:", cohens_d_paired, "\n")## 配对数据的Cohen's d: 2.0686175.4 单样本比例检验
数学基础
检验样本比例是否与已知总体比例有显著差异。
检验统计量（大样本）： \[z = \frac{\hat{p} - p_0}{\sqrt{p_0(1-p_0)/n}}\]
精确检验：使用二项分布。
# 单样本比例检验示例
# 场景：检验某疾病的治愈率是否为70%
# 数据：100名患者中80人治愈
n <- 100
successes <- 80
p_null <- 0.70
# 精确二项检验
binom_result <- binom.test(successes, n, p = p_null)
binom_result##
## Exact binomial test
##
## data: successes and n
## number of successes = 80, number of trials = 100, p-value = 0.02896
## alternative hypothesis: true probability of success is not equal to 0.7
## 95 percent confidence interval:
## 0.7081573 0.8733444
## sample estimates:
## probability of success
## 0.8# 正态近似检验
prop_result <- prop.test(successes, n, p = p_null, correct = FALSE)
prop_result##
## 1-sample proportions test without continuity correction
##
## data: successes out of n, null probability p_null
## X-squared = 4.7619, df = 1, p-value = 0.0291
## alternative hypothesis: true p is not equal to 0.7
## 95 percent confidence interval:
## 0.7111708 0.8666331
## sample estimates:
## p
## 0.8# 带连续性校正
prop_corrected <- prop.test(successes, n, p = p_null, correct = TRUE)
prop_corrected##
## 1-sample proportions test with continuity correction
##
## data: successes out of n, null probability p_null
## X-squared = 4.2976, df = 1, p-value = 0.03817
## alternative hypothesis: true p is not equal to 0.7
## 95 percent confidence interval:
## 0.7056770 0.8707518
## sample estimates:
## p
## 0.8# 比较结果
cat("\n三种方法比较:\n")##
## 三种方法比较:cat("精确二项检验 p值:", binom_result$p.value, "\n")## 精确二项检验 p值: 0.02896126cat("正态近似 p值:", prop_result$p.value, "\n")## 正态近似 p值: 0.02909633cat("连续性校正 p值:", prop_corrected$p.value, "\n")## 连续性校正 p值: 0.03816577# 单侧检验（治愈率是否高于70%）
prop_greater <- prop.test(successes, n, p = p_null,
alternative = "greater", correct = FALSE)
prop_greater##
## 1-sample proportions test without continuity correction
##
## data: successes out of n, null probability p_null
## X-squared = 4.7619, df = 1, p-value = 0.01455
## alternative hypothesis: true p is greater than 0.7
## 95 percent confidence interval:
## 0.7266962 1.0000000
## sample estimates:
## p
## 0.85.5 两样本比例检验
数学基础
检验两个独立样本的比例是否有显著差异。
检验统计量： \[z = \frac{\hat{p}_1 - \hat{p}_2}{\sqrt{\hat{p}(1-\hat{p})(\frac{1}{n_1} + \frac{1}{n_2})}}\]
其中 \(\hat{p} = \frac{x_1 + x_2}{n_1 + n_2}\) 是合并比例。
# 两样本比例检验示例
# 场景：比较两种治疗方法的治愈率
# 治疗组A：120人中95人治愈
# 治疗组B：110人中78人治愈
successes <- c(95, 78)
totals <- c(120, 110)
# 卡方检验（比例差异）
prop_comparison <- prop.test(successes, totals, correct = FALSE)
prop_comparison##
## 2-sample test for equality of proportions without continuity correction
##
## data: successes out of totals
## X-squared = 2.0994, df = 1, p-value = 0.1474
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.02915428 0.19430579
## sample estimates:
## prop 1 prop 2
## 0.7916667 0.7090909# 带Yates连续性校正
prop_yates <- prop.test(successes, totals, correct = TRUE)
prop_yates##
## 2-sample test for equality of proportions with continuity correction
##
## data: successes out of totals
## X-squared = 1.6797, df = 1, p-value = 0.195
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.0378664 0.2030179
## sample estimates:
## prop 1 prop 2
## 0.7916667 0.7090909# Fisher精确检验（小样本时推荐）
fisher_result <- fisher.test(matrix(c(95, 25, 78, 32), nrow = 2))
fisher_result##
## Fisher's Exact Test for Count Data
##
## data: matrix(c(95, 25, 78, 32), nrow = 2)
## p-value = 0.1699
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.8173419 2.9886143
## sample estimates:
## odds ratio
## 1.555915# 计算比例和置信区间
p1 <- 95/120
p2 <- 78/110
cat("治疗组A治愈率:", round(p1*100, 1), "%\n")## 治疗组A治愈率: 79.2 %cat("治疗组B治愈率:", round(p2*100, 1), "%\n")## 治疗组B治愈率: 70.9 %cat("比例差:", round((p1-p2)*100, 1), "%\n")## 比例差: 8.3 %# 效应量（Odds Ratio）
or <- (95/25) / (78/32)
cat("Odds Ratio:", or, "\n")## Odds Ratio: 1.5589745.6 单样本方差检验
数学基础
检验总体方差是否等于某个特定值。
检验统计量： \[\chi^2 = \frac{(n-1)s^2}{\sigma_0^2}\]
服从自由度为n-1的卡方分布。
# 单样本方差检验示例
# 场景：检验某生产过程的方差是否为100（标准差为10）
set.seed(44444)
sample_data <- rnorm(30, mean = 50, sd = 12)
# 假设检验
# H₀: σ² = 100 vs H₁: σ² ≠ 100
n <- length(sample_data)
sample_var <- var(sample_data)
null_var <- 100
# 卡方统计量
chi_sq <- (n - 1) * sample_var / null_var
chi_sq## [1] 30.80413# p值（双侧）
p_value <- 2 * min(pchisq(chi_sq, df = n-1), 1 - pchisq(chi_sq, df = n-1))
p_value## [1] 0.7494319# 使用DescTools包
DescTools::VarTest(sample_data, sigma.squared = 100)##
## One Sample Chi-Square test on variance
##
## data: sample_data
## X-squared = 30.804, df = 29, p-value = 0.6226
## alternative hypothesis: true variance is not equal to 100
## 95 percent confidence interval:
## 67.37224 191.96107
## sample estimates:
## variance of x
## 106.22115.7 两样本方差比较（F检验）
数学基础
检验两个独立样本的方差是否相等。
检验统计量： \[F = \frac{s_1^2}{s_2^2}\]
服从自由度为\((n_1-1, n_2-1)\)的F分布。
# 两样本方差比较示例
set.seed(55555)
group1 <- rnorm(25, mean = 100, sd = 15)
group2 <- rnorm(30, mean = 100, sd = 10)
# F检验
var_test <- var.test(group1, group2)
var_test##
## F test to compare two variances
##
## data: group1 and group2
## F = 1.771, num df = 24, denom df = 29, p-value = 0.1425
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.8222047 3.9271662
## sample estimates:
## ratio of variances
## 1.771034# Levene检验（更稳健，对正态性假设不那么敏感）
car::leveneTest(c(group1, group2),
factor(rep(c(1, 2), c(25, 30))))## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 2.6082 0.1123
## 53# Bartlett检验（对正态性敏感）
bartlett.test(list(group1, group2))##
## Bartlett test of homogeneity of variances
##
## data: list(group1, group2)
## Bartlett's K-squared = 2.1142, df = 1, p-value = 0.1459# 结果解读
cat("F检验结果:\n")## F检验结果:cat("F统计量:", var_test$statistic, "\n")## F统计量: 1.771034cat("p值:", var_test$p.value, "\n")## p值: 0.1425442cat("结论:", ifelse(var_test$p.value < 0.05,
"方差不等",
"可以认为方差相等"), "\n")## 结论: 可以认为方差相等5.8 两总体均值比较的Z检验（大样本）
数学基础
当样本量足够大时（通常n ≥ 30），可以使用Z检验代替t检验。
检验统计量： \[z = \frac{\bar{X}_1 - \bar{X}_2}{\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}}\]
# 大样本Z检验示例
set.seed(66666)
# 大样本数据
group1 <- rnorm(100, mean = 100, sd = 15)
group2 <- rnorm(100, mean = 105, sd = 15)
# 手动计算Z检验
n1 <- length(group1)
n2 <- length(group2)
mean1 <- mean(group1)
mean2 <- mean(group2)
sd1 <- sd(group1)
sd2 <- sd(group2)
# Z统计量
z_stat <- (mean1 - mean2) / sqrt(sd1^2/n1 + sd2^2/n2)
z_stat## [1] -0.8940959# p值（双侧）
p_value_z <- 2 * pnorm(-abs(z_stat))
p_value_z## [1] 0.3712706# 与t检验比较
t_result <- t.test(group1, group2)
t_result##
## Welch Two Sample t-test
##
## data: group1 and group2
## t = -0.8941, df = 194.69, p-value = 0.3724
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -6.348764 2.388003
## sample estimates:
## mean of x mean of y
## 102.5070 104.4874cat("\nZ检验与t检验比较:\n")##
## Z检验与t检验比较:cat("Z统计量:", z_stat, "\n")## Z统计量: -0.8940959cat("t统计量:", t_result$statistic, "\n")## t统计量: -0.8940959cat("Z检验p值:", p_value_z, "\n")## Z检验p值: 0.3712706cat("t检验p值:", t_result$p.value, "\n")## t检验p值: 0.3723745# 大样本时两者结果非常接近5.9 符号检验与Wilcoxon符号秩检验
虽然这些是非参数检验，但常作为参数检验的替代方法，在此简要介绍。
数学基础
符号检验：仅考虑差值的符号，不考虑大小。
Wilcoxon符号秩检验：考虑差值的符号和相对大小。
# 符号检验与Wilcoxon符号秩检验示例
set.seed(77777)
before <- rnorm(20, mean = 100, sd = 15)
after <- before + rnorm(20, mean = 5, sd = 8)
# 差值
differences <- after - before
# 符号检验
# 计算正差值的数量
n_positive <- sum(differences > 0)
n_negative <- sum(differences < 0)
n_total <- n_positive + n_negative
# 二项检验
binom.test(n_positive, n_total, p = 0.5)##
## Exact binomial test
##
## data: n_positive and n_total
## number of successes = 15, number of trials = 20, p-value = 0.04139
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.5089541 0.9134285
## sample estimates:
## probability of success
## 0.75# Wilcoxon符号秩检验
wilcox_signed <- wilcox.test(before, after, paired = TRUE)
wilcox_signed##
## Wilcoxon signed rank exact test
##
## data: before and after
## V = 41, p-value = 0.01531
## alternative hypothesis: true location shift is not equal to 0# 与配对t检验比较
t_paired <- t.test(before, after, paired = TRUE)
cat("\n三种方法比较:\n")##
## 三种方法比较:cat("配对t检验 p值:", t_paired$p.value, "\n")## 配对t检验 p值: 0.02080963cat("Wilcoxon符号秩检验 p值:", wilcox_signed$p.value, "\n")## Wilcoxon符号秩检验 p值: 0.01531219第五章小结
| 检验类型 | 适用场景 | R函数 | 假设条件 |
|---|---|---|---|
| 单样本t检验 | 样本均值与已知值比较 | t.test(x, mu=) | 正态分布 |
| 独立样本t检验 | 两组独立样本均值比较 | t.test(formula) | 正态、独立、方差齐 |
| 配对t检验 | 配对数据均值比较 | t.test(x, y, paired=TRUE) | 差值正态 |
| 单样本比例检验 | 比例与已知值比较 | prop.test(), binom.test() | 大样本或精确 |
| 两样本比例检验 | 两组比例比较 | prop.test() | 大样本 |
| 单样本方差检验 | 方差与已知值比较 | DescTools::VarTest() | 正态分布 |
| 两样本方差比较 | 两组方差比较 | var.test() | 正态分布 |
| Z检验 | 大样本均值比较 | 手动计算 | 大样本 |
第六章：非参数检验
非参数检验不依赖于数据的具体分布形式，适用于数据不满足正态分布假设、样本量较小、或数据为等级/有序类型的情况。本章将介绍常用的非参数检验方法。
非参数检验的优势： 1. 不要求总体分布形式 2. 对异常值不敏感 3. 适用于等级数据和有序数据 4. 小样本情况下更稳健
非参数检验的局限： 1. 当数据满足参数检验假设时，非参数检验功效较低 2. 信息利用不如参数检验充分
6.1 单样本分布检验
6.1.1 Kolmogorov-Smirnov检验
数学基础：
检验样本是否来自某个理论分布。
检验统计量： \[D = \max|F_n(x) - F_0(x)|\]
其中\(F_n(x)\)是经验分布函数，\(F_0(x)\)是理论分布函数。
# Kolmogorov-Smirnov检验示例
set.seed(88888)
# 生成数据
normal_data <- rnorm(100, mean = 50, sd = 10)
exp_data <- rexp(100, rate = 0.1)
# 检验是否来自正态分布
# 需要先标准化数据
ks_normal <- ks.test(scale(normal_data), "pnorm", mean = 0, sd = 1)
ks_normal##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: scale(normal_data)
## D = 0.04866, p-value = 0.9719
## alternative hypothesis: two-sided# 检验是否来自指数分布
ks_exp <- ks.test(exp_data, "pexp", rate = 0.1)
ks_exp##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: exp_data
## D = 0.10507, p-value = 0.2196
## alternative hypothesis: two-sided# 可视化比较
par(mfrow = c(1, 2))
# 正态数据
hist(normal_data, breaks = 20, probability = TRUE,
main = "正态数据分布", xlab = "值", col = "lightblue")
curve(dnorm(x, mean = mean(normal_data), sd = sd(normal_data)),
add = TRUE, col = "red", lwd = 2)
# 指数数据
hist(exp_data, breaks = 20, probability = TRUE,
main = "指数数据分布", xlab = "值", col = "lightgreen")
curve(dexp(x, rate = 0.1), add = TRUE, col = "red", lwd = 2)
par(mfrow = c(1, 1))6.1.2 卡方拟合优度检验
数学基础：
检验观测频数与期望频数是否一致。
检验统计量： \[\chi^2 = \sum\frac{(O_i - E_i)^2}{E_i}\]
# 卡方拟合优度检验示例
# 场景：检验骰子是否公平
# 观测频数（掷骰子120次）
observed <- c(18, 22, 15, 20, 25, 20)
names(observed) <- 1:6
# 期望频数（公平骰子，每个面期望20次）
expected <- rep(20, 6)
# 卡方拟合优度检验
chisq_goodness <- chisq.test(observed, p = rep(1/6, 6))
chisq_goodness##
## Chi-squared test for given probabilities
##
## data: observed
## X-squared = 2.9, df = 5, p-value = 0.7154# 查看期望频数
chisq_goodness$expected## 1 2 3 4 5 6
## 20 20 20 20 20 20# 残差分析
chisq_goodness$residuals## 1 2 3 4 5 6
## -0.4472136 0.4472136 -1.1180340 0.0000000 1.1180340 0.0000000# 可视化
barplot(rbind(observed, expected), beside = TRUE,
names.arg = 1:6,
main = "骰子各面出现频数",
xlab = "骰子面", ylab = "频数",
col = c("steelblue", "coral"))
legend("topright", legend = c("观测值", "期望值"),
fill = c("steelblue", "coral"))
6.2 两独立样本非参数检验
Mann-Whitney U检验（Wilcoxon秩和检验）
数学基础：
检验两个独立样本是否来自同一分布。
检验统计量： \[U = n_1 n_2 + \frac{n_1(n_1+1)}{2} - R_1\]
其中\(R_1\)是第一组样本的秩和。
# Mann-Whitney U检验示例
# 场景：比较两种药物的疗效评分（有序数据）
set.seed(99999)
drug_A <- sample(1:5, 30, replace = TRUE, prob = c(0.1, 0.15, 0.2, 0.3, 0.25))
drug_B <- sample(1:5, 30, replace = TRUE, prob = c(0.2, 0.25, 0.25, 0.2, 0.1))
# 创建数据框
score_data <- data.frame(
score = c(drug_A, drug_B),
drug = factor(rep(c("A", "B"), each = 30))
)
# 描述性统计
table(score_data$score, score_data$drug)##
## A B
## 1 3 7
## 2 7 6
## 3 2 6
## 4 8 7
## 5 10 4# Wilcoxon秩和检验（Mann-Whitney U）
wilcox_test <- wilcox.test(score ~ drug, data = score_data)## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with tieswilcox_test##
## Wilcoxon rank sum test with continuity correction
##
## data: score by drug
## W = 572.5, p-value = 0.0649
## alternative hypothesis: true location shift is not equal to 0# 使用coin包进行精确检验
coin_test <- coin::wilcox_test(score ~ drug, data = score_data,
distribution = "exact")
coin_test##
## Exact Wilcoxon-Mann-Whitney Test
##
## data: score by drug (A, B)
## Z = 1.8535, p-value = 0.06261
## alternative hypothesis: true mu is not equal to 0# 可视化
ggplot2::ggplot(score_data, ggplot2::aes(x = drug, y = score, fill = drug)) +
ggplot2::geom_boxplot() +
ggplot2::geom_jitter(width = 0.2, alpha = 0.5) +
ggplot2::labs(title = "两种药物疗效评分比较",
x = "药物", y = "疗效评分") +
ggplot2::theme_minimal()
6.3 两配对样本非参数检验
Wilcoxon符号秩检验
数学基础：
检验配对差值的中位数是否为0。
# Wilcoxon符号秩检验示例
set.seed(111111)
before <- rnorm(20, mean = 50, sd = 10)
after <- before + rnorm(20, mean = 3, sd = 5)
# Wilcoxon符号秩检验
wilcox_paired <- wilcox.test(before, after, paired = TRUE)
wilcox_paired##
## Wilcoxon signed rank exact test
##
## data: before and after
## V = 48, p-value = 0.03277
## alternative hypothesis: true location shift is not equal to 0# 与配对t检验比较
t_paired <- t.test(before, after, paired = TRUE)
cat("Wilcoxon符号秩检验 p值:", wilcox_paired$p.value, "\n")## Wilcoxon符号秩检验 p值: 0.03276825cat("配对t检验 p值:", t_paired$p.value, "\n")## 配对t检验 p值: 0.06696449# 符号检验（更简单，但功效较低）
diff <- after - before
n_pos <- sum(diff > 0)
n_neg <- sum(diff < 0)
binom.test(n_pos, n_pos + n_neg, p = 0.5)##
## Exact binomial test
##
## data: n_pos and n_pos + n_neg
## number of successes = 14, number of trials = 20, p-value = 0.1153
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.4572108 0.8810684
## sample estimates:
## probability of success
## 0.76.4 多组独立样本非参数检验
Kruskal-Wallis检验
数学基础：
检验多组独立样本是否来自同一分布（非参数方差分析）。
检验统计量： \[H = \frac{12}{N(N+1)}\sum\frac{R_i^2}{n_i} - 3(N+1)\]
# Kruskal-Wallis检验示例
# 场景：比较三种治疗方法的疗效评分
set.seed(222222)
treatment_A <- rnorm(20, mean = 60, sd = 15)
treatment_B <- rnorm(20, mean = 70, sd = 15)
treatment_C <- rnorm(20, mean = 55, sd = 15)
# 创建数据框
kw_data <- data.frame(
score = c(treatment_A, treatment_B, treatment_C),
group = factor(rep(c("A", "B", "C"), each = 20))
)
# Kruskal-Wallis检验
kw_test <- kruskal.test(score ~ group, data = kw_data)
kw_test##
## Kruskal-Wallis rank sum test
##
## data: score by group
## Kruskal-Wallis chi-squared = 15.181, df = 2, p-value = 0.0005051# 事后多重比较（Dunn检验）
# 使用pairwise.wilcox.test
pairwise.wilcox.test(kw_data$score, kw_data$group,
p.adjust.method = "bonferroni")##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: kw_data$score and kw_data$group
##
## A B
## B 0.24293 -
## C 0.03958 0.00036
##
## P value adjustment method: bonferroni# 可视化
ggplot2::ggplot(kw_data, ggplot2::aes(x = group, y = score, fill = group)) +
ggplot2::geom_boxplot() +
ggplot2::geom_jitter(width = 0.2, alpha = 0.5) +
ggplot2::labs(title = "三种治疗方法疗效比较",
x = "治疗组", y = "疗效评分") +
ggplot2::theme_minimal()
6.5 多组相关样本非参数检验
Friedman检验
数学基础：
检验多组相关样本（如重复测量）是否来自同一分布。
# Friedman检验示例
# 场景：同一患者接受三种不同治疗的评分
set.seed(333333)
n_subjects <- 15
# 生成数据（每个受试者三种治疗）
treatment1 <- rnorm(n_subjects, mean = 60, sd = 10)
treatment2 <- treatment1 + rnorm(n_subjects, mean = 10, sd = 5)
treatment3 <- treatment1 + rnorm(n_subjects, mean = 5, sd = 8)
# 创建数据矩阵
friedman_matrix <- cbind(treatment1, treatment2, treatment3)
# Friedman检验
friedman_test <- friedman.test(friedman_matrix)
friedman_test##
## Friedman rank sum test
##
## data: friedman_matrix
## Friedman chi-squared = 13.733, df = 2, p-value = 0.001042# 使用rstatix包（更现代的方法）
friedman_data <- data.frame(
subject = factor(rep(1:n_subjects, 3)),
treatment = factor(rep(c("T1", "T2", "T3"), each = n_subjects)),
score = c(treatment1, treatment2, treatment3)
)
rstatix::friedman_test(friedman_data, score ~ treatment | subject)## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <dbl> <dbl> <chr>
## 1 score 15 13.7 2 0.00104 Friedman test# 事后比较
pairwise.wilcox.test(friedman_data$score, friedman_data$treatment,
paired = TRUE, p.adjust.method = "bonferroni")##
## Pairwise comparisons using Wilcoxon signed rank exact test
##
## data: friedman_data$score and friedman_data$treatment
##
## T1 T2
## T2 0.00037 -
## T3 0.04523 0.36163
##
## P value adjustment method: bonferroni6.6 分类变量关联性检验
6.6.1 卡方独立性检验
数学基础：
检验两个分类变量是否独立。
# 卡方独立性检验示例
# 场景：检验吸烟与肺癌的关系
# 创建列联表
smoking_cancer <- matrix(c(60, 40, 30, 70), nrow = 2)
rownames(smoking_cancer) <- c("吸烟", "不吸烟")
colnames(smoking_cancer) <- c("肺癌", "无肺癌")
smoking_cancer## 肺癌 无肺癌
## 吸烟 60 30
## 不吸烟 40 70# 卡方独立性检验
chisq_ind <- chisq.test(smoking_cancer)
chisq_ind##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: smoking_cancer
## X-squared = 16.99, df = 1, p-value = 3.758e-05# 查看期望频数
chisq_ind$expected## 肺癌 无肺癌
## 吸烟 45 45
## 不吸烟 55 55# 查看残差
chisq_ind$residuals## 肺癌 无肺癌
## 吸烟 2.236068 -2.236068
## 不吸烟 -2.022600 2.022600# 效应量（Cramér's V）
vcd::assocstats(smoking_cancer)## X^2 df P(> X^2)
## Likelihood Ratio 18.480 1 1.7167e-05
## Pearson 18.182 1 2.0079e-05
##
## Phi-Coefficient : 0.302
## Contingency Coeff.: 0.289
## Cramer's V : 0.3026.6.2 Fisher精确检验
数学基础：
适用于小样本或期望频数较小的情况。
# Fisher精确检验示例
# 场景：小样本的药物治疗效果
# 创建2x2列联表
treatment_outcome <- matrix(c(8, 2, 3, 7), nrow = 2)
rownames(treatment_outcome) <- c("治疗组", "对照组")
colnames(treatment_outcome) <- c("有效", "无效")
treatment_outcome## 有效 无效
## 治疗组 8 3
## 对照组 2 7# Fisher精确检验
fisher_test <- fisher.test(treatment_outcome)
fisher_test##
## Fisher's Exact Test for Count Data
##
## data: treatment_outcome
## p-value = 0.06978
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.8821175 127.0558418
## sample estimates:
## odds ratio
## 8.153063# 单侧检验
fisher.test(treatment_outcome, alternative = "greater")##
## Fisher's Exact Test for Count Data
##
## data: treatment_outcome
## p-value = 0.03489
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
## 1.155327 Inf
## sample estimates:
## odds ratio
## 8.153063# 使用exact2x2包进行更精确的计算
exact2x2::exact2x2(treatment_outcome)##
## Two-sided Fisher's Exact Test (usual method using minimum likelihood)
##
## data: treatment_outcome
## p-value = 0.06978
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.0000 84.0382
## sample estimates:
## odds ratio
## 8.1530636.7 一致性检验
McNemar检验
数学基础：
检验配对分类变量的变化是否显著。
# McNemar检验示例
# 场景：治疗前后症状的变化
# 创建配对列联表
# 行：治疗前（有症状/无症状）
# 列：治疗后（有症状/无症状）
mcnemar_table <- matrix(c(20, 30, 10, 40), nrow = 2)
rownames(mcnemar_table) <- c("治疗前有症状", "治疗前无症状")
colnames(mcnemar_table) <- c("治疗后有症状", "治疗后无症状")
mcnemar_table## 治疗后有症状 治疗后无症状
## 治疗前有症状 20 10
## 治疗前无症状 30 40# McNemar检验
mcnemar_test <- mcnemar.test(mcnemar_table)
mcnemar_test##
## McNemar's Chi-squared test with continuity correction
##
## data: mcnemar_table
## McNemar's chi-squared = 9.025, df = 1, p-value = 0.002663# 精确McNemar检验（小样本）
# 使用exact2x2包
exact2x2::mcnemar.exact(mcnemar_table)##
## Exact McNemar test (with central confidence intervals)
##
## data: mcnemar_table
## b = 10, c = 30, p-value = 0.002221
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.1453636 0.7005703
## sample estimates:
## odds ratio
## 0.33333336.8 游程检验（随机性检验）
数学基础：
检验数据是否随机分布。
# 游程检验示例
set.seed(444444)
# 生成随机序列
random_seq <- sample(c(0, 1), 50, replace = TRUE)
# 生成有趋势的序列
trend_seq <- c(rep(0, 25), rep(1, 25))
# 使用DescTools包的RunsTest
DescTools::RunsTest(random_seq)##
## Runs Test for Randomness
##
## data: random_seq
## z = -0.9904, runs = 22, m = 26, n = 24, p-value = 0.322
## alternative hypothesis: true number of runs is not equal the expected numberDescTools::RunsTest(trend_seq)##
## Runs Test for Randomness
##
## data: trend_seq
## z = -6.7157, runs = 2, m = 25, n = 25, p-value = 1.872e-11
## alternative hypothesis: true number of runs is not equal the expected number# 可视化
par(mfrow = c(1, 2))
plot(random_seq, type = "s", main = "随机序列", ylab = "值")
plot(trend_seq, type = "s", main = "有趋势序列", ylab = "值")
par(mfrow = c(1, 1))6.9 相关性非参数检验
Spearman秩相关与Kendall’s tau
# 非参数相关检验示例
set.seed(555555)
x <- rnorm(30, mean = 50, sd = 10)
y <- x + rnorm(30, mean = 0, sd = 15)
# Spearman秩相关
spearman_cor <- cor.test(x, y, method = "spearman")
spearman_cor##
## Spearman's rank correlation rho
##
## data: x and y
## S = 1968, p-value = 0.001468
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.5621802# Kendall's tau
kendall_cor <- cor.test(x, y, method = "kendall")
kendall_cor##
## Kendall's rank correlation tau
##
## data: x and y
## T = 302, p-value = 0.002238
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.3885057# 与Pearson相关比较
pearson_cor <- cor.test(x, y, method = "pearson")
cat("Pearson相关系数:", pearson_cor$estimate, "p值:", pearson_cor$p.value, "\n")## Pearson相关系数: 0.6425019 p值: 0.0001290706cat("Spearman相关系数:", spearman_cor$estimate, "p值:", spearman_cor$p.value, "\n")## Spearman相关系数: 0.5621802 p值: 0.001467922cat("Kendall's tau:", kendall_cor$estimate, "p值:", kendall_cor$p.value, "\n")## Kendall's tau: 0.3885057 p值: 0.0022378# 可视化
plot(x, y, main = "散点图", xlab = "X", ylab = "Y", pch = 16, col = "steelblue")
abline(lm(y ~ x), col = "red", lwd = 2)
第六章小结
| 检验类型 | 适用场景 | R函数 | 参数检验对应 |
|---|---|---|---|
| K-S检验 | 分布拟合检验 | ks.test() | - |
| 卡方拟合优度 | 频数分布检验 | chisq.test() | - |
| Mann-Whitney U | 两独立样本比较 | wilcox.test() | 独立样本t检验 |
| Wilcoxon符号秩 | 两配对样本比较 | wilcox.test(paired=TRUE) | 配对t检验 |
| Kruskal-Wallis | 多组独立样本 | kruskal.test() | 单因素ANOVA |
| Friedman | 多组相关样本 | friedman.test() | 重复测量ANOVA |
| 卡方独立性 | 分类变量关联 | chisq.test() | - |
| Fisher精确检验 | 小样本列联表 | fisher.test() | - |
| McNemar检验 | 配对分类变量 | mcnemar.test() | - |
| Spearman相关 | 等级相关 | cor.test(method=“spearman”) | Pearson相关 |
第七章：方差分析（ANOVA）
方差分析是用于比较多组均值差异的统计方法，是t检验的扩展。本章将介绍各种方差分析方法及其在R中的实现。
方差分析的基本思想：
将总变异分解为组间变异和组内变异，通过比较两者的大小来判断组间均值是否有显著差异。
7.1 方差分析的基本思想
数学基础
总变异分解： \[SS_{total} = SS_{between} + SS_{within}\]
F统计量： \[F = \frac{MS_{between}}{MS_{within}} = \frac{SS_{between}/(k-1)}{SS_{within}/(N-k)}\]
其中： - \(SS_{between}\)：组间平方和 - \(SS_{within}\)：组内平方和 - k：组数 - N：总样本量
假设条件： 1. 各组数据来自正态分布 2. 各组方差相等（方差齐性） 3. 观测值独立
# 方差分析基本概念演示
set.seed(666666)
# 生成三组数据
group1 <- rnorm(20, mean = 50, sd = 10)
group2 <- rnorm(20, mean = 55, sd = 10)
group3 <- rnorm(20, mean = 60, sd = 10)
# 创建数据框
anova_data <- data.frame(
value = c(group1, group2, group3),
group = factor(rep(c("A", "B", "C"), each = 20))
)
# 计算各部分变异
grand_mean <- mean(anova_data$value)
group_means <- tapply(anova_data$value, anova_data$group, mean)
# 总平方和
SS_total <- sum((anova_data$value - grand_mean)^2)
# 组间平方和
SS_between <- sum(20 * (group_means - grand_mean)^2)
# 组内平方和
SS_within <- sum(sapply(1:3, function(i) {
sum((anova_data$value[anova_data$group == unique(anova_data$group)[i]] -
group_means[i])^2)
}))
# 均方
MS_between <- SS_between / 2 # df = k-1 = 2
MS_within <- SS_within / 57 # df = N-k = 57
# F统计量
F_stat <- MS_between / MS_within
cat("总平方和 SS_total:", SS_total, "\n")## 总平方和 SS_total: 6434.028cat("组间平方和 SS_between:", SS_between, "\n")## 组间平方和 SS_between: 1673.201cat("组内平方和 SS_within:", SS_within, "\n")## 组内平方和 SS_within: 4760.827cat("验证: SS_between + SS_within =", SS_between + SS_within, "\n")## 验证: SS_between + SS_within = 6434.028cat("F统计量:", F_stat, "\n")## F统计量: 10.01638# p值
p_value <- 1 - pf(F_stat, 2, 57)
cat("p值:", p_value, "\n")## p值: 0.00018714737.2 单因素方差分析
数学基础
检验一个因素的多个水平之间均值是否有显著差异。
假设： - H₀: \(\mu_1 = \mu_2 = ... = \mu_k\) - H₁: 至少有两个均值不相等
# 单因素方差分析示例
# 场景：比较三种药物的降压效果
set.seed(777777)
drug_A <- rnorm(25, mean = 10, sd = 3)
drug_B <- rnorm(25, mean = 12, sd = 3)
drug_C <- rnorm(25, mean = 8, sd = 3)
drug_data <- data.frame(
reduction = c(drug_A, drug_B, drug_C),
drug = factor(rep(c("A", "B", "C"), each = 25))
)
# 描述性统计
drug_data %>%
dplyr::group_by(drug) %>%
dplyr::summarise(
n = dplyr::n(),
mean = mean(reduction),
sd = sd(reduction)
)## # A tibble: 3 × 4
## drug n mean sd
## <fct> <int> <dbl> <dbl>
## 1 A 25 9.27 2.84
## 2 B 25 12.1 3.88
## 3 C 25 9.02 2.65# 可视化
ggplot2::ggplot(drug_data, ggplot2::aes(x = drug, y = reduction, fill = drug)) +
ggplot2::geom_boxplot() +
ggplot2::geom_jitter(width = 0.2, alpha = 0.5) +
ggplot2::labs(title = "三种药物降压效果比较",
x = "药物", y = "血压降低值 (mmHg)") +
ggplot2::theme_minimal()
# 单因素ANOVA
anova_result <- aov(reduction ~ drug, data = drug_data)
summary(anova_result)## Df Sum Sq Mean Sq F value Pr(>F)
## drug 2 150.6 75.29 7.506 0.00109 **
## Residuals 72 722.1 10.03
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 使用rstatix包（更现代的方法）
rstatix::anova_test(drug_data, reduction ~ drug)## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 drug 2 72 7.506 0.001 * 0.1737.3 事后检验
当ANOVA结果显著时，需要进行事后检验确定哪些组之间存在差异。
常用方法
| 方法 | 特点 | 适用场景 |
|---|---|---|
| Tukey HSD | 最常用，控制整体错误率 | 各组样本量相等 |
| Bonferroni | 保守，简单 | 任意样本量 |
| Dunnett | 与对照组比较 | 有对照组的设计 |
| Games-Howell | 不要求方差齐性 | 方差不齐时 |
# 事后检验示例
# Tukey HSD检验
tukey_result <- TukeyHSD(anova_result)
tukey_result## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = reduction ~ drug, data = drug_data)
##
## $drug
## diff lwr upr p adj
## B-A 2.8760696 0.7324225 5.0197167 0.0055523
## C-A -0.2443531 -2.3880001 1.8992940 0.9598319
## C-B -3.1204227 -5.2640697 -0.9767756 0.0024103# 可视化Tukey结果
plot(tukey_result)
# 使用multcomp包
library(multcomp)
tukey_mc <- glht(anova_result, linfct = mcp(drug = "Tukey"))
summary(tukey_mc)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = reduction ~ drug, data = drug_data)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B - A == 0 2.8761 0.8958 3.211 0.00551 **
## C - A == 0 -0.2444 0.8958 -0.273 0.95983
## C - B == 0 -3.1204 0.8958 -3.484 0.00243 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)# Dunnett检验（与对照比较，假设A为对照组）
dunnett_mc <- glht(anova_result, linfct = mcp(drug = "Dunnett"))
summary(dunnett_mc)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Dunnett Contrasts
##
##
## Fit: aov(formula = reduction ~ drug, data = drug_data)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B - A == 0 2.8761 0.8958 3.211 0.00383 **
## C - A == 0 -0.2444 0.8958 -0.273 0.94710
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)# Bonferroni校正的成对t检验
pairwise.t.test(drug_data$reduction, drug_data$drug,
p.adjust.method = "bonferroni")##
## Pairwise comparisons using t tests with pooled SD
##
## data: drug_data$reduction and drug_data$drug
##
## A B
## B 0.0059 -
## C 1.0000 0.0025
##
## P value adjustment method: bonferroni# Games-Howell检验（方差不齐时）
rstatix::games_howell_test(drug_data, reduction ~ drug)## # A tibble: 3 × 8
## .y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 reduction A B 2.88 0.546 5.21 0.012 *
## 2 reduction A C -0.244 -2.12 1.63 0.947 ns
## 3 reduction B C -3.12 -5.40 -0.840 0.005 **7.4 多因素方差分析
数学基础
检验多个因素及其交互作用对结果的影响。
# 多因素方差分析示例
# 场景：研究药物类型和剂量对血压的影响
set.seed(888888)
# 2x3设计：2种药物 x 3种剂量
n_per_group <- 10
factorial_data <- data.frame(
drug = factor(rep(rep(c("A", "B"), each = 3), each = n_per_group)),
dose = factor(rep(rep(c("Low", "Medium", "High"), 2), each = n_per_group)),
reduction = c(
rnorm(n_per_group, 5, 2), # A, Low
rnorm(n_per_group, 8, 2), # A, Medium
rnorm(n_per_group, 10, 2), # A, High
rnorm(n_per_group, 6, 2), # B, Low
rnorm(n_per_group, 12, 2), # B, Medium
rnorm(n_per_group, 15, 2) # B, High
)
)
# 描述性统计
factorial_data %>%
dplyr::group_by(drug, dose) %>%
dplyr::summarise(
n = dplyr::n(),
mean = mean(reduction),
sd = sd(reduction)
)## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by drug and dose.
## ℹ Output is grouped by drug.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(drug, dose))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.## # A tibble: 6 × 5
## # Groups: drug [2]
## drug dose n mean sd
## <fct> <fct> <int> <dbl> <dbl>
## 1 A High 10 9.88 0.858
## 2 A Low 10 5.42 1.46
## 3 A Medium 10 6.53 1.30
## 4 B High 10 14.5 1.71
## 5 B Low 10 5.81 1.55
## 6 B Medium 10 12.2 1.67# 可视化交互效应
ggplot2::ggplot(factorial_data, ggplot2::aes(x = dose, y = reduction,
color = drug, group = drug)) +
ggplot2::stat_summary(fun = mean, geom = "point", size = 3) +
ggplot2::stat_summary(fun = mean, geom = "line") +
ggplot2::stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
ggplot2::labs(title = "药物与剂量的交互效应",
x = "剂量", y = "血压降低值") +
ggplot2::theme_minimal()
# 两因素ANOVA
factorial_anova <- aov(reduction ~ drug * dose, data = factorial_data)
summary(factorial_anova)## Df Sum Sq Mean Sq F value Pr(>F)
## drug 1 191.0 191.03 90.5 3.88e-13 ***
## dose 2 438.1 219.03 103.8 < 2e-16 ***
## drug:dose 2 78.1 39.05 18.5 7.59e-07 ***
## Residuals 54 114.0 2.11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 使用rstatix
rstatix::anova_test(factorial_data, reduction ~ drug * dose)## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 drug 1 54 90.499 3.88e-13 * 0.626
## 2 dose 2 54 103.764 3.17e-19 * 0.794
## 3 drug:dose 2 54 18.502 7.59e-07 * 0.4077.5 重复测量方差分析
数学基础
同一受试者在不同时间点或条件下的多次测量。
# 重复测量方差分析示例
# 场景：患者治疗前、中、后的血压变化
set.seed(999999)
n_subjects <- 20
# 生成重复测量数据
baseline <- rnorm(n_subjects, mean = 140, sd = 10)
midterm <- baseline - rnorm(n_subjects, mean = 5, sd = 5)
final <- baseline - rnorm(n_subjects, mean = 10, sd = 6)
rm_data <- data.frame(
subject = factor(rep(1:n_subjects, 3)),
time = factor(rep(c("Baseline", "Midterm", "Final"), each = n_subjects),
levels = c("Baseline", "Midterm", "Final")),
bp = c(baseline, midterm, final)
)
# 描述性统计
rm_data %>%
dplyr::group_by(time) %>%
dplyr::summarise(
n = dplyr::n(),
mean = mean(bp),
sd = sd(bp)
)## # A tibble: 3 × 4
## time n mean sd
## <fct> <int> <dbl> <dbl>
## 1 Baseline 20 135. 10.2
## 2 Midterm 20 131. 11.7
## 3 Final 20 125. 13.2# 可视化
ggplot2::ggplot(rm_data, ggplot2::aes(x = time, y = bp, group = subject)) +
ggplot2::geom_line(alpha = 0.3) +
ggplot2::geom_point(alpha = 0.5) +
ggplot2::stat_summary(fun = mean, geom = "point", size = 3, color = "red") +
ggplot2::stat_summary(fun = mean, geom = "line", color = "red", size = 1) +
ggplot2::labs(title = "患者血压随时间变化",
x = "时间点", y = "血压 (mmHg)") +
ggplot2::theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# 使用rstatix进行重复测量ANOVA
rstatix::anova_test(rm_data, bp ~ time, within = subject)## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 time 2 57 3.564 0.035 * 0.111# 使用afex包
afex::aov_ez(data = rm_data, dv = "bp", within = "time", id = "subject")## Anova Table (Type 3 tests)
##
## Response: bp
## Effect df MSE F ges p.value
## 1 time 1.45, 27.59 35.61 19.15 *** .111 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG# 事后检验
rstatix::pairwise_t_test(rm_data, bp ~ time, paired = TRUE,
p.adjust.method = "bonferroni")## # A tibble: 3 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 bp Baseline Midte… 20 20 4.10 19 6.11e-4 2 e-3 **
## 2 bp Baseline Final 20 20 6.53 19 2.95e-6 8.85e-6 ****
## 3 bp Midterm Final 20 20 2.56 19 1.9 e-2 5.8 e-2 ns7.6 方差分析的假设检验
正态性检验与方差齐性检验
# 假设检验示例
# 正态性检验（各组）
by(drug_data$reduction, drug_data$drug, shapiro.test)## drug_data$drug: A
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.97157, p-value = 0.6852
##
## ------------------------------------------------------------
## drug_data$drug: B
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.96863, p-value = 0.6107
##
## ------------------------------------------------------------
## drug_data$drug: C
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.95611, p-value = 0.3425# 方差齐性检验
# Bartlett检验（对正态性敏感）
bartlett.test(reduction ~ drug, data = drug_data)##
## Bartlett test of homogeneity of variances
##
## data: reduction by drug
## Bartlett's K-squared = 4.1012, df = 2, p-value = 0.1287# Levene检验（更稳健）
car::leveneTest(reduction ~ drug, data = drug_data)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.9096 0.4073
## 72# Fligner-Killeen检验（非参数）
fligner.test(reduction ~ drug, data = drug_data)##
## Fligner-Killeen test of homogeneity of variances
##
## data: reduction by drug
## Fligner-Killeen:med chi-squared = 1.7173, df = 2, p-value = 0.42377.7 稳健方差分析
当假设条件不满足时的替代方法。
# 稳健方差分析示例
# Welch's ANOVA（方差不齐时）
oneway.test(reduction ~ drug, data = drug_data, var.equal = FALSE)##
## One-way analysis of means (not assuming equal variances)
##
## data: reduction and drug
## F = 5.981, num df = 2.000, denom df = 46.999, p-value = 0.004851# 使用rstatix
rstatix::welch_anova_test(drug_data, reduction ~ drug)## # A tibble: 1 × 7
## .y. n statistic DFn DFd p method
## * <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 reduction 75 5.98 2 47.0 0.005 Welch ANOVA# Games-Howell事后检验
rstatix::games_howell_test(drug_data, reduction ~ drug)## # A tibble: 3 × 8
## .y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 reduction A B 2.88 0.546 5.21 0.012 *
## 2 reduction A C -0.244 -2.12 1.63 0.947 ns
## 3 reduction B C -3.12 -5.40 -0.840 0.005 **# Kruskal-Wallis检验（非参数替代）
kruskal.test(reduction ~ drug, data = drug_data)##
## Kruskal-Wallis rank sum test
##
## data: reduction by drug
## Kruskal-Wallis chi-squared = 13.135, df = 2, p-value = 0.0014057.8 效应量计算
# 效应量计算示例
# Eta-squared (η²)
eta_sq <- summary(anova_result)[[1]]["Sum Sq"][1,1] /
sum(summary(anova_result)[[1]]["Sum Sq"][,1])
cat("Eta-squared:", eta_sq, "\n")## Eta-squared: 0.1725334# Partial eta-squared
# 使用rstatix
rstatix::eta_squared(anova_result)## drug
## 0.1725334# Omega-squared（更准确的估计）
DescTools::EtaSq(anova_result, type = 2)## eta.sq eta.sq.part
## drug 0.1725334 0.1725334# Cohen's f
cohens_f <- sqrt(eta_sq / (1 - eta_sq))
cat("Cohen's f:", cohens_f, "\n")## Cohen's f: 0.4566267第八章：相关分析与基础回归
本章介绍变量间关系的分析方法，包括相关分析和线性回归。相关分析用于度量变量间的关联程度，回归分析则用于建立变量间的数学模型。
8.1 Pearson相关分析
数学基础
Pearson相关系数：
\[r = \frac{\sum(X_i - \bar{X})(Y_i - \bar{Y})}{\sqrt{\sum(X_i - \bar{X})^2 \sum(Y_i - \bar{Y})^2}}\]
性质： - 取值范围：[-1, 1] - r = 0：无线性相关 - r > 0：正相关 - r < 0：负相关 - |r| 接近1：强相关
假设条件： 1. 两个变量都是连续变量 2. 两个变量服从二元正态分布 3. 关系是线性的 4. 无显著异常值
# Pearson相关分析示例
set.seed(111111)
# 生成相关数据
n <- 50
x <- rnorm(n, mean = 100, sd = 15)
y <- 0.7 * x + rnorm(n, mean = 0, sd = 10)
# 散点图
plot(x, y, main = "散点图", xlab = "X", ylab = "Y", pch = 16, col = "steelblue")
abline(lm(y ~ x), col = "red", lwd = 2)
# 计算Pearson相关系数
cor_xy <- cor(x, y)
cat("Pearson相关系数:", cor_xy, "\n")## Pearson相关系数: 0.7931909# 相关系数的假设检验
cor_test <- cor.test(x, y)
cor_test##
## Pearson's product-moment correlation
##
## data: x and y
## t = 9.024, df = 48, p-value = 6.537e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6607185 0.8777471
## sample estimates:
## cor
## 0.7931909# 提取结果
cat("相关系数:", cor_test$estimate, "\n")## 相关系数: 0.7931909cat("95%置信区间:", cor_test$conf.int, "\n")## 95%置信区间: 0.6607185 0.8777471cat("t统计量:", cor_test$statistic, "\n")## t统计量: 9.024023cat("p值:", cor_test$p.value, "\n")## p值: 6.537419e-12# 相关系数矩阵
set.seed(222222)
data_matrix <- data.frame(
x1 = rnorm(100),
x2 = rnorm(100),
x3 = rnorm(100)
)
data_matrix$x4 <- 0.5 * data_matrix$x1 + 0.3 * data_matrix$x2 + rnorm(100, sd = 0.5)
# 相关系数矩阵
cor_matrix <- cor(data_matrix)
cor_matrix## x1 x2 x3 x4
## x1 1.00000000 0.19367575 0.04872119 0.67136898
## x2 0.19367575 1.00000000 -0.03092974 0.51238135
## x3 0.04872119 -0.03092974 1.00000000 -0.04532995
## x4 0.67136898 0.51238135 -0.04532995 1.00000000# 可视化相关矩阵
library(corrplot)## corrplot 0.95 loadedcorrplot(cor_matrix, method = "color", type = "upper",
tl.col = "black", addCoef.col = "black")
8.2 相关系数的假设检验与置信区间
数学基础
检验统计量： \[t = r\sqrt{\frac{n-2}{1-r^2}}\]
服从自由度为n-2的t分布。
# 相关系数检验示例
# 不同相关强度的比较
set.seed(333333)
# 强正相关
x1 <- rnorm(50)
y1 <- 0.9 * x1 + rnorm(50, sd = 0.3)
# 弱正相关
x2 <- rnorm(50)
y2 <- 0.3 * x2 + rnorm(50, sd = 0.8)
# 无相关
x3 <- rnorm(50)
y3 <- rnorm(50)
# 比较检验结果
cor_test1 <- cor.test(x1, y1)
cor_test2 <- cor.test(x2, y2)
cor_test3 <- cor.test(x3, y3)
cat("强正相关: r =", cor_test1$estimate, "p =", cor_test1$p.value, "\n")## 强正相关: r = 0.9537358 p = 1.060319e-26cat("弱正相关: r =", cor_test2$estimate, "p =", cor_test2$p.value, "\n")## 弱正相关: r = 0.45306 p = 0.000953545cat("无相关: r =", cor_test3$estimate, "p =", cor_test3$p.value, "\n")## 无相关: r = 0.06427032 p = 0.657459# 可视化比较
par(mfrow = c(1, 3))
plot(x1, y1, main = paste("r =", round(cor_test1$estimate, 3)), pch = 16)
plot(x2, y2, main = paste("r =", round(cor_test2$estimate, 3)), pch = 16)
plot(x3, y3, main = paste("r =", round(cor_test3$estimate, 3)), pch = 16)
par(mfrow = c(1, 1))
# Fisher Z变换计算置信区间
fisher_z <- atanh(cor_xy)
se_z <- 1 / sqrt(n - 3)
z_lower <- fisher_z - 1.96 * se_z
z_upper <- fisher_z + 1.96 * se_z
ci_lower <- tanh(z_lower)
ci_upper <- tanh(z_upper)
cat("Fisher Z变换置信区间: [", ci_lower, ",", ci_upper, "]\n")## Fisher Z变换置信区间: [ 0.6607155 , 0.8777483 ]8.3 偏相关分析
数学基础
偏相关系数：在控制其他变量影响后，两个变量之间的相关程度。
\[r_{xy.z} = \frac{r_{xy} - r_{xz}r_{yz}}{\sqrt{(1-r_{xz}^2)(1-r_{yz}^2)}}\]
# 偏相关分析示例
set.seed(444444)
n <- 100
# 生成三个相关变量
z <- rnorm(n, mean = 50, sd = 10)
x <- 0.7 * z + rnorm(n, sd = 5)
y <- 0.6 * z + rnorm(n, sd = 5)
# 创建数据框
partial_data <- data.frame(x = x, y = y, z = z)
# 简单相关
cor(partial_data)## x y z
## x 1.0000000 0.5364063 0.8320560
## y 0.5364063 1.0000000 0.7420707
## z 0.8320560 0.7420707 1.0000000# x和y的简单相关
cor_xy <- cor(x, y)
cat("x和y的简单相关:", cor_xy, "\n")## x和y的简单相关: 0.5364063# 偏相关（控制z）
# 使用ppcor包
library(ppcor)
pcor_result <- pcor(partial_data)
pcor_result$estimate## x y z
## x 1.0000000 -0.2179486 0.7671659
## y -0.2179486 1.0000000 0.6317601
## z 0.7671659 0.6317601 1.0000000# 提取偏相关系数
cat("控制z后，x和y的偏相关:", pcor_result$estimate["x", "y"], "\n")## 控制z后，x和y的偏相关: -0.2179486# 偏相关检验
pcor.test(x, y, z)## estimate p.value statistic n gp Method
## 1 -0.2179486 0.03022388 -2.199418 100 1 pearson# 可视化
pairs(partial_data, main = "变量间散点图矩阵")
8.4 简单线性回归
数学基础
回归模型： \[Y = \beta_0 + \beta_1 X + \epsilon\]
最小二乘估计： \[\hat{\beta}_1 = \frac{\sum(X_i - \bar{X})(Y_i - \bar{Y})}{\sum(X_i - \bar{X})^2}\]
\[\hat{\beta}_0 = \bar{Y} - \hat{\beta}_1\bar{X}\]
# 简单线性回归示例
# 场景：研究身高与体重的关系
set.seed(555555)
height <- rnorm(50, mean = 170, sd = 8)
weight <- 0.8 * height - 60 + rnorm(50, sd = 5)
# 创建数据框
body_data <- data.frame(height = height, weight = weight)
# 散点图
plot(height, weight, main = "身高与体重的关系",
xlab = "身高", ylab = "体重", pch = 16, col = "steelblue")
# 拟合线性回归模型
lm_model <- lm(weight ~ height, data = body_data)
lm_model##
## Call:
## lm(formula = weight ~ height, data = body_data)
##
## Coefficients:
## (Intercept) height
## -56.038 0.778# 回归直线
abline(lm_model, col = "red", lwd = 2)
# 详细结果
summary(lm_model)##
## Call:
## lm(formula = weight ~ height, data = body_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.775 -3.748 1.351 4.143 8.071
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -56.03843 17.01130 -3.294 0.00186 **
## height 0.77799 0.09877 7.877 3.38e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.415 on 48 degrees of freedom
## Multiple R-squared: 0.5638, Adjusted R-squared: 0.5547
## F-statistic: 62.04 on 1 and 48 DF, p-value: 3.383e-10# 提取系数
coefficients(lm_model)## (Intercept) height
## -56.0384349 0.7779929# 置信区间
confint(lm_model)## 2.5 % 97.5 %
## (Intercept) -90.2419427 -21.8349271
## height 0.5794004 0.9765855# 回归系数的解释
cat("截距:", coefficients(lm_model)[1], "\n")## 截距: -56.03843cat("斜率:", coefficients(lm_model)[2], "\n")## 斜率: 0.7779929cat("解释：身高每增加1cm，体重平均增加",
round(coefficients(lm_model)[2], 2), "kg\n")## 解释：身高每增加1cm，体重平均增加 0.78 kg8.5 回归系数的假设检验
数学基础
t检验： \[t = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)}\]
决定系数R²： \[R^2 = \frac{SS_{regression}}{SS_{total}} = 1 - \frac{SS_{residual}}{SS_{total}}\]
# 回归系数检验示例
# 模型摘要
model_summary <- summary(lm_model)
# F检验（整体显著性）
model_summary$fstatistic## value numdf dendf
## 62.04286 1.00000 48.00000# 提取p值
pf(model_summary$fstatistic[1],
model_summary$fstatistic[2],
model_summary$fstatistic[3],
lower.tail = FALSE)## value
## 3.383059e-10# R²和调整R²
cat("R²:", model_summary$r.squared, "\n")## R²: 0.5638063cat("调整R²:", model_summary$adj.r.squared, "\n")## 调整R²: 0.554719# 残差标准误
cat("残差标准误:", model_summary$sigma, "\n")## 残差标准误: 5.414723# 系数表
model_summary$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -56.0384349 17.01129838 -3.294189 1.859111e-03
## height 0.7779929 0.09877106 7.876729 3.383059e-10# 使用broom包整理结果
library(broom)
tidy(lm_model)## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -56.0 17.0 -3.29 1.86e- 3
## 2 height 0.778 0.0988 7.88 3.38e-10glance(lm_model)## # A tibble: 1 × 12
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.564 0.555 5.41 62.0 3.38e-10 1 -154. 315. 321.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>8.6 回归诊断基础
残差分析与异常点识别
# 回归诊断示例
# 残差分析
par(mfrow = c(2, 2))
plot(lm_model)
par(mfrow = c(1, 1))
# 提取诊断统计量
diagnostics <- augment(lm_model)
head(diagnostics)## # A tibble: 6 × 8
## weight height .fitted .resid .hat .sigma .cooksd .std.resid
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 83.2 169. 75.2 8.07 0.0239 5.34 0.0279 1.51
## 2 76.1 174. 79.2 -3.11 0.0210 5.45 0.00361 -0.580
## 3 88.2 180. 83.8 4.35 0.0400 5.43 0.0140 0.821
## 4 85.1 173. 78.7 6.35 0.0205 5.39 0.0147 1.18
## 5 85.0 183. 86.0 -0.984 0.0566 5.47 0.00105 -0.187
## 6 85.7 177. 81.7 4.02 0.0282 5.44 0.00823 0.753# 残差
residuals(lm_model)## 1 2 3 4 5 6
## 8.07105975 -3.10841787 4.35494657 6.34796049 -0.98433690 4.01932230
## 7 8 9 10 11 12
## 0.77936225 0.65445868 0.45311388 -5.33170000 1.41843633 -10.86676289
## 13 14 15 16 17 18
## 5.96328697 -5.12153466 5.95673597 -2.06528622 5.36220703 2.70322455
## 19 20 21 22 23 24
## 3.04767940 4.17806352 -9.24528044 -8.23728182 -3.96093892 -7.65683211
## 25 26 27 28 29 30
## 4.16058354 -7.71815307 3.94227627 5.02499606 7.09754238 6.35363091
## 31 32 33 34 35 36
## 1.74261734 1.28441220 -10.16583885 2.37686092 -4.71091345 -9.77836835
## 37 38 39 40 41 42
## 3.17719223 0.57784113 4.59433949 -5.64918821 -11.77487940 3.04018748
## 43 44 45 46 47 48
## 2.90806488 -2.36558877 2.01520955 0.06689361 -2.36764488 4.72141066
## 49 50
## 4.08959830 0.62543217# 标准化残差
rstandard(lm_model)## 1 2 3 4 5 6
## 1.50871289 -0.58018808 0.82084654 1.18453383 -0.18716361 0.75299828
## 7 8 9 10 11 12
## 0.14644224 0.12455545 0.08517070 -1.01552521 0.31036406 -2.02826074
## 13 14 15 16 17 18
## 1.12240406 -0.95754002 1.12840310 -0.38670764 1.00917510 0.50517043
## 19 20 21 22 23 24
## 0.57136552 0.78367883 -1.72555620 -1.53700156 -0.73999903 -1.45430761
## 25 26 27 28 29 30
## 0.77619507 -1.45685078 0.74523615 0.96946295 1.38423615 1.18532460
## 31 32 33 34 35 36
## 0.32582715 0.24745478 -1.89956664 0.44913897 -0.87910444 -1.83001057
## 37 38 39 40 41 42
## 0.59306217 0.10972123 0.86389336 -1.05413631 -2.19668240 0.57226561
## 43 44 45 46 47 48
## 0.54432397 -0.44188486 0.37845540 0.01247989 -0.44262163 0.88280819
## 49 50
## 0.76356334 0.11796197# 学生化残差
rstudent(lm_model)## 1 2 3 4 5 6
## 1.52962312 -0.57613639 0.81801268 1.18964663 -0.18527134 0.74955352
## 7 8 9 10 11 12
## 0.14494115 0.12327109 0.08428520 -1.01586343 0.30742271 -2.09898305
## 13 14 15 16 17 18
## 1.12551901 -0.95669447 1.13169818 -0.38325571 1.00937307 0.50121470
## 19 20 21 22 23 24
## 0.56731498 0.78048169 -1.76304520 -1.55977414 -0.73646308 -1.47187147
## 25 26 27 28 29 30
## 0.77293326 -1.47456432 0.74173597 0.96884328 1.39792718 1.19046476
## 31 32 33 34 35 36
## 0.32277240 0.24501990 -1.95457724 0.44537267 -0.87698752 -1.87752892
## 37 38 39 40 41 42
## 0.58901391 0.10858590 0.86157127 -1.05538559 -2.29193400 0.56821482
## 43 44 45 46 47 48
## 0.54029420 -0.43814977 0.37505238 0.01234922 -0.43888329 0.88074321
## 49 50
## 0.76019865 0.11674365# 杠杆值（hat values）
hatvalues(lm_model)## 1 2 3 4 5 6 7
## 0.02389756 0.02098654 0.03996030 0.02046272 0.05660797 0.02822624 0.03396364
## 8 9 10 11 12 13 14
## 0.05835658 0.03465663 0.05984931 0.28759935 0.02096034 0.03723470 0.02426062
## 15 16 17 18 19 20 21
## 0.04953609 0.02715768 0.03705479 0.02335694 0.02958395 0.03055915 0.02089463
## 22 23 24 25 26 27 28
## 0.02035956 0.02280449 0.05456174 0.02002679 0.04270871 0.04554874 0.08366032
## 29 30 31 32 33 34 35
## 0.10330843 0.02002084 0.02438831 0.08110770 0.02315683 0.04480171 0.02056238
## 36 37 38 39 40 41 42
## 0.02619262 0.02110972 0.05401760 0.03534201 0.02045134 0.02000047 0.03738437
## 43 44 45 46 47 48 49
## 0.02648927 0.02252035 0.03292870 0.02006929 0.02407748 0.02443028 0.02159506
## 50
## 0.04120911# Cook距离
cooks.distance(lm_model)## 1 2 3 4 5 6
## 2.786387e-02 3.607945e-03 1.402276e-02 1.465573e-02 1.050989e-03 8.234663e-03
## 7 8 9 10 11 12
## 3.769845e-04 4.807274e-04 1.302132e-04 3.282563e-02 1.944359e-02 4.403679e-02
## 13 14 15 16 17 18
## 2.436104e-02 1.139861e-02 3.318063e-02 2.087304e-03 1.959503e-02 3.051588e-03
## 19 20 21 22 23 24
## 4.976181e-03 9.679797e-03 3.177118e-02 2.454824e-02 6.389564e-03 6.102919e-02
## 25 26 27 28 29 30
## 6.156146e-03 4.734482e-02 1.325197e-02 4.290377e-02 1.103781e-01 1.435192e-02
## 31 32 33 34 35 36
## 1.326934e-03 2.702460e-03 4.276941e-02 4.730778e-03 8.112367e-03 4.503841e-02
## 37 38 39 40 41 42
## 3.792442e-03 3.437190e-04 1.367125e-02 1.160003e-02 4.924010e-02 6.359199e-03
## 43 44 45 46 47 48
## 4.031013e-03 2.249343e-03 2.438459e-03 1.594880e-06 2.416746e-03 9.758270e-03
## 49 50
## 6.434220e-03 2.990359e-04# 识别异常点
# 高杠杆点（大于2p/n）
n <- nrow(body_data)
p <- 2 # 参数个数
leverage_threshold <- 2 * p / n
high_leverage <- which(hatvalues(lm_model) > leverage_threshold)
cat("高杠杆点:", high_leverage, "\n")## 高杠杆点: 11 28 29 32# 影响点（Cook距离 > 4/n）
cook_threshold <- 4 / n
influential <- which(cooks.distance(lm_model) > cook_threshold)
cat("影响点:", influential, "\n")## 影响点: 29# 可视化Cook距离
plot(cooks.distance(lm_model), type = "h",
main = "Cook距离", ylab = "Cook距离")
abline(h = cook_threshold, col = "red", lty = 2)
8.7 多元线性回归基础
数学基础
多元回归模型： \[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_p X_p + \epsilon\]
# 多元线性回归示例
# 场景：预测血压（基于年龄、体重、运动时间）
set.seed(666666)
n <- 100
age <- rnorm(n, mean = 50, sd = 10)
weight <- rnorm(n, mean = 70, sd = 12)
exercise <- rnorm(n, mean = 3, sd = 1.5) # 每周运动小时数
# 血压模型
bp <- 80 + 0.5 * age + 0.3 * weight - 2 * exercise + rnorm(n, sd = 5)
bp_data <- data.frame(
bp = bp,
age = age,
weight = weight,
exercise = exercise
)
# 多元回归
lm_multi <- lm(bp ~ age + weight + exercise, data = bp_data)
summary(lm_multi)##
## Call:
## lm(formula = bp ~ age + weight + exercise, data = bp_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.7378 -3.5813 -0.0857 3.9048 12.0688
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 77.29303 4.22257 18.305 < 2e-16 ***
## age 0.49361 0.05248 9.406 2.82e-15 ***
## weight 0.32083 0.04129 7.771 8.65e-12 ***
## exercise -1.67331 0.34683 -4.825 5.27e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.251 on 96 degrees of freedom
## Multiple R-squared: 0.6259, Adjusted R-squared: 0.6142
## F-statistic: 53.55 on 3 and 96 DF, p-value: < 2.2e-16# 标准化回归系数
library(lm.beta)
lm.beta(lm_multi)##
## Call:
## lm(formula = bp ~ age + weight + exercise, data = bp_data)
##
## Standardized Coefficients::
## (Intercept) age weight exercise
## NA 0.5898026 0.4872280 -0.3012182# 模型比较（使用anova）
lm_age <- lm(bp ~ age, data = bp_data)
lm_age_weight <- lm(bp ~ age + weight, data = bp_data)
lm_full <- lm(bp ~ age + weight + exercise, data = bp_data)
anova(lm_age, lm_age_weight, lm_full)## Analysis of Variance Table
##
## Model 1: bp ~ age
## Model 2: bp ~ age + weight
## Model 3: bp ~ age + weight + exercise
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 98 4938.3
## 2 97 3288.6 1 1649.69 59.833 1.028e-11 ***
## 3 96 2646.9 1 641.77 23.277 5.274e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# AIC比较
AIC(lm_age, lm_age_weight, lm_full)## df AIC
## lm_age 3 679.7489
## lm_age_weight 4 641.0935
## lm_full 5 621.38378.8 回归中分类变量的处理
虚拟变量
# 分类变量回归示例
set.seed(777777)
n <- 60
# 生成数据
treatment <- factor(rep(c("A", "B", "C"), each = 20))
effect <- ifelse(treatment == "A", 10,
ifelse(treatment == "B", 15, 12)) + rnorm(n, sd = 2)
cat_data <- data.frame(
treatment = treatment,
effect = effect
)
# 使用虚拟变量的回归
lm_cat <- lm(effect ~ treatment, data = cat_data)
summary(lm_cat)##
## Call:
## lm(formula = effect ~ treatment, data = cat_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.7255 -1.1945 -0.1065 1.3797 6.4063
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.5484 0.5055 18.887 < 2e-16 ***
## treatmentB 5.2298 0.7149 7.315 9.50e-10 ***
## treatmentC 3.0903 0.7149 4.322 6.26e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.261 on 57 degrees of freedom
## Multiple R-squared: 0.4869, Adjusted R-squared: 0.4689
## F-statistic: 27.05 on 2 and 57 DF, p-value: 5.494e-09# 对比矩阵
contrasts(cat_data$treatment)## B C
## A 0 0
## B 1 0
## C 0 1# 更改参照组
cat_data$treatment <- relevel(cat_data$treatment, ref = "B")
lm_cat2 <- lm(effect ~ treatment, data = cat_data)
summary(lm_cat2)##
## Call:
## lm(formula = effect ~ treatment, data = cat_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.7255 -1.1945 -0.1065 1.3797 6.4063
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.7782 0.5055 29.232 < 2e-16 ***
## treatmentA -5.2298 0.7149 -7.315 9.5e-10 ***
## treatmentC -2.1395 0.7149 -2.993 0.00408 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.261 on 57 degrees of freedom
## Multiple R-squared: 0.4869, Adjusted R-squared: 0.4689
## F-statistic: 27.05 on 2 and 57 DF, p-value: 5.494e-09# 使用emmeans进行边际均值比较
library(emmeans)
emmeans(lm_cat, pairwise ~ treatment)## $emmeans
## treatment emmean SE df lower.CL upper.CL
## A 9.55 0.506 57 8.54 10.6
## B 14.78 0.506 57 13.77 15.8
## C 12.64 0.506 57 11.63 13.7
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## A - B -5.23 0.715 57 -7.315 <0.0001
## A - C -3.09 0.715 57 -4.322 0.0002
## B - C 2.14 0.715 57 2.993 0.0112
##
## P value adjustment: tukey method for comparing a family of 3 estimates8.9 稳健相关方法
# 稳健相关示例
set.seed(888888)
# 生成含异常值的数据
x <- rnorm(50)
y <- 0.7 * x + rnorm(50, sd = 0.5)
# 添加异常值
x[1:3] <- c(5, -5, 6)
y[1:3] <- c(-3, 4, -2)
# Pearson相关（受异常值影响）
cor_pearson <- cor.test(x, y)
# Spearman相关（稳健）
cor_spearman <- cor.test(x, y, method = "spearman")
# Kendall's tau（稳健）
cor_kendall <- cor.test(x, y, method = "kendall")
# 比较
cat("Pearson r:", cor_pearson$estimate, "p =", cor_pearson$p.value, "\n")## Pearson r: -0.3652727 p = 0.009097938cat("Spearman rho:", cor_spearman$estimate, "p =", cor_spearman$p.value, "\n")## Spearman rho: 0.4287635 p = 0.002067728cat("Kendall tau:", cor_kendall$estimate, "p =", cor_kendall$p.value, "\n")## Kendall tau: 0.3714286 p = 0.0001412265# 可视化
plot(x, y, main = "含异常值的数据", pch = 16)
points(x[1:3], y[1:3], col = "red", pch = 16, cex = 2)
legend("topleft", legend = "异常值", col = "red", pch = 16)
第九章：样本量与功效分析
功效分析是研究设计的重要组成部分，帮助确定研究所需的样本量，并评估检测到真实效应的能力。
9.1 功效分析基本概念
数学基础
关键概念：
- 效应量（Effect Size）：效应的实际大小，独立于样本量
- 显著性水平（α）：第I类错误概率，通常设为0.05
- 功效（Power）：正确拒绝错误H₀的概率，1-β
- 样本量（n）：研究中观测值的数量
四者关系：已知任意三个，可以计算第四个。
# 功效分析基本概念演示
library(pwr)
# 功效曲线
effect_sizes <- seq(0.1, 1, by = 0.1)
n_values <- c(20, 50, 100, 200)
power_curves <- data.frame()
for (n in n_values) {
for (d in effect_sizes) {
power <- pwr.t.test(d = d, n = n, sig.level = 0.05,
type = "two.sample")$power
power_curves <- rbind(power_curves,
data.frame(n = n, d = d, power = power))
}
}
# 可视化功效曲线
ggplot2::ggplot(power_curves, ggplot2::aes(x = d, y = power, color = factor(n))) +
ggplot2::geom_line(size = 1) +
ggplot2::geom_hline(yintercept = 0.8, linetype = "dashed", color = "red") +
ggplot2::labs(title = "功效曲线",
x = "效应量",
y = "功效",
color = "样本量") +
ggplot2::theme_minimal()
9.2 单样本t检验的功效与样本量
# 单样本t检验功效分析示例
# 场景：检验血压是否不同于标准值
# 已知效应量，计算所需样本量
# 假设效应量d = 0.5（中等效应）
pwr_one_sample <- pwr.t.test(
d = 0.5, # 效应量
sig.level = 0.05, # 显著性水平
power = 0.8, # 期望功效
type = "one.sample"
)
pwr_one_sample##
## One-sample t test power calculation
##
## n = 33.36713
## d = 0.5
## sig.level = 0.05
## power = 0.8
## alternative = two.sided# 已知样本量，计算功效
pwr_one_sample2 <- pwr.t.test(
d = 0.5,
n = 30,
sig.level = 0.05,
type = "one.sample"
)
pwr_one_sample2##
## One-sample t test power calculation
##
## n = 30
## d = 0.5
## sig.level = 0.05
## power = 0.7539647
## alternative = two.sided# 已知样本量和功效，计算可检测的最小效应量
pwr_one_sample3 <- pwr.t.test(
n = 50,
power = 0.8,
sig.level = 0.05,
type = "one.sample"
)
pwr_one_sample3##
## One-sample t test power calculation
##
## n = 50
## d = 0.4041852
## sig.level = 0.05
## power = 0.8
## alternative = two.sided# 可视化
plot(pwr_one_sample)
9.3 两样本t检验的功效与样本量
# 两样本t检验功效分析示例
# 场景：比较两种药物的疗效
# 计算所需样本量
pwr_two_sample <- pwr.t.test(
d = 0.5, # 中等效应量
sig.level = 0.05,
power = 0.8,
type = "two.sample"
)
pwr_two_sample##
## Two-sample t test power calculation
##
## n = 63.76561
## d = 0.5
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* group# 配对样本t检验
pwr_paired <- pwr.t.test(
d = 0.5,
sig.level = 0.05,
power = 0.8,
type = "paired"
)
pwr_paired##
## Paired t test power calculation
##
## n = 33.36713
## d = 0.5
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number of *pairs*# 比较独立样本与配对样本所需样本量
cat("独立样本每组需要:", ceiling(pwr_two_sample$n), "人\n")## 独立样本每组需要: 64 人cat("配对样本需要:", ceiling(pwr_paired$n), "对\n")## 配对样本需要: 34 对# 不等样本量（已知n1，计算n2）
pwr_unequal <- pwr.t2n.test(
n1 = 50,
n2 = NULL, # 需要计算n2
d = 0.5,
sig.level = 0.05,
power = 0.8
)
pwr_unequal##
## t test power calculation
##
## n1 = 50
## n2 = 87.70891
## d = 0.5
## sig.level = 0.05
## power = 0.8
## alternative = two.sided9.4 方差分析的功效与样本量
# ANOVA功效分析示例
# 场景：比较三种治疗方法
# 计算所需样本量
# 效应量f = 0.25（中等效应）
pwr_anova <- pwr.anova.test(
k = 3, # 组数
f = 0.25, # 效应量
sig.level = 0.05,
power = 0.8
)
pwr_anova##
## Balanced one-way analysis of variance power calculation
##
## k = 3
## n = 52.3966
## f = 0.25
## sig.level = 0.05
## power = 0.8
##
## NOTE: n is number in each group# 计算给定样本量的功效
pwr_anova2 <- pwr.anova.test(
k = 3,
n = 20, # 每组样本量
f = 0.25,
sig.level = 0.05
)
pwr_anova2##
## Balanced one-way analysis of variance power calculation
##
## k = 3
## n = 20
## f = 0.25
## sig.level = 0.05
## power = 0.3744311
##
## NOTE: n is number in each group# Cohen's f的解释
# f = 0.1: 小效应
# f = 0.25: 中等效应
# f = 0.4: 大效应
# 从η²计算f
eta_sq <- 0.06 # 假设η² = 0.06
f_from_eta <- sqrt(eta_sq / (1 - eta_sq))
cat("η² =", eta_sq, "对应的f =", f_from_eta, "\n")## η² = 0.06 对应的f = 0.25264569.5 比例检验的样本量计算
# 比例检验功效分析示例
# 场景：比较两组治愈率
# 单样本比例检验
# 检验比例是否不同于0.5
pwr_prop_single <- pwr.p.test(
h = 0.3, # 效应量h
sig.level = 0.05,
power = 0.8
)
pwr_prop_single##
## proportion power calculation for binomial distribution (arcsine transformation)
##
## h = 0.3
## n = 87.20956
## sig.level = 0.05
## power = 0.8
## alternative = two.sided# 两样本比例检验
pwr_prop_two <- pwr.2p.test(
h = 0.3,
sig.level = 0.05,
power = 0.8
)
pwr_prop_two##
## Difference of proportion power calculation for binomial distribution (arcsine transformation)
##
## h = 0.3
## n = 174.4191
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: same sample sizes# 从比例计算效应量h
p1 <- 0.6
p2 <- 0.4
h <- ES.h(p1, p2)
cat("p1 =", p1, "p2 =", p2, "效应量h =", h, "\n")## p1 = 0.6 p2 = 0.4 效应量h = 0.4027158# 使用计算出的h计算样本量
pwr.2p.test(h = h, sig.level = 0.05, power = 0.8)##
## Difference of proportion power calculation for binomial distribution (arcsine transformation)
##
## h = 0.4027158
## n = 96.79194
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: same sample sizes9.6 相关性检验的样本量
# 相关性检验功效分析示例
# 场景：检验两个变量的相关性
# 计算所需样本量
pwr_cor <- pwr.r.test(
r = 0.3, # 期望相关系数
sig.level = 0.05,
power = 0.8
)
pwr_cor##
## approximate correlation power calculation (arctangh transformation)
##
## n = 84.07364
## r = 0.3
## sig.level = 0.05
## power = 0.8
## alternative = two.sided# 计算给定样本量的功效
pwr_cor2 <- pwr.r.test(
r = 0.3,
n = 50,
sig.level = 0.05
)
pwr_cor2##
## approximate correlation power calculation (arctangh transformation)
##
## n = 50
## r = 0.3
## sig.level = 0.05
## power = 0.5715558
## alternative = two.sided# 可检测的最小相关系数
pwr_cor3 <- pwr.r.test(
n = 100,
sig.level = 0.05,
power = 0.8
)
pwr_cor3##
## approximate correlation power calculation (arctangh transformation)
##
## n = 100
## r = 0.275866
## sig.level = 0.05
## power = 0.8
## alternative = two.sided# 相关系数效应量解释
# r = 0.1: 小效应
# r = 0.3: 中等效应
# r = 0.5: 大效应9.7 线性回归的样本量考虑
# 线性回归功效分析示例
# 场景：多元回归
# 计算所需样本量
# f² = R²/(1-R²)
# 假设R² = 0.15，则f² = 0.15/0.85 ≈ 0.176
pwr_reg <- pwr.f2.test(
u = 3, # 预测变量个数
v = NULL, # 残差自由度（待计算）
f2 = 0.176, # 效应量f²
sig.level = 0.05,
power = 0.8
)
pwr_reg##
## Multiple regression power calculation
##
## u = 3
## v = 61.98755
## f2 = 0.176
## sig.level = 0.05
## power = 0.8# 样本量 = u + v + 1
n_needed <- 3 + pwr_reg$v + 1
cat("所需样本量:", ceiling(n_needed), "\n")## 所需样本量: 66# 经验法则
# 每个预测变量至少10-20个观测值
# n ≥ 10k 或 n ≥ 20k（k为预测变量数）
# Green公式
# n ≥ 50 + 8k（用于检验整体回归）
# n ≥ 104 + k（用于检验单个系数）
k <- 3
n_green1 <- 50 + 8 * k
n_green2 <- 104 + k
cat("Green公式（整体）:", n_green1, "\n")## Green公式（整体）: 74cat("Green公式（单个系数）:", n_green2, "\n")## Green公式（单个系数）: 1079.8 使用pwr包进行功效分析总结
# pwr包函数汇总
# 创建汇总表
pwr_functions <- data.frame(
函数 = c("pwr.t.test()", "pwr.t2n.test()", "pwr.anova.test()",
"pwr.p.test()", "pwr.2p.test()", "pwr.r.test()", "pwr.f2.test()"),
用途 = c("单样本/两样本/配对t检验",
"不等样本量两样本t检验",
"单因素ANOVA",
"单样本比例检验",
"两样本比例检验",
"相关性检验",
"线性回归"),
效应量 = c("d (Cohen's d)", "d", "f", "h", "h", "r", "f²")
)
pwr_functions## 函数 用途 效应量
## 1 pwr.t.test() 单样本/两样本/配对t检验 d (Cohen's d)
## 2 pwr.t2n.test() 不等样本量两样本t检验 d
## 3 pwr.anova.test() 单因素ANOVA f
## 4 pwr.p.test() 单样本比例检验 h
## 5 pwr.2p.test() 两样本比例检验 h
## 6 pwr.r.test() 相关性检验 r
## 7 pwr.f2.test() 线性回归 f²# 效应量解释标准
effect_sizes <- data.frame(
效应量 = c("d", "d", "d", "f", "f", "f", "r", "r", "r"),
大小 = rep(c("小", "中", "大"), 3),
值 = c(0.2, 0.5, 0.8, 0.1, 0.25, 0.4, 0.1, 0.3, 0.5)
)
effect_sizes## 效应量 大小 值
## 1 d 小 0.20
## 2 d 中 0.50
## 3 d 大 0.80
## 4 f 小 0.10
## 5 f 中 0.25
## 6 f 大 0.40
## 7 r 小 0.10
## 8 r 中 0.30
## 9 r 大 0.50第十章：统计模拟与重采样方法
统计模拟和重采样方法是现代统计学的重要工具，它们通过计算机模拟来研究统计性质、估计参数和进行假设检验，特别适用于理论推导困难的情况。
10.1 蒙特卡洛模拟基础
数学基础
蒙特卡洛方法：通过大量随机模拟来估计统计量或验证统计性质。
基本步骤： 1. 定义统计模型 2. 生成随机样本 3. 计算统计量 4. 重复多次 5. 分析结果分布
# 蒙特卡洛模拟基础示例
# 验证样本均值的性质
set.seed(12345)
# 模拟参数
n_sim <- 10000 # 模拟次数
n_sample <- 30 # 每次抽样的样本量
true_mean <- 100 # 真实均值
true_sd <- 15 # 真实标准差
# 存储结果
sample_means <- numeric(n_sim)
# 进行模拟
for (i in 1:n_sim) {
sample_data <- rnorm(n_sample, mean = true_mean, sd = true_sd)
sample_means[i] <- mean(sample_data)
}
# 分析结果
cat("模拟次数:", n_sim, "\n")## 模拟次数: 10000cat("样本均值的理论期望:", true_mean, "\n")## 样本均值的理论期望: 100cat("模拟均值的平均值:", mean(sample_means), "\n")## 模拟均值的平均值: 100.03cat("样本均值的理论标准误:", true_sd / sqrt(n_sample), "\n")## 样本均值的理论标准误: 2.738613cat("模拟均值的标准差:", sd(sample_means), "\n")## 模拟均值的标准差: 2.734596# 可视化
hist(sample_means, breaks = 50, probability = TRUE,
main = "样本均值的抽样分布",
xlab = "样本均值", col = "lightblue")
curve(dnorm(x, mean = true_mean, sd = true_sd/sqrt(n_sample)),
add = TRUE, col = "red", lwd = 2)
legend("topright", legend = c("模拟分布", "理论分布"),
col = c("lightblue", "red"), lwd = c(10, 2))
10.2 模拟验证统计性质
验证中心极限定理
# 验证中心极限定理
set.seed(23456)
# 从指数分布（非正态）抽样
n_sim <- 10000
n_values <- c(5, 10, 30, 100)
par(mfrow = c(2, 2))
for (n in n_values) {
sample_means <- numeric(n_sim)
for (i in 1:n_sim) {
sample_means[i] <- mean(rexp(n, rate = 1))
}
hist(sample_means, breaks = 50, probability = TRUE,
main = paste("n =", n),
xlab = "样本均值", col = "lightgreen")
# 理论正态分布
curve(dnorm(x, mean = 1, sd = 1/sqrt(n)),
add = TRUE, col = "red", lwd = 2)
}
par(mfrow = c(1, 1))验证t检验的性质
# 验证t检验的第I类错误率
set.seed(34567)
n_sim <- 10000
n <- 20
alpha <- 0.05
type1_errors <- numeric(n_sim)
for (i in 1:n_sim) {
# 从正态分布生成数据（H₀为真）
data <- rnorm(n, mean = 0, sd = 1)
# 进行t检验
test_result <- t.test(data, mu = 0)
# 记录是否拒绝H₀
type1_errors[i] <- test_result$p.value < alpha
}
cat("理论第I类错误率:", alpha, "\n")## 理论第I类错误率: 0.05cat("模拟第I类错误率:", mean(type1_errors), "\n")## 模拟第I类错误率: 0.0547cat("95%置信区间:",
binom.test(sum(type1_errors), n_sim, p = alpha)$conf.int, "\n")## 95%置信区间: 0.05032338 0.0593381610.3 自举法基础
数学基础
自举法（Bootstrap）：通过有放回重抽样来估计统计量的分布。
基本步骤： 1. 从原始样本中有放回抽取n个观测值 2. 计算统计量 3. 重复B次 4. 使用统计量分布进行推断
# 自举法基础示例
# 估计中位数的标准误
set.seed(45678)
# 原始数据
original_data <- rnorm(50, mean = 100, sd = 20)
# 自举过程
n_boot <- 1000
boot_medians <- numeric(n_boot)
for (i in 1:n_boot) {
boot_sample <- sample(original_data, size = length(original_data),
replace = TRUE)
boot_medians[i] <- median(boot_sample)
}
# 结果
cat("原始样本中位数:", median(original_data), "\n")## 原始样本中位数: 103.7049cat("自举中位数的均值:", mean(boot_medians), "\n")## 自举中位数的均值: 100.7825cat("自举估计的标准误:", sd(boot_medians), "\n")## 自举估计的标准误: 4.569661# 可视化
hist(boot_medians, breaks = 30, probability = TRUE,
main = "中位数的自举分布",
xlab = "中位数", col = "lightcoral")
abline(v = median(original_data), col = "blue", lwd = 2, lty = 2)
legend("topright", legend = "原始中位数", col = "blue", lty = 2, lwd = 2)
10.4 使用boot包进行自举
# 使用boot包示例
library(boot)
# 定义统计量函数
median_fun <- function(data, indices) {
d <- data[indices]
return(median(d))
}
# 进行自举
boot_result <- boot(data = original_data,
statistic = median_fun,
R = 1000)
boot_result##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = original_data, statistic = median_fun, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 103.7049 -2.820728 4.424422# 查看自举分布
plot(boot_result)
# 计算置信区间
# 正态近似
boot.ci(boot_result, type = "norm")## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_result, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 97.9, 115.2 )
## Calculations and Intervals on Original Scale# 百分位数法
boot.ci(boot_result, type = "perc")## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_result, type = "perc")
##
## Intervals :
## Level Percentile
## 95% ( 91.1, 106.0 )
## Calculations and Intervals on Original Scale# BCa法（偏差校正加速）
boot.ci(boot_result, type = "bca")## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_result, type = "bca")
##
## Intervals :
## Level BCa
## 95% ( 91.1, 105.8 )
## Calculations and Intervals on Original Scale10.5 自举法置信区间
不同置信区间方法比较
# 自举置信区间方法比较
# 估计相关系数的置信区间
set.seed(56789)
n <- 50
x <- rnorm(n)
y <- 0.6 * x + rnorm(n, sd = 0.8)
# 定义相关系数函数
cor_fun <- function(data, indices) {
d <- data[indices, ]
cor(d[, 1], d[, 2])
}
# 创建数据框
cor_data <- data.frame(x = x, y = y)
# 自举
boot_cor <- boot(data = cor_data, statistic = cor_fun, R = 2000)
boot_cor##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = cor_data, statistic = cor_fun, R = 2000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.6545179 -0.002059037 0.06285708# 不同方法的置信区间
cat("样本相关系数:", cor(x, y), "\n")## 样本相关系数: 0.6545179cat("\n各种置信区间方法:\n")##
## 各种置信区间方法:cat("正态近似:\n")## 正态近似:print(boot.ci(boot_cor, type = "norm"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_cor, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.5334, 0.7798 )
## Calculations and Intervals on Original Scalecat("\n百分位数法:\n")##
## 百分位数法:print(boot.ci(boot_cor, type = "perc"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_cor, type = "perc")
##
## Intervals :
## Level Percentile
## 95% ( 0.5082, 0.7585 )
## Calculations and Intervals on Original Scalecat("\nBCa法:\n")##
## BCa法:print(boot.ci(boot_cor, type = "bca"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_cor, type = "bca")
##
## Intervals :
## Level BCa
## 95% ( 0.4944, 0.7518 )
## Calculations and Intervals on Original Scale# 与传统方法比较
cat("\n传统Fisher Z变换置信区间:\n")##
## 传统Fisher Z变换置信区间:cor_test <- cor.test(x, y)
print(cor_test$conf.int)## [1] 0.4599689 0.7891038
## attr(,"conf.level")
## [1] 0.9510.6 参数自举与非参数自举
非参数自举
# 非参数自举示例
# 从数据本身重抽样
set.seed(67890)
sample_data <- rnorm(30, mean = 50, sd = 10)
# 非参数自举
nonparam_boot <- function(data, n_boot = 1000) {
n <- length(data)
boot_stats <- numeric(n_boot)
for (i in 1:n_boot) {
boot_sample <- sample(data, size = n, replace = TRUE)
boot_stats[i] <- mean(boot_sample)
}
return(boot_stats)
}
np_result <- nonparam_boot(sample_data)
cat("非参数自举均值标准误:", sd(np_result), "\n")## 非参数自举均值标准误: 1.899314参数自举
# 参数自举示例
# 从拟合的分布中生成数据
param_boot <- function(data, n_boot = 1000) {
n <- length(data)
# 估计参数
mu_hat <- mean(data)
sigma_hat <- sd(data)
boot_stats <- numeric(n_boot)
for (i in 1:n_boot) {
# 从估计的分布中生成数据
boot_sample <- rnorm(n, mean = mu_hat, sd = sigma_hat)
boot_stats[i] <- mean(boot_sample)
}
return(boot_stats)
}
p_result <- param_boot(sample_data)
cat("参数自举均值标准误:", sd(p_result), "\n")## 参数自举均值标准误: 1.862321# 比较
cat("理论标准误:", sd(sample_data) / sqrt(length(sample_data)), "\n")## 理论标准误: 1.84682610.7 置换检验基础
数学基础
置换检验（Permutation Test）：通过随机排列数据来构建检验统计量的零分布。
基本步骤： 1. 计算原始数据的检验统计量 2. 随机排列组标签 3. 计算排列后的统计量 4. 重复多次 5. 比较原始统计量与置换分布
# 置换检验示例
# 比较两组均值
set.seed(78901)
# 生成数据
group1 <- rnorm(20, mean = 50, sd = 10)
group2 <- rnorm(20, mean = 55, sd = 10)
all_data <- c(group1, group2)
group_labels <- rep(c(1, 2), each = 20)
# 原始差异
obs_diff <- mean(group2) - mean(group1)
cat("观测到的组间差异:", obs_diff, "\n")## 观测到的组间差异: -0.9269954# 置换检验
n_perm <- 10000
perm_diffs <- numeric(n_perm)
for (i in 1:n_perm) {
perm_labels <- sample(group_labels)
perm_diffs[i] <- mean(all_data[perm_labels == 2]) -
mean(all_data[perm_labels == 1])
}
# 计算p值
p_value <- mean(abs(perm_diffs) >= abs(obs_diff))
cat("置换检验p值:", p_value, "\n")## 置换检验p值: 0.7505# 与t检验比较
t_result <- t.test(group1, group2)
cat("t检验p值:", t_result$p.value, "\n")## t检验p值: 0.747373# 可视化置换分布
hist(perm_diffs, breaks = 50, probability = TRUE,
main = "置换分布",
xlab = "均值差异", col = "lightyellow")
abline(v = obs_diff, col = "red", lwd = 2)
abline(v = -obs_diff, col = "red", lwd = 2)
legend("topright", legend = "观测值", col = "red", lwd = 2)
10.8 使用coin包进行置换检验
# 使用coin包进行置换检验
# 创建数据框
perm_data <- data.frame(
value = c(group1, group2),
group = factor(rep(c("A", "B"), each = 20))
)
# 两独立样本置换检验
coin_result <- coin::oneway_test(value ~ group, data = perm_data,
distribution = approximate(nresample = 10000))
coin_result##
## Approximative Two-Sample Fisher-Pitman Permutation Test
##
## data: value by group (A, B)
## Z = 0.32834, p-value = 0.7537
## alternative hypothesis: true mu is not equal to 0# Wilcoxon-Mann-Whitney置换检验
coin_wilcox <- coin::wilcox_test(value ~ group, data = perm_data,
distribution = approximate(nresample = 10000))
coin_wilcox##
## Approximative Wilcoxon-Mann-Whitney Test
##
## data: value by group (A, B)
## Z = 0.2164, p-value = 0.8436
## alternative hypothesis: true mu is not equal to 0# 多组比较置换检验
set.seed(89012)
group_A <- rnorm(15, mean = 50, sd = 8)
group_B <- rnorm(15, mean = 55, sd = 8)
group_C <- rnorm(15, mean = 60, sd = 8)
anova_perm_data <- data.frame(
value = c(group_A, group_B, group_C),
group = factor(rep(c("A", "B", "C"), each = 15))
)
# Kruskal-Wallis置换检验
coin_kw <- coin::kruskal_test(value ~ group, data = anova_perm_data,
distribution = approximate(nresample = 10000))
coin_kw##
## Approximative Kruskal-Wallis Test
##
## data: value by group (A, B, C)
## chi-squared = 19.831, p-value < 1e-0410.9 交叉验证基础
数学基础
交叉验证：评估模型在独立数据上的表现，防止过拟合。
K折交叉验证： 1. 将数据分成K份 2. 用K-1份训练，1份验证 3. 重复K次 4. 平均结果
# 交叉验证示例
# 线性回归的K折交叉验证
set.seed(90123)
n <- 100
cv_data <- data.frame(
x1 = rnorm(n),
x2 = rnorm(n),
y = 2 * rnorm(n) + rnorm(n) + rnorm(n, sd = 0.5)
)
# K折交叉验证函数
k_fold_cv <- function(data, k = 5) {
n <- nrow(data)
fold_size <- n %/% k
# 随机分配折
folds <- sample(rep(1:k, length.out = n))
mse_values <- numeric(k)
for (i in 1:k) {
# 分割数据
test_idx <- which(folds == i)
train_data <- data[-test_idx, ]
test_data <- data[test_idx, ]
# 训练模型
model <- lm(y ~ x1 + x2, data = train_data)
# 预测并计算误差
predictions <- predict(model, newdata = test_data)
mse_values[i] <- mean((test_data$y - predictions)^2)
}
return(list(
mse = mse_values,
mean_mse = mean(mse_values),
sd_mse = sd(mse_values)
))
}
cv_result <- k_fold_cv(cv_data, k = 5)
cat("各折MSE:", cv_result$mse, "\n")## 各折MSE: 5.632498 2.723833 6.574188 9.723657 4.461611cat("平均MSE:", cv_result$mean_mse, "\n")## 平均MSE: 5.823157cat("MSE标准差:", cv_result$sd_mse, "\n")## MSE标准差: 2.611222# 使用boot包的cv.glm
library(boot)
glm_model <- glm(y ~ x1 + x2, data = cv_data)
cv_glm_result <- cv.glm(cv_data, glm_model, K = 5)
cat("cv.glm估计的预测误差:", cv_glm_result$delta[1], "\n")## cv.glm估计的预测误差: 6.13596610.10 模拟在功效分析中的应用
# 模拟功效分析示例
# 检验复杂设计的功效
set.seed(12345)
# 模拟函数：计算给定设计下的功效
simulate_power <- function(n_per_group, effect_size, n_sim = 1000, alpha = 0.05) {
significant_count <- 0
for (i in 1:n_sim) {
# 生成数据
group1 <- rnorm(n_per_group, mean = 0, sd = 1)
group2 <- rnorm(n_per_group, mean = effect_size, sd = 1)
# 进行检验
test_result <- t.test(group1, group2)
# 记录结果
if (test_result$p.value < alpha) {
significant_count <- significant_count + 1
}
}
return(significant_count / n_sim)
}
# 不同样本量下的功效
sample_sizes <- c(20, 30, 50, 80, 100)
effect_size <- 0.5
powers <- sapply(sample_sizes, function(n) {
simulate_power(n, effect_size, n_sim = 1000)
})
# 结果
power_df <- data.frame(
n = sample_sizes,
power = powers
)
print(power_df)## n power
## 1 20 0.340
## 2 30 0.492
## 3 50 0.696
## 4 80 0.889
## 5 100 0.940# 可视化
plot(sample_sizes, powers, type = "b", pch = 16, col = "steelblue",
main = "功效随样本量变化",
xlab = "每组样本量", ylab = "功效")
abline(h = 0.8, col = "red", lty = 2)
legend("bottomright", legend = "目标功效 0.8", col = "red", lty = 2)
10.11 自举回归分析
# 自举回归分析示例
set.seed(23456)
# 生成数据
n <- 50
x <- rnorm(n, mean = 50, sd = 10)
y <- 2 + 0.5 * x + rnorm(n, sd = 3)
reg_data <- data.frame(x = x, y = y)
# 定义回归系数函数
coef_fun <- function(data, indices) {
d <- data[indices, ]
model <- lm(y ~ x, data = d)
return(coef(model))
}
# 自举
boot_reg <- boot(data = reg_data, statistic = coef_fun, R = 1000)
boot_reg##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = reg_data, statistic = coef_fun, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 1.3620514 0.100526667 1.81270720
## t2* 0.5015121 -0.001924008 0.03349922# 置信区间
cat("截距置信区间:\n")## 截距置信区间:boot.ci(boot_reg, type = "bca", index = 1)## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_reg, type = "bca", index = 1)
##
## Intervals :
## Level BCa
## 95% (-1.475, 6.123 )
## Calculations and Intervals on Original Scalecat("\n斜率置信区间:\n")##
## 斜率置信区间:boot.ci(boot_reg, type = "bca", index = 2)## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_reg, type = "bca", index = 2)
##
## Intervals :
## Level BCa
## 95% ( 0.4224, 0.5547 )
## Calculations and Intervals on Original Scale# 与传统方法比较
lm_model <- lm(y ~ x, data = reg_data)
confint(lm_model)## 2.5 % 97.5 %
## (Intercept) -2.1132767 4.8373794
## x 0.4345668 0.5684574全书总结
本书系统介绍了R语言统计分析的核心内容，涵盖了从描述性统计到高级统计方法的完整知识体系。
内容回顾
| 章节 | 主题 | 核心内容 |
|---|---|---|
| 第一章 | 描述性统计 | 集中趋势、离散程度、分布形状 |
| 第二章 | 概率分布 | 常见分布、随机数生成、正态性检验 |
| 第三章 | 参数估计 | 点估计、置信区间、自举法 |
| 第四章 | 假设检验 | 基本概念、p值、效应量 |
| 第五章 | 参数检验 | t检验、比例检验、方差检验 |
| 第六章 | 非参数检验 | 秩检验、卡方检验、Fisher检验 |
| 第七章 | 方差分析 | 单因素、多因素、重复测量ANOVA |
| 第八章 | 相关与回归 | 相关系数、线性回归、回归诊断 |
| 第九章 | 功效分析 | 样本量计算、功效曲线 |
| 第十章 | 模拟方法 | 蒙特卡洛、自举、置换检验 |
学习建议
- 理论与实践结合：每学完一个概念，用R代码实践
- 理解假设条件：了解每种方法的适用条件
- 关注效应量：不仅仅关注p值，更要关注实际意义
- 掌握多种方法：参数与非参数方法都要掌握
- 善用可视化：图形化展示数据和分析结果
常用R包速查
# 常用统计包汇总
packages_summary <- data.frame(
包名 = c("stats", "dplyr", "ggplot2", "psych", "car",
"pwr", "boot", "coin", "rstatix", "DescTools"),
主要功能 = c("基础统计函数", "数据处理", "可视化", "心理统计", "回归诊断",
"功效分析", "自举法", "置换检验", "现代统计检验", "描述统计工具")
)
packages_summary## 包名 主要功能
## 1 stats 基础统计函数
## 2 dplyr 数据处理
## 3 ggplot2 可视化
## 4 psych 心理统计
## 5 car 回归诊断
## 6 pwr 功效分析
## 7 boot 自举法
## 8 coin 置换检验
## 9 rstatix 现代统计检验
## 10 DescTools 描述统计工具希望本书能帮助您建立扎实的统计分析基础，并在实际研究中灵活运用！