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

• 2019 年 12 月 5 日
• 筆記

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与数量遗传学 方差组分估计之约束最大似然