R语言求解混合线性方程组（无系谱）

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