思考问题的熊
很高兴和你相遇
这里正在记录我的所思所学
订阅免费邮件通讯接收最新内容
首页 归档 想法 工具 通讯 播客 简历 主页

统计学基础与 R-5

2017-10-01, 2533 words, 11 min read

写在前面

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

单双样本均值分析

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

假设检验一般步骤

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

  1. 确定假设

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

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

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

真实情况H0H_0 真实情况H1H_1
接受H0H_0 H0H_0为真且被接受 H1H_1是真，H0H_0被接受
拒绝H0H_0 H0H_0为真但是被拒绝 H1H_1是真，H0H_0被拒绝

其中H0H_0为真但是被拒绝的错误概率称为** I 型错误**，通常用α\alpha表示，也称为显著性水平，H1H_1为真但是接受H0H_0的错误概率称为** II 型错误**，常用β\beta表示。

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

  1. 选择检验统计量

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

  1. 决定拒绝域

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

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

  1. 求出 p 值

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

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

  2. 做出决策

单样本 t-test

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

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

统计量计算公式为

t=X‾−μ0s/nt=\frac{\overline{X}-\mu_0}{s/\sqrt{n}}

其中X‾\overline X表示样本均值，n 表示样本量，s 是样本标准差（总体方差未知），μ0\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=X‾1−X‾2S2n1+S2n2t = \frac{\overline{X}_1 - \overline{X}_2}{\sqrt{ \frac{S^2}{n_1} + \frac{S^2}{n_2} }}

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

s2=(n1−1)s12+(n2−1)s22n1+n2−2s^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) 进行许可。

如果你对这篇文章感兴趣，欢迎通过邮箱或者微信订阅我的 「熊言熊语」会员通讯，我将第一时间与你分享肿瘤生物医药领域最新行业研究进展和我的所思所学所想，点此链接即可进行免费订阅。


· 分享链接 https://kaopubear.top/blog/2017-10-01-RandStatistics5/

R 数据科学 统计
Powered byGridea&Theme OnePlus
△