我们之前讲过相关性分析,生物信息学常见的相关性分析是pearson相关和spearman相关。
(生物信息学)R语言与统计学入门(六)—— Pearson和Spearman相关性分析_Lijingxian教你学生信的博客
但是相关性分析也有它的的缺点。
相关分析只能得出两个变量之间是否相关, 但却不能回答在两个变量之间
存在相关关系时, 它们之间是如何联系的, 即无法找出刻画它们之间因果关系
的函数关系. 回归分析就可以解决这一问题, 先从一元线性回归讲起。
设变量x和y之间存在一定的相关关系, 回归分析方法即找出Y 的值是如
何随X的值的变化而变化的规律, 我们称Y 为因变量(或响应变量), X为自变
量(或解释变量), 现通过例子说明如何来确定Y 与X之间的关系。
下面我们来看一看,比如两个基因的表达,他们之间的相关性可以用pearson和spearman分析,但是他们之间的线性关系,则得用一元线性回归分析。
下面先看看代码,我们构建两个基因的表达量read count:
GeneA <- c(318, 910, 200, 409, 425, 502, 314, 1210, 1022, 1225)
GeneB <- c(524, 1019, 638, 815, 913, 928, 605, 1516, 1219, 1624)
GeneA <- log2(GeneA)
GeneB <- log2(GeneB)
这里取log2是希望两个数据能够正态化。
plot(GeneA, GeneB)
可以看到两个数据的关系如图所示。
在R中, 由函数lm( )可以非常方便地求出回归方程, 函数confint( )可求 出参数的置信区间. 与回归分析有关的函数还有summary( ), anova( ) 和predict( )等. 函数lm( )的的调用格式为:
lm(formula, data, subset, weights, na.action,
method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
singular.ok = TRUE, contrasts = NULL, offset, ...)
说明: formula是显示回归模型, data是数据框, subset是样本观察的子 集, weights是用于拟合的加权向量, na.action显示数据是否包含缺失值, method是指出用于拟合的方法, model, x, y, qr是逻辑表达, 如果是TRUE, 应 返回其值. 除了第一个选项formula是必选项, 其它都是可选项.
函数confint( )的调用格式为:
confint(object, parm, level = 0.95, ...)
说明: object是指回归模型, parm要求指出所求区间估计的参数, 默认值为所 有的回归参数, level是指置信水平.
下面我们来算回归方程:
GeneA <- c(318, 910, 200, 409, 425, 502, 314, 1210, 1022, 1225)
GeneB <- c(524, 1019, 638, 815, 913, 928, 605, 1516, 1219, 1624)
GeneA <- log2(GeneA)
GeneB <- log2(GeneB)
lm.reg<-lm(GeneA~1+GeneB)
summary(lm.reg)
这里的lm(GeneA~1+GeneB) 也可以写成lm(GeneA~GeneB)
结果如下所示:
> summary(lm.reg)
Call:
lm(formula = GeneA ~ 1 + GeneB)
Residuals:
Min 1Q Median 3Q Max
-0.6350 -0.1507 -0.0466 0.2339 0.5012
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.1998 2.2703 -2.731 0.025814 *
GeneB 1.5539 0.2303 6.747 0.000146 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3794 on 8 degrees of freedom
Multiple R-squared: 0.8505, Adjusted R-squared: 0.8318
F-statistic: 45.52 on 1 and 8 DF, p-value: 0.0001456
结果解读:
-6.1998是结局,1.5539是系数,因此公式为
GeneA = -6.1998+1.5539*GeneB
下面计算置信区间:
confint(lm.reg, level=0.95)
> confint(lm.reg, level=0.95)
2.5 % 97.5 %
(Intercept) -11.435212 -0.964442
GeneB 1.022792 2.085079
得到了回归方程, 还可以对误差项独立同正态分布的假设进行检验. 在R中只需再执行一个plot命令:
op<-par(mfrow=c(2, 2))
plot(lm.reg)
par(op)
上面的命令plot(lm.reg)实际上使用了四次plot(x, y), 产生四个图形, 它们分别为: 1) Residual vs fitted为拟合值yˆ对残差的图形, 可以看出, 数据点都基本均匀 地分布在直线y 0的两侧, 无明显趋势; 2) Normal QQ-plot图中数据点分布趋于一条直线, 说明残差是服从正态分 布的; 3) Scale-Location 图显示了标准化残差(standardized residuals)的平方根 的分布情况. 最高点为残差最大值点; 4) Cook距离(Cook’s distance)图显示了对回归的影响点。
了解了一元线性回归以后,我们便可以将拟合的模型绘制在散点图中: