[生物医学统计]R语言-t检验与Wilcoxon符号秩检验

2021-04-05 本文已影响0人戴钢盔的熊

这篇文章第一次使用markdown格式，如果排版感觉不太舒服，还望读者们提出。本文以Bernard Rosner的《Fundamentals of Biostatistics》的一道习题为例，整理总结了t-检验和wilcoxon符号秩检验的相关R代码。如有改进之处，欢迎交流！

高血压——红花油对SBD的效果评价报告

背景：
饮食中多不饱和脂肪酸对心血管病的一些危险因素有有利影响，其中多不饱和脂肪酸主要是亚油酸。为了检验饮食中补充亚油酸对血压的影响，选17例成人连续4周每日消耗23g红花油（亚油酸含量高）。在基线（摄入红花油以前）及1个月后测量血压，每次随访时取几次读数的平均值，数据见表1。

要检验亚油酸对血压的影响，用什么参数检验方法？
用1的方法进行检验，给出p值。
要检验亚油酸对血压的影响，用什么非参数方法？
用3的上述方法进行检验，给出p值。
5.比较2和4的结果，讨论你认为哪一种方法适合。

P1:

用配对双样本t检验

P2:

首先使用levene test进行方差齐性检验，p=0.9587>0.05，故接受备择假设：基线SBP和1个月后的SBP的方差没有显著差异；
然后使用Shapiro-Wilk normality test对差值进行差值正态性检验，p=0.975>0.05，说明差值数据服从正态整体分布；
利用配对双样本t检验来检验亚油酸对血压的影响：p-value = 0.0009303，说明亚油酸能显著降低SBP，双尾检验的结果为p=0.001861
t-test代码：

library(carData)
library(car)
library(ggpubr)
library(ggplot2)
library(magrittr)

SBP0 = c(119.67,100.00,123.56,109.89,96.22,133.33,115.78,126.39,122.78
         ,117.44,111.33,117.33,120.67,131.67,92.39,134.44,108.67
)
SBP1 = c(117.33,98.78,123.83,107.67,95.67,128.89,113.22,121.56,126.33
         ,110.39,107.00,108.44,117.00,126.89,93.06,126.67,108.67
)
#方差齐性检验
y=c(SBP0,SBP1)
p=rep(1,17)
group = as.factor(c(rep(1,17),rep(2,17)))
leveneTest(y,group)
#差值正态性检验
dSBP = c(2.34,1.22,-0.27,2.22,0.55,4.44,2.56,4.83,-3.55
         ,7.05,4.33,8.89,3.67,4.78,-0.67,7.77,0.00
)
shapiro.test(dSBP)
#t检验
t1=t.test(SBP0, SBP1, paired=T, alternative="two.sided", cond.lvel=0.95)    #双尾检验
t2=t.test(SBP0, SBP1, paired=T, alternative="greater", cond.lvel=0.95)        #单尾检验

下面绘制箱型图，

mydata = data.frame(group = rep(c("Conrtrast", "Trial"), each = 17),
                    SBP = c(SBP0,SBP1)
)
mydata$group=as.factor(mydata$group)
my_comparisons = list(c("Conrtrast", "Trial"))
pdf(file="t.test.boxplot.pdf", width=6, height=5)
ggboxplot(mydata
          ,x="group"
          ,y = "SBP"
          ,fill = "group"
          ,palette = "npg"
          ,linetype = "solid"
          ,bxp.errorbar=T
          ,bxp.errorbar.width=0.1
          ,add = "point"
          ,short.panel.labs = FALSE
)+ 
  stat_compare_means(comparisons=my_comparisons
                     ,aes(label = ..p.format..)
                     ,method = "t.test"
                     ,paired=T
                     ,label.x = 1.5)+
  theme(
    axis.text.x   = element_text(color = 'black', size = 16, angle = 0)
    ,axis.text.y   = element_text(color = 'black', size = 16, angle = 0)
    ,axis.title.x  = element_text(color = 'black', size = 16, angle = 0)
    ,axis.title.y  = element_text(color = 'black', size = 16, angle = 90)
    ,legend.title  = element_text(color = 'black', size  = 16)
    ,legend.text   = element_text(color = 'black', size   = 16)
    ,axis.line.y = element_line(color = 'black', linetype = 'solid')
    ,axis.line.x = element_line (color = 'black',linetype = 'solid') 
    ,panel.border = element_rect(linetype = 'solid', size = 1.2,fill = NA) # 图四周框起来
  )

如下图所示：

配对样本t检验结果

P3:

用Wilcoxon signed-rank检验。

P4:

单尾检验的p-value=0.002053
双尾检验的p-value= 0.004107
代码跟t-test类似，因为是非参数检验，所以少了对方差的检验：

library(ggpubr)
library(ggplot2)
library(magrittr)

SBP0 = c(119.67,100.00,123.56,109.89,96.22,133.33,115.78,126.39,122.78
         ,117.44,111.33,117.33,120.67,131.67,92.39,134.44,108.67
)
SBP1 = c(117.33,98.78,123.83,107.67,95.67,128.89,113.22,121.56,126.33
         ,110.39,107.00,108.44,117.00,126.89,93.06,126.67,108.67
)
#Wilcoxon符号秩检验
wiltest1=wilcox.test(SBP0, SBP1, paired=T, alternative="two.sided", cond.lvel=0.95)
wiltest2=wilcox.test(SBP0, SBP1, paired=T, alternative="greater", cond.lvel=0.95)

下面绘制箱型图：

mydata = data.frame(group = rep(c("Conrtrast", "Trial"), each = 17),
                    SBP = c(SBP0,SBP1)
)
mydata$group=as.factor(mydata$group)
my_comparisons = list(c("Conrtrast", "Trial"))
pdf(file="wilcox.test.boxplot.pdf", width=6, height=5)
ggboxplot(mydata
          ,x="group"
          ,y = "SBP"
          ,fill = "group"
          ,palette = "npg"
          ,linetype = "solid"
          ,bxp.errorbar=T
          ,bxp.errorbar.width=0.1
          ,add = "point"
          ,short.panel.labs = FALSE
)+ 
  stat_compare_means(comparisons=my_comparisons         #设置比较组
                     ,aes(label = ..p.format..)
                     ,method = "wilcox.test"         #默认方法
                     ,paired=T
                     ,label.x = 1.5)+
  theme(
    axis.text.x   = element_text(color = 'black', size = 16, angle = 0)
    ,axis.text.y   = element_text(color = 'black', size = 16, angle = 0)
    ,axis.title.x  = element_text(color = 'black', size = 16, angle = 0)
    ,axis.title.y  = element_text(color = 'black', size = 16, angle = 90)
    ,legend.title  = element_text(color = 'black', size  = 16)
    ,legend.text   = element_text(color = 'black', size   = 16)
    ,axis.line.y = element_line(color = 'black', linetype = 'solid')
    ,axis.line.x = element_line (color = 'black',linetype = 'solid') 
    ,panel.border = element_rect(linetype = 'solid', size = 1.2,fill = NA) # 图四周框起来
  )

如下图所示：

Wilcoxon signed-rank检验结果

P5:

t检验的p值更小，所以t检验更好。(这个回答可能有误，还望指正。)

在医学写作中，Wilcoxon检验相对t检验使用得更多，因为数据的可能不能通过方差齐性检验，尤其是对照实验，也就是配对t检验；但一旦能够满足t检验的条件，参数检验将会被我们优先考虑使用，绘制箱型图时使用了ggpubr这个包，相对只用ggplot2来绘制方便了许多，并且能够通过参数设定不同期刊的绘图风格，另外这次在写代码时有意识地将逗号写在前面，养成良好的编写习惯，方便后期修改，之后会整理总结利用循环对features进行分组wilcoxon符号秩检验。

[生物医学统计]R语言-t检验与Wilcoxon符号秩检验

高血压——红花油对SBD的效果评价报告

P1:

P2:

P3:

P4:

P5:

猜你喜欢

热点阅读