博文
AFEchidna包示例22之管道操作符%＞%的妙用
||
使用管道操作符`%>%`，可以直接把数据传递给下一个函数调用或表达式。目的有二：一是减少代码开发时间，提高代码的可读性和维护性；二是让代码更简短。
`%>%`是最常用的操作符，把左侧的数据或表达式，传递给右侧的函数调用或表达式运行，可以像个链条一样连续操作。
在使用AFEchidna包进行多个随机效应模型的分析时，管道操作符`%>%`的优势立马体现出来。
简单示例如下：
```R
// multiple G-structure models
library(AFEchidna)
res12<-echidna(fixed=cbind(h3,h4)~Trait+Trait:Rep,
random=~us(Trait):Fam,
residual=~units:us(Trait),
predict='Fam',mulT=TRUE,
qualifier = '!filter Spacing !select 1',
es0.file="fm.es0")
res12b <- update(res12, random=c(G1~us(Trait):Fam,
G2~diag(Trait):Fam,
G0~Trait:Fam),
batch.G=TRUE )
```
有时上述3个模型的随机效应差异较大，直接对`res12b`使用Var()来获取各模型的方差分量会出错：
```R
> Var(res12b)
Error in Var.esR(res12b) : 找不到对象'var1'
```
此时，可以考虑使用管道操作符`%>%`：
```R
> res12b %>% b2s %>% lapply(.,Var)
$G1
Term Sigma SE Z.ratio
1 Residual 1.00 0.000 Inf
2 us(Trait):Fam;us(Trait) 132.43 56.428 2.346885
3 us(Trait):Fam;us(Trait) 168.81 75.851 2.225547
4 us(Trait):Fam;us(Trait) 243.39 115.710 2.103448
5 units:us(Trait);us(Trait) 1589.10 100.300 15.843470
6 units:us(Trait);us(Trait) 1964.20 138.080 14.225087
7 units:us(Trait);us(Trait) 3596.00 226.990 15.842108
$G2
Term Sigma SE Z.ratio
1 Residual 1.000 0.000 Inf
2 diag(Trait):Fam;diag(Trait) 17.918 15.986 1.1208557
3 diag(Trait):Fam;diag(Trait) 9.185 32.960 0.2786711
4 units:us(Trait);us(Trait) 1692.500 103.450 16.3605607
5 units:us(Trait);us(Trait) 2121.000 141.610 14.9777558
6 units:us(Trait);us(Trait) 3826.100 234.020 16.3494573
$G0
Term Sigma SE Z.ratio
1 Residual 1.000 0.000 Inf
2 Trait:Fam 15.695 10.719 1.464222
3 units:us(Trait);us(Trait) 1694.500 103.260 16.410033
4 units:us(Trait);us(Trait) 2119.800 141.520 14.978802
5 units:us(Trait);us(Trait) 3817.500 231.060 16.521683
```
通过管道操作符`%>%`来进行模型的比较：
```R
> mc<-function(x,LRT0=T){
+ G1<-x[[1]];G2<-x[[2]];G0<-x[[3]]
+ model.comp(G1,G2,G0,LRT=LRT0)
+ }
> res12b %>% b2s %>% mc
Model comparison results as following:
parNO LogL AIC BIC AIC.State BIC.State
G0 4 -5562.18 11132.37 11152.40
G2 5 -5562.23 11134.46 11159.49
G1 6 -5552.98 11117.96 11148.01 better better
=====================================
Likelihood ratio test (LRT) results:
note:left model before "/" is full model,right is reduced.
Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)
df LR-statistic Pr(Chisq)
G2/G0 1 -0.1 0.5
G1/G2 1 18.5 8.495e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
```
同样道理，不能直接对`res12b`使用pin()来获取各模型的遗传参数，可考虑通过管道操作符`%>%`依次进行模型参数的计算，例如第一个模型：
```R
> getM<-function(x) x[[1]]
> res12b %>% b2s %>% getM %>% pin(.,mulp=c(gcor~V3/sqrt(V2*V4)))
variance components are as following:
Term Sigma SE vcS
1 Residual 1.00 0.000 V1
2 us(Trait):Fam;us(Trait) 132.43 56.428 V2
3 us(Trait):Fam;us(Trait) 168.81 75.851 V3
4 us(Trait):Fam;us(Trait) 243.39 115.710 V4
5 units:us(Trait);us(Trait) 1589.10 100.300 V5
6 units:us(Trait);us(Trait) 1964.20 138.080 V6
7 units:us(Trait);us(Trait) 3596.00 226.990 V7
pin formula:
gcor ~ V3/sqrt(V2 * V4)
<environment: 0x000001ce5c8e82e8>
Term Estimate SE
1 gcor 0.940 0.0612
```
通过管道操作符`%>%`来展示各模型的固定效应值：
```R
> solR<-coef(res12b)
> solR %>% lapply(., function(x) x$fixed) %>%
+ purrr::map(head)
$G1
Term Level Effect SE
1 Trait:Rep 1.r1 0.000 0.000
2 Trait:Rep 1.r2 -11.753 5.641
3 Trait:Rep 1.r3 46.756 5.534
4 Trait:Rep 1.r4 23.774 5.521
5 Trait:Rep 1.r5 4.891 5.510
6 Trait:Rep 2.r1 0.000 0.000
$G2
Term Level Effect SE
1 Trait:Rep 1.r1 0.000 0.000
2 Trait:Rep 1.r2 -13.246 5.775
3 Trait:Rep 1.r3 45.737 5.647
4 Trait:Rep 1.r4 22.641 5.627
5 Trait:Rep 1.r5 3.702 5.625
6 Trait:Rep 2.r1 0.000 0.000
$G0
Term Level Effect SE
1 Trait:Rep 1.r1 0.000 0.000
2 Trait:Rep 1.r2 -13.296 5.778
3 Trait:Rep 1.r3 45.691 5.649
4 Trait:Rep 1.r4 22.587 5.629
5 Trait:Rep 1.r5 3.676 5.627
6 Trait:Rep 2.r1 0.000 0.000
```
类似地，各模型随机效应值的展示：
```R
> solR %>% lapply(.,function(x) x$random) %>%
+ purrr::map(head)
$G1
Term Level Effect SE
1 us(Trait):Fam 1.70048 3.991 7.776
2 us(Trait):Fam 1.70017 -0.412 8.569
3 us(Trait):Fam 1.70002 2.503 8.394
4 us(Trait):Fam 1.70010 3.319 7.781
5 us(Trait):Fam 1.70041 2.792 8.754
6 us(Trait):Fam 1.70031 6.732 8.388
$G2
Term Level Effect SE
1 diag(Trait):Fam 1.70048 -0.947 3.497
2 diag(Trait):Fam 1.70017 -3.595 3.684
3 diag(Trait):Fam 1.70002 0.761 3.645
4 diag(Trait):Fam 1.70010 1.969 3.498
5 diag(Trait):Fam 1.70041 1.498 3.725
6 diag(Trait):Fam 1.70031 0.423 3.643
$G0
Term Level Effect SE
1 Trait:Fam 1.70048 -0.779 3.362
2 Trait:Fam 1.70017 -3.125 3.514
3 Trait:Fam 1.70002 0.674 3.482
4 Trait:Fam 1.70010 1.729 3.363
5 Trait:Fam 1.70041 1.313 3.548
6 Trait:Fam 1.70031 0.412 3.481
```
参考文献：
Zhang WH, Wei RY, Liu Y, Lin YZ. AFEchidna is a R package for genetic evaluation of plant and animal breeding datasets. BioRxiv. DOI: 10.1101/2021.06.24.449740.
首发博文：https://blog.csdn.net/yzhlinscau/article/details/124561257?spm=1001.2014.3001.5501
https://blog.sciencenet.cn/blog-1114360-1336951.html
上一篇：《动植物育种遗传数据分析》之AFEchidna篇--第5章
下一篇：性状的pearson相关与表型相关的差异
