# Fitting Mixed Model
data(tpod)
S = seq(1,350,50)
test1=gibbs(y,X=~factor(fam),Z=gen[,S],S=0.5)
plot(test1$Fit.mean,y)
# Fitting GBLUP
K = tcrossprod(gen)
K = K/mean(diag(K))
iK = chol2inv(K)
test2=gibbs(y,iK=iK)
plot(test2$Fit.mean,y)
Run the code above in your browser using DataLab