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

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包。

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 213,558评论 6赞 492
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 91,002评论 3赞 387
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 159,036评论 0赞 349
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 57,024评论 1赞 285
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 66,144评论 6赞 385
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 50,255评论 1赞 292
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 39,295评论 3赞 412
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 38,068评论 0赞 268
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 44,478评论 1赞 305
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 36,789评论 2赞 327
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 38,965评论 1赞 341
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 34,649评论 4赞 336
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 40,267评论 3赞 318
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 30,982评论 0赞 21
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,223评论 1赞 267
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 46,800评论 2赞 365
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 43,847评论 2赞 351

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)

推荐阅读更多精彩内容