# A sire Model 
# Linear Models for the prediction of Animal Breeding values.
# Julius Mugambe PhD

# variance components
sigma_s <- 5.0
sigma_e <- 55
alpha   <- sigma_e/sigma_s
# AINV
Ainv = matrix(data = c(1.333,0.000,-0.667,0.000,1.000,0.000,-0.667,0.000,1.333), byrow= T, ncol = 3)
X = matrix(data = c (1,0,0,1,0,1,1,0,1,0), byrow=T,ncol=2)
Z = matrix(data = c (1,0,0,0,1,0,1,0,0,0,0,1,0,1,0), byrow = TRUE, ncol = 3)
y = matrix(data = c (4.5,2.9,3.9,3.5,5.0))
#Mixed model equations
XtX <- crossprod(X)
XtZ <- crossprod(X, Z)
ZtX <- crossprod(Z, X)
ZtZ <- crossprod(Z)
# The addition of Ainv_alpha to ZtZ in the LSE yeilds the MME
ZtZA <- crossprod(Z) + (Ainv) * alpha
Xty <- crossprod(X, y)
Zty <- crossprod(Z, y)
LHS1 <- cbind(XtX, XtZ)
LHS2 <- cbind(ZtX, ZtZA)
LHS <- rbind(LHS1, LHS2)
RHS <- rbind(Xty, Zty)
# solutions
sol <- solve(LHS) %*% RHS
round(sol, digits=3)
##        [,1]
## [1,]  4.336
## [2,]  3.382
## [3,]  0.022
## [4,]  0.014
## [5,] -0.043