BGLR with Reps
Published:
To fit a model that has reps or any unbalanced design in BGLR, you need to first design the Z matrix, you can use model.matrix function for that. Then, in the ETA object you need to provide a ZGZ’ or ZM instead of standard X or K matrices. More on it here. Example below:
library(BGLR)
data(wheat)
nIter=2000
burnIn=1000
X=scale(wheat.X)/sqrt(ncol(wheat.X))
y=wheat.Y[,2]
fm1=BGLR( y=y,ETA=list(mrk=list(X=X,model='BRR')),
nIter=nIter,burnIn=burnIn,saveAt='brr_'
)
varE=scan('brr_varE.dat')
varU=scan('brr_ETA_mrk_varB.dat')
h2_1=varU/(varU+varE)
G=tcrossprod(X)
fm2=BGLR( y=y,ETA=list(G=list(K=G,model='RKHS')),
nIter=nIter,burnIn=burnIn,saveAt='eig_'
)
varE=scan( 'eig_varE.dat')
varU=scan('eig_ETA_G_varU.dat')
h2_2=varU/(varU+varE)
## With reps
y2=c(wheat.Y[,2],wheat.Y[,3])
Z = rbind( diag(nrow(wheat.Y)), diag(nrow(wheat.Y)))
ZGZ = Z %*% G %*% t(Z)
fm3=BGLR( y=y2,ETA=list(G=list(K=ZGZ,model='RKHS')),
nIter=nIter,burnIn=burnIn,saveAt='eig3_'
)
varE=scan( 'eig3_varE.dat')
varU=scan('eig3_ETA_G_varU.dat')
h2_3=varU/(varU+varE)
## With reps BRR
y2=c(wheat.Y[,2],wheat.Y[,3])
Z = rbind( diag(nrow(wheat.Y)), diag(nrow(wheat.Y)))
ZX = Z %*% X
fm4=BGLR( y=y2,ETA=list(mrk=list(X=ZX,model='BRR')),
nIter=nIter,burnIn=burnIn,saveAt='brr_'
)
varE=scan('brr_varE.dat')
varU=scan('brr_ETA_mrk_varB.dat')
h2_4=varU/(varU+varE)
plot(fm4$yHat,fm3$yHat)
## GEBVs will be avg per line (or get the first value matching the line, will give similar result)
gebv4 = tapply(fm4$yHat,rep(1:length(y),2),"mean")
gebv3 = tapply(fm3$yHat,rep(1:length(y),2),"mean")