#  Iteration on Data - Solve MME - Summer Course  2012

#  Years and animals are confounded, thus if year effects are included in
#  the model of analysis they take some of the genetic effects 
#  If mu replaces year effects, then all genetic effects go into the
#   animal solutions, but also some of the year effects go into the
#   animal solutions.
#  So the model must have random year effects and an overall mean.

load(file="Larry.Rdata")
names(MYdat2)
nrec = nrow(MYdat2)
obs = MYdat2$obs
yrfac = MYdat2$year
cgfac = MYdat2$cgrp
anwr  = MYdat2$anwr
bi    = 1/MYdat2$bi
aid   = MYdat2$anwr
sid   = MYdat2$sid
did   = MYdat2$did
atbv  = MYdat2$atbv

# Pedigree Preparation
mmm = aid[1] - 1
NP = aid[length(aid)]+1
sanm = rep(0,NP)
syr = rep(0,NP)
scg = rep(0,NP)
smu=0
aaid = c(1:mmm, aid, NP)
ssid = c(rep(NP,mmm), sid, NP)
ddid = c(rep(NP,mmm), did, NP)
bbi  = aaid*0
mwr  = mmm + 1
nam  = length(aaid)-1
bbi[1:mmm]=1
bbi[mwr:nam]=bi
bbi[NP]=0

#  More pedigree prep
mtotal = length(aaid)
adiag = rep(0,mtotal)
for(i in 1:nam){
  kis = ssid[i]; kid = ddid[i]
  d = bbi[i];  x = 0.25*d
  adiag[i]   = adiag[i] + d
  adiag[kis] = adiag[kis] + x
  adiag[kid] = adiag[kid] + x
 }
adiag[NP]= 0
zz = c(rep(0,mmm), rep(1,nrec), 0)

#  Initialize
iters = 0; nmax=1000; ccc = 100
vare = 64; alpha = 1.6; alphc = 4; alphy=200

#  Loop
while( ccc > 0.00000001){
   iters = iters + 1
   rhs = (obs - sanm[aid] - scg[cgfac] - syr[yrfac])
   smu = mean(rhs)
# Year effects, fixed - changed to an overall mean
   rhs = (obs - sanm[aid] - scg[cgfac] - smu)
   dyr = tapply(rep(1,nrec),yrfac,sum) + alphy
   syr = tapply(rhs,yrfac,sum)/dyr
#   CG effects, random
   rhs = (obs - smu - sanm[aid] - syr[yrfac])
   dcg = tapply(rep(1,nrec),cgfac,sum) + alphc
   scg = tapply(rhs,cgfac,sum)/dcg
#  Animals
  rhs = (obs - smu - syr[yrfac] - scg[cgfac])
  rhs = c( rep(0,mmm), tapply(rhs,aid,sum), 0)
#  take care of A-inverse elements
ccd = 0; ccs=0
for(i in nam:1){
  kis = ssid[i]; kid = ddid[i]; d = alpha*bbi[i]
  x = zz[i] + alpha*adiag[i]
  rhs[i]=rhs[i]+ 0.5*d*(sanm[kis]+sanm[kid])
  anew = rhs[i]/x
  diff = anew - sanm[i]
  ccd = ccd + diff*diff
  ccs = ccs + anew*anew
  sanm[i] = anew
  rhs[kis]=rhs[kis] + d*(0.5*sanm[i] - 0.25*sanm[kid])
  rhs[kid]=rhs[kid] + d*(0.5*sanm[i] - 0.25*sanm[kis])
 }
  ccc = ccd/ccs
  if(iters > nmax)ccc=0
 }
#  end of while loop
ccc
iters
ebv = sanm[aid]
gtrnd = tapply(ebv,yrfac,mean)
truet = tapply(atbv,yrfac,mean)
plot(truet, type="b", lwd=3, col="red",xlab="Year of Birth", 
  ylab="Trait Units",ylim=c(-1,30))
lines(gtrnd, type="b", lwd=3, col="blue")
title(main="Genetic Change By Year of Birth")
cor(ebv,MYdat2$atbv)