使用R语言进行协整关系检验

Python019

使用R语言进行协整关系检验,第1张

使用R语言进行协整关系检验

协整检验是为了检验非平稳序列的因果关系,协整检验是解决伪回归为问题的重要方法。首先回归伪回归例子:

伪回归Spurious regression伪回归方程的拟合优度、显著性水平等指标都很好,但是其残差序列是一个非平稳序列,拟合一个伪回归:

#调用相关R包

library(lmtest)

library(tseries)

#模拟序列

set.seed(123456)

e1=rnorm(500)

e2=rnorm(500)

trd=1:500

y1=0.8*trd+cumsum(e1)

y2=0.6*trd+cumsum(e2)

sr.reg=lm(y1~y2)

#提取回归残差

error=residuals(sr.reg)

#作残差散点图

plot(error, main="Plot of error")

#对残差进行单位根检验

adf.test(error)

## Dickey-Fuller = -2.548, Lag order = 7, p-value = 0.3463

## alternative hypothesis: stationary

#伪回归结果,相关参数都显著

summary(sr.reg)

## Residuals:

## Min 1Q Median 3Q Max

## -30.654 -11.526 0.359 11.142 31.006

## Coefficients:

## Estimate Std. Error t value Pr(>|t|)

## (Intercept) -29.32697 1.36716 -21.4 <2e-16 ***

## y2 1.44079 0.00752 191.6 <2e-16 ***

## Residual standard error: 13.7 on 498 degrees of freedom

## Multiple R-squared: 0.987, Adjusted R-squared: 0.987

## F-statistic: 3.67e+04 on 1 and 498 DF, p-value: <2e-16

dwtest(sr.reg)

## DW = 0.0172, p-value <2.2e-16

恩格尔-格兰杰检验Engle-Granger第一步:建立两变量(y1,y2)的回归方程,第二部:对该回归方程的残差(resid)进行单位根检验其中,原假设两变量不存在协整关系,备择假设是两变量存在协整关系。利用最小二乘法对回归方程进行估计,从回归方程中提取残差进行检验。

set.seed(123456)

e1=rnorm(100)

e2=rnorm(100)

y1=cumsum(e1)

y2=0.6*y1+e2

# (伪)回归模型

lr.reg=lm(y2~y1)

error=residuals(lr.reg)

adf.test(error)

## Dickey-Fuller = -3.988, Lag order = 4, p-value = 0.01262

## alternative hypothesis: stationary

error.lagged=error[-c(99,100)]

#建立误差修正模型ECM.REG

dy1=diff(y1)

dy2=diff(y2)

diff.dat=data.frame(embed(cbind(dy1, dy2),2))#emed表示嵌入时间序列dy1,dy2到diff.dat

colnames(diff.dat)=c("dy1","dy2","dy1.1","dy2.1")

ecm.reg=lm(dy2~error.lagged+dy1.1+dy2.1, data=diff.dat)

summary(ecm.reg)

## Residuals:

## Min 1Q Median 3Q Max

## -2.959 -0.544 0.137 0.711 2.307

## Coefficients:

## Estimate Std. Error t value Pr(>|t|)

## (Intercept) 0.0034 0.1036 0.03 0.97

## error.lagged -0.9688 0.1585 -6.11 2.2e-08 ***

## dy1.1 0.8086 0.1120 7.22 1.4e-10 ***

## dy2.1 -1.0589 0.1084 -9.77 5.6e-16 ***

## Residual standard error: 1.03 on 94 degrees of freedom

## Multiple R-squared: 0.546, Adjusted R-squared: 0.532

## F-statistic: 37.7 on 3 and 94 DF, p-value: 4.24e-16

par(mfrow=c(2,2))

plot(ecm.reg)

Johansen-Juselius(JJ)协整检验法,该方法是一种用向量自回归(VAR)模型进行检验的方法,适用于对多重一阶单整I(1)序列进行协整检验。JJ检验有两种:特征值轨迹检验和最大特征值检验。我们可以调用urca包中的ca.jo命令完成这两种检验。其语法:

ca.jo(x, type = c("eigen", "trace"), ecdet = c("none", "const", "trend"), K = 2,spec=c("longrun", "transitory"), season = NULL, dumvar = NULL)

其中:x为矩阵形式数据框;type用来设置检验方法;ecdet用于设置模型形式:none表示不带截距项,const表示带常数截距项,trend表示带趋势项。K表示自回归序列的滞后阶数;spec表示向量误差修正模型反映的序列间的长期或短期关系;season表示季节效应;dumvar表示哑变量设置。

set.seed(12345)e1=rnorm(250,0,0.5)e2=rnorm(250,0,0.5)e3=rnorm(250,0,0.5)#模拟没有移动平均的向量自回归序列;u1.ar1=arima.sim(model=list(ar=0.75), innov=e1, n=250)u2.ar1=arima.sim(model=list(ar=0.3), innov=e2, n=250)y3=cumsum(e3)y1=0.8*y3+u1.ar1y2=-0.3*y3+u2.ar1#合并y1,y2,y3构成进行JJ检验的数据库;y.mat=data.frame(y1, y2, y3)#调用urca包中cajo命令对向量自回归序列进行JJ协整检验vecm=ca.jo(y.mat)jo.results=summary(vecm)#cajorls命令可以得到限制协整阶数的向量误差修正模型的最小二乘法回归结果vecm.r2=cajorls(vecm, r=2)vecm.r2## Call:lm(formula = substitute(form1), data = data.mat)## Coefficients:## y1.d y2.d y3.d## ect1 -0.33129 0.06461 0.01268## ect2 0.09447 -0.70938 -0.00916## constant 0.16837 -0.02702 0.02526## y1.dl1-0.22768 0.02701 0.06816## y2.dl1 0.14445 -0.71561 0.04049## y3.dl1 0.12347 -0.29083 -0.07525## $beta## ect1 ect2## y1.l2 1.000e+00 0.0000## y2.l2 -3.402e-18 1.0000## y3.l2 -7.329e-01 0.2952

数据准备

t检验,亦称student t检验(Student's t test),主要用于样本含量较小(例如n <30),总体标准差σ未知的正态分布。t检验是用t分布理论来推论差异发生的概率,从而比较两个平均数的差异是否显著。

t检验的适用条件为样本分布符合正态分布。

R中检验正态分布的函数:

shapiro.test()

结果p值要是小于0.05,样本分布是非正态分布,如果大于0.05,样本分布是正态分布。

t检验可分为单总体检验和双总体检验,以及配对样本检验。

单总体t检验是检验一个样本平均数与一个已知的总体平均数的差异是否显著。

个人理解的应用实例:已知一个玉米品种的产量是8000 kg/ha,在一个田间试验中测定这个玉米品种的产量,单样本t检验要做的就是检验田间试验测定的产量与已知产量是否相等。

单样本t检验的假设:

H0:样本均值与已知的总体均值相等。

H1:样本均值与已知的总体均值不相等。

t统计量的计算:

m:样本平均值;

:已知总体的均值;

S:样本标准差,自由度df=n-1。

n:样本量。

单样本t检验R调用函数:

t.test(x, mu, alternative = "two.sided")

x:数据向量;

mu:理论平均值。默认为0,可根据自己统计计算需求更改;

alternative:备择假设。允许值为“two.sided”(默认),也可以根据需要设置为“greater”或“less”之一。

结果解释:p值小于0.05,结论是v1的平均值与理论值1.5有显著差异。

检验两个样本平均数与其各自所代表的总体的差异是否显著。

个人理解的应用实例:检验两个玉米品种产量是否存在差异。

t.test(y ~ x, data)

其中的y是一个数值型变量,x是一个二分变量。

t.test(y1, y2)

其中的y1和y2为数值型向量(即各组的结果变量)。可选参数data的取值为一个包含了这些变量的矩阵或数据框。

t检验默认假定方差不相等,并使用Welsh的修正自由度。你可以添加一个参数var.equal=TRUE以假定方差相等,并使用合并方差估计。默认的备择假设是双侧的(即均值不相等,但大小的方向不确定)。你可以添加一个参数alternative="less"或alternative="greater"来进行有方向的检验。

结果解读:得到结果中P值小于0.05,说明要拒绝原假设(两品种v1值无差异),接受备择假设,即两品种v1值差异显著。

非独立样本的t检验假定组间的差异呈正态分布。

个人理解的应用实例:一个玉米品种接受两个施氮处理,两个施氮处理下玉米的产量是否存在差异。

t.test(y1, y2, paired=TRUE)

其中的y1和y2为两个非独立组的数值向量。

结果解读:不同氮素水平的比较显示p值小于0.05,说明v1值在两个氮水平间差异显著;而两个年份下v1值无显著差异。

如果想在多于两个的组之间进行比较,应该怎么做?如果能够假设数据是从正态总体中独立抽样而得的,那么你可以使用方差分析(ANOVA)。ANOVA是一套覆盖了许多实验设计和准实验设计的综合方法。

参考资料:

在单细胞数据分析的过程中,寻找差异基因的过程需要用到对基因统计的假设检验(例如函数FindAllMarkers中的test.use参数),我们这里来深入了解一下假设检验的方法和应用环境。

秩和检验适用于广泛的统计学环境,秩和检验是检验总体分布位置是否相同,因而称为非参数检验(Nonparametric test)。秩和检验(rank sum test)是一类常用的非参数检验。秩和检验首先将数据按从小到大或等级从弱到强转换成秩(也就是顺序),然后求秩和并计算秩和统计量,最后做出统计推断。本文简单介绍秩和检验的原理并基于R语言进行秩和检验的操作。

假设我们从总体A和总体B中分别采样n_a和n_b个样本构成样本集合a和b。通过将样本集a和b中的所有样本按从小到大顺序转化为秩之后我们可以通过绘图的方式对转换的结果进行展示,在图中我们使用“•”代表来自样本集a,使用“o”代表数据来自样本集b。

如果总体A和总体B总体分布位置分布相同(H_0:A=B),那么转换的结果如下图所示:

首先是python(范例),借助于python模块scipy来实现。

其次是R的实现:(wilcox.test的函数)

这里可以发现,秩和检验仅仅和数据的总体分布有关,适用于一般的环境 ,在单细胞数据中寻找markergene 的过程中,大部分默认就是采用此方法,当然,这种检验只是一种很常规的检验,离我们真正的数据分析还很遥远。

t检验,亦称student t检验(Student's t test),主要用于样本含量较小(例如n <30),总体标准差σ未知的正态分布。 [1] t检验是用t分布理论来推论差异发生的 概率 ,从而比较两个平均数的差异是否显著。它与 f检验 、 卡方检验 并列。

这里我们需要注意一下:

(1)t检验对于大样本分布需要转换,而我们单细胞的数据分布属于大样本分布。

(2)数据分布为正态分布,单细胞数据是否为正态分布,在我的文章 单细胞数据分析之PCA再认识与ScaleData函数 做了详细的介绍,大家可以看一下。

t检验最常见的四个用途:

1、 单样本均值检验(One-sample t-test)

用于检验 总体方差未知、正态数据或近似正态的单样本的均值是否与已知的总体均值相等

2、两独立样本均值检验(Independent two-sample t-test)

用于检验两对独立的正态数据或近似正态的样本的均值是否相等,这里可根据总体方差是否相等分类讨论

3、配对样本均值检验(Dependent t-test for paired samples)

用于检验 一对配对样本的均值的差是否等于某一个值

4、回归系数的显著性检验(t-test for regression coefficient significance)

用于检验回归模型的解释变量对被解释变量是否有显著影响。

单样本T检验用于比较一组数据与一个特定数值之间的差异情况。

应用场景:

某个医生检测40名从事铅作业工人的血红蛋白含量,其均数为130.83g/L,标准差为25.74g/L,试分析从事铅作业的工人血红蛋白含量是否不同于正常成年人平均值140g/L?

我们来看一下这个结果,以p=0.2696,以p=0.05为界,说明没有统计意义。

两独立样本t检验的目的是利用来自两个总体的独立样本,推断两个总体的均值是否存在显著差异。

2、使用的前提条件

(1)两个样本应该是相互独立的;

(2)样本来自的两个总体应该服从正态分布。

显然单细胞使用的就是两独立样本均值检验。

用于分析配对定量数据之间的差异对比关系。与独立样本t检验相比,配对样本T检验要求样本是配对的。两个样本的样本量要相同;样本先后的顺序是一一对应的。

配对样本t检验用于样品的两个相关组之间的比较手段。在这种情况下,同一样本有两个值(即一对值)。

举个例子,在1个月内有20只小鼠接受了治疗X。我们想知道处理X是否会对小鼠的体重产生影响。

为了回答这个问题,在治疗之前和之后测量了20只小鼠的体重。通过测量相同小鼠体重的两次,我们得到了治疗前的20组值和治疗后的20组值。

在这种情况下,可以使用配对t检验比较治疗前后的平均体重。

似然比(likelihood ratio, LR) 是反映真实性的一种指标,属于同时反映灵敏度和特异度的复合指标。即有病者中得出某一筛检试验结果的概率与无病者得出这一概率的比值。该指标全面反映筛检试验的诊断价值,且非常稳定。似然比的计算只涉及到灵敏度与特异度,不受患病率的影响。因检验结果有阳性与阴性之分,似然比可相应地区分为阳性似然比(positive likelihood ratio, +LR)和阴性似然比(negative likelihood ratio, -LR)。阳性似然比是筛检结果的真阳性率与假阳性率之比。说明筛检试验正确判断阳性的可能性是错误判断阳性可能性的倍数。比值越大,试验结果阳性时为真阳性的概率越大。阴性似然比是筛检结果的假阴性率与真阴性率之比。表示错误判断阴性的可能性是正确判断阴性可能性的倍数。其比值越小,试验结果阴性时为真阴阳性的可能性越大。

似然比检验(likelihood ratio test, LRT) 是一种检验 参数能否反映真实约束 的方法(分布或模型的某参数 θ 等于 θ 0 是否为真实约束)。似然比检验的思想是:“如果参数约束是有效的,那么加上这样的约束不应该引起似然函数最大值的大幅度降低。也就是说似然比检验的实质是在 比较有约束条件下的似然函数最大值与无约束条件下似然函数最大值 。” 可以看出,似然比检验是一种通用的检验方法(比 t 检验、 Κ 2 检验等具有更广的适用范围)。

这个有点难,我们不展开讨论了,主要就是检验分群结果结束以后,基因的表达分布是否是受到约束的

Identifies 'markers' of gene expression using ROC analysis. For each gene, evaluates (using AUC) a classifier built on that gene alone, to classify between two groups of cells. An AUC value of 1 means that expression values for this gene alone can perfectly classify the two groupings (i.e. Each of the cells in cells.1 exhibit a higher level than each of the cells in cells.2). An AUC value of 0 also means there is perfect classification, but in the other direction. A value of 0.5 implies that the gene has no predictive power to classify the two groups. Returns a 'predictive power' (abs(AUC-0.5) * 2) ranked matrix of putative differentially expressed genes.

关于roc的讲解在我的文章里 深入理解R包AUcell对于分析单细胞的作用 详细提到过,大家可以看一下。