R语言实战(12)——重抽样与自助法

引言:在第7到9章中,我们学习了假设检验和总体参数的置信区间估计等统计方法,但这些都是基于观测数据抽样是来自正态分布或者其他性质较好的分布数据下的本章,我们将学习当以上假设不能满足的条件下,比如数据抽样于未知或混合分布、样本量过小等等情况,我们常用的两种应用广泛的依据随机化思想的统计方法:置换检验和自助法。

12.1 置换检验

也称随机化检验或重随机化检验。置换方法和参数方法都计算了相同的t统计量。但置换方法并不是将统计量与理论分布进行比较,而是将其与置换观测后获得的经验分布进行比较,根据统计量值的极端性判断是否有足够理由拒绝零假设。这种逻辑可以延伸到大部分经典统计检验和线性模型上来。

常用的置换检验包有两个:coin包和lmPerm包。coin包对于独立性问题提供了一个非常全面的置换检验框架,而lmPerm包则专门用来做方差分析和回归分析的置换检验。

要安装coin包,可以使用:

install.packages(c("coin")

要安装lmPerm,可以使用:

(1) http://cran.r-project.org/src/contrib/Archive/lmPerm下载文件lmPerm_1.1-2.tar.gz,并将其保存在硬盘上。
(2) MS Windows用户:从http://cran.r-project.org/bin/windows/RTools安装RToolsMacLinux用户可以跳过此步骤。
(3) R内部 执行函数:

install.packages(file.choose(), repos=NULL, type="source")

当弹出对话框,找到并选择lmPerm_1.1-2.tar.gz文件。把这个包安在我们的计算机上.

#注意:置换检验都是使用伪随机数来从所有可能的排列组合中进行抽样(当作近似检验时)。因此,每次检验的结果都有所不同。为了重现结果,设定随机数种子为 1234(即set.seed(1234) )。

12.2 用 coin 包做置换检验


1 、注意 在coin包中, 类别型变量和序数变量必须分别转化为因子和有序因子。另外, 数据要以数据框形式存储。
2、coin包提供了一个置换检验的一般性框架, 可以分析一组变量相对于其他任意变量, 是否与第二组变量( 可根据一个区组变量分层) 相互独立.。包中提供的函数(见表12-1

12-1 相对于传统检验,提供可选置换检验的 coin函数

检 验  coin函数
两样本和K样本置换检验 oneway_test(y ~ A)
含一个分层(区组)因子的两样本和K样本置换检验 oneway_test(y ~  A | C)
Wilcoxon-Mann-Whitney秩和检验 wilcox_test(y ~ A)
Kruskal-Wallis检验 kruskal_ test(y ~ A)
Pearson卡方检验 chisq_test(A ~ B)
Cochran-Mantel-Haenszel检验 cmh_test(A ~ B  | C)
线性关联检验 lbl_test(D ~ E)
Spearman检验 spearman_test(y  ~ x)
Friedman检验 friedman_test(y ~ A | C)
Wilcoxon符号秩检验 wilcoxsign_test(y1  ~ y2)

注:在coin函数中, yx是数值变量, AB是分类因子, C是类别型区组变量, DE是有序因子, y1y2 是相匹配的数值变量。

表中的函数形式:

f u n c t i o n _ n a m e( f o r m u l a, d a t a, distribution= )

  • formula描述的是要检验变量间的关系,示例可参见表12-1;

  •  data是一个数据框;

  • distribution指定经验分布在零假设条件下的形式,可能值有exact、 asymptotic和approximate。

 

12.2.1 独立两样本和 K 样本检验
首先,根据表12-2的虚拟数据,我们对独立样本t检验和单因素精确检验进行比较,

#虚拟数据中的t检验与单因素置换检验> library(coin)> score <- c(40, 57, 45, 55, 58, 57, 64, 55, 62, 65)> treatment <- factor(c(rep("A",5), rep("B",5)))> mydata <- data.frame(treatment, score)> t.test(score~treatment, data=mydata, var.equal=TRUE)#传统t检验表明存在显著性差异( p<0.05), Two Sample t-testdata: score by treatmentt = -2.3, df = 8, p-value = 0.04705alternative hypothesis: true difference in means is not equal to 095 percent confidence interval: -19.04 -0.16sample estimates:mean in group A mean in group B 51 61> oneway_test(score~treatment, data=mydata, distribution="exact")#精确检验却表明差异并不显著( p>0.072) Exact 2-Sample Permutation Testdata: score by treatment (A, B)Z = -1.9, p-value = 0.07143alternative hypothesis: true mu is not equal to 0

结论:传统t检验表明存在显著性差异( p<0.05),而精确检验却表明差异并不显著( p>0.072)。由于只有10个观测,我更倾向于相信置换检验的结果,在做出最后结论之前,还要多收集些数据。

7章中,我们用wilcox.test() 函数检验了美国南部监禁概率与非南部间的差异。现使用Wilcoxon秩和检验,如下:

> library(MASS)> UScrime <- transform(UScrime, So = factor(So))#数值变量So被转化为因子> wilcox_test(Prob ~ So, data=UScrime, distribution="exact") Exact Wilcoxon Mann-Whitney Rank Sum Testdata: Prob by So (0, 1)Z = -3.7, p-value = 8.488e-05alternative hypothesis: true mu is not equal to 0

结论:此处结果与第7wilcox.test() 计算结果一样,这是因为wilcox.test() 默认计算的也是精确分布。

 

探究K样本检验问题:
在第9章,对于50个患者的样本,我们使用了单因素方差分析来评价五种药物疗法对降低胆固醇的效果。下面代码对其做了近似的K样本置换检验:

> library(multcomp)> set.seed(1234)#设定随机数种子可让结果重现> oneway_test(response~trt, data=cholesterol,distribution=approximate(B=9999)) Approximative K-Sample Permutation Testdata: response by trt (1time, 2times, 4times, drugD, drugE)maxT = 4.7623, p-value < 2.2e-16

结论:结果表明各组间病人的响应值显著不同。

1 2.2.2 列联表中的独立性
置换检验判断两类别型变量的独立性, 可用chisq_test() cmh_test() 函数**
若是大于3个类别型变量, 只可以用cmh_test() 函数**
若变量都是有序型, 可使用lbl_test() 函数来检验是否存在线性趋势**

 

例如,第7章中,我们用卡方检验评价了关节炎的治疗( treatment)与效果( improvement)间的关系,若想实施卡方检验的置换版本,可用如下代码:

> library(coin)> library(vcd)> Arthritis <- transform(Arthritis, Improved=as.factor(as.numeric(Improved)))> set.seed(1234)> chisq_test(Treatment~Improved, data=Arthritis, distribution=approximate(B=9999)) Approximative Pearson's Chi-Squared Testdata: Treatment by Improved (1, 2, 3)chi-squared = 13.055, p-value = 0.0018

12.2.3 数值变量间的独立性
spearman_test() 函数提供了两数值变量的独立性置换检验。

 

7章中,我们检验了美国文盲率与谋杀率间的相关性。如下代码可进行相关性的置换检验:

> states <- as.data.frame(state.x77)> set.seed(1234)> spearman_test(Illiteracy~Murder, data=states,distribution=approximate(B=9999)) Approximative Spearman Correlation Testdata: Illiteracy by MurderZ = 4.7065, p-value < 2.2e-16alternative hypothesis: true mu is not equal to 0

12.2.4两样本和 K 样本相关性检验
当处于不同组的观测已经被分配得当,或者使用了重复测量时,对 于 两 配 对 组 的 置 换 检 验 , 可 使 用 wilcoxsign_test() 函 数 ;多 于 两 组 时 , 使 用friedman_test() 函数。

 

7章中,我们比较了城市男性中1424年龄段( U1)与3539年龄段( U2 )间的失业率差异。此次我们可使用精确Wilcoxon符号秩检验来判断两个年龄段间的失业率是否相等(state是匹配变量):

> library(coin)> library(MASS)> wilcoxsign_test(U1~U2, data=UScrime, distribution="exact")#结果表明两者的失业率是不同的 Exact Wilcoxon-Signed-Rank Testdata: y by x (neg, pos) stratified by blockZ = 5.9691, p-value = 1.421e-14alternative hypothesis: true mu is not equal to 0

 

12.3 lmPerm 包的置换检验(做线性模型)

lmPerm包可做线性模型的置换检验。比如lmp()和aovp()lm()aov()函数的修改版,额外添加了perm=参数)能够进行置换检验,而非正态理论检验。

perm=ProbProb从所有可能的排列中不断抽样,直至估计的标准差在估计的p0.1之下。

perm=SPRSPR使用贯序概率比检验来判断何时停止抽样。

注意,若观测数大于10perm=”Exact” 自动默认转为perm=”Prob” ,因为精确检验只适用于小样本问题。

 

1 2.3.1 简单回归和多项式回归

8章中,我们使用线性回归研究了15名女性的身高和体重间的关系。用lmp() 代替lm()可获得如下代码中的置换检验的结果。

简单线性回归的置换检验

> library(lmPerm)> set.seed(1234)> fit <- lmp(weight~height, data=women, perm="Prob")[1] "Settings: unique SS : numeric variables centered"> summary(fit)Call:lmp(formula = weight ~ height, data = women, perm = "Prob")Residuals: Min     1Q   Median   3Q  Max-1.733 -1.133 -0.383  0.742 3.117Coefficients: Estimate Iter Pr(Prob)Height 3.45 5000 <2e-16 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 1.5 on 13 degrees of freedomMultiple R-Squared: 0.991, Adjusted R-squared: 0.99F-statistic: 1.43e+03 on 1 and 13 DF, p-value: 1.09e-14

多项式回归的置换检验

> library(lmPerm)> set.seed(1234)> fit <- lmp(weight~height + I(height^2), data=women, perm="Prob")[1] "Settings: unique SS : numeric variables centered"> summary(fit)Call:lmp(formula = weight ~ height + I(height^2), data = women, perm = "Prob")Residuals: Min 1Q Median 3Q Max-0.5094 -0.2961 -0.0094 0.2862 0.5971Coefficients: Estimate Iter Pr(Prob)height -7.3483 5000 <2e-16 ***I(height^2) 0.0831 5000 <2e-16 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.38 on 12 degrees of freedomMultiple R-Squared: 0.999, Adjusted R-squared: 0.999F-statistic: 1.14e+04 on 2 and 12 DF, p-value: <2e-16

1 2.3.2 多元回归

在第8章,多元回归被用来通过美国50个州的人口数、文盲率、收入水平和结霜天数预测犯

罪率。将lmp() 函数应用到此问题,如下代码:

多元回归的置换检验

> library(lmPerm)> set.seed(1234)> states <- as.data.frame(state.x77)> fit <- lmp(Murder~Population + Illiteracy+Income+Frost, data=states, perm="Prob")[1] "Settings: unique SS : numeric variables centered"> summary(fit)Call:lmp(formula = Murder ~ Population + Illiteracy + Income + Frost, data = states, perm = "Prob")Residuals: Min 1Q Median 3Q Max-4.79597 -1.64946 -0.08112 1.48150 7.62104Coefficients: Estimate Iter Pr(Prob)Population 2.237e-04 51 1.0000Illiteracy 4.143e+00 5000 0.0004 ***Income 6.442e-05 51 1.0000Frost 5.813e-04 51 0.8627---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '. ' 0.1 ' ' 1Residual standard error: 2.535 on 45 degrees of freedomMultiple R-Squared: 0.567, Adjusted R-squared: 0.5285F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08##第8章,正态理论中Population和Illiteracy均显著( p<0.05)。而该置换检验中,Population不再显著

1 2.3.3 单因素方差分析和协方差分析

9章中任意一种方差分析设计都可进行置换检验。首先看9.1节中的单因素方差分析问题——各种疗法对降低胆固醇的影响。

单因素方差分析的置换检验

> library(lmPerm)> library(multcomp)> set.seed(1234)> fit <- aovp(response~trt, data=cholesterol, perm="Prob")[1] "Settings: unique SS "> anova(fit)Component 1 : Df R Sum Sq R Mean Sq Iter Pr(Prob)trt 4 1351.37 337.84 5000 < 2.2e-16 ***Residuals 45 468.75 10.42---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '. ' 0.1 ' ' 1#结果表明各疗法的效果不全相同

单因素协方差分析的置换检验

协方差分析的置换检验取自第9章的问题:当控制妊娠期时间相同时,观测四种药物剂量对
鼠崽体重的影响。

> library(lmPerm)> set.seed(1234)> fit <- aovp(weight ~ gesttime + dose, data=litter, perm="Prob")[1] "Settings: unique SS : numeric variables centered"> anova(fit)Component 1 : Df R Sum Sq R Mean Sq Iter Pr(Prob)gesttime 1 161.49 161.493 5000 0.0006 ***dose 3 137.12 45.708 5000 0.0392 *Residuals 69 1151.27 16.685---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#依据p值可知,当控制妊娠期时间相同时,四种药物剂量对鼠崽的体重影响不相同

1 2.3.4 双因素方差分析

以第9章中的维生素C对豚鼠牙齿生长的影响为例,该实验两个可操作的因子是剂量(三水平)和喂食方式(两水平)。10只豚鼠分别被分配到每种处理组合中,形成一个3×2的析因实验设计,置换检验结果如下:

> library(lmPerm)> set.seed(1234)> fit <- aovp(len~supp*dose, data=ToothGrowth, perm="Prob")[1] "Settings: unique SS : numeric variables centered"> anova(fit)Component 1 :          Df R Sum Sq R Mean Sq Iter Pr(Prob)supp      1 205.35 205.35 5000 < 2e-16 ***dose      1 2224.30 2224.30 5000 < 2e-16 ***supp:dose 1 88.92 88.92 2032 0.04724 *Residuals 56 933.63 16.67---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

12.4置换检验点
1、除coinlmPerm包外, R还提供了其他可做置换检验的包。

   perm包能实现coin包中的部分功能, 因此可作为coin包所得结果的验证。corrperm包提供了有重复测量的相关性的置换检验。logregperm包提供了Logistic回归的置换检验。另外一个非常重要的包是glmperm,它涵盖了广义线性模型的置换检验。对于广义线性模型,请参见第13章。
2、依靠基础的抽样分布理论知识,置换检验提供了另外一个十分强大的可选检验思路。对于上面描述的每一种置换检验,我们完全可以在做统计假设检验时不理会正态分布、 t分布、 F分布或者卡方分布。
3、置换检验真正发挥功用的地方是处理非正态数据(如分布偏倚很大)、存在离群点、样本很小或无法做参数检验等情况。不过,如果初始样本对感兴趣的总体情况代表性很差,即使是置换检验也无法提高推断效果。
4、置换检验主要用于生成检验零假设的p值,它有助于回答“效应是否存在”这样的问题。不过,置换方法对于获取置信区间和估计测量精度是比较困难的。

12.5 自助法 

自助法即是从初始样本重复随机替换抽样,生成一个或一系列待检验统计量的经验分布。无需假设一个特定的理论分布,便可生成统计量的置信区间,并能检验统计假设。即使潜在分布未知,或出现了离群点,或者样本量过小,再或者是没有可供选择的参数方法,自助法将是生成置信区间和做假设检验的很好的方法。

举个例子:你打算计算一个样本均值95%的置信区间。样本现有10个观测,均值为40,标准差为5。如果假设均值的样本分布为正态分布,那么(1–α/2)%的置信区间计算如下:

其中, t是自由度为n–1 t分布的1–α上界值。对于95%的置信区间,可得40–2.262(5/3.163)<
μ<40+2.262(5/3.162) 或者36.424<μ<43.577。以这种方式创建的95%置信区间将会包含真实的总体均值。倘若你假设均值的样本分布不是正态分布,该怎么办呢?可使用自助法。
(1) 从样本中随机选择10个观测,抽样后再放回。有些观测可能会被选择多次,有些可能一直都不会被选中。

(2) 计算并记录样本均值。
(3) 重复12一千次。
(4) 1000个样本均值从小到大排序。
(5) 找出样本均值2.5%97.5%的分位点。此时即初始位置和最末位置的第25个数,它们就限
定了95%的置信区间。

12.6 boot 包中的自助法

install.packages("boot")

自助法有三步骤

1写一个能返回待研究统计量值的函数。

2、为生成R中自助法所需的有效统计量重复数,使用boot() 函数
b  o o t o b j e c t <­ boot(data=, statistic=, R=, …) #R是自助抽样的次数

3、用boot.ci() 函数获取步骤(2)生成的统计量的置信区间
boot.ci(  b o o t o b j e c t, conf=, type= )

主要的自助法函数是boot() ,它的格式为:

b o o t o b j e c t <- boot(data=, statistic=, R=, ...)
参数描述见表12-2

12-2 boot () 函数的参数

参数 描述
data 向量、矩阵或者数据框
statistic 生成k个统计量以供自举的函数(k=1时对单个统计量进行自助抽样)。函数需包括indices参数,以便
boot() 函数用它从每个重复中选择实例(例子见下文)
R 自助抽样的次数
其他对生成待研究统计量有用的参数,可在函数中传输

boot() 函数调用统计量函数R次,每次都从整数1:nrow(data) 中生成一列有放回的随机指标,这些指标被统计量函数用来选择样本。统计量将根据所选样本进行计算,结果存储在bootobject中。bootobject结构的描述见表12-3

12-3 boot () 函数中返回对象所含的元素

元 素 描 述
t0 从原始数据得到的k个统计量的观测值
t 一个R×k矩阵,每行即k个统计量的自助重复值

你可以如bootobject$t0bootobject$t这样来获取这些元素。
一旦生成了自助样本,可通过print() plot() 来检查结果。如果结果看起来还算合理,
使用boot.ci() 函数获取统计量的置信区间。格式如下:

boot.ci( b o o t o b j e c t, conf=, type= )

参数见表12-4

12-4 boot .ci() 函数的参数

参 数 描 述
bootobject boot() 函数返回的对象
conf 预期的置信区间(默认:conf=0.95
type 返回的置信区间类型。可能值为normbasicstudpercbcaall(默认:type=”all”

type参数设定了获取置信区间的方法。 perc方法(分位数)展示的是样本均值, bca将根据偏差对区间做简单调整。 

 

接下来介绍如何对单个统计量和统计量向量使用自助法。

1 2.6.1 对单个统计量( R平方) 用自助法

       以1974Motor Trend杂志中的mtcars数据集为例,它包含32辆汽车的信息。假设你正使用多元回归根据车重和发动机排量来预测汽车每加仑行驶的英里数,除了标准的回归统计量,同时想获得95%R平方值的置信区间(预测变量对响应变量可解释的方差比),那么便可使用非参数的自助法来获取置信区间。

#先写一个获取R平方值的函数rsq <- function(formula, data, indices) { d <- data[indices,] fit <- lm(formula, data=d) return(summary(fit)$r.square)
#做大量的自助抽样(比如1000)library(boot)set.seed(1234)results <- boot(data=mtcars, statistic=rsq, R=1000, formula=mpg~wt+disp)  #可以输出boot的对象 > print(results)ORDINARY NONPARAMETRIC BOOTSTRAPCall:boot(data = mtcars, statistic = rsq, R = 1000, formula = mpg ~ wt + disp)Bootstrap Statistics : original bias std. errort1* 0.7809306 0.01333670 0.05068926
#绘制结果,图形见图12-1plot(results) 

12-1 自助法所得R平方值的均值

在图12-1中可以看到,自助的R平方值不呈正态分布。它的95%的置信区间可以通过如下代
码获得:

> boot.ci(results, type=c("perc""bca"))BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 1000 bootstrap replicatesCALL :boot.ci(boot.out = results, type = c("perc", "bca"))Intervals :Level Percentile BCa95% ( 0.6838, 0.8833 ) ( 0.6344, 0.8549 )Calculations and Intervals on Original ScaleSome BCa intervals may be unstable

结论:生成置信区间的不同方法将会导致获得不同的区间。本例中的依偏差调整区间方法与分位数方法稍有不同。两例中,由于0都在置信区间外,零假设H0R平方值=0都被拒绝。

1 2.6.2 多个统计量的自助法

继续该例,首先获取一个统计量向量——三个回归系数(截距项、车重和发动机排量) ——95%的置信区间。

#首先,创建一个返回回归系数向量的函数:bs <- function(formula, data, indices) {d <- data[indices,]fit <- lm(formula, data=d)return(coef(fit))}#使用该函数自助抽样1000次:library(boot)set.seed(1234)results <- boot(data=mtcars, statistic=bs, R=1000, formula=mpg~wt+disp)> print(results)ORDINARY NONPARAMETRIC BOOTSTRAPCall:boot(data = mtcars, statistic = bs, R = 1000, formula = mpg ~ wt + disp)Bootstrap Statistics : original bias std. errort1* 34.9606 0.137873 2.48576t2* -3.3508 -0.053904 1.17043t3* -0.0177 -0.000121 0.00879
#绘制车重结果,图形结果见图12-2,plot(results, index=2)#对多个统计量自助抽样时,添加一个索引参数,指明plot() 和boot.ci() 函数所分析bootobject$t的列。在本例中,索引 1 指截距项,索引 2指车重,索引 3指发动机排量。
#获得车重和发动机排量95%的置信区间> boot.ci(results, type="bca", index=2)#索引 2指车重BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 1000 bootstrap replicatesCALL :boot.ci(boot.out = results, type = "bca", index = 2)Intervals :Level BCa95% (-5.66, -1.19 )Calculations and Intervals on Original Scale
> boot.ci(results, type="bca", index=3)#索引 3指发动机排量BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 1000 bootstrap replicatesCALL :boot.ci(boot.out = results, type = "bca", index = 3)Intervals :Level BCa95% (-0.0331, 0.0010 )Calculations and Intervals on Original Scale

12-2 自助法所得车重回归系数的分布

关于自助法,也许有同学会问:

1、初始样本需要多大?

2、应该重复多少次?

第一个问题没固定的答案,只要样本能够较好地代表总体,初始样本大小为2030也可以了;第二个问题,一般情况下1000次重复在大部分情况下都可满足要求。但如果你愿意,也可以增加重复的次数。

 

小结:本章我们介绍了一系列基于随机化和重抽样的计算机密集型方法,在无需理论分布的知识便能够进行假设检验,获得置信区间。当数据来自未知分布,或者存在严重的离群点,或者样本量过小,又或者没有参数方法可以解决我们的假设问题时,当标准的数据假设不满足,或者你对于解决这些问题毫无头绪时,利用它们可以另辟蹊径。但是,置换检验和自助法并不是万能的,它们无法将烂数据转化为好数据。如果初始样本对于总体情况的代表性不佳,或者样本量过小而无法准确地反映总体情况,这些方法也是爱莫能助。在下一章,我们将学习一些数据模型,它们的变量服从已知但不必为正态的分布。

生物信息学

R语言实战(11)——中级绘图

2019-12-17 21:46:03

生物信息学

R语言实战(13)——广义线性模型

2019-12-17 21:49:00

0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
购物车
优惠劵
今日签到
有新私信 私信列表
搜索