BLUP

#1、混合线性模型一般形式及含义
$$y=X\beta +Zu+e$$

n为观测值的个数，q为个体数目,一个个体可以有多个观测值
$y为观测值向量（n\times 1）$
$\beta 为固定效应向量(p\times 1)，X为固定效应的设计矩阵(n\times p)$
$u为个体育种值（随机效应）向量(q\times 1),Z为育种值得设计矩阵(n\times q)$
$e为残差向量（n\times 1）$

$$E(y)=X\beta ,E(u)=0,E(e)=0$$

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

$$Var(u)=G,Var(e)=R,Cov(u,e^{{}^{\prime }})=0$$

$混合模型中，u和e服从m\times p维正态分布，即$
$$u\sim N(0,G),e\sim N(0,R)$$

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

$$\left[\begin{array}{c}{y}_{11}\\ y_{12}\\ y_{13}\\ y_{21}\\ y_{22}\\ y_{31}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ 0& 1\\ 0& 1\\ 1& 0\\ 1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\end{array}\right]+\left[\begin{array}{ccc}1& 0& 0\\ 1& 0& 0\\ 1& 0& 0\\ 0& 1& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right]+\left[\begin{array}{c}{e}_{11}\\ e_{12}\\ e_{13}\\ e_{21}\\ e_{22}\\ e_{31}\end{array}\right]$$

#2、观测值和育种值的联合概率密度
多维正态分布的一般形式：

$$y\sim N_p(\mathbf{u},\mathbf{S})$$
$$f(y)=\frac{1}{(2\pi {)}^{\frac{p}{2}}|\mathbf{S}{|}^{\frac{1}{2}}}\cdot exp\left(-\frac{1}{2}(y-\mathbf{u}{)}^{{}^{\prime }}{\mathbf{S}}^{-1}(y-\mathbf{u})\right)$$
$$Y=\left[\begin{array}{cc}{y}_1& y_2\end{array}\right]\sim N_2\left(\left[\begin{array}{cc}{\mu }_1& {\mu }_2\end{array}\right],\left[\begin{array}{cc}{\sigma }_1^2& \rho \\ \rho & {\sigma }_2^2\end{array}\right]\right)$$
$$f(y_1,y_2|\mu ,\Sigma )=(2\pi \sqrt{1-{\rho }^2}{)}^{-1}exp\left[\frac{1}{2(1-{\rho }^2)}\left\{{\left(\frac{y_1-{\mu }_1}{{\sigma }_1}\right)}^2-2\rho \left(\frac{y_1-{\mu }_1}{{\sigma }_1}\right)\left(\frac{y_2-{\mu }_2}{{\sigma }_2}\right)+{\left(\frac{y_2-{\mu }_2}{{\sigma }_2}\right)}^2\right\}\right]$$
$y_1,y_2为随机变量，{\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2,\rho 均为某一个固定的数值$

产生多维正态分布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)=f_1(y|u)\cdot f_2(u)$$
根据多维正态分布的概率密度公式,$f_1(y|u)$如下：

$$\begin{array}{rlr}& & y|u=y-X\beta -Zu=e\sim N(0,R)\\ & & f_1(y|u)=C_1\cdot exp\left\{-\frac{1}{2}(y-X\beta -Zu{)}^{{}^{\prime }}{R}^{-1}(y-X\beta -Zu)\right\}\\ & & C_1=(2\pi {)}^{-\frac{m}{2}}|R{|}^{-\frac{1}{2}}\end{array}$$

根据多维正态分布的概率密度公式,$f_2(u)$如下
$$\begin{array}{rlr}& & u\sim N(0,G)\\ & & f_2(u)=C_2\cdot exp\left\{-\frac{1}{2}{u}^{{}^{\prime }}{G}^{-1}u\right\}\\ & & C_2=(2\pi {)}^{-\frac{m}{2}}|G{|}^{-\frac{1}{2}}\end{array}$$

所以联合概率密度$f(y,u)$为：
$$f(y,u)=C\cdot exp\left\{-\frac{1}{2}(y-X\beta -Zu{)}^{{}^{\prime }}{R}^{-1}(y-X\beta -Zu)-\frac{1}{2}{u}^{{}^{\prime }}{G}^{-1}u\right\}$$

3、对联合概率密度求解

$求此函数关于\beta 和u的极大值，因此令偏导等于0，可得：$
∂ f ( y , u ) ∂ β = C ⋅ e x p { a } ( X ′ R − 1 y − X ′ R − 1 X β − X ′ R − 1 Z u ) = 0 ∂ f ( y , u ) ∂ u = C ⋅ e x p { a } ( Z ′ R − 1 y − Z ′ R − 1 X β − Z ′ R − 1 Z u − G − 1 u ) = 0 \begin{aligned} \frac{\partial f(y,u)}{\partial \beta} &=& C\cdot exp\{a\}(X^{'}R^{-1}y-X^{'}R^{-1}X\beta-X^{'}R^{-1}Zu) =0 \\ \frac{\partial f(y,u)}{\partial u} &=& C\cdot exp\{a\}(Z^{'}R^{-1}y-Z^{'}R^{-1}X\beta-Z^{'}R^{-1}Zu-G^{-1}u)=0 \end{aligned} ∂β∂f(y,u)​∂u∂f(y,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)中的指数函数的幂。上式整理后得$
$$\begin{array}{r}{X}^{{}^{\prime }}{R}^{-1}X\hat{\beta }+X^{{}^{\prime }}{R}^{-1}Z\hat{u}=X^{{}^{\prime }}{R}^{-1}y\\ Z^{{}^{\prime }}{R}^{-1}X\hat{\beta }+(Z^{{}^{\prime }}{R}^{-1}Z+G^{-1})\hat{u}=Z^{{}^{\prime }}{R}^{-1}y\end{array}$$

其矩阵形式为：
$$\left[\begin{array}{cc}{X}^{{}^{\prime }}{R}^{-1}X& X^{{}^{\prime }}{R}^{-1}Z\\ Z^{{}^{\prime }}{R}^{-1}X& Z^{{}^{\prime }}{R}^{-1}Z+G^{-1}\end{array}\right]\left[\begin{array}{c}\hat{\beta }\\ \hat{u}\end{array}\right]=\left[\begin{array}{c}{X}^{{}^{\prime }}{R}^{-1}y\\ Z^{{}^{\prime }}{R}^{-1}y\end{array}\right]$$
此方程组称为混合模型方程组（mixed model equations，MME）。

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