1，直接计算
数据
试验数据
模型
Model
- X是固定效应矩阵，b是Herd
- Z是随机效应矩阵，u是Sire
- e是残差
矩阵的写法：
Matrix
计算公式及结果
R语言代码
> # data
> y=c(110,100,110,100,100,110,110,100,100)
> X = matrix(c(1,1,0,0,0,0,0,0,0,
+ 0,0,1,1,1,0,0,0,0,
+ 0,0,0,0,0,1,1,1,1),9, byrow=F)
> X
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 1 0 0
[3,] 0 1 0
[4,] 0 1 0
[5,] 0 1 0
[6,] 0 0 1
[7,] 0 0 1
[8,] 0 0 1
[9,] 0 0 1
> Z = matrix(c(1,0,0,0,0,0,0,0,0,
+ 0,0,1,0,0,0,0,0,0,
+ 0,0,0,0,0,1,1,0,0,
+ 0,1,0,1,1,0,0,1,1),9, byrow=F)
> Z
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 0 0 1
[3,] 0 1 0 0
[4,] 0 0 0 1
[5,] 0 0 0 1
[6,] 0 0 1 0
[7,] 0 0 1 0
[8,] 0 0 0 1
[9,] 0 0 0 1
> I1=diag(4)
> I2=diag(9)
> # 方差组分固定值
> se=1
> su=0.1
> G=I1*su
> R=I2*se
>
> V = Z%*%G%*%t(Z) + R
> V
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
[2,] 0.0 1.1 0.0 0.1 0.1 0.0 0.0 0.1 0.1
[3,] 0.0 0.0 1.1 0.0 0.0 0.0 0.0 0.0 0.0
[4,] 0.0 0.1 0.0 1.1 0.1 0.0 0.0 0.1 0.1
[5,] 0.0 0.1 0.0 0.1 1.1 0.0 0.0 0.1 0.1
[6,] 0.0 0.0 0.0 0.0 0.0 1.1 0.1 0.0 0.0
[7,] 0.0 0.0 0.0 0.0 0.0 0.1 1.1 0.0 0.0
[8,] 0.0 0.1 0.0 0.1 0.1 0.0 0.0 1.1 0.1
[9,] 0.0 0.1 0.0 0.1 0.1 0.0 0.0 0.1 1.1
>
> Vinv <- solve(V)
>
> blue <- solve(t(X)%*%Vinv%*%X)%*%t(X)%*%Vinv%*%y
> blue
[,1]
[1,] 105.6387
[2,] 104.2757
[3,] 105.4584
>
> blup <- G%*%t(Z)%*%Vinv%*%(y - X%*%blue)
> blup
[,1]
[1,] 0.3964857
[2,] 0.5203875
[3,] 0.7569272
[4,] -1.6738004
2. 用混合线性方程组计算
公式计算
转化为：
result
BLUE值和BLUP值的计算
R语言解决方案
左边公式计算：
LHS
alpha <- se/su
XpX=crossprod(X) #X’X
XpZ=crossprod(X,Z) #X’Z
ZpX=crossprod(Z,X) #Z’X
ZpZ=crossprod(Z) #Z’Z
Xpy=crossprod(X,y) #X’y
Zpy=crossprod(Z,y) #Z'y
LHS=rbind(cbind(XpX,XpZ),cbind(ZpX,ZpZ+diag(4)*alpha)) #LHS
> LHS
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 2 0 0 1 0 0 1
[2,] 0 3 0 0 1 0 2
[3,] 0 0 4 0 0 2 2
[4,] 1 0 0 11 0 0 0
[5,] 0 1 0 0 11 0 0
[6,] 0 0 2 0 0 12 0
[7,] 1 2 2 0 0 0 15
右边计算
RHS
> RHS=rbind(Xpy,Zpy) #RHS
> RHS
[,1]
[1,] 210
[2,] 310
[3,] 420
[4,] 110
[5,] 110
[6,] 220
[7,] 500
计算结果
> sol=solve(LHS)%*%RHS #
> sol
[,1]
[1,] 105.6386574
[2,] 104.2757378
[3,] 105.4584366
[4,] 0.3964857
[5,] 0.5203875
[6,] 0.7569272
[7,] -1.6738004
asreml 的计算结果
代码如下：
> Ve = 1; Va = 0.1
> names(Ve) <- c("F")
> names(Va) <- c("F")
> dat$Herd <- as.factor(dat$Herd)
> model <- asreml(y ~ Herd-1, random = ~ idv(Sire,init=Va),
+ family=asreml.gaussian(dispersion = Ve), data=dat)
ASReml: Wed May 24 17:15:56 2017
LogLik S2 DF wall cpu
-85.4750 1.0000 6 17:15:56 0.0
-85.4750 1.0000 6 17:15:56 0.0
-85.4750 1.0000 6 17:15:56 0.0
-85.4750 1.0000 6 17:15:56 0.0
Finished on: Wed May 24 17:15:56 2017
LogLikelihood Converged
> summary(model)$varcomp
gamma component std.error z.ratio constraint
Sire!Sire.var 0.1 0.1 NA NA Fixed
R!variance 1.0 1.0 NA NA Fixed
> coef(model)
$fixed
effect
Herd_1 105.6387
Herd_2 104.2757
Herd_3 105.4584
$random
effect
Sire_AD -1.6738004
Sire_BB 0.5203875
Sire_CC 0.7569272
Sire_ZA 0.3964857
$sparse
NULL
> sol
[,1]
[1,] 105.6386574
[2,] 104.2757378
[3,] 105.4584366
[4,] 0.3964857
[5,] 0.5203875
[6,] 0.7569272
[7,] -1.6738004
