ASReml-R Added functions
Attension: The version AAfun on github is very old and might not install online. Please note that I do not pay the fee for ASreml any more, therefore AAfun package will not be updated anymore.
The AAfun is builded on the base of package 'asreml'
for some additional functions, such as calculating standard error (se), running batch analysis, getting cov/var/corr matrix for FA models, etc.
install.packages(c('amap',"agricolae","devtools","desplot","plyr","reshape2"))
devtools::install_github('yzhlinscau/AAfun')
- pin() to calculate heritability or corr with se;
- asreml.batch() to run batch analysis;
- model.comp() to run model comparisons;
- met.corr() to get cov/var/corr matrix for FA models;
- met.plot() to plot MET data;
- met.biplot() to run biplot for MET factor analytic results;
- etc...
library(asreml)
library(AAfun)
demo('pin')
data(PrSpa,package='AAfun')
df<-PrSpa
data(dfm2,package='AAfun')
df2<-dfm2
fm<-asreml(h5~1+Rep, random=~Fam, subset=Spacing=='3',data=df,trace=F)
summary(fm)$varcomp[,1:3]
calculate heritability:
AAfun::pin(fm, h2 ~4*V1/(V1+V2),signif=T)
##
## Estimate SE
## h2 0.306 0.124
fm2<-asreml(cbind(dj,h5)~ trait+trait:Rep,
random=~ us(trait):Fam, rcov=~units:us(trait),
subset=Spacing=='3',data=df,maxit=40,trace=F)
summary(fm2)$varcomp[,1:3]
calculate heritability for both traits:
pin(fm2, h2_A ~ 4 * V1/(V1+V5))
##
## Estimate SE
## h2_A 0.455 0.145
pin(fm2, h2_B ~ 4 * V3/(V3+V7))
calculate genetic and phenotypic corr between traits:
pin(fm2, gCORR ~ V2/sqrt(V1*V3),signif=TRUE)
##
## Estimate SE sig.level
## gCORR 0.442 0.257 *
## ---------------
## Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
pin(fm2, pCORR ~ (V2+V6)/sqrt((V1+V5)*(V3+V7)),signif=TRUE)
fm3<-asreml(cbind(dj,h3,h5)~ trait+trait:Rep,
random=~ corgh(trait):Fam, rcov=~units:us(trait),
subset=Spacing=='3',data=df,maxit=40,trace=F)
summary(fm3)$varcomp[,1:3]
return corr results:
pin(fm3,corN=3)
##
## Estimate SE sig.level
## trait:Fam!trait.h3:!trait.dj.cor 0.751 0.233 ***
## trait:Fam!trait.h5:!trait.dj.cor 0.448 0.257 *
## trait:Fam!trait.h5:!trait.h3.cor 0.798 0.123 ***
## ---------------
## Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
just return the first corr:
pin(fm3)
pin() also works for data with pedigree files.
df1=subset(df,Spacing==3)
AAfun::asreml.batch(data=df1,factorN=1:5,traitN=c(6:13),
FMod=y~1+Rep+Plot,RMod=~Fam,
pformula=h2 ~ 4 * V1/(V1+V2))
##
## ASReml-R batch analysis results:
## Fixed Factors -- Rep
## Randomed Factors -- Fam R
## Index formula -- h2 ~ 4 * V1/(V1 + V2)
##
## Variance order: Fam, R
## Trait V1 V2 V1.se V2.se h2 h2.se Converge Maxit
## 1 dj 0.0001 5.00e-04 0.0000 0.0000 0.447 0.144 TRUE 6
## 2 dm 0.0001 2.00e-03 0.0001 0.0001 0.266 0.122 TRUE 6
## 3 wd 0.0001 6.00e-04 0.0000 0.0000 0.475 0.151 TRUE 6
## 4 h1 12.0487 4.31e+01 3.1499 2.7194 0.875 0.187 TRUE 7
## 5 h2 54.2123 6.21e+02 22.4845 39.1922 0.321 0.127 TRUE 6
## 6 h3 132.5869 1.59e+03 56.4592 100.2910 0.308 0.125 TRUE 6
## 7 h4 241.3909 3.60e+03 115.5059 227.5395 0.251 0.116 TRUE 6
## 8 h5 441.9941 5.33e+03 187.1784 337.2136 0.306 0.124 TRUE 6
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),
mulT=TRUE,mulN=2,mulR=TRUE,corMout=T,
pformula=r.g ~ V2/sqrt(V1*V3),
pformula1=h2.A ~ 4*V1/(V1+V5),
pformula2=h2.B ~ 4*V3/sqrt(V3+V7))
asreml.batch(data=df1,factorN=1:5,traitN=c(10:13),
FMod=cbind(y1,y2,y3)~trait+trait:Rep,
RMod=~corgh(trait):Fam,
EMod=~units:us(trait), maxit=30,
mulT=TRUE,mulN=3,mulR=TRUE,corM=TRUE)
ped<-df2[,1:3]
pedinv<-asreml.Ainverse(ped)$ginv
df2a=subset(df2,Spacing==3)
asreml.batch(data=df2a,factorN=1:6,traitN=c(7:14),
FMod=y~1+Rep,RMod=~ped(TreeID),
ped=T,pedinv=pedinv,
ginverse=list(TreeID=pedinv),
pformula=h2 ~ V1/(V1+V2))
asreml.batch(data=df2a,factorN=1:6,traitN=c(10:14),
FMod=cbind(y1,y2)~trait+trait:Rep,
RMod=~us(trait):ped(TreeID),
EMod=~units:us(trait),maxit=40,
mulT=TRUE,mulN=2,mulR=TRUE,corMout=F,
ped=T,pedinv=pedinv,
ginverse=list(TreeID=pedinv),
pformula=r.g ~ V2/sqrt(V1*V3),
pformula1=h2.A ~ V1/(V1+V5),
pformula2=h2.B ~ V3/(V3+V7))
asreml.batch(data=df2a,factorN=1:6,traitN=c(10:14),
FMod=cbind(y1,y2)~trait+trait:Rep,
RMod=~corgh(trait):ped(TreeID),
EMod=~units:us(trait),maxit=40,
mulT=TRUE,mulN=2,corM=TRUE,
ped=T,pedinv=pedinv,
ginverse=list(TreeID=pedinv))
fm1<-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, rcov=~units:us(trait),
subset=Spacing=='3',data=df,maxit=40)
fm1b<-asreml(cbind(dj,h5)~ trait+trait:Rep,
random=~ diag(trait):Fam, rcov=~units:us(trait),
subset=Spacing=='3',data=df,maxit=40)
fm3a<-asreml(cbind(dj,h3,h5)~ trait+trait:Rep,
random=~ diag(trait):Fam, rcov=~units:us(trait),
subset=Spacing=='3',data=df,maxit=40)
fm3b<-asreml(cbind(dj,h3,h5)~ trait+trait:Rep,
random=~ diag(trait):Fam, rcov=~units:us(trait),
subset=Spacing=='3',data=df,maxit=40)
comparison between two models:
AAfun::model.comp(m1=fm1a,m2=fm1b)
model.comp(m1=fm1a,m2=fm1b,LRT=TRUE)
model.comp(m1=fm1a,m2=fm1b,LRT=TRUE,rdDF=TRUE)
##
## Attension:
## Fixed factors should be the same!
##
## Model LogL Npm AIC AIC.State
## 1 fm1b -868 6 1749
## 2 fm1a -867 7 1748 better
## -----------------------------
## Lower AIC is better model.
##
## Attension: Please check every asreml results' length is 43;
## if the length < 43, put the object at the end of Nml.
## In the present, just allow one object's length < 43.
## =====================================
## Likelihood ratio test (LRT) results:
##
## Model compared between fm2a -- fm2b :
## Model LogL Npm AIC Pr(>F) Sig.level
## 1 fm1b -868 6 1749
## 2 fm1a -867 7 1748 0.038 *
## ---------------
## Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
## =====================================
## Attension: Ddf=Ddf-0.5.
## When for corr model, against +/-1.
comparison among more than two models:
model.comp(Nml=c(fm3a,fm3b,fm1a,fm1b,fm1),mulM=TRUE)
model.comp(Nml=c(fm3a,fm3b,fm1a,fm1b,fm1),mulM=TRUE,LRT=TRUE)
## met.plot(): plots MET data
## met.corr(): calculate var/cov/corr from asreml MET factor analytic results
## met.biplot(): biplots MET factor analytic results from asreml
data(MET,package='AAfun')
## plot MET data -- example 1
# variable order: genotype,yield,site,row,col
MET2<-MET[,c(1,9,2,4:5)]
AAfun::met.plot(MET2)
## plot MET data -- example 2
MET3<-MET[,c(1,9,2,4:7)] # add variable order on MET2: Rep, Block
met.plot(MET3,"My met trials")
MET$yield<-0.01*MET$yield
#summary(MET$yield)
met1.asr<-asreml(yield~Loc, random=~ diag(Loc):Rep + Genotype:fa(Loc,2),
rcov=~ at(Loc):units,
data=MET, maxiter=50)
met2.asr<-asreml(yield~Loc, random=~ Genotype:fa(Loc,2),
rcov=~ at(Loc):ar1(Col):ar1(Row),
data=MET, maxiter=50)
met3.asr<-asreml(yield~Loc, random=~ Genotype:fa(Loc,3),
rcov=~ at(Loc):ar1(Col):ar1(Row),
data=MET, maxiter=50,trace=F)
## count var/cov/corr matrix, etc.
AAfun::met.corr(met2.asr, site=MET$Loc, faN=2, kn=2)
##
## Site cluster results:
## S1 S2 S3 S4 S5 S6
## 1 1 1 2 1 2
##
## Cov\Var\Corr matrix
## S1 S2 S3 S4 S5 S6
## S1 5.38 0.816 0.875 0.520 0.979 0.130
## S2 3.90 4.246 0.688 0.472 0.811 0.161
## S3 4.97 3.469 5.988 0.042 0.759 -0.366
## S4 2.84 2.288 0.243 5.535 0.682 0.914
## S5 4.31 3.169 3.525 3.044 3.598 0.328
## S6 0.32 0.353 -0.954 2.291 0.662 1.135
met.corr(met1.asr, site=MET$Loc, faN=2, kn=2)
met.corr(met3.asr, site=MET$Loc, faN=3, kn=2)
## biplot asreml-met results
AAfun::met.biplot(met2.asr, siteN=6, VarietyN=36, faN=2)
met.biplot(met3.asr, siteN=6, VarietyN=36, faN=3)
met.biplot(met2.asr, siteN=6, VarietyN=36, faN=2, dSco.u=1.8, dLam.u=0.8)
met.biplot(met2.asr, siteN=nlevels(MET$Loc), VarietyN=nlevels(MET$Genotype), faN=2)
# dLam.u -- least distance from center
# dSco.u -- least score of Variety breeding value
# if can not draw fig 3, try multiplying or being devided by 10 for aim trait data.