博文
新包breedRPlus的用法
|||
本人最近在breedR包基础上，编写了新包breedRPlus，用于添加部分函数，便于breedR包的使用，本程序包已上传至github，可免费在线安装。欢迎各位网友使用和反馈意见。
breedRPlus官方网站:https://github.com/yzhlinscau/breedRPlus
林元震
BreedRPlus manual
About breedRPlus
The breedRPlus is builded on the base of package 'breedR' for some
additional functions, such as calculating standard error (se), running
batch analysis, etc. At present, breedRPlus only works for 'ai' method
in 'breedR', not well worked for spatial (and competition) model at
present.
Install package
install.packages(c('devtools','dplyr','nadiv','ggplot2',
'plyr','msm','reshape2'))
devtools::install_github(c('yzhlinscau/breedRPlus'))
breedRPlus function
var() to get variance components, similar to asreml;
plot1() to test trait’s norm, similar to asreml;
pin() to calculate heritability or corr with se;
batch1() to run batch analysis for single trait.
00 Example data
Using data 'douglas' from package 'breedR':
data(douglas)
S3 <- subset(douglas, site == 's3' )
S3a <- filter(S3, is.na(dad))
S3a <- transform(S3a, Mum = factor(mum))
S3a <- droplevels(S3a)
01 single trait analysis
01A family model
m1 <- remlf90(fixed = H04 ~ 1 + orig,
random = ~ Mum + block,
data = S3a)
analysis results
test trait’s norm
plot1(m1)
get variance components:
var(m1)
## gamma component std.error z.ratio constraint
## Mum 0.07329389 624.31 303.18 2.059206 Positive # V1
## block 0.05959450 507.62 214.96 2.361463 Positive # V2
## Residual 1.00000000 8517.90 436.16 19.529301 Positive # V3
calculate heritability:
pin(m1, h2 ~ 4 * V1/(V1 + V3))
##
## Estimate SE
## h2 0.27315 0.12501
01B tree(animal) model
im1 <- remlf90(fixed = H04 ~ 1 + orig,
random = ~ block,
genetic = list(model = 'add_animal',
pedigree = S3a[,1:3],
id = 'self'),
data=S3a)
analysis results
test trait’s norm:
plot1(im1)
get variance components:
var(im1)
## gamma component std.error z.ratio constraint
## block 0.07639127 507.62 214.96 2.361463 Positive # V1
## genetic 0.37580135 2497.20 1212.70 2.059207 Positive # V2
## Residual 1.00000000 6645.00 1030.90 6.445824 Positive # V3
calculate heritability:
pin(im1, h2 ~ V2/(V2 + V3))
##
## Estimate SE
## h2 0.27315 0.12501
02 singe trait-batch anlysis
data reshape
names(S3a)[7:8] <- c('x1','y1')
df <- tidyr::gather(S3a, key = Trait, y, c(-1:-8,-12,-14:-16))
02A family model
# for family model
fixed <- y ~ 1 + orig
random1 <- ~ Mum + block
pformula1 <- h2 ~ 4 * V1/(V1 + V3)
mm1 <- plyr::ddply(df,'Trait',
function(dat) batch1(dat,FMod = fixed, RMod = random1, pformula = pformula1)
)
batch analysis result:
mm1
## Trait AIC logLik Mum block Residual Mum.se block.se
## 1 C13 10095.399 -5044.700 2218.20 537.08 21156.0 979.67 356.640
## 2 H02 9056.080 -4525.040 203.66 211.55 3818.3 112.64 94.259
## 3 H03 9399.865 -4696.933 387.74 303.46 5831.9 195.49 137.680
## 4 H04 9790.566 -4892.283 624.31 507.62 8517.9 303.18 214.960
## Residual.se h2 h2.se
## 1 1111.00 0.380 0.154
## 2 196.70 0.203 0.108
## 3 300.24 0.249 0.119
## 4 436.16 0.273 0.125
02B tree model
# for tree model
fixed <- y ~ 1 + orig
random2 <- ~ block
genetic <- list(model = 'add_animal',
pedigree = S3a[,1:3],
id = 'self')
pformula2 <- h2 ~ V2/(V2 + V3)
mm2 <- plyr::ddply(df,'Trait',
function(dat) batch1(dat,FMod = fixed,
RMod = random2,
geneticM = genetic,
pformula = pformula2)
)
batch analysis result:
mm2
## Trait AIC logLik block genetic Residual block.se genetic.se
## 1 C13 10095.399 -5044.700 537.08 8872.90 14501.0 356.640 3918.70
## 2 H02 9056.080 -4525.040 211.55 814.64 3207.3 94.259 450.54
## 3 H03 9399.865 -4696.933 303.46 1551.00 4668.7 137.680 781.97
## 4 H04 9790.566 -4892.283 507.62 2497.20 6645.0 214.960 1212.70
## Residual.se h2 h2.se
## 1 3192.20 0.380 0.154
## 2 402.83 0.203 0.108
## 3 674.84 0.249 0.119
## 4 1030.90 0.273 0.125
03 multi-trait model
03A family model
mt1 <- remlf90(fixed = cbind(H04, C13) ~ orig,
random = ~ Mum + block,
data = S3a)
get variance components:
var(mt1, mulT = TRUE)
## gamma component std.error z.ratio
## Mum.H04 620.71 620.71 301.90 2.056012
## Mum.H04_Mum.C13 962.38 962.38 503.90 1.909863
## Mum.C13 2208.60 2208.60 987.11 2.237441
## block.H04 531.88 531.88 221.03 2.406370
## block.H04_block.C13 362.63 362.63 270.11 1.342527
## block.C13 879.39 879.39 452.53 1.943274
## Residual.H04 8510.30 8510.30 435.58 19.537858
## Residual.H04_Residual.C13 12256.00 12256.00 676.82 18.108212
## Residual.C13 22894.00 22894.00 1194.90 19.159762
## constraint
## Mum.H04 Positive
## Mum.H04_Mum.C13 Positive
## Mum.C13 Positive
## block.H04 Positive
## block.H04_block.C13 Positive
## block.C13 Positive
## Residual.H04 Positive
## Residual.H04_Residual.C13 Positive
## Residual.C13 Positive
calculate heritability for both traits:
pin(mt1, h2.H04~ 4 * V1/(V1 + V7))
##
## Estimate SE
## h2.H04 0.27191 0.12468
pin(mt1, h2.C13~ 4 * V3/(V3 + V9))
##
## Estimate SE
## h2.C13 0.35193 0.14523
calculate genetic, residula, and phenotypic correlations:
pin(mt1, gcorr ~ V2/sqrt(V1 * V3), signif=TRUE)
##
## Estimate SE sig.level
## gcorr 0.82195 0.10289 ***
## ---------------
## Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
pin(mt1, ecorr ~ V8/sqrt(V7 * V9), signif=TRUE)
##
## Estimate SE sig.level
## ecorr 0.87804 0.00848 ***
## ---------------
## Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
pin(mt1, pcorr ~ (V2+V8)/sqrt((V1+V7)*(V3+V9)), signif=TRUE)
##
## Estimate SE sig.level
## pcorr 0.87309 0.01033 ***
## ---------------
## Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
03B tree(animal) model
mt2 <- remlf90(fixed = cbind(H04, C13) ~ orig,
random = ~ block,
genetic = list(model = 'add_animal',
pedigree = S3a[,1:3],
id = 'self'),
data = S3a)
get variance components:
var(mt2, mulT = TRUE )
## gamma component std.error
## block.H04 183970.0 183970.0 4592.1
## block.H04_block.C13 393310.0 393310.0 4916.5
## block.C13 842610.0 842610.0 0.0
## genetic.direct.H04 2435.6 2435.6 1195.5
## genetic.direct.H04_genetic.direct.C13 3698.6 3698.6 1974.4
## genetic.direct.C13 8415.0 8415.0 3835.2
## Residual.H04 6695.7 6695.7 1020.3
## Residual.H04_Residual.C13 9535.2 9535.2 1664.5
## Residual.C13 16714.0 16714.0 3177.0
## z.ratio constraint
## block.H04 40.062281 Positive
## block.H04_block.C13 79.997966 Positive
## block.C13 Inf Positive
## genetic.direct.H04 2.037307 Positive
## genetic.direct.H04_genetic.direct.C13 1.873278 Positive
## genetic.direct.C13 2.194149 Positive
## Residual.H04 6.562482 Positive
## Residual.H04_Residual.C13 5.728567 Positive
## Residual.C13 5.260938 Positive
calculate heritability for both traits:
pin(mt2, h2.H04 ~ V4/(V4 + V7))
##
## Estimate SE
## h2.H04 0.26673 0.12361
pin(mt2, h2.C13 ~ V6/(V6 + V9))
##
## Estimate SE
## h2.C13 0.33487 0.14161
calculate genetic, residula, and phenotypic correlations:
pin(mt2, gcorr ~ V5/sqrt(V4 * V6), signif=TRUE)
##
## Estimate SE sig.level
## gcorr 0.81697 0.10635 ***
## ---------------
## Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
pin(mt2, ecorr ~ V8/sqrt(V7 * V9), signif=TRUE)
##
## Estimate SE sig.level
## ecorr 0.90135 0.03005 ***
## ---------------
## Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
pin(mt2, pcorr ~ (V5+V8)/sqrt((V4+V7)*(V6+V9)), signif=TRUE)
##
## Estimate SE sig.level
## pcorr 0.87364 0.01025 ***
## ---------------
## Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
More details will be updated in the futures.
https://blog.sciencenet.cn/blog-1114360-1126208.html
上一篇：ggplot2绘制带概率或显著性的相关矩阵的热图
下一篇：浅谈统计方法的重要性
扫一扫，分享此博文
全部作者的精选博文
全部作者的其他最新博文
全部精选博文导读
相关博文
- • Semantic Mathematics: Abstraction and Real-World AI（初学者版）
- • The 9-No Problems within the DIKWP Model（初学者版）
- • Beyond Prof. Yucong Duan\´s Paradox of Mathematics（初学者版）
- • Challenging Traditional Mathematical Paradigms for AI （初学者版）
- • Yucong Duan\´s Paradox of Mathematics for AI Semantics（初学者版）
- • Paradox of Mathematics by Yucong Duan in DIKWP（初学者版）