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

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

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

1、确定显著性的因素（固定因素、随机因素），剔除不显著的

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

设置初始值的目的：

1、容易拟合

2、快速拟合

具体代码如下：结论：

> 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语言汇总统计，包括分组数据的汇总，结果友好（代码分享）