# Mrode 1996, pp 74-76

library(pedigree)
id <- 1:15 # individual ID
dam <- c(0,0,0,0,0,2,2,2,4,4,4,4,5,5,5) # Dams; 0 for unknown
sire <- c(0,0,0,0,0,1,1,1,3,3,3,3,1,1,1) # Sires; 0 for unknown
ped <- data.frame(id,dam,sire)
makeA(ped,which = c(rep(FALSE,0),rep(TRUE,15)))
K <- read.table("A.txt")

A <- matrix(NA, nrow=15, ncol=15)
A[upper.tri(A,diag=TRUE)]<-K$V3
A <- forceSymmetric(A)

y<-matrix(c(90,70,65,98,106,60,80,100,85,68),nrow=10)
X<-matrix(c(1,0,0,0,1,0,0,1,0,1,0,1,1,1,0,1,1,0,1,0),nrow=10)
Z<-cbind(matrix(0,10,5),diag(10))
W<-matrix(c(1,1,1,0,0,0,0,0,0,0,
            0,0,0,1,1,1,1,0,0,0,
            0,0,0,0,0,0,0,1,1,1),nrow=10)
vu<-20
vc<-15
ve<-65     
a1<-ve/vu
a2<-ve/vc

# crossproducts
XX<-crossprod(X,X)
XZ<-t(X) %*% Z
XW<-t(X) %*% W
ZX<-t(Z) %*% X
ZZ<-crossprod(Z,Z)+a1*solve(as.matrix(A))
ZW<-t(Z) %*% W
WX<-t(W) %*% X
WZ<-t(W) %*% Z
WW<-crossprod(W,W)+a2*diag(3)

# mixed model equations
# coefficient matrix and right hand side
C<-rbind(cbind(XX,XZ,XW),cbind(ZX,ZZ,ZW),cbind(WX,WZ,WW))
rhs<-rbind(t(X) %*% y,t(Z) %*% y,t(W) %*% y)

#solution
theta.hat <- solve(C) %*% rhs
theta.hat