Mixed models are widely used in animal genetics to deal with genetic and environmental effects (see for example (Mrode 2014)). In practice, these parameters are frequently estimated using restricted maximum likelihood (REML) (Patterson and Thompson 1971). In the following sections, the simulation of an animal population is described then the code to fit the animal model using Asreml (Gilmour et al. 2015) and Asremel-R (Butler et al. 2017) are given for standard mixed model, multivariate mixed model and random regression model.
1 Population definition
The population generated mimics the schema of (González-Diéguez et al. 2020), following (Rohmer, Ricard, and David 2022). The reproducers are chosen at random. That is to say no strategy of selection is considered here.
\[ \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}{\ldots} \quad G_d \]
NMaleG0=12;NFemaleG0=204 # number of males and females in G0
p<-4 # total number of generations (with random selection of the reproducers)
off<-12 # Number of offspring by female
G<-NMale<-NFemale<-rep(0,p)
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 of generation G1 - Gp
print(G) # number of cumulative animals by generations## [1] 0 216 2664 5112 7560 10008###Note that for the moment, only the Animals of generation G1 - Gd are consideredThey are \(N_4=10008\) animal in the considered pedigree.
###Pedigree initialization
pedi<-cbind(1:N,rep(0,N),rep(0,N))
sex<-rep(0,N)
sex[1:NMaleG0]<-1 ##1 stands for male, 0 for 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, with our index, 2 males 10 females by litters
pedi<-as.data.frame(cbind(pedi,sex))
colnames(pedi)<-c("ANIMAL","SIR","DAM","SEX")2 Mendelian sampling
The breeding vectors are simulated by Mendelian sampling. First, the BVs are initialized:
library(mvtnorm) #Multivariate Gaussian simulation using a Cholesky decomposition
#library(copula) #copula simulation
loop=5
set.seed(loop) #Fixe the seed for the simulations.
REPTRA="/travail/trohmer/testR" #Working repertory
REPMOD="/modgen/trohmer/testR"
##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.59,0.59,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)
#BVs for the multitrait model
a.mod1[set,]<-(a.mod1[ii[,1],]+a.mod1[ii[,2],])/2 + a.mod1[set,]/sqrt(2)
#BVs for the RR model
a.mod2[set,]<-(a.mod2[ii[,1],]+a.mod2[ii[,2],])/2 + a.mod2[set,]/sqrt(2)
pedi[set,2:3]<-iiAs mentioned before, the reproducers were chosen at random (particularly 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(p>1)
for(k in 2:p){
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
#i.e. 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.58 \\ 0.58 & 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
cove<-sin(pi/2 * taue)*sqrt(vare1*vare2)
Sigmae<-matrix(c(vare1,cove,cove,vare2),2,2)
E<-rmvnorm(N,rep(0,2),Sigmae)
#For non-Gaussian residual (copula based-residual),
#we can use the following R code using the R package Copula
# 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) X.mod1 = cbind(rep(X.1%*%beta1),rep(X.1%*%beta2))[1:N,1:2]
e.mod1 = E[1:N,]
y.mod1<- X.mod1[1:N,]+a.mod1[1:N,]+e.mod1[1:N,]
data.mod1<-data.frame(ANIMAL=pedi$ANIMAL[1:G[p+2]],identif=pedi$ANIMAL[1:G[p+2]],
fac1=X11[1:N],fac2=X12[1:N],fac3=X13[1:N],y=y.mod1[1:N,])
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 Mar 22 16:15:04 2023print(head(data.mod1))## ANIMAL identif fac1 fac2 fac3 y.1 y.2
## 1 1 1 2 2 1 4.754244 7.679623
## 2 2 2 4 1 1 6.967537 8.844630
## 3 3 3 4 2 2 8.187719 8.764477
## 4 4 4 2 1 3 3.395347 6.351866
## 5 5 5 2 1 1 2.212313 5.288640
## 6 6 6 3 2 3 7.745904 5.504286##A inverse
ainv=ainverse(pedi[1:N,])
##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)## Warning in na.action$y == "include" && any(is.na(data[, lhs, with = FALSE])):
## ‘length(x) = 3 > 1’ dans la conversion automatique vers ‘logical(1)’## Online License checked out Wed Mar 22 16:15:07 2023
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Wed Mar 22 16:15:07 2023
## LogLik Sigma2 DF wall cpu
## 1 -9965.699 1.0 20000 16:15:12 0.1
## 2 -6605.897 1.0 20000 16:15:12 0.1
## 3 -5995.110 1.0 20000 16:15:12 0.1
## 4 -5763.125 1.0 20000 16:15:12 0.1
## 5 -5710.002 1.0 20000 16:15:12 0.1
## 6 -5708.902 1.0 20000 16:15:12 0.1
## 7 -5708.901 1.0 20000 16:15:12 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] -5708.901
##
## $nedf
## [1] 20000
##
## $sigma
## [1] 1
##
## $varcomp
## component std.error z.ratio bound %ch
## vm(ANIMAL, ainv):trait!trait_y.1:y.1 0.6783700 0.05215561 13.00666 P 0
## vm(ANIMAL, ainv):trait!trait_y.2:y.1 0.5942371 0.04908439 12.10644 P 0
## vm(ANIMAL, ainv):trait!trait_y.2:y.2 0.6676145 0.05202522 12.83252 P 0
## units:trait!R 1.0000000 NA NA F 0
## units:trait!trait_y.1:y.1 0.9756715 0.03113778 31.33402 P 0
## units:trait!trait_y.2:y.1 0.8712543 0.02940786 29.62658 P 0
## units:trait!trait_y.2:y.2 0.9884597 0.03123297 31.64795 P 0
##
## $bic
## [1] 11477.22
## attr(,"parameters")
## [1] 6
##
## $aic
## [1] 11429.8
## attr(,"parameters")
## [1] 6
##
## attr(,"class")
## [1] "summary.asreml"gc()#garbage #free the memory## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2087796 111.6 4169358 222.7 2761571 147.5
## Vcells 4464871 34.1 10146329 77.5 8366880 63.9detach("package:asreml")3.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="")
)write.table(data.mod1,file=paste(REPTRA,"/data.dat",sep=""),row.names=FALSE,col.names=FALSE)
write.table(pedi[1:N,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.975675 0.9756715 1.00
## units:trait!trait_y.2:y.1 0.871258 0.8712543 0.89
## units:trait!trait_y.2:y.2 0.988463 0.9884597 1.00
## vm(ANIMAL, ainv):trait!trait_y.1:y.1 0.678363 0.6783700 0.67
## vm(ANIMAL, ainv):trait!trait_y.2:y.1 0.594230 0.5942371 0.59
## vm(ANIMAL, ainv):trait!trait_y.2:y.2 0.667608 0.6676145 0.674 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
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)
#For non-Gaussian residuals,
#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) 4.1 Longitudinal phenotype construction
U2<-P2<-X2<-matrix(rep(0,2*n*N),n*N,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:N,1])+
kronecker(phi[,2],a.mod2[1:N,2])+
kronecker(phi[,3],a.mod2[1:N,3])
U2[,2]=kronecker(phi[,1],a.mod2[1:N,4])+
kronecker(phi[,2],a.mod2[1:N,5])+
kronecker(phi[,3],a.mod2[1:N,6])
P2[,1]=kronecker(phi[,1],e[1:N,1],)+
kronecker(phi[,2],e[1:N,2])+
kronecker(phi[,3],e[1:N,3])
P2[,2]=kronecker(phi[,1],e[1:N,4])+
kronecker(phi[,2],e[1:N,5])+
kronecker(phi[,3],e[1:N,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*N),1],N,n,byrow=TRUE))
X2[,2]<-as.vector(matrix(X[1:(n*N),2],N,n,byrow=TRUE))
y2<-X2[1:(n*N),]+U2[1:(n*N),]+P2[1:(n*N),]+E2[1:(n*N),]
data.mod2<-data.frame(ANIMAL=rep(pedi$ANIMAL[1:G[p+2]],n),identif=rep(pedi$ANIMAL[1:G[p+2]],n),
fac1=rep(X11[1:N],n),fac2=rep(X12[1:N],n),fac3=rep(X13[1:N],n),time=rep(1:n,each=N),y=y2[1:(n*N),])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 2 1 1 1.872537 4.4656825
## 2 2 2 4 1 1 1 13.353741 8.0315186
## 3 3 3 4 2 2 1 9.389222 11.3830031
## 4 4 4 2 1 3 1 -2.952822 -0.7504775
## 5 5 5 2 1 1 1 8.100532 -3.8392779
## 6 6 6 3 2 3 1 14.241269 7.2102702##A inverse
ainv=ainverse(pedi[1:N,])
##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)## Warning in na.action$y == "include" && any(is.na(data[, lhs, with = FALSE])):
## ‘length(x) = 3 > 1’ dans la conversion automatique vers ‘logical(1)’## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Wed Mar 22 16:15:37 2023
## LogLik Sigma2 DF wall cpu
## 1 -968481.3 2.37895 1000792 16:15:42 3.8
## 2 -832314.4 1.78376 1000792 16:15:44 2.6
## 3 -715522.8 1.38602 1000792 16:15:47 2.6
## 4 -645566.8 1.18400 1000792 16:15:50 2.6
## 5 -607441.3 1.07898 1000792 16:15:52 2.6
## 6 -592535.9 1.03448 1000792 16:15:55 2.6
## 7 -587511.6 1.01431 1000792 16:15:57 2.6 (1 restrained)
## 8 -586334.1 1.00562 1000792 16:16:00 2.6
## 9 -586202.0 1.00285 1000792 16:16:03 2.6
## 10 -586199.0 1.00242 1000792 16:16:05 2.6
## 11 -586199.0 1.00241 1000792 16:16:08 2.6#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] -586199
##
## $nedf
## [1] 1000792
##
## $sigma
## [1] 1.001204
##
## $varcomp
## component std.error
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order0:order0 4.1348219 0.39598952
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order0 0.5956330 0.18102314
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order1 4.4553403 0.38399072
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order0 0.4665148 0.06155100
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order1 -0.1837315 0.05723133
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order2 0.7631065 0.04808452
## leg(time, 2):ANIMAL!leg(time, 2)_order0 10.8107629 0.27246514
## leg(time, 2):ANIMAL!leg(time, 2)_order1 8.4748635 0.24346629
## leg(time, 2):ANIMAL!leg(time, 2)_order2 0.6285926 0.02606453
## units!R 1.0024091 0.00143880
## z.ratio bound %ch
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order0:order0 10.441746 P 0
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order0 3.290369 P 0
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order1 11.602729 P 0
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order0 7.579322 P 0
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order1 -3.210330 P 0
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order2 15.870108 P 0
## leg(time, 2):ANIMAL!leg(time, 2)_order0 39.677600 P 0
## leg(time, 2):ANIMAL!leg(time, 2)_order1 34.809186 P 0
## leg(time, 2):ANIMAL!leg(time, 2)_order2 24.116783 P 0
## units!R 696.697927 P 0
##
## $bic
## [1] 1172536
## attr(,"parameters")
## [1] 10
##
## $aic
## [1] 1172418
## attr(,"parameters")
## [1] 10
##
## attr(,"class")
## [1] "summary.asreml"gc()#garbage #free the memory## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 3140878 167.8 5809097 310.3 3253886 173.8
## Vcells 37157014 283.5 55518102 423.6 55410778 422.8detach("package:asreml")4.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="")
)write.table(data.mod2,file=paste(REPTRA,"/data.dat",sep=""),row.names=FALSE,col.names=FALSE)
write.table(pedi[1:N,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="")))The lanceRR.sh file is available at https://genoweb.toulouse.inra.fr/~trohmer/divers/
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 1.002410 1.0012038
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order0:order0 4.124860 4.1348219
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order0 0.594209 0.5956330
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order1:order1 4.444620 4.4553403
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order0 0.465392 0.4665148
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order1 -0.183290 -0.1837315
## leg(time, 2):vm(ANIMAL, ainv)!leg(time, 2)_order2:order2 0.761273 0.7631065
## leg(time, 2):ANIMAL!leg(time, 2)_order0 10.784800 10.8107629
## leg(time, 2):ANIMAL!leg(time, 2)_order1 8.454500 8.4748635
## leg(time, 2):ANIMAL!leg(time, 2)_order2 0.627082 0.6285926
## 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="")
)write.table(data.mod2,file=paste(REPTRA,"/data.dat",sep=""),row.names=FALSE,col.names=FALSE)
write.table(pedi[1:N,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 > paramThe lanceRR2.sh file is available at https://genoweb.toulouse.inra.fr/~trohmer/divers/
library(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 1.002410 1.000
## residual 2 1 0.893543 0.891
## residual 2 2 1.002630 1.000
## var.gen 1 1 4.115820 4.000
## var.gen 2 1 0.578207 0.620
## var.gen 2 2 4.519210 5.000
## var.gen 3 1 0.473105 0.360
## var.gen 3 2 -0.191050 -0.280
## var.gen 3 3 0.775155 0.800
## var.gen 4 1 0.738600 0.588
## var.gen 4 2 0.387242 0.200
## var.gen 4 3 0.206639 0.200
## var.gen 4 4 4.025920 4.000
## var.gen 5 1 0.237009 0.200
## var.gen 5 2 0.662808 0.588
## var.gen 5 3 0.196144 0.200
## var.gen 5 4 0.782055 0.620
## var.gen 5 5 4.859790 5.000
## var.gen 6 1 0.298541 0.200
## var.gen 6 2 0.272039 0.200
## var.gen 6 3 0.547782 0.588
## var.gen 6 4 0.398364 0.360
## var.gen 6 5 -0.306389 -0.280
## var.gen 6 6 0.752900 0.800
## var.per 1 1 10.823400 11.000
## var.per 2 1 0.000000 0.000
## var.per 2 2 8.439840 8.000
## var.per 3 1 0.000000 0.000
## var.per 3 2 0.000000 0.000
## var.per 3 3 0.622737 0.600
## var.per 4 1 0.305157 0.588
## var.per 4 2 0.281400 0.200
## var.per 4 3 0.218194 0.200
## var.per 4 4 10.829800 11.000
## var.per 5 1 0.193315 0.200
## var.per 5 2 0.513664 0.588
## var.per 5 3 0.173511 0.200
## var.per 5 4 0.000000 0.000
## var.per 5 5 8.109110 8.000
## var.per 6 1 0.239838 0.200
## var.per 6 2 0.226219 0.200
## var.per 6 3 0.619409 0.588
## var.per 6 4 0.000000 0.000
## var.per 6 5 0.000000 0.000
## var.per 6 6 0.641482 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.1131224 0.073773044.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")
Acknowledgements
Thanks to H. Chapuis and I. David for their help.
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.
Mrode, Raphael A. 2014. Linear Models for the Prediction of Animal Breeding Values. 3rd ed. Wallingford: CABI.
Patterson, H Desmond, and Robin Thompson. 1971. “Recovery of Inter-Block Information When Block Sizes Are Unequal.” Biometrika 58: 545–54.
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.
GenPhySE, Université de Toulouse, INRAE, ENVT, F-31326 Castanet Tolosan, tom.rohmer@inrae.fr↩︎