R语言统计入门-第六章

Mar 20, 2018 12 分钟阅读时长 R, 统计

第六章 回归于相关性

6.1 简单线性回归

线性回归模型是如下定义的： \[ y_i=\alpha+\beta x_i+\epsilon_i \] 其中，假设\(\epsilon_i\)是独立的，并且来自于\(N(0,\sigma^2)\)。这个方程非随机的部分用一条直线描述\(y_i\)。这条直线的斜率（回归系数）是\(\beta\)，即\(x\)的每变化一单位给\(y\)所带来的增长。这条线与\(y\)轴交于截距点\(\alpha\)。

系数\(\alpha\)，\(\beta\)和\(\sigma^2\)都能用最小二乘法来估计。找到让残差平方和最小的\(\alpha\)，\(\beta\)： \[ SS_{res}=\sum_{i}[y_i-(\alpha+\beta x_i)]^2 \] 对于这些参数，能够推出使得SSres最小的显式表达式： \[ \hat \beta=\frac{\sum(x_i-\bar x)(y_i-\bar y)}{\sum(x_i-\bar x)^2} \]

\[ \hat\alpha =\bar y-\hat \beta\bar x \]

残差的方差可以通过SSres/(n-2)来估计，标准差自然是这个值的平方根。

本节使用thuesen作为例子

library(ISwR)
attach(thuesen)
#使用函数lm（线性模型）进行线性回归分析
#lm里面的参数是模型方程，波浪号~读作通过……来描述
lm(short.velocity ~ blood.glucose)
## 
## Call:
## lm(formula = short.velocity ~ blood.glucose)
## 
## Coefficients:
##   (Intercept)  blood.glucose  
##       1.09781        0.02196

函数lm的原始输出格式非常简单，只有估计出来的截距α和斜率β。可以看到最优拟合直线为short.velocity=1.098+0.0220*blood.glucose，但是没有给出任何像显著性检验之类的其他信息。

lm的输出结果是一个模型对象。能够通过用析取函数来得到想要的结果。

#一个基本的析取函数是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

下面对输出结果进行解剖：

call：

lm(formula = short.velocity ~ blood.glucose)

与t.test中相似，输出的开头本质上在重复一个函数调用，若用户只是在R中将其输出，那么这部分意义不大。但是如果结果被保存在一个变量中，之后查看输出的时候，这部分就很有用了。

Residuals:

min 1Q Median 3Q Max

-0.40141 -0.14760 -0.02202 0.03001 0.43490

这部分简单的描述了残差的分布，可以帮助用户对分不行的假设做快速检查。根据定义，残差的均值是0，所以中位数应该离0不远，然后最大值、最小值的绝对值也应该大致相当。这个例子中，第三分位数明显过于接近0，不过考虑到这里样本量不多，所以无须太过担心。

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 *

这里再次见到了回归系数和截距，不过这次还伴随着标准误， t 检验和 p 值。最右边的符号是显著性水平。这个符号可以通过options(show.signif.stars=FALSE)来关闭。

Residual standard error: 0.2167 on 21 degrees of freedom

(1 observation deleted due to missingness)

这是残差的波动情况，通过观测值在回归线附近的波动情况来估计模型参数σ。因为short.velocity有一个缺失值，所以这个模型并没有对整个数据集进行拟合。

Multiple R-squared: 0.1737， Adjusted R-squared: 0.1343

第一个是R2，在简单线性回归中能将其理解为Person相关参数的平方，即R2=r2。另一个是修正后R2；如果将其乘上100%，它就可以被理解为“方差降低的百分率”。

F-statistic: 4.414 on 1 and 21 DF, p-value: 0.0479

这是对假设回归系数是0而进行的F检验。这个检验在简单线性回归中并不特别，因为它只是对已有信息的重复——它在解释变量不止一个时就变得更有意义了。注意，这里的结果与斜率的 t 检验结果一模一样。实际上，F 检验是 t 检验的平方：4.414=(2.101)2。

#这里展示如何画出残差图
#首先把数据点和回归线画出来
plot(blood.glucose, short.velocity)
abline(lm(short.velocity ~ blood.glucose))

#abline函数就是(a,b)-线段的意思。根据截距α和斜率β画一条直线
#它能够接受数值参数，比如abline(1.1, 0.022)

6.2 残差与回归值

已经用summary从回归分析的结果中提取更多信息。另外两个进一步的析取函数是fitted和resid。它们的用法如下。为了方便，将lm的返回值存入变量lm.velo中。

lm.velo <- lm(short.velocity ~ blood.glucose)
#函数fitted返回的是回归值
#即根据最佳拟合直线与给定的x值计算出来的y值
#对于这个结果就是1.098+0.0220*blood.glucose
fitted(lm.velo)
##        1        2        3        4        5        6        7        8 
## 1.433841 1.335010 1.275711 1.526084 1.255945 1.214216 1.302066 1.341599 
##        9       10       11       12       13       14       15       17 
## 1.262534 1.365758 1.244964 1.212020 1.515103 1.429449 1.244964 1.190057 
##       18       19       20       21       22       23       24 
## 1.324029 1.372346 1.451411 1.389916 1.205431 1.291085 1.306459
#函数resid显示的残差是short.velocity的回归值与观测值之差
resid(lm.velo)
##            1            2            3            4            5 
##  0.326158532  0.004989882 -0.005711308 -0.056084062  0.014054962 
##            6            7            8            9           10 
##  0.275783754  0.007933665 -0.251598875 -0.082533795 -0.145757649 
##           11           12           13           14           15 
##  0.005036223 -0.022019994  0.434897199 -0.149448964  0.275036223 
##           17           18           19           20           21 
## -0.070057471  0.045971143 -0.182346406 -0.401411486 -0.069916424 
##           22           23           24 
## -0.175431237 -0.171085074  0.393541161

注意的是，回归值与残差都带着数据框thuesen的行名。同时，都没有第16个观测值的信息，因为原数据在这里缺少相应变量的值。

在这里有必要讨论缺少数据时出现的一些问题。

虽然用abline(lm.velo)更方便，可能仍会想到用lines将回归线画在图上，不过：

#plot(blood.glucose,short.velocity)
#lines(blood.glucose,fitted(lm.velo))
#Error in xy.coords(x, y) : 'x' and 'y' lengths differ
#Calls: local ... eval -> lines -> lines.default -> plot.xy -> xy.coords

就会出现这个情况，一共有24个观测值，但其中只有23个回归值，因为short.velocity的值是NA。碰巧的是，这个错误在一串交织在一起的函数调用中出现，它们都在报错信息中标出来以便于理解。

我们需要的是blood.glucose，但是需要病人的short.velocity也被记录了才行。

#函数is.na能够产生一个向量，在参数为NA的对应位置上标记为TRUE
#这个方法的好处之一就是回归线不会超过数据的范围
#lines(blood.glucose [!is.na(short.velocity)],fitted(lm.velo))
#但是在许多个变量中都有缺失值的时候就变得很繁琐
#...blood.glucose[!is.na(short.velocity) & !is.na(blood.glucose)]...
#用函数complete.cases更简单
#它能够筛选出在若干个变量甚至是整个数据框中都没有缺失值的观测
cc <- complete.cases(thuesen)
#然后就能用thuesen[cc,]来进行分析了
#不过有更好的方法
#可以用na.exclude方法处理缺失值
#它既可以作为lm的一个参数，也可以作为一个选项
options(na.action=na.exclude)
lm.velo <- lm(short.velocity ~ blood.glucose)
fitted(lm.velo)
##        1        2        3        4        5        6        7        8 
## 1.433841 1.335010 1.275711 1.526084 1.255945 1.214216 1.302066 1.341599 
##        9       10       11       12       13       14       15       16 
## 1.262534 1.365758 1.244964 1.212020 1.515103 1.429449 1.244964       NA 
##       17       18       19       20       21       22       23       24 
## 1.190057 1.324029 1.372346 1.451411 1.389916 1.205431 1.291085 1.306459
#注意的是，第16个观测带有一个缺失值出现了
#在改变了选项后，有必要重新计算lm.velo
#为了在一幅图中通过将观测值与对应的回归值连起来而显示出残差
#函数segments用来画线段，它的参数是两端断点的坐标（x1,y1,x2,y2）
#但需要注意，segments函数不能单独使用，需要先plot一个图
#以下的plot与abline前文中已经有了，但是在RMD中最好还是重新调用下
plot(blood.glucose,short.velocity)
abline(lm(short.velocity ~ blood.glucose))
segments(blood.glucose,fitted(lm.velo),
        blood.glucose,short.velocity)

#也能通过Q-Q图的线性性来检验残差的正态性
qqnorm(resid(lm.velo))

6.3 预测与置信带

回归线通常与不确切的边界带一起展示。一般有两种边界带，通常被称为“窄”边界和“宽”边界。 窄边界，又叫置信带，反映了这条线本身的不确定性，就像SEM反映了一个已知均值的准确度。如果观测数量很多的话，这个边界会很窄，这意味着这是一个比较准确的线。这个边界通常有明显的弧度，因为回归线在点阵的中心通常更准确。 宽边界，又叫预测带，包含了未来观测值的不确定性。这些边界关心大部分数据，同时不会随着观测数量的增加而缩成一条线。这个边界是由回归线±2倍的标准差（95%的水平）而得到的。在小样本情况下，这个边界也是有弧度的，因为它们包含了这条直线本身的不确定性，只不过这个弧度没有置信带的那么明显。很显然，这些边界十分依赖残差正态性以及方差齐性的假设，所以如果你的数据不太满足这些性质，那么最好不要用这个边界。 无论是否计算了置信带与预测带，我们都能用函数predict析取出预测值。不加其他参数的话，它就只会输出回归值。

predict(lm.velo)
##        1        2        3        4        5        6        7        8 
## 1.433841 1.335010 1.275711 1.526084 1.255945 1.214216 1.302066 1.341599 
##        9       10       11       12       13       14       15       16 
## 1.262534 1.365758 1.244964 1.212020 1.515103 1.429449 1.244964       NA 
##       17       18       19       20       21       22       23       24 
## 1.190057 1.324029 1.372346 1.451411 1.389916 1.205431 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
## 4  1.526084 1.306561 1.745607
## 5  1.255945 1.139367 1.372523
## 6  1.214216 1.069315 1.359118
## 7  1.302066 1.205244 1.398889
## 8  1.341599 1.246317 1.436881
## 9  1.262534 1.149694 1.375374
## 10 1.365758 1.263750 1.467765
## 11 1.244964 1.121641 1.368287
## 12 1.212020 1.065457 1.358583
## 13 1.515103 1.305352 1.724854
## 14 1.429449 1.290217 1.568681
## 15 1.244964 1.121641 1.368287
## 16       NA       NA       NA
## 17 1.190057 1.026217 1.353898
## 18 1.324029 1.230050 1.418008
## 19 1.372346 1.267629 1.477064
## 20 1.451411 1.295446 1.607377
## 21 1.389916 1.276444 1.503389
## 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")
## Warning in predict.lm(lm.velo, int = "p"): predictions on current data refer to _future_ responses
##         fit       lwr      upr
## 1  1.433841 0.9612137 1.906469
## 2  1.335010 0.8745815 1.795439
## 3  1.275711 0.8127292 1.738693
## 4  1.526084 1.0248161 2.027352
## 5  1.255945 0.7904672 1.721423
## 6  1.214216 0.7408499 1.687583
## 7  1.302066 0.8411393 1.762993
## 8  1.341599 0.8809929 1.802205
## 9  1.262534 0.7979780 1.727090
## 10 1.365758 0.9037136 1.827802
## 11 1.244964 0.7777510 1.712177
## 12 1.212020 0.7381424 1.685898
## 13 1.515103 1.0180367 2.012169
## 14 1.429449 0.9577873 1.901111
## 15 1.244964 0.7777510 1.712177
## 16       NA        NA       NA
## 17 1.190057 0.7105546 1.669560
## 18 1.324029 0.8636906 1.784367
## 19 1.372346 0.9096964 1.834996
## 20 1.451411 0.9745421 1.928281
## 21 1.389916 0.9252067 1.854626
## 22 1.205431 0.7299634 1.680899
## 23 1.291085 0.8294798 1.752690
## 24 1.306459 0.8457315 1.767186

结果中fit变量表示了期望得到的值，在这里就等于回归值，而lwr和upr就分别是下界和上界，即对那些blood.glucose取此值的病人预测short.velocity时的边界。警告信息并不是说有什么做错了，而是提醒这里有一个陷阱：这个边界不能被用来考量我们用来做回归的已观测数据。因为在x的极值处，数据影响力更大，所以这会导致这些地方的边界离回归线更近，也就是说，预测带弯向了错误的方向。 将置信带与预测带加到散点图上的最好方法是通过函数matlines，它能将矩阵的每一列以某一个向量作为x轴画出来。 不过，有几个小障碍： + blood.glucose的值是随机排列的，不希望置信曲线上的线段杂乱无章的排列； + 下方的预测区间超出了画图区域； + matlines的命令需要放置不停更迭的线段样式和颜色。 解决方法是用合适的x（这里是blood.glucose）生成一个新数据框，然后在新数据框上进行预测：

#用希望进行预测的blood.glucose值生成了一个新的数据框
pred.frame <- data.frame(blood.glucose=4:20)
#pp和pc用来记录predict函数在pred.frame上的结果，并且保留了预测带和置信带
pp <- predict(lm.velo, int="p", newdata = pred.frame)
pc <- predict(lm.velo, int="c", newdata = pred.frame)
#先plot画一个标准的散点图，并为预测带预留了足够的空间，即参数ylim
#函数range返回一个长度为2的向量，其中是传入参数的最大值和最小值
#并需要使用na.rm=T在计算中忽略缺失值
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 = "green")

6.4 相关性

相关系数是一个对称并且不随尺度变化的量，用于衡量两个随机变量之间的关联程度。它的值域是-1到1，这两个极端表示完美的相关，0则表示没有相关性。一个变量的较大值与另一个变量的较小值有关联时，相关系数是负的；如果两个变量有同时变大或减小的趋势，那么相关系数是负的。

6.4.1 皮尔逊相关系数

皮尔逊相关系数扎根于二维正态分布中，其中理论上的相关性描述了密度函数的椭圆等高线。如果两个变量的方差都变换成了1，那么相关系数为0就对应圆形的等高线，其他情况下，椭圆变得越来越窄，最后在相关系数等于±1的时候坍缩成一条直线。 经验相关系数是： \[ r=\frac{\sum(x_i-\bar x)(y_i-\bar y)}{\sqrt{\sum(x_i-\bar x)^2\sum(y_i-\bar y)}} \] 可以证明，在\(x_i\)和\(y_i\)没有完美线性相关的情况下，\(\left| r \right|\)恒小于1，皮尔逊相关系数也因此又是被称为“线性相关系数”。 函数cor能计算两个或多个向量之间的相关系数。但如果对thuesen中的两个向量也这样操作，就会出现下面情况：

cor(blood.glucose,short.velocity)
## [1] NA
#Error in cor(blood.glucose,short.velocity):
#           missing observations in cov/cor

R中所有的基本统计函数都要求输入的参数没有缺失值，或者要明确指定如何处理缺失值。对于函数mean，var，sd以及类似的单向量函数，可以传递na.rm=T这个参数告诉它们在计算之前应该移除缺失值。对于函数cor，可以写成：

cor(blood.glucose,short.velocity,use = "complete.obs")
## [1] 0.4167546

函数cor不使用na.rm=T，因为在移除缺失值与计算出错之外还有很多其他可能性。如果要对超过两个变量进行计算，也能用所有非缺失值对的信息来进行计算。 可以通过如下代码得到一个数据框中所有变量的相关系数矩阵：

cor(thuesen,use = "complete.obs")
##                blood.glucose short.velocity
## blood.glucose      1.0000000      0.4167546
## short.velocity     0.4167546      1.0000000

但是，上面的计算并没有告诉我们这个相关系数是不是显著不为0的。为做到这点，需要使用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

6.4.2 斯皮尔曼相关系数

同样也有非参数的方法。即斯皮尔曼秩相关系数\(\rho\)。它将观测值替换为它们的秩，再计算相关系数。在两变量独立的零假设下，可以计算出\(\rho\)的精确分布。 与之前一个函数对应一个检验不同，相关性检验的集中方法都打包进了cor.test中。

cor.test(blood.glucose,short.velocity,method = "spearman")
## Warning in cor.test.default(blood.glucose, short.velocity, method =
## "spearman"): Cannot compute exact p-value with ties
## 
##  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

6.4.3 肯德尔等级相关系数τ

第三个可供选择的算法是肯德尔等级相关系数，这个算法基于统计一致对和不一致对的数量。如果x坐标之差的符号与y坐标之差是一致的，那么这一对点是个一致对。在一个完美的单调关系中，所有的点对要么都是一致的；要么都是不一致的。在独立的情况下，一致对和不一致对的数量应该一样多。 相关系数\(\tau\)的一大优点是其意义比斯皮尔曼更好解释，但是除此之外，就没有什么优于对方的特点了。

cor.test(blood.glucose,short.velocity,method = "kendall")
## Warning in cor.test.default(blood.glucose, short.velocity, method =
## "kendall"): Cannot compute exact p-value with ties
## 
##  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

练习题

  1. 在rmr数据集中，画出代谢率关于体重的散点图。对这个关系拟合一个曲线。这个模型预测70 kg的体重对应的代谢率是多少？对这条线的斜率给出一个95%置信区间。
  2. 在juul数据集中，拟合一个IGF-I集中度的平方根关于25岁以上人群年龄的线性回归模型。
  3. 在malaria数据集中，分析对数变换后的antibody水平关于年龄的关系。画出一个关系图。
  4. 一个人能够这样对指定\(\rho\)的二维正态分布中生成模拟数据：（a）以均值为0，标准差为1的正态分布生成\(X\);（b）以均值为\(\rho X\)，标准差为\(\sqrt{1-\rho^2}\)的正态分布生成\(Y\)。对于指定的相关系数，用这个方法画出模拟数据的散点图。计算部分数据的斯皮尔曼相关系数与肯德尔等级相关系数。
R 统计
Avatar
Dr.二哈
在读苦逼科研狗

研究方向：脂质营养，业余时间自学R。

相关