BLUP

1、混合线性模型一般形式及含义

y=Xβ+Zu+e

n为观测值的个数，q为个体数目
y为观测值向量（n×1）
β为固定效应向量(p×1)，X为固定效应的设计矩阵(n×p)
u为个体育种值（随机效应）向量(q×1),Z为育种值得设计矩阵(n×q)
e为残差向量（n×1）

E(y)=Xβ,E(u)=0,E(e)=0

?E(y)=Xβ,为n×1向量，每个个体观测值分布的期望为每个个体的固定效应的值
E(u)=0,为q×1 0向量，即每个个体育种值分布的期望为0
?E(e)=0,为n×1 0向量，即每个个体的残差分布的期望为0
? 每个个体虽然只有一两个观测值，但观测值服从某一个分布，同样育种值也服从一个分布，最后会产生一个残差分布。

Var(u)=G,Var(e)=R,Cov(u,e′)=0

混合模型中，u和e服从m×p维正态分布，即

u∼N(0,G),e∼N(0,R)

u∼N(0,G),0为q×1向量；G为q×q的方阵，对角线为个体育种值分布的方差，非对角线不同个体育种值的协方差
e∼N(0,R),0为n×1向量；R为n×n的方阵，对角线为方差，非对角线为协方差

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢y11y12y13y21y22y31⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢000110111001⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥[b1b2]+⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢111000000110000001⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢a1a2a3⎤⎦⎥+⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢e11e12e13e21e22e31⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥

2、观测值和育种值的联合概率密度

多维正态分布的一般形式：

y∼Np(u,S)

f(y)=1(2π)p2|S|12⋅exp(−12(y−u)′S−1(y−u))

Y=[y1y2]∼N2([μ1μ2],[σ21ρρσ22])

f(y1,y2|μ,Σ)=(2π1−ρ2−−−−−√)−1exp[12(1−ρ2){(y1−μ1σ1)2−2ρ(y1−μ1σ1)(y2−μ2σ2)+(y2−μ2σ2)2}]

y1,y2为随机变量，μ1,μ2,σ1,σ2,ρ均为某一个固定的数值

产生多维正态分布R语言代码：

library(MASS) Sigma <- matrix(c(10,3,3,2),2,2) Sigma x=mvrnorm(n=1000, rep(0, 2), Sigma)#x为二维正态分布

y和u的联合密度函数为：

f(y,u)=f1(y|u)⋅f2(u)

根据多维正态分布的概率密度公式, f1(y|u)如下：

y|u=y−Xβ−Zu=e∼N(0,R)f1(y|u)=C1⋅exp{−12(y−Xβ−Zu)′ R−1 (y−Xβ−Zu)}C1=(2π)−m2|R|−12

根据多维正态分布的概率密度公式, f2(u)如下

u∼N(0,G)f2(u)=C2⋅exp{−12u′G−1u}C2=(2π)−m2|G|−12

所以联合概率密度f(y,u)为：

f(y,u)=C⋅exp{−12(y−Xβ−Zu)′ R−1 (y−Xβ−Zu)−12u′G−1u}

3、对联合概率密度求解

求此函数关于β和u的极大值，因此令偏导等于0，可得：

∂f(y,u)∂β∂f(y,u)∂u==C⋅exp{a}(X′R−1y−X′R−1Xβ−X′R−1Zu)=0C⋅exp{a}(Z′R−1y−Z′R−1Xβ−Z′R−1Zu−G−1u)=0

其中的a为f(y,u)中的指数函数的幂。上式整理后得

X′R−1Xβ^+X′R−1Zu^=X′R−1yZ′R−1Xβ^+(Z′R−1Z+G−1)u^=Z′R−1y

其矩阵形式为：

[X′R−1XZ′R−1XX′R−1ZZ′R−1Z+G−1][β^u^]=[X′R−1yZ′R−1y]

此方程组称为混合模型方程组（mixed model equations，MME）。

http://wenku.baidu.com/view/ce677909f78a6529647d539b