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

博文

ASReml-R多性状分析,怎么设置初始值?

已有 3045 次阅读 2016-6-23 18:22 |个人分类:农学统计|系统分类:科研笔记

建立了一个数据分析500人群,感兴趣的老师、同学

添加微信号: joinVSNC 为好友,拉你进入数据分析群


先进行单变量分析的原因:

1、确定显著性的因素(固定因素、随机因素),剔除不显著的

2、确定方差组分,为多变量分析时提供拟合初始值。


设置初始值的目的:

1、容易拟合

2、快速拟合

具体代码如下:

library(asreml)

#载入数据

data(harvey)

head(harvey)

#查看数据

str(harvey)

#模型1,性状y1

my1 <- asreml(fixed = y1~1,random = ~ Sire,data=harvey)summary(my1)$varcomp

#模型2,性状y2

my2 <- asreml(fixed = y2~1,random = ~ Sire,data=harvey)summary(my2)$varcomp



##两性状合并,默认初始值mula <- asreml(cbind(y1,y2)~trait,random = ~ us(trait):Sire,rcov=~ units:us(trait),data = harvey)summary(mula)$varcomp

####两性状合并,设置初始值,第一个是y1的方差分量,第二个是协方差,默认0.1

####第三个是y2的方差分量,残差也是这样。

mulb <- asreml(cbind(y1,y2)~trait,    

              random = ~ us(trait,init=c(88.78,0.1,84.87)):Sire,  

              rcov=~ units:us(trait,init=c(132.27,0.1,1638.40)),

              data = harvey)summary(mulb)$varcomp


结论:

> library(asreml)

> #载入数据

> data(harvey)

> head(harvey)

 Calf   Sire Dam Line ageOfDam  y1  y2  y3

1  101 Sire_1   0    1        3 192 390 224

2  102 Sire_1   0    1        3 154 403 265

3  103 Sire_1   0    1        4 185 432 241

4  104 Sire_1   0    1        4 183 457 225

5  105 Sire_1   0    1        5 186 483 258

6  106 Sire_1   0    1        5 177 469 267

> #查看数据

> str(harvey)

'data.frame':65 obs. of  8 variables:

$ Calf    : Factor w/ 65 levels "101","102","103",..: 1 2 3 4 5 6 7 8 9 10 ...

$ Sire    : Factor w/ 9 levels "Sire_1","Sire_2",..: 1 1 1 1 1 1 1 1 2 2 ...

$ Dam     : Factor w/ 1 level "0": 1 1 1 1 1 1 1 1 1 1 ...

$ Line    : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...

$ ageOfDam: int  3 3 4 4 5 5 5 5 4 4 ...

$ y1      : int  192 154 185 183 186 177 177 163 188 178 ...

$ y2      : int  390 403 432 457 483 469 428 439 439 407 ...

$ y3      : int  224 265 241 225 258 267 271 247 229 226 ...

> #模型1,性状y1

> my1 <- asreml(fixed = y1~1,random = ~ Sire,data=harvey)

ASReml: Fri Jul 29 16:58:33 2016


    LogLik         S2      DF      wall     cpu

  -201.0882    172.6344    64  16:58:33     0.0

  -199.5907    159.6531    64  16:58:33     0.0

  -198.2345    146.5045    64  16:58:33     0.0

  -197.5085    136.0169    64  16:58:33     0.0

  -197.4312    132.7488    64  16:58:33     0.0

  -197.4297    132.2748    64  16:58:33     0.0

  -197.4297    132.2688    64  16:58:33     0.0

  -197.4297    132.2689    64  16:58:33     0.0


Finished on: Fri Jul 29 16:58:33 2016

LogLikelihood Converged

> summary(my1)$varcomp

                 gamma component std.error  z.ratio constraint

Sire!Sire.var 0.6712362  88.78366  53.52084 1.658862   Positive

R!variance    1.0000000 132.26887  24.97546 5.295952   Positive

> #模型2,性状y2

> my2 <- asreml(fixed = y2~1,random = ~ Sire,data=harvey)

ASReml: Fri Jul 29 16:58:36 2016


    LogLik         S2      DF      wall     cpu

  -272.2863   1597.3735    64  16:58:36     0.0

  -272.2354   1611.3459    64  16:58:36     0.0

  -272.2050   1627.9973    64  16:58:36     0.0

  -272.2001   1639.1024    64  16:58:36     0.0

  -272.2001   1638.3948    64  16:58:36     0.0

  -272.2001   1638.4026    64  16:58:36     0.0


Finished on: Fri Jul 29 16:58:36 2016

LogLikelihood Converged

> summary(my2)$varcomp

                  gamma  component std.error   z.ratio constraint

Sire!Sire.var 0.05179991   84.86911  160.4663 0.5288907   Positive

R!variance    1.00000000 1638.40263  309.1192 5.3002291   Positive

> ##两性状合并,默认初始值

> mula <- asreml(cbind(y1,y2)~trait,random = ~ us(trait):Sire,rcov=~ units:us(trait),data = harvey)

ASReml: Fri Jul 29 16:58:41 2016


US matrix updates modified 1 times to remain positive definite.

    LogLik         S2      DF      wall     cpu

  -516.9617      1.0000   128  16:58:41     0.0 (6 restrained)

  -473.1624      1.0000   128  16:58:41     0.0 (3 restrained)

  -472.2445      1.0000   128  16:58:41     0.0 (1 restrained)

  -472.0819      1.0000   128  16:58:41     0.0

  -466.0507      1.0000   128  16:58:41     0.0

  -464.1416      1.0000   128  16:58:41     0.0

  -463.7376      1.0000   128  16:58:41     0.0

  -463.7018      1.0000   128  16:58:41     0.0

  -463.7014      1.0000   128  16:58:41     0.0

  -463.7014      1.0000   128  16:58:41     0.0

US variance structures were modified in 1 instances to make them positive definite


Finished on: Fri Jul 29 16:58:41 2016

LogLikelihood Converged

> summary(mula)$varcomp

                           gamma  component std.error   z.ratio constraint

trait:Sire!trait.y1:y1   90.26179   90.26179  54.25120 1.6637749   Positive

trait:Sire!trait.y2:y1   50.00361   50.00361  69.09275 0.7237172   Positive

trait:Sire!trait.y2:y2   73.28462   73.28462 155.42184 0.4715208   Positive

R!variance                1.00000    1.00000        NA        NA      Fixed

R!trait.y1:y1           132.06585  132.06585  24.91871 5.2998665   Positive

R!trait.y2:y1           191.28116  191.28116  67.30377 2.8420572   Positive

R!trait.y2:y2          1646.72905 1646.72905 311.40526 5.2880579   Positive

> ####两性状合并,设置初始值,第一个是y1的方差分量,第二个是协方差,默认0.1

> ####第三个是y2的方差分量,残差也是这样。

> mulb <- asreml(cbind(y1,y2)~trait,

+ random = ~ us(trait,init=c(88.78,0.1,84.87)):Sire,

+ rcov=~ units:us(trait,init=c(132.27,0.1,1638.40)),

+ data = harvey)

ASReml: Fri Jul 29 16:58:50 2016


    LogLik         S2      DF      wall     cpu

  -469.6227      1.0000   128  16:58:50     0.0 (2 restrained)

  -469.5193      1.0000   128  16:58:50     0.0 (2 restrained)

  -467.9921      1.0000   128  16:58:50     0.0

  -466.4986      1.0000   128  16:58:50     0.0

  -463.9073      1.0000   128  16:58:50     0.0

  -463.7038      1.0000   128  16:58:50     0.0

  -463.7014      1.0000   128  16:58:50     0.0

  -463.7014      1.0000   128  16:58:50     0.0

  -463.7014      1.0000   128  16:58:50     0.0


Finished on: Fri Jul 29 16:58:50 2016

LogLikelihood Converged

> summary(mulb)$varcomp

                           gamma  component std.error   z.ratio constraint

trait:Sire!trait.y1:y1   90.26179   90.26179  54.25132 1.6637713   Positive

trait:Sire!trait.y2:y1   50.00350   50.00350  69.09328 0.7237101   Positive

trait:Sire!trait.y2:y2   73.28490   73.28490 155.42495 0.4715131   Positive

R!variance                1.00000    1.00000        NA        NA      Fixed

R!trait.y1:y1           132.06585  132.06585  24.91871 5.2998670   Positive

R!trait.y2:y1           191.28117  191.28117  67.30366 2.8420618   Positive

R!trait.y2:y2          1646.72876 1646.72876 311.40406 5.2880773   Positive




https://blog.sciencenet.cn/blog-2577109-986408.html

上一篇:SAS 处理裂区试验代码
下一篇:R语言汇总统计,包括分组数据的汇总,结果友好(代码分享)
收藏 IP: 113.53.149.*| 热度|

0

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

数据加载中...

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

GMT+8, 2024-9-19 09:18

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部