MATERNAL EFFECTS MODELS

Dams provide a maternal environment to their offspring.
The maternal ability is a genetically controlled trait of the dam.

Weight at age x = Age of dam x gender of lamb(calf,piglet)
                  + Year-Month of birth
                  + CG contemporary group
                  + animal (lamb,calf,piglet)
                  + maternal effect of dam
                  + mat. permanent environmental effect
                  + residual




Pedigree         Data
aid  sid  did   animal  dam  age  gender  Weight
  1   -    -     10      2    4     1=F    36
  2   -    -     11      4    4     2=M    41
  3   -    -     12      6    3     1      32
  4   -    -     13      8    3     2      38
  5   -    -     14      2    5     1      31
  6   -    -     15      8    4     2      35
  7   -    -     16     10    3     1      34
  8   -    -     17      6    4     2      37
  9   -    -     18     12    3     1      33
 10   1    2     19     14    3     2      40
 11   3    4
 12   3    6
 13   5    8
 14   7    2
 15   3    8
 16   5   10
 17   9    6
 18   9   12
 19  11   14

Priors  
G = 3.3  0.0
    0.0  1.6   vare= 6.1  varp=0.9  
degrees 0f belief = 20

Set up MME
One Gibbs sampling
One REML  DF or EM methods



library(MASS)
zlrs=file.choose()  # lrscrips.R
source(zlrs)

#  Enter known information
aid=c(1:19)
sid=c(rep(0,9),1,3,3,5,7,3,5,9,9,11)
did=c(rep(0,9),2,4,6,8,2,8,10,6,12,14)

G=matrix(data=c(3.3,0,0,1.6),byrow=TRUE,nrow=2)
vare=6.1
varp=0.9

y=matrix(data=c(36,41,32,38,31,35,34,37,33,40),ncol=1)
gend=c(1,2,1,2,1,2,1,2,1,2)
age=c(4,4,3,3,5,4,3,4,3,3) - 2

dam=c(2,4,6,8,2,8,10,6,12,14)  # maternal effects and maternal PE
anwr=c(10:19)

# Create design matrices
Za=desgn(anwr,0)
Xa=desgn(age,0)
Xg=desgn(gend,0)
Zm=desgn(dam,0)
Zmp=Zm[ ,dam]
Z=cbind(Za,Zm,Zmp)
dim(Z)
X=cbind(Xg,Xa)

#  Set up covariance matrices
RI = id(length(y))
bi=c(rep(1,9),rep(0.5,10))
AI=AINV(sid,did,bi)
GI = solve(G)*vare
Gam = GI %x% AI
Gmp = id(length(dam))*(vare/varp)
HI = block(Gam,Gmp)
dim(HI)

T=MME(Xa,Z,HI,RI,y)
bh=T$SOLNS
SSE = sum(y*y) - T$SSR

# REML EM
ka=c(6:24)
km=c(25:43)
kp=c(44:50)
ahat=bh[ka, ]
mhat=bh[km, ]
ghat = cbind(ahat,mhat)
C = T$C
Caa=C[ka,ka]
Cam=C[ka,km]
Cmm=C[km,km]
phat=bh[kp, ]
Cpp=C[kp,kp]

evhat=SSE/(length(y) - 4)
evhat

aAa = t(ghat) %*% AI %*% ghat
trCa = sum(diag(AI %*% Caa))
trCm = sum(diag(AI %*% Cmm))
trCam = sum(diag(AI %*% Cam))
q=length(ahat)
trs=c(trCa,trCam,trCam,trCm)
Cii=matrix(data=trs,byrow=TRUE,ncol=2)
Ghat=(aAa + Cii*evhat)/q
Ghat

Phat = (t(phat)%*%phat + sum(diag(Cpp))*evhat)/length(phat)
Phat