方差组分估计之约束最大似然

推荐三本书。第一本是Linear Models for the Prediction of Animal Breeding Values，这本书可以帮助了解线性模型框架，不过这本书没有详述很多公式的具体来源，方差组分估计一块也介绍的过于简略。第二本是Improved Iterative Schemesfor REML Estimation of Variance Parameters in Linear Mixed Models，严格来说，这是一个博士生的毕业论文，方差组分一块阐述的比较详细。第三本是MATRIX ANALYSIS AND APPLICATIONS，矩阵性质很好的工具书，老教授写的，有中文版，可能想扩大影响力，最近出了英文版。回复Books可以获得三本书，有条件的可以支持一下原版。

-------------------------------------

约束最大似然（Restricted Maximum Likelihood, REML），实质还是最大似然，约束的是固定效应部分，避免了由于模型拟合固定效应而造成的自由度损失。说实话，自由度真是个很复杂的问题，具体怎么损失的，损失后的影响，没细究过。不过结论还是非常明确的，约束最大似然是方差组分估计的标准方法。

模型

其中，y是表型向量，b是固定效应向量，u是随机效应向量，e是随机残差。X和Z是设计矩阵。

假设u和e服从多元正态分布，即

需要注意的是，与以前更新的文章不同，在这里，我们用G表示随机效应的方差协方差矩阵，K表示基因组关系矩阵(kinship)，I是单位阵，和分别是加性方差和随机残差。

约束最大似然

定义L矩阵（nx(n-p)维），满足，这样

其中，n是个体数，p为矩阵X的秩。

原模型转化为

L’y的分布如下

L’y的密度函数

对数似然函数

经过一系列的推导（这句话说的比较敷衍，主要是推导方法不止一种，而且比较复杂，等我都理清楚了专门更新一篇）

结合《方差组分估计之最大似然》和公式1可以看出，约束最大似然相对于最大似然，多了一项。

需要注意的是，这两个公式是线性混合模型约束最大似然函数的通用公式，也就是说，动物模型、重复力模型、多性状模型、随机回归模型等等，公式都是他们，只是方差结构的形式不同而已。

方差组分估计的MME-based REML和directREML

通过最小化公式1或公式2（也就是最大化约束似然函数），我们即可求出方差组分。那么问题来了，这两个公式有什么区别？这就涉及到MME-based REML和direct REML问题了。最早提出这个问题的是Lee and van der Werf (2006)，文章作者之一，Julius van der Werf目前是GSE杂志的chief editor。

最直接的公式，应该是公式1，但是最早应用的公式，却是公式2。两公式区别

a.公式1含V，公式2含C

b.公式1利用的是K，不需要对K求逆，而公式2利用的是K逆

c.传统线性混合模型软件，例如DMU、ASreml、blupf90等等都是利用的公式2，而目前流行的GWAS软件，GCTA、GEMMA、FaST-LMM、EMMAX、TASSEL等等，都是用的公式1。

具体利用哪个公式，主要是考虑矩阵的稀疏性和维数。传统遗传评估、目前流行的一步法基因组选择，利用到系谱关系，A逆和H逆比较稀疏，构建的C有很好的稀疏性，所以利用公式2速度快。而GWAS一般只利用SNP信息构建基因组关系矩阵（K），K是密集矩阵，构建出C后也是密集矩阵，并且C矩阵的维数要比V大，反而会降低速度。导致流行的GWAS软件，普遍利用公式1。

direct REML的Python实现

与前面更新《方差组分估计之最大似然》类似，我们先利用Python自带的求极值函数估计方差组分，常用的迭代估计方差组分的方式，后续介绍。

# -*- coding: UTF-8 -*-

from scipy.optimize import minimize

import numpy as np

from pysnptools.snpreader import Bed

from scipy import linalg

import pandas as pd

###构建G矩阵###

bed_file = 'plink'

M = Bed(bed_file,count_A1=False).read()

freq = np.sum(M.val, axis=0) / (2 *M.iid_count)

scale = np.sum(2 * freq * (1 -freq))

Z = M.val - 2 * freq

K = np.dot(Z, Z.T) / scale

# 读取数据文件，生成设计矩阵

data_file = 'phe'

data = pd.read_table(data_file,header=0, sep='\s+')

print data.head()

y = np.array(data['phe'],dtype=np.float)

X = np.array(data.loc[:,'mean':'treat'], dtype=np.float)

num = len(y)

Z = np.eye(num)

# 定义似然函数，var存储的加性方差和残差的列表

def llf(var):

V = K * var[0] + np.eye(num) * var[1]

VX = np.dot(Vinv, X)

XVX = np.dot(X.T, VX)

XVXinv = linalg.inv(XVX)

P = Vinv - np.dot(np.dot(VX, XVXinv), VX.T)

yPy = np.dot(np.dot(y.T, P), y)

ll = logdetV + yPy + logdetXVX

return ll

# 调用函数求极小值

result = minimize(llf, [1.0, 1.0],method='L-BFGS-B', bounds=[[10e-10, 10e10], [10e-10, 10e10]])

print result

结果展示

------------------------------------------

回复Qgenetics获取本次更新的数据和代码

参考文献

张沅，张勤. 畜禽育种中的线性模型

陈希孺，倪国熙.数量统计学教程

EmmaKnight. Improved Iterative Schemes for REML Estimation of Variance Parametersin Linear Mixed Models

Lee, S.H. and vander Werf, J.H. An efficient variance component approach implementing an averageinformation REML suitable for combined LD and linkage mapping with a generalcomplex pedigree.Genetics, selection,evolution : GSE2006;38(1):25-43.