# 使用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 901 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)?
进行单样本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.6251.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 31.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.598361F=1.6,p大于0.05,不拒绝原假设,可以认为方差齐。