如何求解单边z变换_简单高效的求解线性混合模型（多层线性模型）

本文面向教育统计，心理统计等社会科学统计入门者。
本文阐述了如何用极大似然估计（ML）和约束极大似然估计（REML）对线性混合模型（多层线性模型）进行参数估计，给出了对数似然函数，梯度和黑塞矩阵的计算过程，以及优化矩阵求逆的Woodbury恒等式，并且用python实现了参数估计、标准误、P值、AIC和BIC的计算。
本文有全网唯一正确的线性混合模型的黑塞矩阵解析式

起因

写这篇文章是因为发现github上拥有3K星星的著名计量经济学库statsmodels的线性混合模型的黑塞矩阵计算竟然是错的，而且没人提issue和pull request，其参考的论文《Newton Raphson and EM algorithms for linear mixed effects models for repeated measures data》在梯度和黑塞矩阵上的解析式亦存在大坑，故写此篇。本文主要参考《Generalized, Linear, and Mixed Models》这本书和上面的论文。

模型

线性混合模型，也称多层线性模型(以下简称lmm)，是社科常用的统计模型之一。最简单的lmm如下

构造

，并构造

线性混合模型可以改写为

可以看做随机误差

构造

，所以

所以

令

所以

多元正态分布极大似然估计的对数似然函数直接套公式得

这里

是稀疏矩阵，极其不利于计算，所以需要对

做处理

利用分块矩阵的性质

，

并令

我们把对数似然函数改写为

基于极大似然估计的对数似然公式对

求导

先求微分

是

组的样本量

所以

令

得到

计算

的期望，注意到

是标量，下面的计算参考了《The Matrix Cookbook》第六章的内容，有兴趣也可以看看《linear regression analysis》关于求线性模型期望的内容

由于

，

是

的列数，也即

的元素数，所以这个估计不是无偏估计。普通线性模型不在意这个细节，是因为普通线性模型不在意

估计，而线性混合模型在意

的估计，这个问题就很严重了（当

，这个问题无伤大雅）。

有种情况可以让方差估计为无偏估计，一是模型不包含固定效应参数，即

,

，二是通过矩阵变换，将固定效应参数干掉，即找到一个矩阵

使得

，也即

， 例如

（注意到

），但是这个矩阵的性质不利于后续的计算，所以需要构造一个更利于计算的矩阵，例如

，

，经过计算可以得到一个省略参数

的似然函数（详细的计算篇幅比较大，得专门讲解）。REML的对数似然函数一般为

注意到

对

求导并进一步计算得

式是

的无偏估计。

把

式代入

式，把

代入

式，可以得到不含

的ML对数似然函数和REML对数似然函数，分别为

其中

其中

需要估计的参数

,

和

的下三角或上三角元素。

若

有

列，

有

列，则极大似然估计的参数估计个数是

lmm的参数估计过程

第一种方法：选取

的初值

（一般是

），把

带入

式计算

的广义最小二乘估计

，把

带入

式或

式迭代求解

，把估计得到的

带入

式再求解

的广义最小二乘估计

，最后把

和

带入

式或

式求解

。

第二种方法：选取

的初值

、

的初值

，依据

式或

式迭代求解

和

，最后把估计得到的

和

带入

式或

式求解

。

其他方法，诸如EM，MCMC，变分等方法不在此文讨论范围中。

广义最小二乘

由于lmm参数估计的第一种方法设计涉及广义最小二乘法，这里稍微介绍下。

我们知道，传统的线性模型形式如下

其中

，由最小二乘法得

如果

，将

式左乘

（

的选取方法详情可以参考王松桂的线性模型书）

所以

，所以套用最小二乘法得

上式便是广义最小二乘法

一阶导数

迭代求解

或

需要求

或

的偏导数，为了方便用numpy,r或eigen写代码，我们的求导采取矩阵求导的方式，关于矩阵求导可以参考《The Matrix Cookbook》第二章和第十章的内容。

先写一个接下来常用的微分等式

对

求微分

对

求导

求

把

式带入

式，再把

式带入

式

并令

，

所以

又

所以

求

把

式带入

式，再把

式带入

式，得到

又

所以

求

把

带入

式，得到

令

，

，把

式带入

式

所以

所以

对固定系数

求导

把

代入

，把

和

所以

令

所以

黑塞矩阵（二阶导数）

黑塞矩阵可以用于牛顿迭代的优化，也可以用于构造参数的标准误、P值等信息，黑塞矩阵即观察信息矩阵。

求

对

式求微分

先求

， 运用

式结果

依据

式得

再求

，运用

式结果得到

依据

式

所以

依据

式

依据

式

最后求

把

式代入

式

由

式得

所以

综合

式得到

由上可知，计算二阶导数很冗长，为节省篇幅，以下只给出计算最终结果

Woodbury恒等式

由于对数似然方程和score方程大量存在

的形式， 设

的形状为

，则

为样本量，

为特征数量，一般来说，

远小于

，当

很大时，求

相当浪费时间，为了追求效率，我们可以用Woodbury恒等式简化求逆（Woodbury等式的内容参考自《The Matrix Cookbook》第三章），

在上式中，需要计算逆的矩阵的大小为

，计算量大大小于

我们为上式加一个右乘，可得

进一步简化了我们的计算

代码实现

我们用python实现上面的内容，具体代码如下

import numpy as np
from scipy.optimize import minimize


def gls(group_y, group_x, group_z, group_size, d):
    # 广义最小二乘法
    xvx = 0
    xvy = 0
    d_inv = np.linalg.inv(d)
    for i in range(group_size):
        x = group_x[i]
        y = group_y[i]
        z = group_z[i]
        wdv = woodbury(z, d_inv)
        xvx += x.T.dot(x - wdv.dot(x))
        xvy += x.T.dot(y - wdv.dot(y))
    return np.linalg.solve(xvx, xvy)


def woodbury(z, d_inv):
    # woodbury恒等式
    dzzz = np.linalg.solve(d_inv + z.T.dot(z), z.T)
    zdzzz = z.dot(dzzz)
    return zdzzz


def tril_to_square(tril, square_size):
    # 向量转换为方阵
    square = np.zeros((square_size, square_size))
    square[np.tril_indices_from(square)] = tril
    square[np.triu_indices_from(square)] = tril
    return square


def lmm(y, x, z, groups, is_reml=True):
    """
    线性混合模型(多层线性模型)
    :param y: 因变量
    :param x: 固定效应自变量
    :param z: 随机效应自变量
    :param groups: 组
    :param is_reml: 是否reml
    """
    group_dt = {}
    n = 0
    for i, group in enumerate(groups):
        n += 1
        if group in group_dt:
            group_dt[group].append(i)
        else:
            group_dt[group] = [i]
    group_size = 0
    group_y = []
    group_x = []
    group_z = []
    for val in group_dt.values():
        group_size += 1
        group_y.append(y[val])
        group_x.append(x[val])
        group_z.append(z[val])
    p = group_x[0].shape[1]
    q = group_z[0].shape[1]
    d = np.eye(q)
    tril_d = d[np.tril_indices_from(d)]
    beta = gls(group_y, group_x, group_z, group_size, d)
    res = minimize(fun=lambda *x: -loglik(*x) / n,
                   x0=tril_d,
                   args=(group_y, group_x, group_z, group_size, beta, n, is_reml),
                   jac=jac,
                   method='BFGS')
    d = tril_to_square(res.x, q)
    beta = gls(group_y, group_x, group_z, group_size, d)
    sigma_square = _get_sigma_square(beta, d, group_size, group_x, group_y, group_z, n, is_reml)
    negative_hess = negative_hessian(d, group_y, group_x, group_z, group_size, beta, n, is_reml)
    std_error = np.diag(np.linalg.inv(negative_hess)) ** 0.5
    print('==固定系数==')
    print(beta)
    print('==随机系数方差协方差矩阵==')
    print(sigma_square * d)
    print('==固定系数标准误')
    print(std_error[:p])
    print('==随机系数方差协方差矩阵标准误==')
    print(std_error[p:] * sigma_square ** 0.5)


def _get_sigma_square(beta, d, group_size, group_x, group_y, group_z, n, is_reml):
    # 计算sigma^2参数
    rvr = 0
    d_inv = np.linalg.inv(d)
    for i in range(group_size):
        x = group_x[i]
        r = group_y[i] - x.dot(beta)
        z = group_z[i]
        wdv = woodbury(z, d_inv)
        rvr += r.T.dot(r - wdv.dot(r))
    if is_reml:
        p = group_x[0].shape[1]
        n = n - p
    sigma_square = rvr / n
    return sigma_square


def loglik(tril_d, group_y, group_x, group_z, group_size, beta, n, is_reml):
    # 对数似然函数
    q = group_z[0].shape[1]
    d = tril_to_square(tril_d, q)
    d_inv = np.linalg.inv(d)
    loglik_val = 0
    logv = 0
    rtv_r = 0
    xtv_x = 0
    for i in range(group_size):
        x = group_x[i]
        r = group_y[i] - x.dot(beta)
        z = group_z[i]
        wdv = woodbury(z, d_inv)
        rtv_r += r.T.dot(r - wdv.dot(r))
        v = z.dot(d).dot(z.T)
        v += np.eye(len(v))
        logv += np.log(np.linalg.det(v))
        if is_reml:
            xtv_x += x.T.dot(x - wdv.dot(x))
    if is_reml:
        p = group_x[0].shape[1]
        n = n - p
        loglik_val -= np.log(np.linalg.det(xtv_x))
    loglik_val -= logv + n * np.log(rtv_r)
    loglik_val += n * np.log(n) - n - n * np.log(2 * np.pi)
    loglik_val /= 2
    return loglik_val[0][0]


def jac(tril_d, group_y, group_x, group_z, group_size, beta, n, is_reml):
    # 雅克比矩阵
    q = group_z[0].shape[1]
    d = tril_to_square(tril_d, q)
    d_inv = np.linalg.inv(d)
    zvz = chc = kk = rvr = xvx = 0
    zvx_list = []
    for i in range(group_size):
        x = group_x[i]
        r = group_y[i] - x.dot(beta)
        z = group_z[i]
        wdv = woodbury(z, d_inv)
        vr = r - wdv.dot(r)
        rvr += r.T.dot(vr)
        ztv_r = z.T.dot(vr)
        zvz += z.T.dot(z - wdv.dot(z))
        vx = x - wdv.dot(x)
        zvx = z.T.dot(vx)
        xvx += x.T.dot(vx)
        kk -= ztv_r.dot(ztv_r.T)
        if is_reml:
            zvx_list.append(zvx)

    if is_reml:
        for zvx in zvx_list:
            chc -= zvx.dot(np.linalg.solve(xvx, zvx.T))

    _n = n
    _jac_val = zvz
    if is_reml:
        p = group_x[0].shape[1]
        _n = _n - p
        _jac_val += chc
    _jac_val += (_n * kk / rvr)
    _jac_val /= n
    diag_idx = range(q)
    _jac_val[diag_idx, diag_idx] /= 2
    return _jac_val[np.tril_indices_from(_jac_val)]

黑塞矩阵计算代码如下

def negative_hessian(d, group_y, group_x, group_z, group_size, beta, n, is_reml):
    # 负黑塞矩阵
    d_inv = np.linalg.inv(d)
    rvr = b = kk_kk = zvz_zvz = h = rvz = xvx = ck_kc = 0
    zvr_list = []
    zvx_list = []
    xvr_list = []
    zvz_list = []
    for i in range(group_size):
        x = group_x[i]
        r = group_y[i] - x.dot(beta)
        z = group_z[i]
        wdv = woodbury(z, d_inv)
        vr = r - wdv.dot(r)
        rvr += r.T.dot(vr)
        vz = z - wdv.dot(z)
        rvz += r.T.dot(vz)
        zvz = z.T.dot(vz)
        zvz_list.append(zvz)
        zvr = z.T.dot(vr)
        zvr_list.append(zvr)
        zvr_zvr = zvr.dot(zvr.T)
        b += np.kron(zvz, zvr_zvr) + np.kron(zvr_zvr, zvz)
        zvz_zvz -= np.kron(zvz, zvz.T)
        vx = x - wdv.dot(x)
        h += x.T.dot(vx)
        xvx += x.T.dot(vx)
        zvx = z.T.dot(vx)
        zvx_list.append(zvx)
        ck_kc += np.kron(zvr.T, zvx.T) + np.kron(zvx.T, zvr.T)
        xvr_list.append(x.T.dot(vr))

    logxvx_dd = 0
    kgk_gkk = 0
    gg_gg = 0
    for i in range(group_size):
        zvr = zvr_list[i]
        zvx = zvx_list[i]
        zvz = zvz_list[i]
        zvxhzvx = zvx.dot(np.linalg.solve(h, zvx.T))
        xvr = xvr_list[i]
        if is_reml:
            logxvx_dd += np.kron(zvxhzvx, zvz)
            logxvx_dd += np.kron(zvz, zvxhzvx)
        for j in range(group_size):
            zvr_j = zvr_list[j]
            zvx_j = zvx_list[j]
            xvr_j = xvr_list[j]
            zvr_zvr = zvr.dot(zvr_j.T)
            kk_kk -= np.kron(zvr_zvr, zvr_zvr)
            hzvx = np.linalg.solve(h, zvx_j.T)
            xvr_zvr = xvr_j.dot(zvr.T)
            kgk_gkk -= np.kron(zvr.T, xvr_zvr) + np.kron(xvr_zvr, zvr.T)
            gg_gg -= np.kron(xvr_j, xvr.T) + xvr_j.dot(xvr.T)
            if is_reml:
                logxvx_dd -= np.kron(zvx.dot(hzvx), zvx.dot(hzvx))

    p = group_x[0].shape[1]
    q = group_z[0].shape[1]
    rand_fix_hess = 0.5 * (n - p) * (ck_kc / rvr + kgk_gkk / rvr ** 2)
    rand_hess = 0.5 * zvz_zvz + 0.5 * (n - p) * (b / rvr + kk_kk / (rvr ** 2))
    if is_reml:
        rand_hess += 0.5 * logxvx_dd
    fix_hess = 0.5 * (n - p) * (2 * xvx / rvr + 2 * gg_gg /rvr ** 2)
    negative_hess = np.zeros((p + q * (q + 1) // 2, p + q * (q + 1) // 2))
    delete_rows = []
    delete_cols = []
    for i in range(q):
        for j in range(i):
            rand_hess[i * q + j, :] += rand_hess[j * q + i, :]
            rand_hess[:, i * q + j] += rand_hess[:, j * q + i]
            rand_fix_hess[:, i * q + j] += rand_fix_hess[:, j * q + i]
            delete_rows.append(j * q + i)
            delete_cols.append(j * q + i)
    rand_hess = np.delete(rand_hess, delete_rows, axis=0)
    rand_hess = np.delete(rand_hess, delete_cols, axis=1)
    rand_fix_hess = np.delete(rand_fix_hess, delete_cols, axis=1)
    negative_hess[:p, :p] = fix_hess
    negative_hess[:p, p:] = rand_fix_hess
    negative_hess[p:, :p] = rand_fix_hess.T
    negative_hess[p:, p:] = rand_hess
    return negative_hess

为了确保梯度和黑塞矩阵计算的正确性，我们还用数值解验证了解析解的结果，例如依据导数的定义

，我们以求

数值导数为例，示例代码如下

def f(x):
    return x ** 2


def der(x, delta=1e-8):
    # 求导
    return (f(x + delta) - f(x)) / (delta)


print(der(8))

写计算代码的时候，由于公式的复杂度和命名的复杂度，往往会一不留神犯错，所以用数值方法检验解析方法的正确性是非常好的习惯。

运行我们的lmm参数估计程序

我们以lme4的sleepstudy数据验证我们的程序，数据格式形如下图，其中第一列是因变量，第二列是固定系数和随机系数的协变量，第三列是组数据

sleep = np.loadtxt('sleepstudy.csv', delimiter=',')
# 因变量
y_data = sleep[:, 0, np.newaxis]
# 固定系数的协变量，并加截距
fixed_data = sleep[:, 1]
fixed_data = np.concatenate((np.ones((len(y_data), 1)), fixed_data[:, np.newaxis]), axis=1)
# 随机系数的协变量
random_data = fixed_data
# 组数据
groups_data = sleep[:, 2]
# 求解
lmm(y_data, fixed_data, random_data, groups_data, is_reml=False)

打印结果如下

==固定系数==
[[251.40510485]
 [ 10.46728596]]
==随机系数方差协方差矩阵==
[[565.51604604  11.05536609]
 [ 11.05536609  32.68214786]]
==固定系数标准误
[6.66943562 1.51065183]
==随机系数方差协方差矩阵标准误==
[[11.07071819  1.69007625  0.5679701 ]]

扩展

多随机效应模型（多水平多层线性模型）

形如

，

，

，则

，其他求解过程与上面内容一致

广义线性混合模型

广义线性混合模型的因变量可以是正态分布，也可以是二项分布、泊松分布等等，求解时会遇到难以处理的积分，一般用变分、拉普拉斯，自适应积分等方法处理。

全文结束