本文首发于 ”百味科研芝士“ 微信公众号，转载请注明：百味科研芝士，Focus科研人的百味需求
点击蓝字关注我们
这一系列里前面的三个部分都是用于比较组间差异的各种方法，详情点击：
R语言系列第四期：①R语言单样本双样本差异性检验
、R语言系列第四期：②R语言多组样本方差分析与KW检验
、R语言系列第四期：③R语言表格数据率的比较
在这个部分里，我们会为大家介绍如何使用R进行基础回归和相关分析，以及模型作图、置信区间的预估和展示。
A. 简单线性回归
我们使用数据集thuesen作为这一部分的例子，如下导入：
> library(ISwR)
> attach(thuesen)
我们使用函数lm( linear model,线性模型 )进行线性分析：
> lm(short.velocity~blood.glucose)
Call:
lm(formula = short.velocity ~ blood.glucose)
Coefficients:
(Intercept) blood.glucose
1.09781 0.02196
#Tips：函数lm()里面的参数是模型方程，波浪号(~)可以认为是“通过...来描述”。它在之前出现过几次，比如图形展示部分箱式图boxplot(),t检验，anova检验里等等。
#Tips：lm()函数的原始输出格式非常简单。你能看见的只有估计出来的截距α与斜率β。可以看出最优拟合直线为：short.velocity=1.098+0.0220*blood.glucose，但是没有给出其他任何像显著性检验之类的信息。
#Tips：其实，函数lm()可以处理比简单线性回归复杂很多的模型。除了一个解释变量与一个因变量之外，模型方程还能描述很多其他的情况。比如，要在y上通过x1，x2，x3进行多元线性回归分析（后文会介绍），可以通过y~x1+x2+x3来完成。
如果要获悉模型的其他信息，可以使用summary()：
> summary(lm(short.velocity~blood.glucose))
Call:
lm(formula = short.velocity ~ blood.glucose)
Residuals:
Min 1Q Median 3Q Max
-0.40141 -0.14760 -0.02202 0.03001 0.43490
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.09781 0.11748 9.345 6.26e-09 ***
blood.glucose 0.02196 0.01045 2.101 0.0479 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2167 on 21 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.1737, Adjusted R-squared: 0.1343
F-statistic: 4.414 on 1 and 21 DF, p-value: 0.0479
##Tips：首先，是
Call:
lm(formula = short.velocity ~ blood.glucose)
输出的开头本质上在重复一个函数的调用。当然如果你一开始就把模型赋值给了一个变量，那么调用summary()之后这个部分就是那个变量了。
Residuals:
Min 1Q Median 3Q Max
-0.40141 -0.14760 -0.02202 0.03001 0.43490
这部分简单地描述了残差的分布，可以帮助用户对分布的假设做快速检查。
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.09781 0.11748 9.345 6.26e-09 ***
blood.glucose 0.02196 0.01045 2.101 0.0479 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
这里我们再次见到了回归系数和截距，不过这次还伴有标准误，t检验和p值，以及最右侧的统计检验表示显著性水平的图像化标志。当然，如果你不喜欢这些星号，可以使用options(show.signif.stars=F)关闭它们。
Residual standard error: 0.2167 on 21 degrees of freedom
(1 observation deleted due to missingness)
残差波动情况。
Multiple R-squared: 0.1737, Adjusted R-squared: 0.1343
F-statistic: 4.414 on 1 and 21 DF, p-value: 0.0479
上式第一项是R2，在简单线性回归里可以被理解为Pearson相关系数的平方，另一个是修正后的R2；第二行是对假设回归系数是0进行的F检验，对整体模型的检验。实际上，F检验是t检验的平方：4.414=（2.101）2。
我们先把数据点和回归线画出来：
> plot(blood.glucose,short.velocity)
> abline(lm(short.velocity~blood.glucose))
#Tips：abline()函数根据截距和斜率画一条直线。它能够接受数值参数，比如abline(1.1,0.022)；不过更方便的是，它也能够从一个用lm拟合的线性回归中直接提取相关信息。
B. 预测和置信带
无论是否计算了置信带和预测带，我们都能够用函数predict析取出预测值，不加其他参数，它就只会输出回归值。我们为了方便用变量lm.velo代替模型：
> lm.velo<-lm(short.velocity~blood.glucose)
> predict(lm.velo)
1 2 3 4 5 6 7
1.433841 1.335010 1.275711 1.526084 1.255945 1.214216 1.302066
8 9 10 11 12 13 14
1.341599 1.262534 1.365758 1.244964 1.212020 1.515103 1.429449
15 17 18 19 20 21 22
1.244964 1.190057 1.324029 1.372346 1.451411 1.389916 1.205431
23 24
1.291085 1.306459
如果你加上了interval=”confidence”或者interval=”prediction”，那么你就能在预测值向量的基础上得到各类估计边界的值。同时也可以缩写这两个单词：
> predict(lm.velo,int="c")
fit lwr upr
1 1.433841 1.291371 1.576312
2 1.335010 1.240589 1.429431
3 1.275711 1.169536 1.381887
.....
22 1.205431 1.053805 1.357057
23 1.291085 1.191084 1.391086
24 1.306459 1.210592 1.402326
> predict(lm.velo,int="p")
fit lwr upr
1 1.433841 0.9612137 1.906469
2 1.335010 0.8745815 1.795439
3 1.275711 0.8127292 1.738693
...
22 1.205431 0.7299634 1.680899
23 1.291085 0.8294798 1.752690
24 1.306459 0.8457315 1.767186
Warning message:
In predict.lm(lm.velo, int = "p") :
用当前数据得到的预测结果对_未来_响应有用
#Tips：前一个是置信带，后一个是预测带。lwr和upr分别是下界和上界。Warning信息里提醒我们：这个预测边界不能用来考察我们做回归线所使用的已观测数据。
我们可以利用外部数据根据已有的模型做出它的预测情况图，图形的展示可以通过下面的组合函数来实现：
> pred.frame<-data.frame(blood.glucose=4:20)
> pp<-predict(lm.velo,int="p",newdata=pred.frame)
> pc<-predict(lm.velo,int="c",newdata=pred.frame)
> plot(blood.glucose,short.velocity,ylim=range(short.velocity,pp,na.rm=T))
> pred.gluc<-pred.frame$blood.glucose
> matlines(pred.gluc,pc,lty=c(1,2,2),col="black")
> matlines(pred.gluc,pp,lty=c(1,3,3),col="black")
#Tips：我们采用由4到20这几个数据作为新的x来做出预测情况的图。Predict()函数里newdata=的参数就是调用新数据的参数；plot()函数里的ylim参数使用range()函数来保证图形全部在范围内；matlines()函数里的lty是设置线型。
A. 皮尔逊相关系数
相关系数的计算可以使用cor()函数，但是如果对thuesen中的两个向量也进行这样简单的操作，就会发生下面状况：
> cor(blood.glucose,short.velocity)
[1] NA
R中所有的基本统计函数都要求输入的参数没有缺失值，或者你明确指定如何处理缺失值。对于函数mean（），var（），sd（）以及类似的单向量函数，你可以传递na.rm=T这个参数告诉它们在计算之前应该移除缺失值。对于cor()函数，use=“complete.obs”代表提取所有变量全部非空的观测，你也可以这样写：
> cor(blood.glucose,short.velocity,use="complete.obs")
[1] 0.4167546
我们还可以通过如下的代码得到一个数据框中多种变量的相关系数矩阵：
> cor(thuesen,use="complete.obs")
blood.glucose short.velocity
blood.glucose 1.0000000 0.4167546
short.velocity 0.4167546 1.0000000
#Tips：当然，数据框中变量超过两个结果会更有意思。
但是，上面的结果只是跟我们展示了关联系数，但是没有做出检验，因此我们需要cor.test()函数，可以通过指定两个变量来使用：
> cor.test(blood.glucose,short.velocity)
Pearson's product-moment correlation
data: blood.glucose and short.velocity
t = 2.101, df = 21, p-value = 0.0479
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.005496682 0.707429479
sample estimates:
cor
0.4167546
#Tips：我们也能得到一个真实相关系数的置信区间。注意，这里的p值和之前回归分析的p值是一样的。同样与之前回归模型的anova表里的p值是一样的。
B. 斯皮尔曼相关系数和肯德尔等级相关系数
与前面的部分所讲的单样本和双样本问题一样，相关问题也有非参数的方法，这些方法的优点在于不需要假设数据的正态分布性，而且结果也不会受到单调变换的影响。
相关性检验的几个方法都打包进了cor.test中，没有额外提供专门的spearman.test()函数。所以可以在cor.test()中指明：
> cor.test(blood.glucose,short.velocity,method="spearman")
Spearman's rank correlation rho
data: blood.glucose and short.velocity
S = 1380.4, p-value = 0.1392
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.318002
Warning message:
In cor.test.default(blood.glucose, short.velocity, method = "spearman") :
无法给连结计算精確p值
而等级相关同理，只需要在method参数中做出修改就可以实现：
> cor.test(blood.glucose,short.velocity,method="kendall")
Kendall's rank correlation tau
data: blood.glucose and short.velocity
z = 1.5604, p-value = 0.1187
alternative hypothesis: true tau is not equal to 0
sample estimates:
tau
0.2350616
Warning message:
In cor.test.default(blood.glucose, short.velocity, method = "kendall") :
无法给连结计算精確p值
#Tips：注意，这两个非参数方法的相关系数在5%置信水平下不显著，而pearson相关系数也仅仅是刚过线的显著。
到目前为止，简单的统计分析方法就告一段落了，而进一步的统计分析探讨我们会在下面一个系列为大家介绍，敬请期待。
参考资料：
1.《R语言统计入门（第二版）》 人民邮电出版社 Peter Dalgaard著
2.《R语言初学者指南》 人民邮电出版社 Brian Dennis著
3.Vicky的小笔记本《blooming for you》 by Vicky
网友评论