# 1，直接计算

• X是固定效应矩阵，b是Herd
• Z是随机效应矩阵，u是Sire
• e是残差
矩阵的写法：

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. 用混合线性方程组计算

BLUE值和BLUP值的计算

R语言解决方案

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=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

https://wap.sciencenet.cn/blog-2577109-1132837.html

## 相关博文

GMT+8, 2021-12-9 15:31