AAfun4包的部分更新

最新推荐文章于 2022-02-25 22:10:15 发布

yzhlinscau

最新推荐文章于 2022-02-25 22:10:15 发布

阅读量1k

点赞数

文章标签： R ASReml

本文链接：https://blog.csdn.net/yzhlinscau/article/details/96477198

版权

注：AAfun4包最新版为2.0，同时该版起，不再免费提供。读者如想使用，需付费使用。感兴趣的，可发邮件到yzhlinscau@163.com。

1–pin.batch: pin()的升级

pin.batch可用于批量运行pin函数，简单示例如下：


// An simple example
fm2b.asr<-asreml(cbind(h3,h5)~ trait+trait:Rep, 
                 random=~ us(trait):Fam, 
                 residual=~units:us(trait),
                 subset=Spacing=='3',data=df,maxit=40)
Var(fm2b.asr)
bp<-c(h2_A ~ 4 * V1/(V1+V5), 
      h2_B ~ 4 * V3/(V3+V7),
      gCORR ~ V2/sqrt(V1*V3),
      eCORR ~ V6/sqrt(V5*V7),
      pCORR ~ (V2+V6)/sqrt((V1+V5)*(V3+V7)))
pin.batch(fm2b.asr,bp,digit=5)
pin.batch(fm2b.asr,bp,digit=5,signif=T)

运行结果如下：

// result
> Var(fm2b.asr)
                        component std.error   z.ratio bound %ch
trait:Fam!trait_h3:h3    131.9156  56.33195  2.341754     P   0
trait:Fam!trait_h5:h3    190.1109  89.88093  2.115141     P   0
trait:Fam!trait_h5:h5    440.1269 186.57337  2.359002     P   0
units:trait!R              1.0000        NA        NA     F   0
units:trait!trait_h3:h3 1589.3402 100.32251 15.842308     P   0
units:trait!trait_h5:h3 2006.0292 157.85827 12.707786     P   0
units:trait!trait_h5:h5 5332.9624 336.94459 15.827417     P   0
> pin.batch(fm2b.asr,bp,digit=5)
      Estimate      SE
h2_A   0.30656 0.12515
h2_B   0.30495 0.12360
gCORR  0.78899 0.12806
eCORR  0.68904 0.02348
pCORR  0.69668 0.02241
> pin.batch(fm2b.asr,bp,digit=5,signif=T)
      Estimate      SE siglevel
h2_A   0.30656 0.12515       **
h2_B   0.30495 0.12360       **
gCORR  0.78899 0.12806      ***
eCORR  0.68904 0.02348      ***
pCORR  0.69668 0.02241      ***
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif'

2–asreml.batch的更新

asreml.batch更新了遗传参数计算公式，原来只能计算5个参数，且受参数赋值的限制，现在统一使用参数mulp来赋值，理论上不受计算公式数目的限制了。简单示例如下：

// An example
mulp1=c(h2 ~ 4 * V1/(V1+V2),
       H2~V1/(V1+V2/4))
 > asreml.batch(data=df1,factorN=1:5,traitN=c(9:13),
   +              FMod=y~1+Rep+Plot,
   +              RMod=~Fam,
   +              asrV=4,
   +              #family=asr_gaussian(),
   +              mulp=mulp1
   +              )
Program starts running batch analysis ------     

ASReml-R batch analysis results:

Fixed Factors -- Rep, Plot 
Randomed Factors -- Fam, units 
Index formula -- 
	h2 ~ 4 * V1/(V1 + V2)
	H2 ~ V1/(V1 + V2/4)

Variance order: V1-Fam, V2-units

  Trait       V1        V2    V1.se    V2.se     h2     H2  h2.se  H2.se Maxit
1    h1  11.9614   43.0035   3.1371   2.7240 0.8705 0.5266 0.1872 0.0685     6
2    h2  53.7874  623.9830  22.4836  39.4926 0.3174 0.2564 0.1266 0.0826     4
3    h3 132.2286 1591.7955  56.4882 100.7596 0.3068 0.2494 0.1253 0.0828     4
4    h4 239.4825 3618.3452 115.6495 229.4253 0.2483 0.2093 0.1161 0.0825     5
5    h5 445.0040 5345.0410 188.3244 338.9637 0.3074 0.2498 0.1243 0.0821     4
  Converge
1     TRUE
2     TRUE
3     TRUE
4     TRUE
5     TRUE	

Total running times: 0.37 sec elapsed

多性状的示例如下：

// An highlighted block
> mulp2=c(gcorr ~ V2/sqrt(V1*V3),
+        h2.A ~ 4*V1/sqrt(V1+V5),
+        h2.B ~ 4*V3/sqrt(V3+V7))
> asreml.batch(data=df1,factorN=1:5,traitN=c(10:13),
+              FMod=cbind(y1,y2)~trait+trait:Rep,
+              RMod=~us(trait):Fam,
+              EMod=~units:us(trait),
+              asrV=4,
+              #family=list(asr_gaussian(),asr_gaussian()), 
+              mulT=TRUE,mulN=2,mulR=TRUE,corMout=F,
+              mulp=mulp2
+              ) 
Program starts running batch analysis ------ 

ASReml-R batch analysis results:

Fixed Factors -- Rep 
Randomed Factors -- Fam, units 
Index formula -- 
	gcorr ~ V2/sqrt(V1 * V3)
	h2.A ~ 4 * V1/sqrt(V1 + V5)
	h2.B ~ 4 * V3/sqrt(V3 + V7)

Variance order: V1-Fam_y1:y1, V2-Fam_y2:y1, V3-Fam_y2:y2, V4-Rvar.F, V5-units_y1:y1, V6-units_y2:y1, V7-units_y2:y2

$`Varcomp`
  Trait       V1       V2       V3 V4        V5        V6       V7    V1.se    V2.se
1 h2-h3  54.2228  78.0998 132.0715  1  621.0410  801.1607 1589.248  22.4840  33.2004
2 h2-h4  54.5503  78.9329 241.3618  1  620.8553 1127.1223 3589.885  22.5492  44.4371
3 h2-h5  53.6143 119.8578 438.4297  1  621.3692 1102.6926 5326.062  22.3684  55.2260
4 h3-h4 132.4268 168.8040 243.3841  1 1589.0726 1964.1480 3595.996  56.4251  75.8468
5 h3-h5 131.9156 190.1109 440.1269  1 1589.3402 2006.0292 5332.962  56.3319  89.8809
6 h4-h5 237.9995 303.2088 441.9947  1 3601.0223 3659.1288 5333.342 114.9186 137.7661
     V3.se V4.se    V5.se    V6.se    V7.se  gcorr    h2.A    h2.B gcorr.se h2.A.se
1  56.3534    NA  39.1910  56.9579 100.3130 0.9229  8.3465 12.7332   0.0643  3.3703
2 115.2239    NA  39.1748  83.4551 226.6322 0.6879  8.3960 15.5976   0.1588  3.3789
3 186.0647    NA  39.2206  94.9475 336.4598 0.7818  8.2546 23.0983   0.1428  3.3553
4 115.7026    NA 100.2978 138.0756 226.9823 0.9403 12.7668 15.7116   0.0612  5.3060
5 186.5734    NA 100.3225 157.8583 336.9446 0.7890 12.7184 23.1704   0.1281  5.2985
6 187.1775    NA 227.7240 255.2506 337.2139 0.9349 15.3648 23.2642   0.0620  7.2860
  h2.B.se Maxit Converge
1  5.3001     7     TRUE
2  7.3090     7     TRUE
3  9.5644     7     TRUE
4  7.3300     7     TRUE
5  9.5824     7     TRUE
6  9.6103     7     TRUE

$Corr.erro.matrix
       h2     h3     h4     h5
h2 1.0000 0.9229 0.6879 0.7818
h3 0.0643 1.0000 0.9403 0.7890
h4 0.1588 0.0612 1.0000 0.9349
h5 0.1428 0.1281 0.0620 1.0000

$Corr.sig.matrix
    h2     h3     h4     h5
h2   1 0.9229 0.6879 0.7818
h3 ***      1 0.9403  0.789
h4 ***    ***      1 0.9349
h5 ***    ***    ***      1
=================
upper is corr and lower is error (or sig.level) for corr matrix.
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1

Total running times: 0.57 sec elapsed

3 model.comp的更新

统一使用mulM参数来输入asreml模型结果，简化函数参数，并让结果更直观些。简单示例如下：

// An highlighted block
fm<-asreml(h5~ 1+Rep,random=~ Fam, 
           subset=Spacing=='3',data=df,maxit=40)
fm1a<-asreml(cbind(dj,h5)~ trait+trait:Rep, 
             random=~ us(trait):Fam, 
             residual=~units:us(trait),
             subset=Spacing=='3',data=df,maxit=40)
fm1b<-asreml(cbind(dj,h5)~ trait+trait:Rep, 
             random=~ diag(trait):Fam, 
             residual=~units:us(trait),
             subset=Spacing=='3',data=df,maxit=40)
fm3a<-asreml(cbind(dj,h3,h5)~ trait+trait:Rep, 
             random=~ diag(trait):Fam, 
             residual=~units:us(trait),
             subset=Spacing=='3',data=df,maxit=40)
fm3b<-asreml(cbind(dj,h3,h5)~ trait+trait:Rep, 
             random=~ diag(trait):Fam, 
             residual=~units:us(trait),
             subset=Spacing=='3',data=df,maxit=40)

> #####   model comparison    #####
> model.comp(mulM=c(fm1a,fm1b),asrV=4)
Loading required package: plyr
Attension:
Fixed factors should be the same!

  Model      LogL Npm      AIC      BIC AIC.State BIC.State
1  fm1b -868.4305   6 1748.861 1778.907              better
2  fm1a -867.0246   7 1748.049 1783.102    better          
-----------------------------
Lower AIC and BIC is better model.

> model.comp(mulM=c(fm1a,fm1b),LRT=TRUE,asrV=4)
Attension:
Fixed factors should be the same!

  Model      LogL Npm      AIC      BIC AIC.State BIC.State
1  fm1b -868.4305   6 1748.861 1778.907              better
2  fm1a -867.0246   7 1748.049 1783.102    better          
-----------------------------
Lower AIC and BIC is better model.

Attension: Please check every asreml result's length is 43(V3) or 36(V4);
if the length < 43 or 36(V4), put the object at the end of Nml.
In the present, just allow one object's length < 43 or 36(V4).
=====================================
Model compared between  fm1a -- fm1b :
  Model      LogL Npm      AIC Pr(>F)  Sig.level
1  fm1b -868.4305   6 1748.861                  
2  fm1a -867.0246   7 1748.049  0.094 Not signif
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
=====================================
Attension: Ddf did not minus 0.5. 
When for corr model, against 0.

> model.comp(mulM=c(fm3a,fm3b,fm1a,fm1b),asrV=4)
Attension:
Fixed factors should be the same!

  Model       LogL Npm      AIC      BIC AIC.State BIC.State
1  fm1b  -868.4305   6 1748.861 1778.907              better
2  fm1a  -867.0246   7 1748.049 1783.102    better          
3  fm3a -3038.1773  10 6096.355 6150.494                    
4  fm3b -3038.1773  10 6096.355 6150.494                    
-----------------------------
Lower AIC and BIC is better model.

> model.comp(mulM=c(fm3a,fm3b,fm1a,fm1b),LRT=TRUE,asrV=4)
Attension:
Fixed factors should be the same!
Likelihood ratio test (LRT) results:

  Model       LogL Npm      AIC      BIC AIC.State BIC.State
1  fm1b  -868.4305   6 1748.861 1778.907              better
2  fm1a  -867.0246   7 1748.049 1783.102    better          
3  fm3a -3038.1773  10 6096.355 6150.494                    
4  fm3b -3038.1773  10 6096.355 6150.494                    
-----------------------------
Lower AIC and BIC is better model.

Attension: Please check every asreml result's length is 43(V3) or 36(V4);
if the length < 43 or 36(V4), put the object at the end of Nml.
In the present, just allow one object's length < 43 or 36(V4).
=====================================
Likelihood ratio test (LRT) results:

Model compared between  fm1a -- fm1b :
  Model      LogL Npm      AIC Pr(>F)  Sig.level
1  fm1b -868.4305   6 1748.861                  
2  fm1a -867.0246   7 1748.049  0.094 Not signif
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1

Model compared between  fm1b -- fm3a :
  Model       LogL Npm      AIC Pr(>F)  Sig.level
1  fm1b  -868.4305   6 1748.861                  
2  fm3a -3038.1773  10 6096.355      1 Not signif
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1

Model compared between  fm1b -- fm3b :
  Model       LogL Npm      AIC Pr(>F)  Sig.level
1  fm1b  -868.4305   6 1748.861                  
2  fm3b -3038.1773  10 6096.355      1 Not signif
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1

Model compared between  fm1a -- fm3a :
  Model       LogL Npm      AIC Pr(>F)  Sig.level
1  fm1a  -867.0246   7 1748.049                  
2  fm3a -3038.1773  10 6096.355      1 Not signif
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1

Model compared between  fm1a -- fm3b :
  Model       LogL Npm      AIC Pr(>F)  Sig.level
1  fm1a  -867.0246   7 1748.049                  
2  fm3b -3038.1773  10 6096.355      1 Not signif
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1

Model compared between  fm3a -- fm3b :
  Model      LogL Npm      AIC Pr(>F)  Sig.level
1  fm3a -3038.177  10 6096.355                  
2  fm3b -3038.177  10 6096.355      1 Not signif
---------------
Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
=====================================
Attension: Ddf did not minus 0.5. 
When for corr model, against 0.

yzhlinscau

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
AAfun4包的部分更新

1–pin.batch: pin()的升级pin.batch可用于批量运行pin函数，简单示例如下：// An simple examplefm2b.asr<-asreml(cbind(h3,h5)~ trait+trait:Rep, random=~ us(trait):Fam, residual=~units...
复制链接

扫一扫