1  t检验

t检验主要适用于1组或2组的均数的比较,要求数据符合正态性和方差齐性。关于t检验的本质问题本文不做过多探讨,只是学习如何用R语言实现t检验!

在R中进行t检验非常简单,就是t.test()函数,不管是单样本、两样本都是这一个函数。使用R语言进行统计学比SPSS好的一点是,当你需要使用某种检验时,直接打出这个检验的名字一般就会跳出来相关的函数,而我在学习使用SPSS进行各种统计检验时,可能更多的是关注如何点点点(鼠标点击菜单),不利于记住哪种情况使用什么方法。

1.1 单样本t检验

孙振球《医学统计学》第4版例3-5

某医生测量了36名从事铅作业男性工人的血红蛋白含量,算得其均数为130.83g/L,标准差为25.74g/L。问,从事铅作业男性工人的血红蛋白含量均数(μ)是否不等于正常成年男性的均数140 g/L(μ0)?

# 使用foreign包读取SPSS数据
library(foreign)
data3_5 <- read.spss('datasets/例03-05.sav',to.data.frame = T)

head(data3_5)
##   no  hb
## 1  1 112
## 2  2 137
## 3  3 129
## 4  4 126
## 5  5  88
## 6  6  90

进行单样本t检验,与140进行比较:

t.test(data3_5$hb,mu=140,alternative = 'two.sided') # 双侧检验
## 
##  One Sample t-test
## 
## data:  data3_5$hb
## t = -2.1367, df = 35, p-value = 0.03969
## alternative hypothesis: true mean is not equal to 140
## 95 percent confidence interval:
##  122.1238 139.5428
## sample estimates:
## mean of x 
##  130.8333

结果显示t=-2.1367,自由度df=35,p=0.03969<0.05,拒绝H0,接受H1,差异具有统计学意义,可认为从事铅作业的男性工人平均血红蛋白含量低于正常成年男性。

1.2 配对样本t检验

孙振球《医学统计学》第4版例3-6

为比较两种方法对乳酸饮料中脂肪含量测定结果是否不同,随机抽取了10份乳酸饮料制品,分别用脂肪酸水解法和哥特里-罗紫法测定。问两法测定结果是否不同?

library(foreign)
data3_6 <- read.spss('datasets/例03-06.sav',to.data.frame = T)

head(data3_6)
##   no    x1    x2
## 1  1 0.840 0.580
## 2  2 0.591 0.509
## 3  3 0.674 0.500
## 4  4 0.632 0.316
## 5  5 0.687 0.337
## 6  6 0.978 0.517

数据一共3列10行,第1列是样本编号,第2列和第3列是要比较的值。

进行配对样本t检验:

t.test(data3_6$x1,data3_6$x2,paired = T,var.equal = T)
## 
##  Paired t-test
## 
## data:  data3_6$x1 and data3_6$x2
## t = 7.926, df = 9, p-value = 2.384e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.1946542 0.3501458
## sample estimates:
## mean difference 
##          0.2724

结果显示t=7.926,自由度df=9,p<0.001,拒绝H0,接受H1,可认为两种方法对脂肪含量的测定结果不同,哥特里-罗紫法测定结果较高。

如果数据格式是两列数据,R语言也支持formula的形式。

警告

t.test在R4.4.0及以后不再支持在公式形式中使用paired参数,所以如果你的R版本在4.4.0及以后的,在公式形式中使用paired参数会报错:cannot use 'paired' in formula method 非常变态的更新,

首先我们要把数据变成formula形式需要的格式,也就是长数据:

library(tidyr)

data3_6_long <- data3_6 %>% 
  pivot_longer(2:3,names_to = "group",values_to = "x")

# 数据一列是分组,另一列是数值,还有一列id没什么用
data3_6_long
## # A tibble: 20 × 3
##       no group     x
##    <dbl> <chr> <dbl>
##  1     1 x1    0.84 
##  2     1 x2    0.58 
##  3     2 x1    0.591
##  4     2 x2    0.509
##  5     3 x1    0.674
##  6     3 x2    0.5  
##  7     4 x1    0.632
##  8     4 x2    0.316
##  9     5 x1    0.687
## 10     5 x2    0.337
## 11     6 x1    0.978
## 12     6 x2    0.517
## 13     7 x1    0.75 
## 14     7 x2    0.454
## 15     8 x1    0.73 
## 16     8 x2    0.512
## 17     9 x1    1.2  
## 18     9 x2    0.997
## 19    10 x1    0.87 
## 20    10 x2    0.506

这时可以用formula形式(R4.4.0及以后的版本都不能这样用了):

t.test(x ~ group, data = data3_6_long, paired = T, var.equal = T)

    Two Sample t-test

data:  x by group
t = 3.2894, df = 18, p-value = 0.004076
alternative hypothesis: true difference in means between group x1 and group x2 is not equal to 0
95 percent confidence interval:
 0.09841834 0.44638166
sample estimates:
mean in group x1 mean in group x2 
          0.7952           0.5228 

结果也是一模一样的。

在R语言中,绝大多数统计检验函数都是支持多种形式的输入样式的,需要自己注意。

1.3 两样本t检验

孙振球《医学统计学》第4版例3-7

为研究阿卡波糖胶囊的降血糖效果,某医院用40名2型糖尿病患者进行同期随机对照试验。研究者将这些患者随机等分到试验组(用阿卡波糖胶囊)和对照组(用拜唐苹胶囊),分别测得试验开始前和试验8周时的空腹血糖,算得空腹血糖下降值如下,能否认为阿卡波糖胶囊与拜唐苹胶囊对空腹血糖的降糖效果不同?

library(foreign)
data3_7 <- read.spss('datasets/例03-07.sav',to.data.frame = T)
head(data3_7)
##   no    x        group
## 1  1 -0.7 阿卡波糖胶囊
## 2  2 -5.6 阿卡波糖胶囊
## 3  3  2.0 阿卡波糖胶囊
## 4  4  2.8 阿卡波糖胶囊
## 5  5  0.7 阿卡波糖胶囊
## 6  6  3.5 阿卡波糖胶囊

一共有3列40行,第1列是编号,第2列是血糖下降值,第3列是组别(阿卡波糖组和拜糖平组,每组20个人)。

进行两样本t检验:

t.test(x ~ group, data = data3_7, var.equal = T)
## 
##  Two Sample t-test
## 
## data:  x by group
## t = -0.64187, df = 38, p-value = 0.5248
## alternative hypothesis: true difference in means between group 阿卡波糖胶囊 and group 拜唐苹胶囊 is not equal to 0
## 95 percent confidence interval:
##  -2.326179  1.206179
## sample estimates:
## mean in group 阿卡波糖胶囊   mean in group 拜唐苹胶囊 
##                      2.065                      2.625

结果显示t=-0.64187,自由度df=38,p=0.5248。差异无统计学意义,尚不能认为阿卡波糖胶囊与拜唐苹胶囊对空腹血糖的降糖效果不同。

假如两样本方差不等,可以使用近似t检验,或者使用Mann-Whitney-U检验。孙振球《医学统计学》第4版给出了3种近似t检验的方法(例3-8),而t.test只能给出Welch校正的结果。

孙振球《医学统计学》第4版例3-8只给出了均值和方差,没给原始数据,所以没法复现其结果。下面使用例3-7的数据演示下如何实现近似t检验:

t.test(x ~ group, data = data3_7, var.equal = F)#这里指定方差不等即可
## 
##  Welch Two Sample t-test
## 
## data:  x by group
## t = -0.64187, df = 36.086, p-value = 0.525
## alternative hypothesis: true difference in means between group 阿卡波糖胶囊 and group 拜唐苹胶囊 is not equal to 0
## 95 percent confidence interval:
##  -2.32926  1.20926
## sample estimates:
## mean in group 阿卡波糖胶囊   mean in group 拜唐苹胶囊 
##                      2.065                      2.625

1.4 正态性检验和两样本方差比较的F检验

1.4.1 正态性检验

图示法暂不演示,直接展示计算法。

孙振球《医学统计学》第4版例3-9:对例3-1中计算机模拟抽样所得的100个样本均数进行正态性检验。

# 使用foreign包读取SPSS数据
library(foreign)
data3_9 <- read.spss('datasets/例03-01.sav',to.data.frame = T)
head(data3_9)
##   no   mean   sd
## 1  1 167.41 2.74
## 2  2 165.56 6.57
## 3  3 168.20 5.36
## 4  4 166.67 4.81
## 5  5 164.89 5.41
## 6  6 166.36 4.50

进行Shapiro-Wilk正态性检验:

shapiro.test(data3_9$mean)
## 
##  Shapiro-Wilk normality test
## 
## data:  data3_9$mean
## W = 0.99409, p-value = 0.9444

这个结果和SPSS的结果是一样的。

计算偏度和峰度:

library(moments) # psych

skewness(data3_9$mean)
## [1] 0.1423707
kurtosis(data3_9$mean)
## [1] 3.045566

# 偏度和峰度检验
agostino.test(data3_9$mean)
## 
##  D'Agostino skewness test
## 
## data:  data3_9$mean
## skew = 0.14237, z = 0.61614, p-value = 0.5378
## alternative hypothesis: data have a skewness
anscombe.test(data3_9$mean)
## 
##  Anscombe-Glynn kurtosis test
## 
## data:  data3_9$mean
## kurt = 3.04557, z = 0.41992, p-value = 0.6745
## alternative hypothesis: kurtosis is not equal to 3

1.4.2 方差齐性检验

孙振球《医学统计学》第4版中的例3-11没给数据(使用了例3-8的数据,但是只给了均值和方差,没给原始数据),所以我们用例3-7的数据演示F检验。

library(foreign)
data3_7 <- read.spss('datasets/例03-07.sav',to.data.frame = T)
head(data3_7)
##   no    x        group
## 1  1 -0.7 阿卡波糖胶囊
## 2  2 -5.6 阿卡波糖胶囊
## 3  3  2.0 阿卡波糖胶囊
## 4  4  2.8 阿卡波糖胶囊
## 5  5  0.7 阿卡波糖胶囊
## 6  6  3.5 阿卡波糖胶囊

进行F检验:

var.test(x ~ group, data = data3_7)
## 
##  F test to compare two variances
## 
## data:  x by group
## F = 1.5984, num df = 19, denom df = 19, p-value = 0.3153
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.6326505 4.0381795
## sample estimates:
## ratio of variances 
##           1.598361

F=1.6,p大于0.05,不拒绝原假设,可以认为方差齐。