注: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.