Phenotypes simulation, multitrait and Random Regression models with Asreml
Tom Rohmer (thanks to H. Chapuis and I. David for their help)
tom.rohmer@inrae.fr18/05/2022
library(mvtnorm) #Multivariate Gaussian simulation using a Cholesky decomposition
#library(copula) #copula simulation
loop=5
set.seed(loop) #Fixe the seed for the simulations.
#system("mkdir /travail/trohmer/testR")
REPTRA="/travail/trohmer/testR" #Working repertory
REPMOD="/modgen/trohmer/testR"1 Population definition
The population generated mimics the schema of (González-Diéguez et al. 2020), following (Rohmer, Ricard, and David 2022):
\[ \begin{matrix} & G_0 & \\ \hline \text{male} & & \text{female}\\ \hline 12 & \times & 204 \end{matrix} \quad \underset{\longrightarrow}{\text{12 off/litt.}} \quad \begin{matrix} & G_1 & \\ \hline \text{male} & & \text{female}\\ \hline 408 & & 2040\\ & \text{Random} & \\ \downarrow &\text{selection} & \downarrow\\ 12 & \times & 204 \end{matrix} \quad \underset{\longrightarrow}{\text{12 off/litt.}} \quad \begin{matrix} & G_2 & \\ \hline \text{male} & & \text{female}\\ \hline 408 & & 2040\\ & \text{Random} & \\ \downarrow &\text{selection} & \downarrow\\ 12 & \times & 204 \end{matrix} \underset{\longrightarrow}{\text{12 off/litt.}} \quad G_3 \]
NMaleG0=12;NFemaleG0=204 # number of males and females in G0
p<-8 # total number of generations
off<-12 # Number of offspring by female
d<-4 # number of generations with random selection
G<-NMale<-NFemale<-rep(0,p+1)
NMale[1]<-NMaleG0
NFemale[1]<-NFemaleG0
G[1] <- NMale[1]+NFemale[1]
for(i in 2:(p+1)){
NFemale[i] <-NFemale[1]
NMale[i] <-NMale[1]
G[i]<- G[i-1]+ NFemale[i-1]*off
}
G<-c(0,G)
N=G[p+2] # Total number of animals
Nd=G[d+2] # Animals of generation G1 - Gd
print(G) # number of cumulative animal by generations## [1] 0 216 2664 5112 7560 10008 12456 14904 17352 19800###Note that for the moment, only the Animals of generation G1 - Gd are considered\(G_4\) to \(G_8\) will not be considered here. That is to say, they are \(N_3=7560\) animal in the pedigree.
###Pedigree initialization
pedi<-cbind(1:N,rep(0,N),rep(0,N))
sex<-rep(0,N)
sex[1:NMaleG0]<-1 ##1 pour male, 0 pour female
index=0.83 ##proportion of phenotyped females by litters
offfem<-round(index*off)
offmal<-round((1-index)*off)
sex[(G[2]+1):N]<-rep(c(rep(1,offmal),rep(0,offfem)),NFemaleG0*p)
##That is 2 males 10 females
pedi<-as.data.frame(cbind(pedi,sex))
colnames(pedi)<-c("ANIMAL","SIR","DAM","SEX")2 Mandelian sampling
The breeding vectors are simulated by Mendelian sampling. First, the BVs are initialized:
##Genetic variance matrix for the RR model
sigma1=matrix(c(4,0.62,0.36,0.62,5,-0.28,0.36,-0.28,0.8),3,3)
sigma2=matrix(c(0.588,0.2,0.2,0.2,0.588,0.2,0.2,0.2,0.588),3,3)
sigma.mod2<-rbind(cbind(sigma1,sigma2),cbind(sigma2,sigma1))
##Genetic variance matrix for the multitrait model
sigma.mod1<-matrix(c(0.67, 0.39,0.39,0.67),2,2)
#Initalisation of the BVs
a.mod1<-rmvnorm(N,rep(0,2),sigma.mod1)
a.mod2<-rmvnorm(N,rep(0,6),sigma.mod2) For the founders, the BVs \({\bf a_i}\), \(i=1,\ldots,N_0\) follow a bivariate Gaussian distribution with covariance matrix \(G\) and are directly obtained from the initialization. For the others generations the vector \(\bf a_i\) of the BVs for the \(i\)th individual is
\[ \bf a_i = 0.5(\bf a_{i_S} + \bf a_{i_D}) + \bf M_{i}, \] where \(\bf M_{i}\) follows the bivariate Gaussian distribution with covariance matrix \(G/2\) and \(i_S\) and \(i_D\) stand for the indexes of the Sire and the Dam of the \(i\)th animal. Obviously the vectors \({\bf a}_i\) and \({\bf M}_i\) are independent with each other.
Thus, for the first generation,
k=1
setpm<-(G[k]+1):G[k+1] ##index for generation G0
set<-(G[k+1]+1):G[k+2] ##index for generation G1
pedpm<-pedi[setpm,]
m<-pedpm$ANIMAL[pedpm$SEX==1] ##male of generation G0
f<-pedpm$ANIMAL[pedpm$SEX==0] ##female of generation G0
ip<-rep(m,each=off)
im<-rep(f,each=off)
ii<-cbind(ip,im)
a.mod1[set,]<-(a.mod1[ii[,1],]+a.mod1[ii[,2],])/2 + a.mod1[set,]/sqrt(2) #BVs for the multitrait model
a.mod2[set,]<-(a.mod2[ii[,1],]+a.mod2[ii[,2],])/2 + a.mod2[set,]/sqrt(2) #BVs for the RR model
pedi[set,2:3]<-iiFor the generation \(G_2\) to \(G_d\), the reproducers were chosen at random (selection was made within the progeny of one male). Male and female breeders were randomly mated to produce the next generation, but in such a way, that the full/half siblings were not mated.
if(d>1)
for(k in 2:d){
setpm<-(G[k]+1):G[k+1] ##index of generation G_{k-1}
set<-(G[k+1]+1):G[k+2] ##index of generation G_{k}
pedpm<-pedi[setpm,]
m<-pedpm$ANIMAL[pedpm$SEX==1] ##males of generation G_{k-1}
f<-pedpm$ANIMAL[pedpm$SEX==0] ##females of generation G_{k-1}
##selection intra-pere at random
sire<-by(m,pedpm$SIR[pedpm$SEX==1],sample,size=1,replace=FALSE)
#1 male at random by family
dame<-by(f,pedpm$SIR[pedpm$SEX==0],sample,size=NFemaleG0/NMaleG0,replace=FALSE)
#17 female at random by family
sire<-c(sire)
###no half-sib are matted
ii<-NULL
it=rep(1:12,12)
for(ss in 1:NMaleG0){
im<-rep(dame[[ss]],each=off)
ip<-rep(c(sire[-ss],sire[it[(ss+1):(ss+floor(NFemaleG0/NMaleG0)-NMaleG0+1)]]),each=off)
#the 17 females are matted at random with the 11 males of the other family (6 males are two times matted, the choice of the 6 is taken by rotation)
ii<-rbind(ii,cbind(ip,im))
}
a.mod1[set,]<-(a.mod1[ii[,1],]+a.mod1[ii[,2],])/2 + a.mod1[set,]/sqrt(2)
a.mod2[set,]<-(a.mod2[ii[,1],]+a.mod2[ii[,2],])/2 + a.mod2[set,]/sqrt(2)
pedi[set,2:3]<-ii
}3 Multitrait model
The considered simulated model is the following
\[ \left\{\begin{matrix} {\bf y}_{1} = X\beta_1 + {\bf a}_{1} + \varepsilon_{1}\\ {\bf y}_{2} = X\beta_2 + {\bf a}_{2} + \varepsilon_{2}\\ \end{matrix}\right. \]
with
\[ \left( \begin{matrix} {\bf a}_{1} \\ {\bf a}_{2} \\ \end{matrix} \right) \sim \mathcal N (0, G\otimes A), \quad \text{with } G= \left( \begin{matrix} 0.67 & 0.39 \\ 0.39 & 0.67 \\ \end{matrix} \right), \] and \(A\) is the genetic relationship matrix (here of dimension \(N_3\times N_3\), that is to say, all the animals of the pedigree were been phenotyped). The distribution of the residuals is \[ \left( \begin{matrix} \varepsilon_{1} \\ \varepsilon_{2} \\ \end{matrix} \right) \sim \mathcal N (0, E\otimes I_{N_3}), \quad \text{with } E= \left( \begin{matrix} 1 & 0.89 \\ 0.89 & 1 \\ \end{matrix} \right). \]
We considered three categorical fixed effects with three, two and two modalities for each trait:
mu=c(2,3)
beta11=c(1,2,3)
beta12=c(2,4)
beta13=c(1,2)
beta21=c(3,2,4)
beta22=c(1,3)
beta23=c(2,1)
beta1<-c(mu[1],beta11,beta12,beta13)
beta2<-c(mu[2],beta21,beta22,beta23)
X11<-as.factor(sample(1:4,N,replace=TRUE))
X12<-as.factor(sample(1:3,N,replace=TRUE))
X13<-as.factor(sample(1:3,N,replace=TRUE))
X.1<-model.matrix(~X11+X12+X13)For a bivariate Gaussian vector \((Z_1,Z_2)\) with Kendall’s correlation \(\tau\), the correlation (Pearson’s sense) is
\[ corr(Z_1,Z_2)=\sin(\tau\times\pi/2 ). \] Hence a Kendall’s correlation of 0.7 leads to a Pearson’s correlation of 0.89.
vare1 <- vare2<-1
taue=0.7 #that corresponds to a Pearson's correlation of 0.89 for the Gaussian case
# parJoe<-iTau(joeCopula(),abs(taue))
#bJ<-mvdc(joeCopula(parJoe), margins = c("norm","norm"),paramMargins = list(list(sd = sqrt(vare1)),list(sd = sqrt(vare2))))
#EJ<-rMvdc(n*N,bJ) #For non-Gaussian residual
cove<-sin(pi/2 * taue)*sqrt(vare1*vare2)
Sigmae<-matrix(c(vare1,cove,cove,vare2),2,2)
E<-rmvnorm(N,rep(0,2),Sigmae)X.mod1 = cbind(rep(X.1%*%beta1),rep(X.1%*%beta2))[1:Nd,1:2]
e.mod1 = E[1:Nd,]
y.mod1<- X.mod1[1:Nd,]+a.mod1[1:Nd,]+e.mod1[1:Nd,]
data.mod1<-data.frame(ANIMAL=pedi$ANIMAL[1:G[d+2]],identif=pedi$ANIMAL[1:G[d+2]],fac1=X11[1:Nd],fac2=X12[1:Nd],fac3=X13[1:Nd],y=y.mod1[1:Nd,])
data.mod1$ANIMAL=as.factor(data.mod1$ANIMAL)
data.mod1$fac1 = as.factor(data.mod1$fac1)
data.mod1$fac2 = as.factor(data.mod1$fac2)
data.mod1$fac3 = as.factor(data.mod1$fac3)Model estimations were carried out using (Gilmour et al. 2015) and (Butler et al. 2017), as described below.
3.1 ASREML-R
library(asreml)## Le chargement a nécessité le package : Matrix## Online License checked out Wed Oct 5 10:09:35 2022print(head(data.mod1))## ANIMAL identif fac1 fac2 fac3 y.1 y.2
## 1 1 1 2 3 1 6.984589 9.287795
## 2 2 2 1 3 1 4.960819 5.963034
## 3 3 3 3 2 2 9.220426 8.554550
## 4 4 4 1 3 1 6.089765 5.777495
## 5 5 5 2 1 2 2.595887 6.702446
## 6 6 6 3 2 3 8.775754 6.565975##A inverse
ainv=ainverse(pedi[1:Nd,])
##Memory allocation
control <-asreml.options(workspace=400e6, pworkspace = 400e6)
##That corresponds to 3.2Gb
TOT1<-asreml(fixed=cbind(y.1,y.2)~trait*fac1 +trait*fac2 + trait*fac3,
random=~vm(ANIMAL,ainv):us(trait), residual = ~ units:us(trait),maxit=20,data=data.mod1,trace=T)## Online License checked out Wed Oct 5 10:09:36 2022
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Wed Oct 5 10:09:36 2022
## LogLik Sigma2 DF wall cpu
## 1 -11396.17 1.0 20000 10:09:38 0.1
## 2 -8494.95 1.0 20000 10:09:38 0.1
## 3 -8180.17 1.0 20000 10:09:39 0.1
## 4 -8148.88 1.0 20000 10:09:39 0.1
## 5 -8147.71 1.0 20000 10:09:39 0.1
## 6 -8147.71 1.0 20000 10:09:39 0.1#comp1=summary(TOT)$varcomp$component
#print(comp1)
print(summary(TOT1))## $call
## asreml(fixed = cbind(y.1, y.2) ~ trait * fac1 + trait * fac2 +
## trait * fac3, random = ~vm(ANIMAL, ainv):us(trait), residual = ~units:us(trait),
## data = data.mod1, maxit = 20, trace = T)
##
## $loglik
## [1] -8147.711
##
## $nedf
## [1] 20000
##
## $sigma
## [1] 1
##
## $varcomp
## component std.error z.ratio bound %ch
## vm(ANIMAL, ainv):trait!trait_y.1:y.1 0.6719329 0.05065792 13.26413 P 0
## vm(ANIMAL, ainv):trait!trait_y.2:y.1 0.3648393 0.04328855 8.42808 P 0
## vm(ANIMAL, ainv):trait!trait_y.2:y.2 0.6361512 0.04997204 12.73014 P 0
## units:trait!R 1.0000000 NA NA F 0
## units:trait!trait_y.1:y.1 0.9917101 0.03059609 32.41296 P 0
## units:trait!trait_y.2:y.1 0.8753418 0.02710326 32.29656 P 0
## units:trait!trait_y.2:y.2 0.9803098 0.03027363 32.38164 P 0
##
## $bic
## [1] 16354.84
## attr(,"parameters")
## [1] 6
##
## $aic
## [1] 16307.42
## attr(,"parameters")
## [1] 6
##
## attr(,"class")
## [1] "summary.asreml"gc()#garbage #free the memory## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2002963 107 4052670 216.5 2568265 137.2
## Vcells 4578868 35 10146329 77.5 8385019 64.0detach("package:asreml")##libere le tiquet3.2 Asreml Stand alone
cat(
{
"
!L !DOPATH $1
testmodele copule
anim !P
identif !P
fac1 !A
fac2 !A
fac3 !A
y1
y2
gen.ped !MAKE
data.dat !WORKSPACE 3200 !MAXIT 20
!PATH 1 ##structural version
y1 y2 ~ Trait.mu Trait.fac1 Trait.fac2 Trait.fac3 ,
!r Trait.anim
1 2 1
0 0 ID
Trait 0 US 1 0.1 1 !GPUP
Trait.anim 2
Trait 0 US 1 0.1 1 !GPUP
anim
!PATH 2 ##fonctionnal version
y1 y2 ~ Trait.mu Trait.fac1 Trait.fac2 Trait.fac3 ,
!r nrm(anim).us(Trait !GPUP)
residual id(units).us(Trait !GPUP)
"
},file=paste(REPTRA,"/t2_asreml.as",sep="")
)The lance.sh file is available at https://genoweb.toulouse.inra.fr/~trohmer/divers/
write.table(data.mod1,file=paste(REPTRA,"/data.dat",sep=""),row.names=FALSE,col.names=FALSE)
write.table(pedi[1:Nd,1:3],file=paste(REPTRA,"/gen.ped",sep=""),row.names=FALSE,col.names=FALSE)
setwd(REPTRA)
system(paste(REPMOD,"/lance.sh",sep=""))
##lance.sh contains (chmod a+x before the first using)
##
###!/bin/ksh
##
##unset DISPLAY
##asreml4 /travail/trohmer/testR/t2_asreml.as 2
##grep -w "^ Trait" t2_asreml.asr > asreml.param
##
## var=$(echo `awk -F' ' '{print $5}' asreml.param`)
##
## echo $var > param
#print(readLines(paste(REPTRA,"/t2_asreml.asr",sep="")))3.3 Comparaison SA, R and theoretical values
#library(stringr)
#RES=readLines(paste(REPTRA,"/t2_asreml.asr",sep=""))[c(69:71,74:76)]
#par_SA=as.numeric(str_sub(RES,51,59))
par_SA<-unlist(read.table(paste(REPTRA,"/param",sep="")))
par_R=summary(TOT1)$varcomp$component[c(5:7,1:3)]
theo<-c(1,round(cove,2),1,sigma.mod1[1,1],sigma.mod1[1,2],sigma.mod1[2,2])
TOT=cbind(par_SA,par_R,theo)
rownames(TOT)<-rownames(summary(TOT1)$varcomp[c(5:7,1:3),])
print(TOT)## par_SA par_R theo
## units:trait!trait_y.1:y.1 0.991711 0.9917101 1.00
## units:trait!trait_y.2:y.1 0.875343 0.8753418 0.89
## units:trait!trait_y.2:y.2 0.980310 0.9803098 1.00
## vm(ANIMAL, ainv):trait!trait_y.1:y.1 0.671931 0.6719329 0.67
## vm(ANIMAL, ainv):trait!trait_y.2:y.1 0.364838 0.3648393 0.39
## vm(ANIMAL, ainv):trait!trait_y.2:y.2 0.636150 0.6361512 0.67e.cora<-par_SA[5]/sqrt(par_SA[4]*par_SA[6])
e.core<-par_SA[2]/sqrt(par_SA[1]*par_SA[3])
e.herit1<-par_SA[4]/(par_SA[4]+par_SA[1])
e.herit2<-par_SA[6]/(par_SA[6]+par_SA[3])
cora<-theo[5]/sqrt(theo[4]*theo[6])
core<-theo[2]/sqrt(theo[1]*theo[3])
herit1<-theo[4]/(theo[4]+theo[1])
herit2<-theo[6]/(theo[6]+theo[3])
theo.gen.par<-c(cora,core,herit1,herit2)
est.gen.par<-c(e.cora,e.core,e.herit1,e.herit2)
par=cbind(est.gen.par,theo.gen.par)
colnames(par)<-c("est","theo");rownames(par)<-c("gen cor","res cor","herit1","hrit2")
print(par)## est theo
## gen cor 0.5580305 0.5820896
## res cor 0.8877772 0.8900000
## herit1 0.4038916 0.4011976
## hrit2 0.3935452 0.40119764 Random Regression model
The model is the following:
\[ \left\{\begin{matrix} {\bf y}_{1,t} = X\beta_1 + \sum_{k=0}^2 {\bf a}_{1,k}\phi^{(k)}(t) + \sum_{k=0}^2 {\bf e}_{1,k}\phi^{(k)}(t) + \varepsilon_{1,t}\\ {\bf y}_{2,t} = X\beta_2 + \sum_{k=0}^2 {\bf a}_{2,k}\phi^{(k)}(t) + \sum_{k=0}^2 {\bf e}_{2,k}\phi^{(k)}(t) + \varepsilon_{2,t}\\ \end{matrix}\right. , \]
where \(\phi^{(k)}\) are Legendre’s polynomials of degree \(k\) defined as
legendre<-function(x,deg,max=100,min=1){
t<-2*(x-min)/(max-min) - 1
if(deg==0) leg=rep(0.7071068,length(x))
if(deg==1) leg=1.224745*t
if(deg==2) leg=-0.7905694+2.371708*t^2
leg
}It’s considered that each animals are measured at time \(t=1,\ldots,100\). For \(k=0,1,2\), \(a_{j,k}\) are the vectors of BVs for the trait \(j=1,2\). Their distributions are \[ \left( \begin{matrix} {\bf a}_{1,1} \\ {\bf a}_{1,2} \\ {\bf a}_{1,3} \\ {\bf a}_{2,1} \\ {\bf a}_{2,2} \\ {\bf a}_{2,3} \\ \end{matrix} \right) \sim \mathcal N (0, G\otimes A), \quad \text{with } G= \left( \begin{matrix} G_1 & G_{12}\\ G_{12} & G_2 \end{matrix} \right), \] \[ G_1 = G_2 = \left( \begin{matrix} 4.00 & 0.62 & 0.36 \\ 0.62 & 5.00 & -0.28 \\ 0.36 & -0.28 & 0.80 \\ \end{matrix} \right) \quad \text{and } G_{12} = \left( \begin{matrix} 0.59 & 0.20 & 0.20 \\ 0.20 & 0.59 & 0.20 \\ 0.20 & 0.20 & 0.59\\ \end{matrix} \right), \] and \(A\) is the genetic relationship matrix.
For \(k=0,1,2\), \(e_{j,k}\) are the vectors of permanent effects for the trait \(j=1,2\). Their distributions are \[ \left( \begin{matrix} {\bf e}_{1,1} \\ {\bf e}_{1,2} \\ {\bf e}_{1,3} \\ {\bf e}_{2,1} \\ {\bf e}_{2,2} \\ {\bf e}_{2,3} \\ \end{matrix} \right) \sim \mathcal N (0, P\otimes I_N), \quad \text{with } P= \left(\begin{matrix} P_1 & G_{12}\\ G_{12} & P_2 \end{matrix}\right) \]
\[ P_1=P_2= \left( \begin{matrix} 11.00 & 0.00 & 0.00 \\ 0.00 & 8.00 & 0.00 \\ 0.00 & 0.00 & 0.60 \\ \end{matrix} \right) \]
# Permanent variance matrix for the RR model
Sigmae=matrix(c(11,0,0,0,8,0,0,0,0.6),3,3)
sigmaP<-rbind(cbind(Sigmae,sigma2),cbind(sigma2,Sigmae))
e<-rmvnorm(N,rep(0,6),sigmaP) ##permanant partFinally the residual vectors \(\varepsilon_{t,j}=(\varepsilon_{t,j,1},\ldots,\varepsilon_{t,j,n})\), \(t=1,\ldots,T\) are independent and their distribution is \[ (\varepsilon_{t,1},\varepsilon_{t,2})\sim N\left(0,\left(\begin{matrix}1 & 0.89\\ 0.89 & 1\end{matrix}\right)\otimes I_{N_3}\right). \]
vare1 <- vare2<-1
taue=0.7 #that corresponds to a Pearson's correlation of 0.85
n<-100 #number of observations by animals
#parJoe<-iTau(joeCopula(),abs(taue))
#bJ<-mvdc(joeCopula(parJoe), margins = c("norm","norm"),paramMargins = list(list(sd = sqrt(vare1)),list(sd = sqrt(vare2))))
#EJ<-rMvdc(n*N,bJ) #For non-Gaussian residual
cove<-sin(pi/2 * taue)*sqrt(vare1*vare2)
Sigmares<-matrix(c(vare1,cove,cove,vare2),2,2)
E2<-rmvnorm(n*N,rep(0,2),Sigmares)4.1 Longitudinal phenotype construction
U2<-P2<-X2<-matrix(rep(0,2*n*Nd),n*Nd,2) #initialization
phi=cbind(legendre(1:n,deg=0,max=n),legendre(1:n,deg=1,max=n),legendre(1:n,deg=2,max=n))
U2[,1]=kronecker(phi[,1],a.mod2[1:Nd,1])+kronecker(phi[,2],a.mod2[1:Nd,2])+kronecker(phi[,3],a.mod2[1:Nd,3])
U2[,2]=kronecker(phi[,1],a.mod2[1:Nd,4])+kronecker(phi[,2],a.mod2[1:Nd,5])+kronecker(phi[,3],a.mod2[1:Nd,6])
P2[,1]=kronecker(phi[,1],e[1:Nd,1],)+kronecker(phi[,2],e[1:Nd,2])+kronecker(phi[,3],e[1:Nd,3])
P2[,2]=kronecker(phi[,1],e[1:Nd,4])+kronecker(phi[,2],e[1:Nd,5])+kronecker(phi[,3],e[1:Nd,6])
#
X.1<-model.matrix(~X11+X12+X13)
X<-cbind(rep(X.1%*%beta1,each=n),rep(X.1%*%beta2,each=n))
X2[,1]<-as.vector(matrix(X[1:(n*Nd),1],Nd,n,byrow=TRUE))
X2[,2]<-as.vector(matrix(X[1:(n*Nd),2],Nd,n,byrow=TRUE))
y2<-X2[1:(n*Nd),]+U2[1:(n*Nd),]+P2[1:(n*Nd),]+E2[1:(n*Nd),]
data.mod2<-data.frame(ANIMAL=rep(pedi$ANIMAL[1:G[d+2]],n),identif=rep(pedi$ANIMAL[1:G[d+2]],n),fac1=rep(X11[1:Nd],n),fac2=rep(X12[1:Nd],n),fac3=rep(X13[1:Nd],n),time=rep(1:n,each=Nd),y=y2[1:(n*Nd),])4.2 Univariate RR model
Consider the univariate RR model
\[ y_{n,1,t} = X\beta_1 + \sum_{k=0}^2 a_{1,k}\phi^{(k)}(t) + \sum_{k=0}^2 e_{1,k}\phi^{(k)}(t) + \varepsilon_{1,t}\\ \]
4.2.1 ASREML-R
library(asreml)
data.mod2$ANIMAL=as.factor(data.mod2$ANIMAL)
data.mod2$fac1 = as.factor(data.mod2$fac1)
data.mod2$fac2 = as.factor(data.mod2$fac2)
data.mod2$fac3 = as.factor(data.mod2$fac3)
print(head(data.mod2))## ANIMAL identif fac1 fac2 fac3 time y.1 y.2
## 1 1 1 2 3 1 1 5.986050 19.337369
## 2 2 2 1 3 1 1 5.843148 11.038674
## 3 3 3 3 2 2 1 13.417274 13.034949
## 4 4 4 1 3 1 1 15.361263 11.193494
## 5 5 5 2 1 2 1 5.575422 25.832604
## 6 6 6 3 2 3 1 10.088774 8.925774##A inverse
ainv=ainverse(pedi[1:Nd,])
##Memory allocation
control <-asreml.options(workspace=400e6, pworkspace = 400e6)
##That corresponds to 3.2Gb
TOTJ1<-asreml(fixed=y.1~fac1 +fac2 + fac3,
random=~us(leg(time,2)):vm(ANIMAL,ainv)+diag(leg(time,2)):id(ANIMAL),
maxit=20,data=data.mod2,trace=T,asreml.options=control)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Wed Oct 5 10:10:06 2022
## LogLik Sigma2 DF wall cpu
## 1 -966252.0 2.36836 1000792 10:10:13 3.8
## 2 -829977.1 1.77540 1000792 10:10:16 2.7
## 3 -713116.6 1.37931 1000792 10:10:18 2.8
## 4 -643129.4 1.17818 1000792 10:10:21 2.7
## 5 -604988.6 1.07362 1000792 10:10:24 2.7
## 6 -590072.3 1.02934 1000792 10:10:27 2.8
## 7 -585038.3 1.00926 1000792 10:10:29 2.8 (1 restrained)
## 8 -583851.2 1.00063 1000792 10:10:32 2.8
## 9 -583715.9 0.99786 1000792 10:10:35 2.8
## 10 -583712.6 0.99742 1000792 10:10:38 2.8
## 11 -583712.6 0.99741 1000792 10:10:41 2.9#comp1=summary(TOTJ1)$varcomp$component
#print(comp1)
summary(TOTJ1)## $call
## asreml(fixed = y.1 ~ fac1 + fac2 + fac3, random = ~us(leg(time,
## 2)):vm(ANIMAL, ainv) + diag(leg(time, 2)):id(ANIMAL), data = data.mod2,
## maxit = 20, trace = T, asreml.options = control)
##
## $loglik
## [1] -583712.6
##
## $nedf
## [1] 1000792
##
## $sigma
## [1] 0.9987053
##
## $varcomp
## component std.error
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order0:order0 3.9990760 0.393068819
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order0 0.4788418 0.181017499
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order1 4.8808825 0.391615466
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order0 0.4208063 0.062419339
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order1 -0.1606951 0.057739362
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order2 0.8395831 0.051245990
## leg(time, 2):ANIMAL!leg(time, 2)_order0 11.0752994 0.274668619
## leg(time, 2):ANIMAL!leg(time, 2)_order1 7.9255929 0.240405509
## leg(time, 2):ANIMAL!leg(time, 2)_order2 0.5891769 0.027091094
## units!R 0.9974123 0.001431628
## z.ratio bound %ch
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order0:order0 10.173984 P 0
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order0 2.645279 P 0
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order1 12.463457 P 0
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order0 6.741601 P 0
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order1 -2.783112 P 0
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order2 16.383391 P 0
## leg(time, 2):ANIMAL!leg(time, 2)_order0 40.322405 P 0
## leg(time, 2):ANIMAL!leg(time, 2)_order1 32.967601 P 0
## leg(time, 2):ANIMAL!leg(time, 2)_order2 21.747991 P 0
## units!R 696.697920 P 0
##
## $bic
## [1] 1167563
## attr(,"parameters")
## [1] 10
##
## $aic
## [1] 1167445
## attr(,"parameters")
## [1] 10
##
## attr(,"class")
## [1] "summary.asreml"gc()#garbage #free the memory## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 3055740 163.2 5686501 303.7 3196433 170.8
## Vcells 41245802 314.7 69997562 534.1 58264601 444.6detach("package:asreml")##libere le tiquet4.2.2 ASREML-stand Alone
cat(
{
"
!L !DOPATH $1
testmodele RR
anim !P
identif !I
fac1 !A
fac2 !A
fac3 !A
time !I
y1
y2
gen.ped !MAKE
data.dat !WORKSPACE 3200 !MAXIT 15
!PATH 2 #fonctionnal, first trait only
y1 ~ mu fac1 fac2 fac3 ,
!r str(leg(time,2).anim us(3 !GPUPUUP).nrm(anim)),
str(leg(time,2).identif us(3 !INIT 1 0 1 0 0 1 !GPFPFFP).id(identif))
residual idv(units)
!PATH 1 #fonctionnal, two traits simultaneously
y1 y2 ~ Trait.mu Trait.fac1 Trait.fac2 Trait.fac3 ,
!r str(Trait.leg(time,2).anim, us(6 !INIT 1 0.1 1 0.1 0.1 1 0 0 0 1 0 0 0 0.1 1 0 0 0 0.1 0.1 1 !GPUPUUPFFFPFFFUPFFFUUP).nrm(anim)) str(Trait.leg(time,2).identif, us(6 !INIT 1 0 1 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1 !GPFPFFPFFFPFFFFPFFFFFP).id(identif))
residual id(units).us(Trait !INIT 1 0 1 !GPFP)
"
},file=paste(REPTRA,"/t2_asreml.as",sep="")
)The lanceRR.sh file is available at https://genoweb.toulouse.inra.fr/~trohmer/divers/
write.table(data.mod2,file=paste(REPTRA,"/data.dat",sep=""),row.names=FALSE,col.names=FALSE)
write.table(pedi[1:Nd,1:3],file=paste(REPTRA,"/gen.ped",sep=""),row.names=FALSE,col.names=FALSE)
setwd(REPTRA)
system(paste(REPMOD,"/lanceRR.sh",sep=""))
##lanceRR.sh contains (chmod a+x before the first using)
##
###!/bin/ksh
##
##unset DISPLAY
##asreml4 /travail/trohmer/testR/t2_asreml.as 2
##grep -w "^ Residual" t2_asreml.asr > asreml.param
##grep -w "^ 3" t2_asreml.asr >> asreml.param
##
##var=$(echo `awk -F' ' '{print $5}' asreml.param`)
##
##echo $var > param
#print(readLines(paste(REPTRA,"/t2_asreml.asr",sep="")))4.2.3 Comparaison SA, R and theoretical values
#library(stringr)
#RES=readLines(paste(REPTRA,"/t2_asreml.asr",sep=""))[86:101]
#par_SA=as.numeric(str_sub(RES[c(4,6:11,14:16)],51,59))
par_SA<-unlist(read.table(paste(REPTRA,"/param",sep="")))
par_SA<-par_SA[par_SA!=0]
par_R=c(summary(TOTJ1)$sigma,summary(TOTJ1)$varcomp$component[1:9])
theo<-c(1,sigma1[1,1],sigma1[1,2],sigma1[2,2],sigma1[3,1],sigma1[3,2],sigma1[3,3],diag(sigmaP)[1:3])
TOT=cbind(par_SA,par_R,theo)
rownames(TOT)<-c("sigma",rownames(summary(TOTJ1)$varcomp[1:9,]))
print(TOT)## par_SA par_R
## sigma 0.997412 0.9987053
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order0:order0 4.009450 3.9990760
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order0 0.480084 0.4788418
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order1 4.893540 4.8808825
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order0 0.421899 0.4208063
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order1 -0.161113 -0.1606951
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order2 0.841760 0.8395831
## leg(time, 2):ANIMAL!leg(time, 2)_order0 11.104000 11.0752994
## leg(time, 2):ANIMAL!leg(time, 2)_order1 7.946160 7.9255929
## leg(time, 2):ANIMAL!leg(time, 2)_order2 0.590706 0.5891769
## theo
## sigma 1.00
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order0:order0 4.00
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order0 0.62
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order1 5.00
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order0 0.36
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order1 -0.28
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order2 0.80
## leg(time, 2):ANIMAL!leg(time, 2)_order0 11.00
## leg(time, 2):ANIMAL!leg(time, 2)_order1 8.00
## leg(time, 2):ANIMAL!leg(time, 2)_order2 0.604.3 Multitrait RR model
4.3.1 Asreml Stand-Alone
cat(
{
"
!L !DOPATH $1
testmodele RR
anim !P
identif !A 7560
fac1 !A
fac2 !A
fac3 !A
time !I
y1
y2
gen.ped !MAKE
data.dat !WORKSPACE 3200 !MAXIT 30
!PATH 1 # poly order 1 structural version
y1 y2 ~ Trait.mu Trait.fac1 Trait.fac2 Trait.fac3 ,
!r Trait.leg(time,1).anim Trait.leg(time,1).identif
0 0 2
Trait.leg(time,1).anim 2
4 0 US 1 0.1 1 0.1 0.1 1 0.1 0.1 0.1 1 !GPUPUUPUUUP
anim 0 ainv
Trait.leg(time,1).identif 2
4 0 US 1 0.1 1 0.1 0.1 1 0.1 0.1 0.1 1 !GPUPUUPUUUP
identif 0 ID
!PATH 3 # poly order 2 structural version
y1 y2 ~ Trait.mu Trait.fac1 Trait.fac2 Trait.fac3 ,
!r Trait.leg(time,2).anim Trait.leg(time,2).identif
0 0 2
Trait.leg(time,2).anim 2
6 0 US 1 0.1 1 0.1 0.1 1 0.1 0.1 0.1 1 0.1 0.1 0.1 0.1 1 0.1 0.1 0.1 0.1 0.1 1 !GPUPUUPUUUPUUUUPUUUUUP
anim 0 ainv
Trait.leg(time,2).identif 2
6 0 US 1 0 1 0 0 1 0.1 0.1 0.1 1 0.1 0.1 0.1 0 1 0.1 0.1 0.1 0 0 1 !GPFPFFPUUUPUUUFPUUUFFP
identif 0 ID
!PATH 2 #fonctionnal version
y1 y2 ~ Trait.mu Trait.fac1 Trait.fac2 Trait.fac3 ,
!r str(Trait.leg(time,2).anim, us(6 !GPUPUUPUUUPUUUUPUUUUUP).nrm(anim)) str(Trait.leg(time,2).identif, us(6 !INIT 1 0 1 0 0 1 0.1 0.1 0.1 1 0.1 0.1 0.1 0 1 0.1 0.1 0.1 0 0 1 !GPFPFFPUUUPUUUFPUUUFFP).id(identif))
residual id(units).us(Trait !GPUP)
"
},file=paste(REPTRA,"/t2_asreml.as",sep="")
)The lanceRR2.sh file is available at https://genoweb.toulouse.inra.fr/~trohmer/divers/
write.table(data.mod2,file=paste(REPTRA,"/data.dat",sep=""),row.names=FALSE,col.names=FALSE)
write.table(pedi[1:Nd,1:3],file=paste(REPTRA,"/gen.ped",sep=""),row.names=FALSE,col.names=FALSE)
setwd(REPTRA)
system(paste(REPMOD,"/lanceRRb.sh",sep=""))
##lanceRR2.sh contains (chmod a+x before the first using)
##
###!/bin/ksh
##
##unset DISPLAY
##asreml4 /travail/trohmer/testR/t2_asreml.as 2
##
##grep -w "^ Trait" t2_asreml.asr > asreml.param
##grep -w "^ 6" t2_asreml.asr >> asreml.param
##
##var=$(echo `awk -F' ' '{print $5}' asreml.param`)
##
##echo $var > paramlibrary(stringr)
#RES=readLines(paste(REPTRA,"/t2_asreml.asr",sep=""))[c(122:124,126:146,149:169)]
#par_SA=as.numeric(str_sub(RES,51,63))
par_SA<-unlist(read.table(paste(REPTRA,"/param",sep="")))
theo<-NULL
for(j in 1:2)
theo<-c(theo,Sigmares[j,1:j]);
for(j in 1:6)
theo<-c(theo,sigma.mod2[j,1:j])
for(j in 1:6)
theo<-c(theo,sigmaP[j,1:j])
TOT=cbind(par_SA,par_theo=round(theo,4))
#rownames(TOT)=str_sub(RES,1,21)
#TOT<-cbind(as.numeric(str_sub(RES,32,32)),as.numeric(str_sub(RES,35,35)),TOT)
c1=c(1,2,2,rep(c(1,rep(2,2),rep(3,3),rep(4,4),rep(5,5),rep(6,6)),2))
c2=c(1,1,2,rep(c(1,1:2,1:3,1:4,1:5,1:6),2))
TOT<-cbind(c1,c2,TOT)
rownames(TOT)<-c(rep("residual",3),rep("var.gen",21),rep("var.per",21))
print(TOT)## c1 c2 par_SA par_theo
## residual 1 1 0.997412 1.000
## residual 2 1 0.887911 0.891
## residual 2 2 0.996295 1.000
## var.gen 1 1 3.917390 4.000
## var.gen 2 1 0.483428 0.620
## var.gen 2 2 4.904660 5.000
## var.gen 3 1 0.419473 0.360
## var.gen 3 2 -0.163037 -0.280
## var.gen 3 3 0.854341 0.800
## var.gen 4 1 0.696031 0.588
## var.gen 4 2 0.413445 0.200
## var.gen 4 3 0.169668 0.200
## var.gen 4 4 3.640230 4.000
## var.gen 5 1 0.156216 0.200
## var.gen 5 2 0.725747 0.588
## var.gen 5 3 0.154820 0.200
## var.gen 5 4 0.404134 0.620
## var.gen 5 5 5.105330 5.000
## var.gen 6 1 0.226872 0.200
## var.gen 6 2 0.260894 0.200
## var.gen 6 3 0.638295 0.588
## var.gen 6 4 0.331169 0.360
## var.gen 6 5 -0.299869 -0.280
## var.gen 6 6 0.839051 0.800
## var.per 1 1 11.124100 11.000
## var.per 2 1 0.000000 0.000
## var.per 2 2 7.913490 8.000
## var.per 3 1 0.000000 0.000
## var.per 3 2 0.000000 0.000
## var.per 3 3 0.581335 0.600
## var.per 4 1 0.355327 0.588
## var.per 4 2 0.153103 0.200
## var.per 4 3 0.213919 0.200
## var.per 4 4 11.433200 11.000
## var.per 5 1 0.432237 0.200
## var.per 5 2 0.534428 0.588
## var.per 5 3 0.206172 0.200
## var.per 5 4 0.000000 0.000
## var.per 5 5 8.028410 8.000
## var.per 6 1 0.199242 0.200
## var.per 6 2 0.195975 0.200
## var.per 6 3 0.573816 0.588
## var.per 6 4 0.000000 0.000
## var.per 6 5 0.000000 0.000
## var.per 6 6 0.590750 0.6004.3.2 plot of the estimated heritabilities
phi=cbind(legendre(1:n,deg=0,max=n),legendre(1:n,deg=1,max=n),legendre(1:n,deg=2,max=n))
#sU=4*legendre(1:n,deg=0,max=n)^2+5*legendre(1:n,deg=1,max=n)^2+0.8*legendre(1:n,deg=2,max=n)^2 + 2*0.62*legendre(1:n,deg=0,max=n)*legendre(1:n,deg=1,max=n) + 2*0.36*legendre(1:n,deg=0,max=n)*legendre(1:n,deg=2,max=n) -2*0.280*legendre(1:n,deg=1,max=n)*legendre(1:n,deg=2,max=n)
sU=diag(phi%*%sigma.mod2[1:3,1:3]%*%t(phi))
#sP=11*legendre(1:n,deg=0,max=n)^2+8*legendre(1:n,deg=1,max=n)^2+0.6*legendre(1:n,deg=2,max=n)^2
sP=diag(phi%*%sigmaP[1:3,1:3]%*%t(phi))
h= (sU)/(1+sU+sP)
est.res<-par_SA[1:3]
est.gen<-par_SA[4:24]
est.per<-par_SA[25:45]
Sigmag.est<-Sigmap.est<-matrix(rep(0,36),6,6)
estg<-est.gen
estp<-est.per
for(j in 1:6){
Sigmag.est[j,1:j]<-estg[1:j]
Sigmag.est[1:j,j]<-estg[1:j]
Sigmap.est[j,1:j]<-estp[1:j]
Sigmap.est[1:j,j]<-estp[1:j]
estg<-estg[-c(1:j)]
estp<-estp[-c(1:j)]
}
sU.est.1<-diag(phi%*%Sigmag.est[1:3,1:3]%*%t(phi))
sU.est.2<-diag(phi%*%Sigmag.est[4:6,4:6]%*%t(phi))
sP.est.1<-diag(phi%*%Sigmap.est[1:3,1:3]%*%t(phi))
sP.est.2<-diag(phi%*%Sigmap.est[4:6,4:6]%*%t(phi))
h.est.1<-sU.est.1/(sP.est.1 + sU.est.1 + est.res[1])
h.est.2<-sU.est.2/(sP.est.2 + sU.est.2 + est.res[3])
plot(sU,type='l',ylab="variance of the genetic part",lty=2,xlab="time")
lines(sU.est.2,col="red")
lines(sU.est.1,col="blue")
plot(sP,type='l',ylab="variance of the permanent part",lty=2,xlab="time")
lines(sP.est.2,col="red")
lines(sP.est.1,col="blue")
plot(h,type='l',ylab="heritability",lty=2,xlab="time")
lines(h.est.2,col="red")
lines(h.est.1,col="blue")
RMSE1=sqrt(sum((h.est.1 - h)^2))
RMSE2=sqrt(sum((h.est.2 - h)^2))
print(cbind(RMSE1,RMSE2))## RMSE1 RMSE2
## [1,] 0.075288 0.13836394.3.3 plot of the estimated correlation
corp_th=diag(phi%*%sigmaP[4:6,1:3]%*%t(phi))/sqrt(diag(phi%*%sigmaP[1:3,1:3]%*%t(phi))*diag(phi%*%sigmaP[4:6,4:6]%*%t(phi)))
corpe=diag(phi%*%Sigmap.est[4:6,1:3]%*%t(phi))/sqrt(diag(phi%*%Sigmap.est[1:3,1:3]%*%t(phi))*diag(phi%*%Sigmap.est[4:6,4:6]%*%t(phi)))
plot(corp_th,type='l',ylab="trait correlation of the permanant part",lty=2,xlab="time",ylim=c(min(corp_th,corpe),max(corp_th,corpe)))
lines(corpe,col="red")
corg_th=diag(phi%*%sigma.mod2[4:6,1:3]%*%t(phi))/sqrt(diag(phi%*%sigma.mod2[1:3,1:3]%*%t(phi))*diag(phi%*%sigma.mod2[4:6,4:6]%*%t(phi)))
corge=diag(phi%*%Sigmag.est[4:6,1:3]%*%t(phi))/sqrt(diag(phi%*%Sigmag.est[1:3,1:3]%*%t(phi))*diag(phi%*%Sigmag.est[4:6,4:6]%*%t(phi)))
plot(corg_th,type='l',ylab="trait correlation of the genetic part",lty=2,xlab="time",ylim=c(min(corg_th,corge),max(corg_th,corge)))
lines(corge,col="red")
References
Butler, DG, BR Cullis, AR Gilmour, BJ Gogel, and R Thompson. 2017. “ASReml-r Reference Manual Version 4.” VSN International Ltd, Hemel Hempstead, Hp1 1es, UK. https://asreml.kb.vsni.co.uk/wp-content/uploads/sites/3/ASReml-R-Reference-Manual-4.pdf.
Gilmour, AR, BJ Gogel, BR Cullis, SJa Welham, and R Thompson. 2015. “ASReml User Guide Release 4.1 Structural Specification.” Hemel Hempstead: VSN International Ltd. https://asreml.kb.vsni.co.uk/wp-content/uploads/sites/3/ASReml-4.2-Functional-Specification.pdf.
González-Diéguez, David, Llibertat Tusell, Alban Bouquet, Andres Legarra, and Zulma G Vitezica. 2020. “Purebred and Crossbred Genomic Evaluation and Mate Allocation Strategies to Exploit Dominance in Pig Crossbreeding Schemes.” G3 (Bethesda) 10: 2829–41.
Rohmer, Tom, Anne Ricard, and Ingrid David. 2022. “Copula Miss-Specification in REML Multivariate Genetic Animal Model Estimation.” Genetics Selection Evolution 54 (1): 1–12.