育种数据分析之放飞自我分享 http://blog.sciencenet.cn/u/yijiaobai 关注:生物统计,数量遗传,混合线性模型,生物信息,R,Perl,Python,GWAS,GS相关方法,文章及代码

博文

R中的线性混合模型介绍(翻译博客)

已有 35358 次阅读 2016-1-12 10:30 |个人分类:翻译博客|系统分类:科研笔记| 混合线性模型, 育种值, BLUP

Linear mixed models in R

Luis


我的工作就是统计相关,但是很大部分的工作都是用线性混合模型,我用线性混合模型做的比较多的工作是分析试验数据的遗传结构,各种遗传参数评估(遗传力、配合力、育种值等)。


A brief history of time

在1992-1995年之间,我刚开始用SAS(首先是proc glm,最后是proc mixed),后来发现当我想用SAS计算动物模型的BLUP(最佳线性无偏预测值)时,过程是很痛苦,而且时间很漫长,于是1995-1996我使用DFREML(Karen Meyer,现在是WOMBAT)和AIREML(现在已经没有了),这些软件是为分析动物育种后代的性状,对于试验设计的数据分析比较麻烦。最后在1996或者1997早期,我开始接触使用ASREML(程序是由Arthur和Robin、Brian编写,商业化运作由VSNi公司)我甚至写了本ASReml的cookbook,那时我还用SAS来清洗数据,直到2003年我还没真正开始使用R(或S)语言,后来我认识到在笔记本电脑里面装个虚拟机来运行SAS,是得不偿失的,于是在2009年我完全放弃了SAS开始使用R语言。

Options

选择是一个很头疼的问题,R中有很多软件包可以做混合线性模型,这里我只介绍nlme、lme4和ASreml(对!ASreml是商业版,但是还有对应的R包),这些都是基于Reml的算法,当然还有一些包是基于贝叶斯的算法,这部分在其他章节介绍。

nlme:

这是一个比较成熟的R包,是R语言安装时默认的包,它除了可以分析分层的线性混合模型,也可以处理非线性模型。在优势方面,个人认为它可以处理相对复杂的线性和非线性模型,可以定义方差协方差结构,可以在广义线性模型中定义几种分布函数和连接函数。

它的短板:

1、随机效应的定义过于呆板

2、数据量大时速度很慢

3、不能处理系谱数据

4、不能处理多变量数据。

lme4

lme4包是由Douglas Bates开发,他也是nlme包的作者之一,相对于nlme包而言,它的运行速度快一点,对于随机效应的结构也可以更复杂一点,但是它的缺点也和nlme一样:

1、不能处理协方差和相关系数结构

2、它可以与构建系数的包连接,比如mmpedigree包,但是结合比较脆弱。

3、不能处理多性状数据


ASReml-R:

ASReml-R是ASReml的R版本,它的优点:

1、可以处理复杂的随机因子结构

2、可以处理多元数据

3、可以处理系谱数据

4、可以处理大批量的数据

主要的缺点:

1、它是收费的,当然它对于不发达国家的科研机构是免费的,不过需要申请和被审核。

它的用户主要是育种公司、科研机构等,它可以在各种平台上运行,包括Windows、Linux、OS X等。


一个简单示例:

我将用一个传统的裂区数据来说明不同软件包的用法,这个数据oats是在MASS包中,是研究大麦品种和N肥处理的裂区试验,其中品种为主区,肥料为裂区。代码如下

1、data


library(MASS)

data(oats)

names(oats) = c('block', 'variety', 'nitrogen', 'yield')

oats$mainplot = oats$variety

oats$subplot = oats$nitrogen

summary(oats)  


2、nlme包

调用这个包很简单,y变量写在左边,然后是固定因子,然后是随机因子,注意1|block/mainplot是裂区试验残差的写法,因为里面有两个残差。代码如下:


library(nlme)

m1.nlme = lme(yield ~ variety*nitrogen,

                     random = ~ 1|block/mainplot,

                     data = oats)

summary(m1.nlme)

它的方差分析结果为:

anova(m1.nlme)

                numDF denDF   F-value p-value

(Intercept)          1    45 245.14333  <.0001

variety              2    10   1.48534  0.2724

nitrogen             3    45  37.68561  <.0001

variety:nitrogen     6    45   0.30282  0.9322

3、lme4包

lme4包的语法也相似,随机效应有着和nlme相同的语法,不同的是lme4包它的结果给出了随机效应的标准差,而不是方差。

代码如下:


library(lme4)

m1.lme4 = lmer(yield ~ variety*nitrogen + (1|block/mainplot),

                      data = oats)

summary(m1.lme4)

anova(m1.lme4)

Analysis of Variance Table

                Df  Sum Sq Mean Sq F value

variety           2   526.1   263.0  1.4853

nitrogen          3 20020.5  6673.5 37.6856

variety:nitrogen  6   321.7    53.6  0.3028


ASReml-R包

它的功能很强大,用在这里有些杀鸡用牛刀的感觉。

代码如下:



library(asreml)

m1.asreml = asreml(yield ~ variety*nitrogen,

                  random = ~ block/mainplot,

                  data = oats)

summary(m1.asreml)$varcomp

wald(m1.asreml)

Wald tests for fixed effects

Response: yield

Terms added sequentially; adjusted for those above

                Df Sum of Sq Wald statistic Pr(Chisq)    

(Intercept)       1     43410        245.141    <2e-16 ***

variety           2       526          2.971    0.2264    

nitrogen          3     20020        113.057    <2e-16 ***

variety:nitrogen  6       322          1.817    0.9357    

residual (MS)             177                            

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

译者注:

翻译的过程也是学习的过程,看着运行的结果,和页面结果上一样,感觉棒棒的。





如果您对于数据分析,对于软件操作,对于数据整理,对于结果理解,有任何问题,欢迎联系我。
r-breeding.png




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

上一篇:GenStat软件介绍
下一篇:R语言中计算blue和blup值的包:lme4、nlme、MCMCglmm、asreml

1 覃伟

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

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2021-1-19 04:55

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部