R语言实现bootstrap和jackknife检验方法

Python020

R语言实现bootstrap和jackknife检验方法,第1张

写在最前面:

首先需要说一下,本文的bootstrap和jackknife都算是蒙特卡罗方法(Monte Carlo method)的一种。应用广泛的的MCMC链(马尔可夫链蒙特卡洛方法Markov chain Monte Carlo)也是蒙特卡罗与马尔可夫链的结合。简单来说,蒙特卡罗方法就是从已知样本的分布中随机抽取新的样本集进行评估,然后放回,再次抽取的方法。根据具体方法的不同,抽取样本集的手段也不同。

bootstrap抽样方法将观测到的样本视为一个有限的总体,是唯一的信息来源,从中有放回的随机抽样来评估总体特征,以及对抽样总体进行推断统计。bootstrap 也分参数bootstrap和非参数bootstrap,前者的分布已完全知道。但在生信领域一般没有这种情况。所以下面讨论的是非参数bootstrap。

直接上例子:

假设现在有bootstrap包中的law数据集如下,

现在我们要计算LSAT成绩(美国法学入学考试)和GPA之间的相关系数。但因为样本量太少了,所以我们使用bootstrap重复抽样评估其标准误。

200次循环抽样后,计算得se.R标准误为0.1474629

得到如下的图:

1e6次循环抽样后,计算得se.R标准误为0.1333802

得到如下的图:

如果用bootstrap包的bootstrap函数会快一些:

bootstrap函数的用法: bootstrap(抽取样本范围,重复次数,进行bootstrap的函数,bootstrap的数据集)

偏差定义为bootstrap结果(多个数值)与原数据统计结果(单个数值)的均值:

得到bias大约为0.001817608,比较小

换一个包,boot包

这里用了三种方法计算置信区间:basic、正态和百分数。样本相关系数分布接近正态,则正态置信区间接近百分数区间。此外还有“Better Bootstrap Confivendence Interval” 更好的bootstrap置信区间,称为BCa区间,使用偏差和偏度对百分数置信区间进行矫正。设置type="bca"即可。

简单的说,bootstrap是从原有真实样本中有放回地抽取n个。jacknife就是每次都抽取n-1个样本,也就是每次只剔除一个原样本。

同样地,如果以bootstrap包中的law数据进行演示:

Jackknife计算的bias为-0.006473623。 这里jackknife的偏差公式相比于bootstrap有一个(n-1)系数,推导就不写了。

标准误se为0.1425186,与bootstrap得出的比较接近。

当统计量不太平滑的时候,Jacknife有很大误差。比如说对中位数进行统计,其变化很大。在进行Jacknife之后最好再跑一次bootstrap,看看是否相差很大。

居然还能这么嵌套着玩,针对每次bootstrap形成的数列向量计算jackknife的标准差,这样可以看出bootstrap若干次取样之间的差异。

算出来分别为0.1344824和0.08545141。后者较小,表面bootstrap取样之间的variance较小。

简单来说就是一种数据分割检验的方法,将数据分割为K份,称为"K-fold"交叉检验,每次第i个子集作为测试集来评估模型,其余的用来构建模型。Admixture使用的就是这个原理。Jackknife也属于Cross Validation的应用之一。

现在我创建一个这样的alignment:

这棵树长这样,符合遗传距离:

进行bootstrap:

phylogeny的bootstrap是对每一个节点都进行bootstrap取样并建树,比如说在9号节点,查看其bootstrap子集建的树符合系统发育关系((human2,human4,human3)(human8,human1,human6,human7,human5))的百分比(不管内部怎么样,先看这个节点)。发现Node1支持率是100(1000次都符合)。而后移到下一个节点,并且只看节点内部的分支支持率是多少。

其实原理都比较简单,计算bootstrap也会有专门的软件。

参考资料:

1)中科大张伟平教授课件

2) https://ecomorph.wordpress.com/2014/10/09/phylogenetic-trees-in-r-4/

Meta分析是一种对同一主题下的多个独立实验(研究)进行综合的统计分析方法。它萌芽于本世纪初[2];1976年由美国教育学家定义为Meta分析,并揭开了它在教育学、心理学及医学中的应用的新篇章[3]。Meta分析在这些学科的应用中取得了极大的成功,发展出了多种分析方法。Mann称其为医学方法学研究中的一次革命,且羽翼渐丰[4]。

直到90年代,此方法才被生态学家发现,虽然目前它在生态学中的应用实例还很少,但已引起了生态学界的高度重视。Gurevitch(1993)出版了第一部生态学中的Meta分析专著[5],并与人合作于1997年发行了MetaWin软件包。

在我国,彭少麟(1988)首次将此方法引入我国生态学界[6],并利用此方法进行生态学分析[7]。

Meta分析目前主要应用于对照实验的综合研究中,目的为判断实验中的处理会对实验对象产生正或负效应;效应是大还是小;同一主题下不同独立实验的结果是否一致,变异程度有多大等问题。

但Meta分析决不仅仅是一个数学分析过程,它本身也是一项研究,需要认真设计。主要步骤如下所述。

提出所要解决的问题并制定搜集、选择文献的标准。搜集文献,这是一项非常繁重且关键的工作。为了能搜集到全面的文献,通过各种途径来最大可能地收集已发表的和未发表文献(包括正式期刊中的论文、会议论文、摘要以及各种私人交换资料等)。

标定各研究的特点,并对其进行分类。根据研究背景特点的不同将所有研究分为几个级别(class),以作比较。

定量测度研究特点。为了避免分析时对质量不等的研究给予相同的结合标准,导致分析结果的不准确,分析家们提出了定性Meta分析,即制定标准,对研究特点进行打分评估;综合研究结果并结合研究特点来分析结果。也有人称这一步为定量Meta分析,以相对于定性Meta分析。

研究特征分析(敏感性分析),分析研究的基本特征(研究对象、研究环境等的特征)和方法学特征对效应值之间的协变关系。

目前已有发展出多种定量Meta分析方法。但它们的基本思想是一致的,那就是先提出假设,构造一个结合统计量,然后计算各研究的结合统计量,并用其在定性Meta分析中所得分数去权重它的结合统计量;计算各级别研究中的加权平均结合统计量(在平均过程中,要根据其各结合统计量的方差进行权重);做各级别研究间统计量的异质性检验。

定量Meta分析方法的不同主要在于结合统计量和统计假设的不同。

2 MetaWin软件的特点

MetaWin是一个主要为生态学工作者设计的定量Meta分析软件,其主要特点如下所述。

2.1 提供了两种假设模型

这两种假设模型为固定效应模型和混合效应模型,具体计算过程见文献[6]。两者的区别主要在于前者假设所综合的研究共享一个真实效应大小,实际测量的效应大小不同是由于随机取样所导致,而后者却假设研究间具有不同的真实效应大小,即所测效应大小的不同是由两部分组成,真实效应的不同,随机取样造成误差。后者更切合实际,区间估计较保守,更受Meta分析家们欢迎。

2.2 提供多种可选择的结合统计量

在生态学领域内的Meta分析中最常用的结合统计量为Hedges’d效应值:d =(Xe -Xc)/(SJ)(其中,Xe、Xc分别为实验组和对照组的测量平均值,S为两组共同标准差,J为小样本较正值),MetaWin还提供了反应比(response ratio):ln(Xe /Xc)(Xe、Xc的意义同上)这是从医学Meta中新引进的一种结合统计量;此外,MetaWin还为对Meta分析较为熟悉的分析者提供了更多的选择机会,如相关系数(correlation coefficient)等。

2.3 提供了两种数据输入方式

对有经验的分析者可直接输入效应值、样本方差等所需数据,其格式称效应数据格式。这种数据输入法的好处在于分析者可根据所收集的文献的实际情况来自己构造结合统计量,也即MetaWin为分析者提供了较大的自由。在文献数据满足前两种结合统计量计算情况下,分析者可以输入原文献中的统计数据,如平均值、样本方差、样本大小来进行计算,比较方便,称原始数据格式。

2.4 提供了一项非参数检验——重取样检验

上述参数模型检验是在假设所有研究中的实验组和对照组观测值均遵循正态分布情况下进行的;许多Meta分析方法基于大样本近似原理,即当实验组和对照组样本大小不小于10时,效应值才趋于正态分布。但如果样本太小,实验组和对照组样本大小太悬殊或效应值太大时,大样本近似原理就变得不准确了[8]。但事实上,许多生态学观察值却违背了上述情况[9]。此外,只有当上述假设被满足时,用于检验研究间效应异质性的Q值才有近似的X 2分布[7]。重取样检验法是取代传统参数和非参数检验的一种好方法。

重取样检验是一种计算机加强(computer intensive)非参数检验方法[10]。MetaWin中提供了随机化检验法(randomization test)和自助法或靴襻法(bootstrap)。前者常被用来决定一个统计量的显著性水平,后者则用于给出统计量的置信区间。

MetaWin中用自助法来计算所有研究总效应值和每一级别加权平均效应值的置信区间,对于样本含量为i的每一级别,我们均以放回式取样选取i个研究并计算其加权效应值,然后重复上述取样方法多次,按大小顺序将效应值排列起来,在两端取2.5%处的值做为5%至信区间的上下限,置信区间包括零在内的级别被认为没有显著不等于零。但当样本含量太小时,会出现区间估计过低,此时,可用偏差较正法[9]。

MetaWin中用随机检验来判断级别间效应大小的差异是否显著。首先用原始数算出QB,然后将j个级别里的所有研究混在一起,再随机将它们分成j个级别,级别含量仍与原来相同,算出QB值,重复此过程多次,得出一个QB值的分布,QB的显著性水平为随机QB值大于等于实际QB值数占重复随机取样数的百分比。

3 MetaWin软件的使用方法

3.1 MetaWin软件构成

运行MetaWin,只需一台装有Windows95、Windows3.1或WindowsNT的IBM兼容机,其中共包括8个文件。(1)MetaWin.exe:在Windows95和WindowsNT下的可执行文件。(2)MetaWin.hlp:在Windows95和WindowsNT下的帮助文件。(3)MetaWin.cnt:在Windows95和WindowsNT下的帮助文件的内容。(4)MetaW16.exe:在Windows3.1下的可执行文件。(5)Meta16.hlp:在Windows3.1下的帮助文件。(6)Raw.dta:以原文献统计数据输入数据的格式示范文件。(7)Effect.dta:以效应大小输入数据的格式示范文件。(8)Gur-hed.dta:作者的示范数据格式文件。

Windows95和Windows3.1版本的不同之处主要在于研究特征类型量、每一特征类型中级别数、每一数据文件中所含研究量及非参数检验中的重复数的最大值的不同,Windows95比Windows3.1范围更广。

3.2 MetaWin软件使用方法

3.2.1 数据输入

打开MetaWin文件,下拉file菜单,点击edit data file,进入数据输入状态,可直接在弹出的窗口中输入数据,也可从file菜单中点击load a file上载已有的文件。原始数据输入格式如下:

sex tree state +/- Nc Ne Xc Xe Sc Se Label

m oak pa + 7 7 78.14 79.71 40.650 40.650 study1

m maple ny + 7 7 18.86 26.00 9.170 9.170 study2

f maple ny - 6 6 -1.80 -2.10 0.490 0.490 study3

其中,第一行为标题行,前3项为级别分类标准,事实上,Windows95版本可允许10个分类标准,Windows3.1为5个;+/-为方向符,如果你所期待的效应值为正值时(即实验中的处理会对实验对象产生正效应),在按所期待趋势应该出现正效应值的研究项中加+,负效应值的研究中加-,它必须紧跟级别组,否则程序运行时不能识别数据文件中共有几项划分级别标准;Ne、Nc分别为实验组和对照组的样本含量;Xe、Xc分别为实验组和对照组的测量平均值;Se、Sc分别为实验组和对照组的标准差;Label为各研究的标记。标题行下面的每一行为一个研究的效应数据。效应数据输入格式为:

sex tree state +/- Nc Ne effect var Label

m oak pa + 7 7 78.0.036 0.286 study1

m maple ny + 7 7 0.565 0.347 study2

f maple ny - 6 6 1.533 0.517 study3

其中,effect一列为效应值;var为效应方差;其它同原始数据格式。

3.2.2 数据分析

在打开MetaWin文件的同时,会自动弹出一个Meta-Analysis窗口。在此窗口上部Type of Input一项中点击raw或 effect(确定数据为原始还是效应格式),原来灰色的Data file就会加亮,点击,从弹出的‘打开’窗口中选定并打开要分析的数据文件。

在Meta-Analysis窗口中部选择固定效应或混合效应模型,结合统计量,也可增加重复检验;窗口下部gourp by中可选择划分级别的标准,并在Refine Analysis中可以在不改变数据文件的情况下去掉一些级别或研究来纯化分析。

所有这些选项选择好后,即可点start键进行运算。运算结束后,会自动弹出一个Meta-Analysis output窗口,显示分析结果。

3.2.3 结果显示

在分析结果中,可看到分析时间,数据来源路径,以及名为Parametric methods和 Meta-Analysis results for groups的两个表。前者为所有研究的效应值表,每一行代表一个研究,包括其名称、所属级别、小样本校正值(J)、对照和实验两组的共有标准差(spool)、效应值(d)、95%的置信区间(95%CI)、各研究的方差(V)、权重(W)。

一般先假设所有研究享有共同的d值进行分析,此时在第二个结果表中可看到所有研究的总平均效应值(d++)、95%的置信区间、同质性(Qwi),自由度(df)、X 2检验的p值。如果其级别内异质性经X 2检验显著,则说明假设不正确,此时按一定的标准将所有研究划分为几个级别,再进行分析。此时的第二个结果表中会显示各级别内所有研究的加权平均效应值(di+)、95%的置信区间、同质性(Qwi),自由度(df)、X 2检验的p值以及级别间同质性(Qb)、级别内总同质性(Qw),总的同质性(Qtotal)。如果级别间同质性(Qb)经检验后显著,则说明级别间差异显著;如果某一级别内同质性(Qwi)经检验显著,说明这一级别内各研究的效应值差异较大,应该进一步划分此级别,再分析,直到Qwi经检验不显著。

如我们在做捕食关系的Meta分析中发现,捕食者导致被捕食者种群数量降低,d++=-0.3855(固定效应模型)d++=-0.4589(混合效应模型);但不同标准划分的级别的效应大小有差异,捕食效应随地带性而变化,热带效应值最大,亚热带、温带、寒带也有效应,其中亚热带最小;按所在生态系统划分级别时,陆生生态系统级别为中效应,淡水生态系统为小效应[10]。

MetaWin是一个操作简单且功能较全的Meta分析软件,以Windows作支持,用户通过界面与机器直接对话,分析过程简单易学,结果输出明了。遗撼的是MetaWin中没有考虑定性Meta分析所得出的研究质量评估值,所以利用原始数据直接输入法不能对效应值进行研究质量权重。同时值得注意的是MetaWin只提供了定量Meta分析方法,而Meta分析本是一项研究,数量分析前需要认真设计,分析后也需对结果进行研究特征分析。一个好的Meta分析不仅要选择好的定量分析方法,而且分析前设计和分析后的特征分析都非常重要,因为统计的目的是为解决问题提供科学依据。

另外,团IDC网上有许多产品团购,便宜有口碑

R语言基本数据分析

本文基于R语言进行基本数据统计分析,包括基本作图,线性拟合,逻辑回归,bootstrap采样和Anova方差分析的实现及应用。

不多说,直接上代码,代码中有注释。

1. 基本作图(盒图,qq图)

#basic plot

boxplot(x)

qqplot(x,y)

2. 线性拟合

#linear regression

n = 10

x1 = rnorm(n)#variable 1

x2 = rnorm(n)#variable 2

y = rnorm(n)*3

mod = lm(y~x1+x2)

model.matrix(mod) #erect the matrix of mod

plot(mod) #plot residual and fitted of the solution, Q-Q plot and cook distance

summary(mod) #get the statistic information of the model

hatvalues(mod) #very important, for abnormal sample detection

3. 逻辑回归

#logistic regression

x <- c(0, 1, 2, 3, 4, 5)

y <- c(0, 9, 21, 47, 60, 63) # the number of successes

n <- 70 #the number of trails

z <- n - y #the number of failures

b <- cbind(y, z) # column bind

fitx <- glm(b~x,family = binomial) # a particular type of generalized linear model

print(fitx)

plot(x,y,xlim=c(0,5),ylim=c(0,65)) #plot the points (x,y)

beta0 <- fitx$coef[1]

beta1 <- fitx$coef[2]

fn <- function(x) n*exp(beta0+beta1*x)/(1+exp(beta0+beta1*x))

par(new=T)

curve(fn,0,5,ylim=c(0,60)) # plot the logistic regression curve

3. Bootstrap采样

# bootstrap

# Application: 随机采样,获取最大eigenvalue占所有eigenvalue和之比,并画图显示distribution

dat = matrix(rnorm(100*5),100,5)

no.samples = 200 #sample 200 times

# theta = matrix(rep(0,no.samples*5),no.samples,5)

theta =rep(0,no.samples*5)

for (i in 1:no.samples)

{

j = sample(1:100,100,replace = TRUE)#get 100 samples each time

datrnd = dat[j,]#select one row each time

lambda = princomp(datrnd)$sdev^2#get eigenvalues

# theta[i,] = lambda

theta[i] = lambda[1]/sum(lambda)#plot the ratio of the biggest eigenvalue

}

# hist(theta[1,]) #plot the histogram of the first(biggest) eigenvalue

hist(theta)#plot the percentage distribution of the biggest eigenvalue

sd(theta)#standard deviation of theta

#上面注释掉的语句,可以全部去掉注释并将其下一条语句注释掉,完成画最大eigenvalue分布的功能

4. ANOVA方差分析

#Application:判断一个自变量是否有影响 (假设我们喂3种维他命给3头猪,想看喂维他命有没有用)

#

y = rnorm(9)#weight gain by pig(Yij, i is the treatment, j is the pig_id), 一般由用户自行输入

#y = matrix(c(1,10,1,2,10,2,1,9,1),9,1)

Treatment <- factor(c(1,2,3,1,2,3,1,2,3)) #each {1,2,3} is a group

mod = lm(y~Treatment) #linear regression

print(anova(mod))

#解释:Df(degree of freedom)

#Sum Sq: deviance (within groups, and residuals) 总偏差和

# Mean Sq: variance (within groups, and residuals) 平均方差和

# compare the contribution given by Treatment and Residual

#F value: Mean Sq(Treatment)/Mean Sq(Residuals)

#Pr(>F): p-value. 根据p-value决定是否接受Hypothesis H0:多个样本总体均数相等(检验水准为0.05)

qqnorm(mod$residual) #plot the residual approximated by mod

#如果qqnorm of residual像一条直线,说明residual符合正态分布,也就是说Treatment带来的contribution很小,也就是说Treatment无法带来收益(多喂维他命少喂维他命没区别)

如下面两图分别是

(左)用 y = matrix(c(1,10,1,2,10,2,1,9,1),9,1)和

(右)y = rnorm(9)

的结果。可见如果给定猪吃维他命2后体重特别突出的数据结果后,qq图种residual不在是一条直线,换句话说residual不再符合正态分布,i.e., 维他命对猪的体重有影响。