线性混合模型系列五：REML实战

1. 混合线性模型

方程：

假定

2. 混合线性模型的似然函数

2.1 混合模型中y的分布

2.2 对V进行公式代换

2.3 写出似然函数

宁超在他的公众号“Pythn与数量遗传学”的“方差组分估计之约束最大似然”文章中，给出了下面两种计算公式，公式一是直接似然函数（direct REML），公式二是间接的似然函数（MME based REML）。

• 公式1：

  常用于GWAS软件中，比如TASSEL、GCTA、GEMMA、FaST-LMM、EMMAX

• 公式2：

  常用于传统遗传评估软件，比如ASREML，DMU，BLUPF90

  
  

3.  一个示例

下面这个示例是张勤老师在2019年统计遗传学暑假学校教材的数据

# zhangqin laoshi PPT codeherd <- factor(c(1,2,2,1,1,2,1,2,2))
sire <- factor(c(1,1,1,2,2,3,4,4,4))
y <- c(240,190,170,180,200,140,170,100,130)
REML_data <- data.frame(herd,sire,y)
X <- model.matrix(y~herd-1)
Z <- model.matrix(y~sire-1)
A <- matrix(c(1,0.25,0,0,0.25,1,0,0,0,0,1,0,0,0,0,1),nr=4)

3.1 公式1

书写似然函数

mixmodel.loglik<-function(theta,y){
  sigmag2<-theta[1]
  sigmae2<-theta[2]

  n = length(y)
  G = A*sigmag2
  R = diag(n)*sigmae2
  V = Z %*% G %*% t(Z) + R
  Vi = solve(V)

  P = Vi - Vi %*% X %*% solve(t(X) %*% Vi %*% X) %*% t(X) %*% Vi
  logl = -0.5*(log(det(V)) + log(det(t(X) %*% Vi %*% X)) + t(y) %*% P %*% y)  return(-logl)
}

optim(c(1,1),mixmodel.loglik,y=y)

3.2 公式2

mixmodel.loglik2 <-function(theta,y){
  sigmag2<-theta[1]
  sigmae2<-theta[2]

  n = length(y)
  G = A*sigmag2
  R = diag(n)*sigmae2
  V = Z %*% G %*% t(Z) + R
  Vi = solve(V)
  P = Vi - Vi %*% X %*% solve(t(X) %*% Vi %*% X) %*% t(X) %*% Vi
  C11 = t(X) %*%  solve(R) %*% X
  C12 = t(X) %*%  solve(R) %*% Z
  C21 = t(Z) %*%  solve(R) %*% X
  C22 = t(Z) %*%  solve(R) %*% Z + solve(G)
  C = rbind(cbind(C11,C12),cbind(C21,C22))
  logl2 = -0.5*(log(det(C)) + log(det(R)) + log(det(G)) + t(y) %*% P %*% y)  return(-logl2)
}
optim(c(1,1),mixmodel.loglik2,y=y)

两者结果一样。

3.3 张勤老师PPT上的EM算法

herd <- factor(c(1,2,2,1,1,2,1,2,2))
sire <- factor(c(1,1,1,2,2,3,4,4,4))
y <- c(240,190,170,180,200,140,170,100,130)
REML_data <- data.frame(herd,sire,y)
X <- model.matrix(y~herd-1)
Z <- model.matrix(y~sire-1)
XX <- crossprod(X)
XZ <- crossprod(X,Z)
ZX <- t(XZ)
ZZ <- crossprod(Z)
LHS <- rbind(cbind(XX,XZ),cbind(ZX,ZZ))
Xy <- crossprod(X,y)
Zy <- crossprod(Z,y)
RHS <- rbind(Xy,Zy)
A <- matrix(c(1,0.25,0,0,0.25,1,0,0,0,0,1,0,0,0,0,1),nr=4)
Ainv <- solve(A)
yy <- crossprod(y)
N <- length(y)
rankX <-qr(X)$rank
q <- length(unique(sire))
nh <- length(unique(herd))
k0 <- 12threshold <- 0.00000001write.table(t(c("sigmaE","sigmaS","k")),file="results.txt",row.names = FALSE, col.names = FALSE)repeat{
    LHS1 <- LHS
    k <- k0
    LHS1[(nh+1):(nh+q), (nh+1):(nh+q)] <- LHS1[(nh+1):(nh+q), (nh+1):(nh+q)]+Ainv*k
    sol <- solve(LHS1,RHS)
  C <- solve(LHS1)
  Czz <- C[(nh+1):(nh+q), (nh+1):(nh+q)]
    sigmaE <- (yy - crossprod(sol,RHS))/(N-rankX)
    sAs <- t(sol[(nh+1):(nh+q)]) %*% Ainv %*% sol[(nh+1):(nh+q)]
    trCA <- sum(diag(Czz %*% Ainv))
    sigmaS <- (sAs + trCA*sigmaE)/q
    k0 <- as.numeric(sigmaE/sigmaS)
    out <- c(sigmaE,sigmaS,k0)
    write.table(t(out), file="results.txt", append = TRUE, row.names = FALSE, col.names = FALSE)
    if(abs(k-k0) < threshold) break}
results <- read.table("results.txt",header = TRUE)
results

最终的结果是，加性方差组分为206.3386，残差的方差组分为848.3219，结果一致

4. 注意

公式中的log，也可以写为ln，是自然对数，在R中log默认的就是自然对数

# 自然对数的3次方exp(3)

20.0855369231877

# 对上面结果求自然对数log(exp(3))

3

在R中，如果想计算10的对数，用函数log10(x)

log(10)

2.30258509299405

log10(10)

1

5. 参考文献

http://people.math.aau.dk/~rw/Undervisning/Mat6F12/Handouts/2.hand.pdf
http://www2.stat.duke.edu/~sayan/Sta613/2018/lec/LMM.pdf
张勤老师 《统计遗传学暑期学校教材》 中国农业大学 2019.07
宁超公众号：Python与数量遗传学 方差组分估计之约束最大似然

之前系列：

线性混合模型系列一：基本定义

线性混合模型系列二：模型假定

线性混合模型系列三：似然函数

线性混合模型系列四：矩阵求解

混合线性模型的其它博文：

混合线性模型如何检测固定因子和随机因子的显著性以及计算R2

混合线性模型如何进行多重比较

如何检测遗传相关的显著性：LRT检验操作方法