DALS016-Linear Models线性模型05共线性

2019-08-27 本文已影响0人 backup备份

title: DALS016-Linear Models线性模型05共线性
date: 2019-08-20 12:0:00
type: "tags"
tags:

矩阵计算
线性模型
categories:
生物统计

前言

这一部分是《Data Analysis for the life sciences》的第5章线性模型的第5小节，这一部分的主要内容涉及矩阵的共线性(Co-linearity)，秩(rank)，QR分解，有关共线性和秩的笔记R markdown文档可以参考作者的Github，QR分解见另外一篇。

共线性(Co-linearity)

如果一个实验设计不合理，那么我们就无法对目标参数进行估计。同样，当我们分析数据时，不合理地使用一个模型时也会得错误的结论。我们使用线性模型时，通过检查设计矩阵的共线性(collinearity)就能从数学角度来考虑我们选择的这个模型是否合理。

方程组案例

我们来看一个方程组，如下所示：
$\begin{align*} a+c &=1\\ b-c &=1\\ a+b &=2 \end{align*}$
这个方程就不止一个解，因为有无数组数据满足a=1-c，b=1+c。例如a=1, b=1, c=0和a=0,b=c,c=1。

矩阵代数方法

现在我们使用矩阵来说明一下上面的情况，上面的方程组我们写成矩阵形式如下所示：
$\, \begin{pmatrix} 1&0&1\\ 0&1&-1\\ 1&1&0\\ \end{pmatrix} \begin{pmatrix} a\\ b\\ c \end{pmatrix} = \begin{pmatrix} 1\\ 1\\ 2 \end{pmatrix}$
从矩阵中我们可以发现，其实矩阵的第3列就是第1列和第2列的线性组合，如下所示：
$$
,
\begin{pmatrix}
1\
0\
1
\end{pmatrix}

-1 \begin{pmatrix}
0\
1\
1
\end{pmatrix}
=
\begin{pmatrix}
1\
-1\
0
\end{pmatrix}
$因此，我们可以说第3列与前2列共线性(collinear)，这也就是说，我们可以将上面的矩阵写成下面的这个形式：$
,
\begin{pmatrix}
1&0&1\
0&1&-1\
1&1&0
\end{pmatrix}
\begin{pmatrix}
a\
b\
c
\end{pmatrix}
=
a
\begin{pmatrix}
1\
0\
1
\end{pmatrix}

b \begin{pmatrix}
0\
1\
1
\end{pmatrix}

c
\begin{pmatrix}
1-0\
0-1\
1-1
\end{pmatrix}
$$

$$
=(a+c)
\begin{pmatrix}
1\
0\
1\
\end{pmatrix}

(b-c)
\begin{pmatrix}
0\
1\
1\
\end{pmatrix}
$$

新的矩阵也就说明了，第3列其实没什么用，我们真正有的是3个方程和2个未知数，即a+c和b+c。一旦我们得到了这两个量（即a+c和b+c），那么a, b, c这个未知数就有无数种解。

共线性和最小平方

现在我们考虑一个设计矩阵 $\mathbf{X}$ ，这个矩阵含有2个共线性列。我们创建一个极端的例子，在这个例子中，有两列互为相反，如下所示：
$\mathbf{X} = \begin{pmatrix} \mathbf{1}&\mathbf{X}_1&\mathbf{X}_2&\mathbf{X}_3\\ \end{pmatrix} \mbox{ with, say, } \mathbf{X}_3 = - \mathbf{X}_2$
这就意味着，我们可以把上面的矩阵写为如下形式：
$\mathbf{Y}- \left\{ \mathbf{1}\beta_0 + \mathbf{X}_1\beta_1 + \mathbf{X}_2\beta_2 + \mathbf{X}_3\beta_3\right\}\\ = \mathbf{Y}- \left\{ \mathbf{1}\beta_0 + \mathbf{X}_1\beta_1 + \mathbf{X}_2\beta_2 - \mathbf{X}_2\beta_3\right\}\\ = \mathbf{Y}- \left\{\mathbf{1}\beta_0 + \mathbf{X}_1 \beta_1 + \mathbf{X}_2(\beta_2 - \beta_3)\right\}$

如果 $\hat{\beta}_1$ , $\hat{\beta}_2$ , $\hat{\beta}_3$ 是一个最小平方解(least squares solution)，那么，例如 $\hat{\beta}_1$ , $\hat{\beta}_2+1$ , $\hat{\beta}_3+1$ 也是一个解。

混淆(confounding)案例分析

现在我们将演示一下共线性是如何帮助我们解决当前实验设计中最常见的错误之一的：混淆(confounding)。为了说明这个问题。现在我们设想一个实验，例如我们感兴趣的是四种处理因素，即A，B，C和D。我们给每种处理因素分配2只小鼠。当我们给A、B处理分配的是雌性小鼠后，我们意识到性别(sex)有可能对实验造成影响。我们决定给C和D两组雄性小鼠，从而估计性别对实验的影响。但是，通过这种方式我们能否估计出性别的影响呢，这个实验我们表示为矩阵形式，则如下所示：
$\, \begin{pmatrix} Sex & A & B & C & D\\ 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & 1\\ \end{pmatrix}$
在这个矩阵里，我们可以看到处理因素与性别产生了混淆(confound)，尤其是，性别这一列可以用C列和D列进行表示（也就是说共线性），如下所示：
$$
,
\begin{pmatrix}
Sex \
0\
0 \
0 \
0 \
1\
1\
1 \
1 \
\end{pmatrix}
=
\begin{pmatrix}
C \
0\
0\
0\
0\
1\
1\
0\
0\
\end{pmatrix}

\begin{pmatrix}
D \
0\
0\
0\
0\
0\
0\
1\
1\
\end{pmatrix}
$$
这就说明了，无法使用最小平方和来对实验结果进行计算。

秩(rank)

一个矩阵的秩(rank)表示的是，矩阵中独立于其他列的数目。如果矩阵的秩小于列的数目，那么LSE则不唯一。在R中，我们可以使用qr()函数来获取一个矩阵的秩，如下所示：

Sex <- c(0,0,0,0,1,1,1,1)
A <- c(1,1,0,0,0,0,0,0)
B <- c(0,0,1,1,0,0,0,0)
C <- c(0,0,0,0,1,1,0,0)
D <- c(0,0,0,0,0,0,1,1)
X <- model.matrix(~Sex+A+B+C+D-1)
cat("ncol=",ncol(X),"rank=", qr(X)$rank,"\n")

结果如下所示：

ncol= 5 rank= 4

移除混淆

前面我们提到的那个实验可以设计得更加完全。使用相同数量的雄性小鼠和雌性小鼠可以很容易设计一个实验来研究性别影响与干预影响。具体来说就是，当我们平衡了性别与干预后，混淆就会被消除，也就是说，通过一定的合理实验设计，就是使秩与列的数目相同，如下所示：

Sex <- c(0,1,0,1,0,1,0,1)
A <- c(1,1,0,0,0,0,0,0)
B <- c(0,0,1,1,0,0,0,0)
C <- c(0,0,0,0,1,1,0,0)
D <- c(0,0,0,0,0,0,1,1)
X <- model.matrix(~Sex+A+B+C+D-1)
cat("ncol=",ncol(X),"rank=", qr(X)$rank,"\n")

计算结果如下所示：

ncol= 5 rank= 5

*练习

P237

矩阵的QR分解

这一部分的R markdown参考文档见Github。

前面我们已经了解到了通过逆转一个矩阵来计算LSE，前面使用了solve()函数。但是，这并非是一种稳定的解决方式方式。当我们进行LSE计算时，可以使用QR分解。

逆转

为了使RSS最小：
$(\mathbf{Y}-\mathbf{X}\boldsymbol{\beta})^\top (\mathbf{Y}-\mathbf{X}\boldsymbol{\beta})$
我们需要解下面方程：
$\mathbf{X}^\top \mathbf{X} \boldsymbol{\hat{\beta}} = \mathbf{X}^\top \mathbf{Y}$
解为：
$\boldsymbol{\hat{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{Y}$
因此，我们需要计算 $(\mathbf{X}^\top \mathbf{X})^{-1}$

这个解在数学上是不稳定的(solve is numerically unstable)。

为了说明什么叫在数学上是不稳定的，我们来看下面的代码：

n <- 50;M <- 500
x <- seq(1,M,len=n)
X <- cbind(1,x,x^2,x^3)
colnames(X) <- c("Intercept","x","x2","x3")
beta <- matrix(c(1,1,1,1),4,1)
set.seed(1)
y <- X%*%beta+rnorm(n,sd=1)
solve(crossprod(X)) %*% crossprod(X,y)

运行后就出错了，如下所示：

> solve(crossprod(X)) %*% crossprod(X,y)
Error in solve.default(crossprod(X)) : 
  system is computationally singular: reciprocal condition number = 2.93617e-17

现在我们来看一下为什么会出现这种错误，如下所示：

options(digits=4)
log10(crossprod(X))

结果如下所示：

> options(digits=4)
> log10(crossprod(X))
          Intercept      x     x2     x3
Intercept     1.699  4.098  6.625  9.203
x             4.098  6.625  9.203 11.810
x2            6.625  9.203 11.810 14.434
x3            9.203 11.810 14.434 17.070

这里我们要注意一下数量级的差异。常用的数字计算机储存的数字范围有限。当我们要要考虑使用非常大的数字时，这会使一些数字看起来像0。这反过来又会导致由0划分导致的错误划分（原文是：This in turn leads to divisions that are practically divisions by 0 errors）。

分解(factorization)

QR分解是基于一种数学结果，它能告诉我们可以分解任意的全秩(any full rank)的 $N \times p$ 矩阵 $\mathbf{X}$ ，如下所示：
$\mathbf{X}=\mathbf{QR}$
其中：

$\mathbf{X}$ 是一个 $N \times p$ 矩阵， $\mathbf{Q}^{T}\mathbf{Q}=I$
$\mathbf{R}$ 是一个 $p \times p$ 矩阵，上三角矩阵。

上三角矩阵在解方程组方面非常有用。

什么是上三角矩阵

在下面的案例中，最左边的就是上三角矩阵：上三角矩阵的特点就是，矩阵对角线的左下方都是0，右上方不为零，根据矩阵的运算规则，这就很容易解方程组，如下所示：
$\, \begin{pmatrix} 1&2&-1\\ 0&1&2\\ 0&0&1\\ \end{pmatrix} \begin{pmatrix} a\\ b\\ c\\ \end{pmatrix} = \begin{pmatrix} 6\\ 4\\ 1\\ \end{pmatrix}$
从上面的形式我们就能马上看出来， $c=1, b + 2 = 4, a = 3$ 。

利用QR来计算LSE

如果我们重新使用 $\mathbf{QR}$ 来替代 $\mathbf{X}$ 来写LSE的公式，则如下所示：

$\mathbf{X}^\top \mathbf{X} \boldsymbol{\beta} = \mathbf{X}^\top \mathbf{Y}$

$(\mathbf{Q}\mathbf{R})^\top (\mathbf{Q}\mathbf{R}) \boldsymbol{\beta} = (\mathbf{Q}\mathbf{R})^\top \mathbf{Y}$

$\mathbf{R}^\top (\mathbf{Q}^\top \mathbf{Q}) \mathbf{R} \boldsymbol{\beta} = \mathbf{R}^\top \mathbf{Q}^\top \mathbf{Y}$

$\mathbf{R}^\top \mathbf{R} \boldsymbol{\beta} = \mathbf{R}^\top \mathbf{Q}^\top \mathbf{Y}$

$(\mathbf{R}^\top)^{-1} \mathbf{R}^\top \mathbf{R} \boldsymbol{\beta} = (\mathbf{R}^\top)^{-1} \mathbf{R}^\top \mathbf{Q}^\top \mathbf{Y}$

$\mathbf{R} \boldsymbol{\beta} = \mathbf{Q}^\top \mathbf{Y}$

其中， $\mathbf{R}$ 是上三角矩阵，它能够更加稳定地求出LSE。此外，由于 $\mathbf{Q}^T\mathbf{Q}=\mathbf{I}$ ，我们知道 $\mathbf{Q}$ 的列能够稳定右侧（这一段不太懂，原文：we know that the columns of Q are in the same scale which stabilizes the right side)。

现在我们就可以通过QR分解来计算LSE了，也就是解下面的这个方程：

$\mathbf{R} \boldsymbol{\beta} = \mathbf{Q}^\top \mathbf{Y}$

在R中我们使用backsolve()函数来角这个方程，如下所示：

QR <- qr(X)
Q <- qr.Q( QR )
R <- qr.R( QR )
(betahat <- backsolve(R, crossprod(Q,y) ) )

结果如下所示：

> QR <- qr(X)
> Q <- qr.Q( QR )
> R <- qr.R( QR )
> (betahat <- backsolve(R, crossprod(Q,y) ) )
          [,1]
[1,] 0.9038372
[2,] 1.0066440
[3,] 0.9999622
[4,] 1.0000001

实际上，R中有内置函数solve.qr()也能实现此功能，如下所示：

QR <- qr(X)
(betahat <- solve.qr(QR, y))

结果如下所示：

> QR <- qr(X)
> (betahat <- solve.qr(QR, y))
               [,1]
Intercept 0.9038372
x         1.0066440
x2        0.9999622
x3        1.0000001

拟合值

这种分解过程也能简化拟合值的计算，如下所示：
$\mathbf{X}\boldsymbol{\hat{\beta}} = (\mathbf{QR})\mathbf{R}^{-1}\mathbf{Q}^\top \mathbf{y}= \mathbf{Q}\mathbf{Q}^\top\mathbf{y}$

在R中，计算如下所示：

library(rafalib)
mypar(1,1)
plot(x,y)
fitted <- tcrossprod(Q)%*%y
lines(x,fitted,col=2)

结果如下所示：

image

标准误

如果要计算LSE的标准误，需要考虑如下公式：

$(\mathbf{X^\top X})^{-1} = (\mathbf{R^\top Q^\top QR})^{-1} = (\mathbf{R^\top R})^{-1}$
使用chol12inv()函数是专门用于计算这个逆矩阵的，如下所示：

df <- length(y) - QR$rank
sigma2 <- sum((y-fitted)^2)/df
varbeta <- sigma2*chol2inv(qr.R(QR))
SE <- sqrt(diag(varbeta))
cbind(betahat,SE)
summary(lm(y~0+X))$coef

结果如下所示：

> cbind(betahat,SE)
                              SE
Intercept 0.9038372 4.508440e-01
x         1.0066440 7.858164e-03
x2        0.9999622 3.661705e-05
x3        1.0000001 4.802429e-08
> summary(lm(y~0+X))$coef
            Estimate   Std. Error      t value      Pr(>|t|)
XIntercept 0.9038372 4.508440e-01 2.004767e+00  5.089408e-02
Xx         1.0066440 7.858164e-03 1.281017e+02  2.171355e-60
Xx2        0.9999622 3.661705e-05 2.730865e+04 1.745091e-167
Xx3        1.0000001 4.802429e-08 2.082280e+07 4.558613e-300

其它线性模型

这一部分的参考R markdown文档见Github。

线性模型的用途使用非常广，并且它还可以扩展到其它方面，后面就介绍一些线性模型的拓展。

稳健线性模型(Robust linear models)

在标准线性模型里，我们在计算它的解以及估计误差时，我们通常会通过计算最小平方误差来实现目标。这就涉及所有点的平方和，同时，这也表明，一旦数据中存在着异常值，那么这些异常值就会对计算结果产生重要的影响。此外，我们在使用线性模型时，我们是假定误差是恒定的（也就是同方差性，或齐方差性，homoskedasticity），不过在实际运用过程中，这种假设并非总成立（如果误差不剂，我们称为异方差性，heteroskedasticity）。因此需要其它一些方法用于产生稳健的解决方案，这些新的方法可以很好地处理异常值，或者是当我们的假设（也是同方差时）不满足时，这些模型也能很好地实现目的。在R的官网上，有介绍这些方差的内容《CRAN Task View: Robust Statistical Methods》。

广义线性模型(Generalized linear models,GLM)

在标准线性模型里，我们对 $\mathbf{Y}$ 的分布没做任何假设，但是在某些情况下，如果我们知道了 $\mathbf{Y}$ ，例如 $\mathbf{Y}$ 仅限于非负整数 $0,1,2,\dots$ ，或者是限于区间 $[0,1]$ 。用于分析这类情况的框架被称为广义线性模型，通常缩写为GLMs。GLM的两个关键组成部分就是链接函数与概率分布。链接函数 $g$ 会通常以下方式将我们熟悉的矩阵乘积 $\mathbf{X} \boldsymbol{\beta}$ 与 $\mathbf{Y}$ 值链接起来：

$\textrm{E}(\mathbf{Y}) = g^{-1}( \mathbf{X} \boldsymbol{\beta} )$

在R中，使用glm()函数来拟合GLMs，并使用与lm()函数类似的形式来进行计算。其它的参数还包括family，它用于指定 $\mathbf{Y}$ 的分布假设。在QuickR网站上，我们可以看到很多GLM的使用案例。

混合效应线性模型(Mixed effects linear models)

在标准线性模型里，我们假设矩阵 $\mathbf{X}$ 是固定的(fixed)，并不随机。例如，我们检测了每对腿的摩擦系数，以及push和pull方向。在 $\mathbf{X}$ 矩阵中，如果在特定的列中含有1，那么这1列就不是随机的，而是由设计矩阵决定的。但是在父子高度的案例中，我们并没有固定父亲身高的数值，而是观察到了这些值（通过一些误差来对抽样数据进行检测）。用于研究 $\mathbf{X}$ 中各个列的随机性影响的框架被称为混合效应模型，这就意味着一些效应是固定的，一些是效应是随机的。R中用于拟合线性模型的一个包是lme4，这个包还发了相应的论文《 Fitting Linear Mixed-Effects Models using lme4》。

贝叶斯线性模型(Bayesian linear models)

这种方法呈现的是假设 $\beta$ 是一个固定参数（非随机）。我们提取了估计此参数的方法，以及估计过程中如何量化这些参数不确定性的标准误差。这种方法是就所谓的频率论方法(frequentist)。另一种方法是假设 $\beta$ 是随机的，它的分布可以对我们先前关于 $\beta$ 应该是什么进行量化。一旦我们观察到数据，那么我们通过计算给定数据的 $\beta$ 条件分布（称为后验分布,posterior)，我们就就更新我们先前的信息，这种方法称贝叶斯方法。例如，一旦我们计算了 $\beta$ 的后验分布，我们就可以计算出某事件最容易发生的一个区间（可信区间）。此外，许多模型可以连接在一起，也就是所谓的分层模型(hierarchical mode)。在后面的内容里，我们会介绍贝叶斯统计和层次模型。贝叶斯分层模型的一本很好的书是《Bayesian Data Analysis》，在CRAN里，也有计算贝叶斯线性模型的包。其它的一些计算贝叶斯模型的软件是stan和BUGS。

惩罚线性模型(Penalized linear models)

如果我们在模型中包含足够的参数，我们就可以实现残差平方和为0。惩罚纯情模型引入了一个惩罚项到最小化的最小二乘方程里。这些惩罚项通常是 $\lambda \sum_{j=1}^p \|\beta_j\|^k$ 这个形式，它们对大的绝对值以及大量的参数进行惩罚。添加这个额外的参数的目的就是避免过拟合。如果要使用这些模型，我们需要选择使用哪种模式来决定我们的惩罚程度。当 $k=2$ 时，这种模型称为岭回归(ridge regression)，Tikhonov正则化(regularization)L2正则化。当 $k=1$ 时，被称为LASSO或L1正则化。有关这些惩罚线性模型的一本书是《Elements of Statistical Learn》（网上有免费的电子版）。与这些模型有关的R函数包括MASS包里的lm.ridge函数，还有lars包，glmnet包。

DALS016-Linear Models线性模型05共线性

前言

共线性(Co-linearity)

方程组案例

矩阵代数方法

共线性和最小平方

混淆(confounding)案例分析

秩(rank)

移除混淆

*练习

矩阵的QR分解

逆转

分解(factorization)

什么是上三角矩阵

利用QR来计算LSE

拟合值

标准误

其它线性模型

稳健线性模型(Robust linear models)

广义线性模型(Generalized linear models,GLM)

混合效应线性模型(Mixed effects linear models)

贝叶斯线性模型(Bayesian linear models)

惩罚线性模型(Penalized linear models)

猜你喜欢

热点阅读