yijiaobai的个人博客分享 http://blog.sciencenet.cn/u/yijiaobai

博文

混合线性模型中BLUE值 VS BLUP值

已有 207 次阅读 2018-5-8 11:35 |个人分类:农学统计|系统分类:科研笔记

混合线性模型中BLUE值 VS BLUP值

最佳线性无偏预测(best linear unbiased prediction, 简称BLUP),又音译为“布拉普”[1],是统计学上用于线性混合模型对随机效应进行预测的一种方法。最佳线性无偏预测由C.R. Henderson提出。随机效应的最佳线性无偏预测(BLUP)等同于固定效应的最佳线性无偏估计(best linear unbiased estimates, BLUE)(参见高斯-马尔可夫定理)。因为对固定效应使用估计一词,而对随机效应使用预测,这两个术语基本是等同的。BLUP被大量使用于动物育种。——wiki

BLUP值,相当于是对混合线性模型中随机因子的预测;

BLUE值,相当于是对混合线性模型中固定因子的估算

predict means:预测均值,固定因子和随机因子都可以预测均值,它的尺度和表型值尺度一致

示例

将处理作为固定因子

setwd("D:\\02 ASReml\\blue VS blup")library(asreml)library(tidyverse)
dat <- read.csv("MaizeRILs.csv",head=T)for (i in 1:4) dat[,i] <- as.factor(dat[,i])
as1 <- asreml(height ~ location/rep + location*RIL,data=dat)
ASReml: Tue May 08 11:07:55 2018

     LogLik         S2      DF      wall     cpu
   -723.8797     64.8862   244  11:07:55     0.1
   -723.8797     64.8862   244  11:07:55     0.0

Finished on: Tue May 08 11:07:55 2018

LogLikelihood Converged 

计算品种的BLUE值

ablue <- coef(as1)$fixed
blue1 <- ablue[grep("^RIL_RIL*",rownames(ablue)),] %>% as.data.frame()
head(blue1)

计算品种的预测均值(predict means)

pv1 <- predict(as1,"RIL")$predictions$pvals
ASReml: Tue May 08 11:13:33 2018

     LogLik         S2      DF      wall     cpu
   -723.8797     64.8862   244  11:13:33     0.0
   -723.8797     64.8862   244  11:13:33     0.0

Finished on: Tue May 08 11:13:33 2018

LogLikelihood Converged 
head(pv1) # 类似SAS中的lsmeans

运行模型:因素作为随机因子

as2 <- asreml(height ~ 1,random = ~location/rep + location*RIL,data=dat)
ASReml: Tue May 08 11:13:34 2018

     LogLik         S2      DF      wall     cpu
  -1646.5302    233.5135   495  11:13:34     0.0
  -1569.0397    137.9186   495  11:13:34     0.0
  -1507.3257     94.6888   495  11:13:34     0.0
  -1471.3354     74.5149   495  11:13:34     0.0
  -1462.9209     67.6142   495  11:13:34     0.0
  -1461.7649     65.3553   495  11:13:34     0.0
  -1461.7228     64.9069   495  11:13:34     0.0
  -1461.7228     64.8863   495  11:13:34     0.0
  -1461.7228     64.8862   495  11:13:34     0.0

Finished on: Tue May 08 11:13:34 2018

LogLikelihood Converged 
blup <- coef(as2)$random
blup2 <- blup[grep("^RIL_RIL-*",rownames(blup)),] %>% as.data.frame()
head(blup2)

预测均值

pv2 <- predict(as2,"RIL")$predictions$pvals
ASReml: Tue May 08 11:13:34 2018

     LogLik         S2      DF      wall     cpu
  -1461.7228     64.8862   495  11:13:34     0.0
  -1461.7228     64.8862   495  11:13:34     0.0
  -1461.7228     64.8862   495  11:13:34     0.0
  -1461.7228     64.8862   495  11:13:34     0.0

Finished on: Tue May 08 11:13:34 2018

LogLikelihood Converged 
head(pv2)

计算遗传力

summary(as2)$varcomp
str(dat)
'data.frame':    496 obs. of  9 variables:
 $ location: Factor w/ 4 levels "ARC","CLY","PPAC",..: 1 1 2 2 3 3 4 4 1 1 ...
 $ rep     : Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2 1 2 ...
 $ block   : Factor w/ 8 levels "1","2","3","4",..: 4 6 5 4 8 5 1 4 1 2 ...
 $ plot    : Factor w/ 122 levels "1","2","3","4",..: 28 47 36 92 64 40 7 27 6 9 ...
 $ RIL     : Factor w/ 62 levels "RIL-1","RIL-11",..: 1 1 1 1 1 1 1 1 2 2 ...
 $ pollen  : int  73 74 71 73 97 95 72 72 69 69 ...
 $ silking : int  77 79 74 77 101 100 78 78 71 72 ...
 $ ASI     : int  4 5 3 4 4 5 6 6 2 3 ...
 $ height  : num  182 169 213 203 156 ...
VSNR::pin(as2,h2 ~ V3/(V3 + V4/4 + V5/(2*4)))

pin

将数据保存到excel中

library(openxlsx)
write.xlsx(blue1,"blue.xlsx")
write.xlsx(blup2,"blup.xlsx")
write.xlsx(pv1,"pm1.xlsx")
write.xlsx(pv2,"pm2.xlsx")

结果解析

图片.png

  • RIL是基因型
  • pm2-random是RIL作为随机因子的预测均值
  • pm1-fixed是RIL作为固定因子时的预测均值
  • blue是RIL作为固定因子的BLUE值
  • blup是RIL作为随机因子的BLUP值
  • pm2-blup 是随机因子的预测均值 减去 随机因子的BLUP值,可以看到得到的是一个常数(均值)
  • pm1-mu-random 是固定因子的预测均值 减去 固定依着你的BLUE值, 可以看到不是一个常数
  • blup/blue_effect=heritibility 是BLUP值 除以 BLUE效应值,得到的是遗传力常数
  • 备注:blue_effect是用固定因子的预测均值 减去 整体均值

结论

  • 1,如果想用固定因子的效应值,可以使用固定因子的predict means(预测均值)-整体均值,因为固定因子的BLUE第一个强制为0,可以用预测均值代替
  • 2,随机因子的预测均值和BLUP值,在数据平衡时,他们之差为常数(整体均值)
  • 3,尽量使用BLUP值(或者随机因子的预测均值)作为排名依据,因为它相对于固定因子的效应值,考虑了遗传力的影响

获取数据、代码及结果

关注公众号,回复blup2,获取下载链接

wechat



http://blog.sciencenet.cn/blog-2577109-1112913.html

上一篇:如何计算遗传相关(How to)
下一篇:GWAS分析中表型值是使用BLUE值还是BLUP值?

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2018-5-21 23:12

Powered by ScienceNet.cn

Copyright © 2007-2017 中国科学报社

返回顶部