|||
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 |
译者注:
翻译的过程也是学习的过程,看着运行的结果,和页面结果上一样,感觉棒棒的。
如果您对于数据分析,对于软件操作,对于数据整理,对于结果理解,有任何问题,欢迎联系我。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-12-5 12:37
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社