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

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

#### 2.3 写出似然函数

• 公式1： 常用于GWAS软件中，比如TASSEL、GCTA、GEMMA、FaST-LMM、EMMAX
• 公式2： 常用于传统遗传评估软件，比如ASREML，DMU，BLUPF90

### 3. 一个示例

```# 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```

sigmaE

sigmaS

k

925.6148

91.78376

10.0847344

897.9372

108.19047

8.2995954

863.6173

129.55950

6.6657968

820.9699

157.69430

5.2060850

768.2845

194.87483

3.9424508

704.4117

243.56957

2.8920351

629.8507

305.53381

2.0614764

548.0643

380.16717

1.4416402

465.8970

462.97884

1.0063030

391.7085

546.12170

0.7172550

331.7683

621.64909

0.5336907

287.8838

684.73480

0.4204312

258.0564

734.15339

0.3515020

238.7261

770.93255

0.3096589

226.5165

797.12077

0.2841684

218.8922

815.07790

0.2685538

214.1497

827.02910

0.2589385

211.2014

834.81106

0.2529931

209.3676

839.80212

0.2493058

208.2261

842.97109

0.2470145

207.5152

844.97003

0.2455888

207.0721

846.22564

0.2447009

206.7960

847.01226

0.2441476

206.6238

847.50423

0.2438027

206.5165

847.81160

0.2435877

206.4496

848.00351

0.2434537

206.4078

848.12328

0.2433701

206.3818

848.19801

0.2433179

206.3655

848.24463

0.2432854

206.3554

848.27371

0.2432651

206.3491

848.29185

0.2432525

206.3452

848.30317

0.2432446

206.3427

848.31022

0.2432397

206.3412

848.31462

0.2432366

206.3402

848.31737

0.2432347

206.3396

848.31908

0.2432335

206.3392

848.32015

0.2432328

206.3390

848.32081

0.2432323

206.3389

848.32123

0.2432320

206.3388

848.32149

0.2432318

206.3387

848.32165

0.2432317

206.3387

848.32175

0.2432316

206.3387

848.32181

0.2432316

206.3387

848.32185

0.2432316

206.3386

848.32188

0.2432315

206.3386

848.32189

0.2432315

206.3386

848.32190

0.2432315

### 4. 注意

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

20.0855369231877

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

3

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

0 条评论

• ### 基因组选择和SNP分析在ASREML-SA中的实现方法

asreml是非常强大的软件, 由于太强大, 很多人不会使用. 基因组选择在育种中的应用, 其基础是常规的系谱动物模型, 动物模型也可以很复杂, 看一下asre...

• ### 怕死秃头工作站？？？

我仔细看了一下，原来Rstudio新建的脚本，默认的是R脚本，如果想要运行Python，需要新建Python脚本。

• ### R语言读取 xlsx 和xls 文件

xlsx文件，是2007,2013,2016版本的Excel文件，R语言中可以通过openxlsx包进行读取。

• ### Docker+Jenkins持续集成环境（1）使用Docker搭建Jenkins+Docker持续集成环境

本文介绍如何通过Jenkins的docker镜像从零开始构建一个基于docker镜像的持续集成环境，包含自动化构建、发布到仓库\并部署上线。 0. 前置条件 服...

• ### libuv源码分析之unix域

unix域是一种基于单主机的进程间通信方式。实现模式类似tcp通信。今天先分析他的实现，后续会分析他的使用。在libuv中，unix域用uv_pipe_t表示。

• ### Seurat 4.0 || 分析scRNA和膜蛋白数据

得益于单细胞多模态技术的发展，允许我们在单细胞水平从不同侧面考察细胞状态，如CITE-seq技术可以同时对单细胞转录组和膜蛋白进行定量。单细胞多模态数据分析面临...

• ### Spark shuffle详细过程

有许多场景下，我们需要进行跨服务器的数据整合，比如两个表之间，通过Id进行join操作，你必须确保所有具有相同id的数据整合到相同的块文件中。那么我们先说一下m...

• ### 常见问题: MongoDB 存储

存储引擎是数据库的一部分，负责管理如何在内存和磁盘上存储数据。许多数据库支持多个存储引擎，其中不同的引擎对特定工作负载的性能会更好。例如，一个存储引擎可能为读取...

• ### 扔掉 Postman 吧，试试 Postwoman 高能神器！

个人觉得，运行在浏览器端这一点比较实用，毕竟我们都是 Web 开发人员，浏览器跨平台的便利性早已深入人心。无需安装，随时随地可以测试接口。按照作者自己的说法，他...

• ### 我决定用 Postwoman 替代 Postman

个人觉得，运行在浏览器端这一点比较实用，毕竟我们都是 Web 开发人员，浏览器跨平台的便利性早已深入人心。无需安装，随时随地可以测试接口。按照作者自己的说法，他...