数据科学19 | 统计推断-t分布置信区间
1. t分布
当样本量足够大，总体标准差已知时，根据中心极限定理可以用标准正态分布估计总体均值；t分布适用于小样本估计呈正态分布的总体均值。
当随机变量X满足 时，服从自由度df为n-1的t分布。
注意，如果用?代替S，则X服从标准正态分布。
t分布的置信区间为 ， 为标准误。
使用manipulate( )观察不同自由度的t分布与标准正态分布：
k <- 1000
xvals <- seq(-5, 5, length = k)
myplot <- function(df){
d <- data.frame(y = c(dnorm(xvals), dt(xvals, df)),
x = xvals,
dist = factor(rep(c("Normal", "T"), c(k,k))))
g <- ggplot(d, aes(x = x, y = y))
g <- g + geom_line(size = 2, aes(colour = dist))
g
}
manipulate(myplot(mu), mu = slider(1, 20, step = 1))
与标准正态分布相比，df为1时t分布的峰值更低，两端的“尾巴”更厚。通过左上角设置图标控制df，df变大，t分布的峰值变高，两端的“尾巴”变低，逐渐接近标准正态分布。
使用manipulate( )观察不同自由度的t分布与标准正态分布的分位数：
pvals <- seq(.5, .99, by = .01)
myplot2 <- function(df){
d <- data.frame(n= qnorm(pvals),t=qt(pvals, df), p = pvals)
g <- ggplot(d, aes(x= n, y = t))
g <- g + geom_abline(size = 2, col = "lightblue")
g <- g + geom_line(size = 2, col = "black")
g <- g + geom_vline(xintercept = qnorm(0.975))
g <- g + geom_hline(yintercept = qt(0.975, df))
g
}
manipulate(myplot2(df), df = slider(1, 20, step = 1))
两个分布对称，零点从第50百分位数开始。
标准正态分布的97.5百分位数约为1.96（蓝色参考线）；自由度为2时，t分布的第97.5分位数大于4（黑色曲线）。自由度越大，t分位数越接近于正态分位数。t分位数（黑色曲线）总是在正态分位数（蓝色参考线）之上，意味着t分布的置信区间总是比正态分布的宽。
2. t分布置信区间
- 当自由度很大时，t分布接近标准正态分布，置信区间收敛于标准正态分布的置信区间。
- 偏态分布的数据不满足t分布置信区间的假设，置信区间的中心落在均值处没有意义，可以考虑使用对数处理数据，或使用其他统计量如中位数。
➢配对样本——配对t检验
例：sleep数据集，10名患者服用2种不同安眠药后睡眠时间增加的数据。
两组样本数据来自于同10名患者，两组样本均值不独立。
data(sleep)
head(sleep)
extra group ID
1 0.7 1 1
2 -1.6 1 2
3 -0.2 1 3
4 -1.2 1 4
5 -0.1 1 5
6 3.4 1 6
#extra为增加的睡眠时间，数据分为两组，10名患者ID记为1-10
ggplot(sleep, aes(x=group,y=extra,color=ID))+
geom_point(color="salmon", size=4, alpha=.3)+
geom_line(aes(group=ID))
计算两组差异的均值的置信区间：
g1 <- sleep$extra[1 : 10]
g2 <- sleep$extra[11 : 20]
difference <- g2 - g1 #计算同一患者对两种药物增加睡眠时间的差值
mn <- mean(difference) #计算差值的均值
s <- sd(difference) #计算样本标准差
n <- 10
#公式计算
mn + c(-1, 1) * qt(.975, n-1) * s / sqrt(n)
[1] 0.7001 2.4599
#单样本t检验
t.test(difference)
One Sample t-test
data: difference
t = 4.1, df = 9, p-value = 0.003
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
0.7001 2.4599
sample estimates:
mean of x
1.58
#配对t检验
t.test(g2, g1, paired = TRUE)
Paired t-test
data: g2 and g1
t = 4.1, df = 9, p-value = 0.003
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.7001 2.4599
sample estimates:
mean of the differences
1.58
#配对t检验
t.test(extra ~ I(relevel(group, 2)), paired = TRUE, data = sleep)
Paired t-test
data: extra by I(relevel(group, 2))
t = 4.1, df = 9, p-value = 0.003
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.7001 2.4599
sample estimates:
mean of the differences
1.58
#组合输出得到的置信区间
rbind(
mn + c(-1, 1) * qt(.975, n-1) * s / sqrt(n),
t.test(difference)$conf,
t.test(g2, g1, paired = TRUE)$conf,
t.test(extra ~ I(relevel(group, 2)), paired = TRUE, data = sleep)$conf)
[,1] [,2]
[1,] 0.7001 2.46
[2,] 0.7001 2.46
[3,] 0.7001 2.46
[4,] 0.7001 2.46以上方法得到的置信区间都为(0.70,2.46)。
➢独立样本，方差齐——成组t检验
对于分组独立且来自正态分布的样本，总体方差相等时，?y-?x的(1-?)×100%置信区间为 。
自由度为nx+ny-2，合并方差为 ， 为合并标准差。
两组的方差相同，需要用两个样本的方差来估计总体方差，这正是合并方差的作用。
例：比较8名口服避孕药及21名空白对照患者的血压。
口服避孕药的患者的平均收缩压 =132.86mmHg，标准差 =15.34mmHg。
对照者的平均收缩压 =127.44mmHg，标准差 =18.23mmHg。
两组受试者之间相互独立且样本数不同时，不能使用配对t检验。
用合并标准差估计均值差异的置信区间：
sp<-sqrt((7*15.34^2+20*18.23^2)/(8+21-2))
132.86-127.44+c(-1,1)*qt(.975,27)*sp*(1/8+1/21)^.5
[1] -9.521 20.361区间包括0，所以不能排除两组之间总体差异性为0的可能性。
例：sleep数据集的错误处理：假设数据集中两组样本不配对且方差齐
n1 <- length(g1); n2 <- length(g2)
sp <- sqrt( ((n1 - 1) * sd(g1)^2 + (n2-1) * sd(g2)^2) / (n1 + n2-2)) #计算合并标准差
md <- mean(g2) - mean(g1)
semd<-sp*sqrt(1/n1+1/n2) #计算均值之差的标准误
rbind(
md+c(-1,1)*qt(.975,n1+n2-2)*semd, #公式计算置信区间
t.test(g2, g1, paired = FALSE, var.equal = TRUE)$conf, #成组t检验
t.test(g2, g1, paired = TRUE)$conf #配对t检验
)
[,1] [,2]
[1,] -0.2039 3.364
[2,] -0.2039 3.364
[3,] 0.7001 2.460如果按配对样本计算，置信区间不包括0，但如果忽视样本配对，置信区间会包括0。
例：ChickWeight数据集，包含4种不同饮食对小鸡生长的影响的数据
library(datasets)
data(ChickWeight)
head(ChickWeight)
weight Time Chick Diet
1 42 0 1 1
2 51 2 1 1
3 59 4 1 1
4 64 6 1 1
5 76 8 1 1
6 93 10 1 1
summary(ChickWeight)
weight Time Chick Diet
Min. : 35 Min. : 0.0 13 : 12 1:220
1st Qu.: 63 1st Qu.: 4.0 9 : 12 2:120
Median :103 Median :10.0 20 : 12 3:120
Mean :122 Mean :10.7 10 : 12 4:118
3rd Qu.:164 3rd Qu.:16.0 17 : 12
Max. :373 Max. :21.0 19 : 12
(Other):506
str(ChickWeight)
Classes ‘nfnGroupedData’, ‘nfGroupedData’, ‘groupedData’ and 'data.frame': 578 obs. of 4 variables:
$ weight: num 42 51 59 64 76 93 106 125 149 171 ...
$ Time : num 0 2 4 6 8 10 12 14 16 18 ...
$ Chick : Ord.factor w/ 50 levels "18"<"16"<"15"<..: 15 15 15 15 15 15 15 15 15 15 ...
$ Diet : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
#weight为每只小鸡从出生开始在不同时间点测的体重
#Time为不同的监测时间
#Chick为每只小鸡的编号
#Diet为4种饮食的编号重组数据：
library(reshape2)
##define weight gain or loss
wideCW <- dcast(ChickWeight, Diet + Chick ~ Time, value.var = "weight")
names(wideCW)[-(1 : 2)] <- paste("time", names(wideCW)[-(1 : 2)], sep = "")
library(dplyr)
wideCW <- mutate(wideCW,
gain = time21 - time0)
head(wideCW)
Diet Chick time0 time2 time4 time6 time8 time10 time12 time14 time16 time18 time20 time21 gain
1 1 18 39 35 NA NA NA NA NA NA NA NA NA NA NA
2 1 16 41 45 49 51 57 51 54 NA NA NA NA NA NA
3 1 15 41 49 56 64 68 68 67 68 NA NA NA NA NA
4 1 13 41 48 53 60 65 67 71 70 71 81 91 96 55
5 1 9 42 51 59 68 85 96 90 92 93 100 100 98 56
6 1 20 41 47 54 58 65 73 77 89 98 107 115 117 76画出原始数据：
meanweight<-ChickWeight %>%
group_by(Time,Diet) %>%
summarise(weight = mean(weight)) #按Time统计weight平均值
meanweight
# A tibble: 48 x 3
# Groups: Time [12]
Time Diet weight
<dbl> <fct> <dbl>
1 0 1 41.4
2 0 2 40.7
3 0 3 40.8
4 0 4 41
5 2 1 47.2
6 2 2 49.4
7 2 3 50.4
8 2 4 51.8
9 4 1 56.5
10 4 2 59.8
# … with 38 more rows
ggplot(ChickWeight, aes(x=Time, y=weight, color=Diet)) +
geom_line(aes(group=Chick))+ #画出不同饮食分组下每只小鸡的体重生长曲线
geom_line(data=meanweight,color="black")+ #画出不同饮食分组的体重生长均值曲线
facet_grid(.~Diet)
第1种饮食的末端变异似乎比第4种饮食的末端变异大得多，但第1种饮食中的鸡比第4种饮食中的鸡数量要多，所以很难真正比较变化。观察每组均值，第1种饮食的平均体重增长似乎确实比第4种饮食的平均体重增长慢。
画出每种饮食的小鸡最终体重增长量：
ggplot(wideCW,aes(x=Diet,y=gain,fill=Diet))+
geom_violin(size=1,color="black")+
labs(x="factor(Diet)",fill="factor(Diet)")
比较第1种饮食和第4种饮食的差异：
wideCW14 <- subset(wideCW, Diet %in% c(1, 4))
rbind(
t.test(gain ~ Diet, paired = FALSE, var.equal = TRUE, data = wideCW14)$conf,
t.test(gain ~ Diet, paired = FALSE, var.equal = FALSE, data = wideCW14)$conf)
[,1] [,2]
[1,] -108.1 -14.81
[2,] -104.7 -18.30两组样本相互独立且样本量不同，不能使用配对t检验，paired = FALSE。方差齐或不齐的情况下，置信区间小于0，表明第1种饮食比第4种饮食的体重增加更少。方差是否一致会影响区间。
➢独立样本，方差不齐——校正t检验
对于分组独立且来自正态分布的样本，若方差不齐性不严重时，可以用校正t检验，
?y-?x的95%置信区间可用 计算，其中tdf用自由度 计算。
实际上，方差不齐的独立样本的相关标准化统计量不服从t分布，当其自由度用这种方式计算下才近似t分布。
例：比较8名口服避孕药及21名空白对照患者的血压。
口服避孕药的患者的平均收缩压 =132.86mmHg，标准差 =15.34mmHg。
对照者的平均收缩压 =127.44mmHg，标准差 =18.23mmHg。
df=15.04，t15.04,0.975=2.13。
计算均值之差的置信区间：
132.86 - 127.44 + c(-1, 1) * 2.13 * (15.34^2/8 + 18.23^2/21)^.5
[1] -8.906 19.746R中可以使用t.test(..., var.equal = FALSE)计算。
编辑：李雪纯 冯文清
校审：张健 罗鹏