||
建立了一个数据分析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
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-20 04:33
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社