统计学基础与R-5

统计学基础与R-5

写在前面

入门生物信息,所有人都绕不开统计基础知识和相关实现方式。本章我们将简要介绍统计学相关基础知识以及如何使用R语言进行简单地计算和分析。

单双样本均值分析

根据数据组数的不同,均值比较可以分为单样本、双样本和多样本。本节首先介绍单样本和双样本。

假设检验一般步骤

谈具体的假设检验之前首先介绍假设检验的一般步骤。

  1. 确定假设

原假设(null hypothesis, $H_0$) 指需要检验的假设,只要当我们有足够的证据时才能否定。在某种意义上与原假设相反的的假设被称为备择假设(alternative hypothesis,$H_1$)。

假设检验是无法给出绝对证明的,我们只能在假定原假设为真的情况下通过假设检验来判断结果是否可信。如果结果极不可能发生,则拒绝原假设,进而接受备择假设。

根据接受假设和真实情况之间的关系,会有四种可能发生。

真实情况$H_0$ 真实情况$H_1$
接受$H_0$ $H_0$为真且被接受 $H_1$是真,$H_0$被接受
拒绝$H_0$ $H_0$为真但是被拒绝 $H_1$是真,$H_0$被拒绝

其中$H_0$为真但是被拒绝的错误概率称为I型错误,通常用$\alpha$表示,也称为显著性水平,$H_1$为真但是接受$H_0$的错误概率称为II型错误,常用$\beta$表示。

检验的功效(power)=1-$\beta$。虽然我们希望$\alpha$和$\beta$都可以尽可能小,但这两者往往是矛盾的,通常我们会先固定$\alpha$在某个水平,再找合适的检验使$\beta$尽可能小。

  1. 选择检验统计量

检验统计量(test statistic)是用于对假设进行检验的统计量,进行检验的过程建立在这个统计量之上。

  1. 决定拒绝域

拒绝域是可以不接受原假设的一组数值,其分界点被称为临界点。为了求出拒绝域首先要确定假设检验的显著性水平,也就是当样本结果的不可能发生程度多大时就可以拒绝原假设。

单尾检验指检验的拒绝域落在可能数据集的一侧(左侧或者右侧),双尾检验的拒绝域则分布在数据集的两侧,对于检验水平是$\alpha$的双尾检测而言,两侧分别是$\alpha/2$。如果备择假设的表示是“不等于”,则应该使用双尾检验。

  1. 求出p值

p值取决于检验统计量和拒绝域。表示得到更加极端结果的概率。

  1. 判断样本结果是否在拒绝域内

  2. 做出决策

单样本t-test

说完假设检验的基本流程之后,下面介绍和均值比较相关的假设检验,首先是单样本t-test

在数据符合正态分布的前提下使用单样本t-test来比较一组样本的均值和已知(理论/总体)均值,所谓的已知均值可能来自于之前的实验数据或者理论值。根据研究问题(原假设)的不同又分为双尾(不等)和单尾检验(大于或者小于)。

统计量计算公式为

$$t=\frac{\overline{X}-\mu_0}{s/\sqrt{n}}$$

其中$\overline X$表示样本均值,n表示样本量,s是样本标准差(总体方差未知),$\mu_0$表示理论值。通过统计量t和自由度,我们可以计算出对应的p值。

使用R进行单样本t-test,这里借用R自带数据集PlantGrowth的10个对照组数据。

weight <- PlantGrowth$weight[PlantGrowth$group=="ctrl"]
shapiro.test(weight)
summary(weight)

# Shapiro-Wilk normality test
#
#data: weight
#W = 0.95668, p-value = 0.7475
#
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 4.170 4.550 5.155 5.032 5.293 6.110

可以发现均值大概是5左右,且数据符合正态分布。下面检验和理论值7相比的情况,原假设是等于7,备择假设是不等于7,所以采用双尾检验。

t.test(weight,mu = 7)

# One Sample t-test
#
#data: weight
#t = -10.673, df = 9, p-value = 2.075e-06
#alternative hypothesis: true mean is not equal to 7
#95 percent confidence interval:
# 4.614882 5.449118
#sample estimates:
#mean of x
# 5.032

通过上述结果看出统计量t为-10.673,自由度是9,p值是2.075e-06<0.05,所以拒绝原假设。

单样本Wilcoxon 符号秩检验

如果样本数据没有通过正态分布检验就要采用单样本wilcoxon符号秩检验进行计算。使用该检验需要满足的条件是样本值均匀地分布在均值两侧。

其R中的函数为wilcox.test()

wilcox.test(weight,mu=7)

通过结果可知拒绝原假设。

独立双样本t-test

两个独立样本是指两个比较的样本之间没有关联,不互相影响。使用双样本独立t-test的前提是两个样本符合正态分布且方差相等

独立双样本t-test其统计量计算公式为

$$t = \frac{\overline{X}_1 - \overline{X}_2}{\sqrt{ \frac{S^2}{n_1} + \frac{S^2}{n_2} }}$$

$\overline{X}$分别表示两个样本的均值,n表示样本量。$s^2$表示两个独立样本的方差合并估计,计算公式为

$$s^2 = \frac{(n_1-1)s^2_1+(n_2-1)s^2_2}{n_1+n_2-2}$$

依旧以PlantGrowth数据进行举例,PlantGrowth共有三组数据,在这里取两组进行分析。

# 查看各组基本信息
library(dplyr)
tmp1<-group_by(PlantGrowth,group) %>% summarise(count=n(),mean=mean(weight),sd=sd(weight))
tmp1

# A tibble: 3 x 4
# group count mean sd
# <fctr> <int> <dbl> <dbl>
#1 ctrl 10 5.032 0.5830914
#2 trt1 10 4.661 0.7936757
#3 trt2 10 5.526 0.4425733

进行独立双样本t-test前,首先验证数据是否符合正态分布(Shapiro-Wilk normality test)以及两个样本方差是否相等(F-test)

关于方差比较

  • 两个方差比较通常使用F-test,在R中的函数为var.test(),进行F-test 检验前一定要首先进行正态分布检验,该统计量对于正态分布这个前提条件非常敏感;

  • 多个方差进行比较还可以使用Bartlett’s test,该方法同样要求满足正态分布前提条件;

  • 如果不能确定数据是否符合正态分布,可以使用Levene’s test 进行检验,其对数据正态分布的要求并非十分严格;

  • Fligner-Killeen test 则是一种非参数检验方法,数据可以不满足正态分布;

  • 关于方差比较部分内容,在本章不再展开。

shapiro.test(PlantGrowth$weight[PlantGrowth$group=="ctrl"])
shapiro.test(PlantGrowth$weight[PlantGrowth$group=="trt1"])
var.test(PlantGrowth$weight[PlantGrowth$group=="ctrl"],
PlantGrowth$weight[PlantGrowth$group=="trt1"])

# Shapiro-Wilk normality test
#
#data: PlantGrowth$weight[PlantGrowth$group == "ctrl"]
#W = 0.95668, p-value = 0.7475
#
#
# Shapiro-Wilk normality test
#
#data: PlantGrowth$weight[PlantGrowth$group == "trt1"]
#W = 0.93041, p-value = 0.4519
#
#
# F test to compare two variances
#
#data: PlantGrowth$weight[PlantGrowth$group == "ctrl"] and PlantGrowth$weight[PlantGrowth$group == "trt1"]
#F = 0.53974, num df = 9, denom df = 9, p-value = 0.3719
#alternative hypothesis: true ratio of variances is not equal to 1
#95 percent confidence interval:
# 0.1340645 2.1730025
#sample estimates:
#ratio of variances
# 0.5397431

满足前提条件,继续方差相等的独立双样本t-test

t.test(PlantGrowth$weight[PlantGrowth$group=="ctrl"],
PlantGrowth$weight[PlantGrowth$group=="trt1"],
var.equal = T )

# Two Sample t-test
#
#data: PlantGrowth$weight[PlantGrowth$group == "ctrl"] and PlantGrowth$weight[PlantGrowth$group == "trt1"]
#t = 1.1913, df = 18, p-value = 0.249
#alternative hypothesis: true difference in means is not equal to 0
#95 percent confidence interval:
# -0.2833003 1.0253003
#sample estimates:
#mean of x mean of y
# 5.032 4.661

独立双样本Wilcoxon test

当两个样本不满足正态分布时,使用Wilcoxon秩和检验进行非参数检验。

wilcox.test(PlantGrowth$weight[PlantGrowth$group=="ctrl"],
PlantGrowth$weight[PlantGrowth$group=="trt1"],
exact = F )

# Wilcoxon rank sum test with continuity correction
#
# data: PlantGrowth$weight[PlantGrowth$group == "ctrl"] and PlantGrowth$weight[PlantGrowth$group == "trt1"]
# W = 67.5, p-value = 0.1986
# alternative hypothesis: true location shift is not equal to 0

非独立双样本t-test

在很多试验中需要比较的两组样本往往有关系的,比如一组病人服药前后的变化,每个人以自己为对照。所谓非独立双样本,就是在第一组样本中的每个数据点都和第二组样本中唯一的数据点对应。

非独立样本和独立样本相比,可以显著降低样本量,提高统计的power。其统计量计算公式和单样本t检验一致。

t.test(PlantGrowth$weight[PlantGrowth$group=="ctrl"],
PlantGrowth$weight[PlantGrowth$group=="trt1"],
paired = T )

# Paired t-test
#
# data: PlantGrowth$weight[PlantGrowth$group == "ctrl"] and PlantGrowth$weight[PlantGrowth$group == "trt1"]
# t = 0.99384, df = 9, p-value = 0.3463
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# -0.4734609 1.2154609
# sample estimates:
# mean of the differences
# 0.371

非独立双样本Wilcoxon test

和单样本类似,如果不符合正态分布则使用非参数检验。

wilcox.test(PlantGrowth$weight[PlantGrowth$group=="ctrl"],
PlantGrowth$weight[PlantGrowth$group=="trt1"],
exact = F , paired = T)

# Wilcoxon signed rank test with continuity correction
#
# data: PlantGrowth$weight[PlantGrowth$group == "ctrl"] and PlantGrowth$weight[PlantGrowth$group == "trt1"]
# V = 37, p-value = 0.359
# alternative hypothesis: true location shift is not equal to 0

本文作者:思考问题的熊

版权声明:本博客所有文章除特别声明外,均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×