# Day 5 - Phantom Groups

# Any particular topics you want me to 
  discuss, write them on a sheet of paper
  and leave them here on the desk - then
  I can make sure to cover them next week

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

y=matrix(data=c(119,130,103,87,96,108,129,
    140,125,93,114,101,117,136,124),ncol=1)
year=c(1,1,1,1,1,2,2,2,2,2,2,3,3,3,3)
gend=c(2,1,1,2,1,1,2,2,1,1,2,1,2,1,1)
anim=c(1:30)
sire=c(rep(29,9),
   1,26,1,2,26,27,3,4,5,27,1,10,3,4,5,
  0,0,0,0,0,0)
dam=c(rep(30,9),
  25,6,7,25,9,7,28,28,6,8,28,9,12,15,18,
  0,0,0,0,0,0)
bi=c(1,1,1,1,1,1,1,1,1,.75,.75,.5,.75,.75,
  .75,.75,.75,.5,.75,.75,.5,.5,.5,.5,
   0,0,0,0,0,0)
cbind(anim,sire,dam,bi)
anwr=c(10:24)
AI=AINV(sire,dam,bi)
Xyr=desgn(year,0)
Xg=desgn(gend,0)
h2=0.2
vratio = (1 - h2)/h2
X=cbind(Xyr,Xg)
Z=desgn(anwr,30)
N=length(y)
RI=id(N)
GI=AI*vratio
qqq=MME(X,Z,GI,RI,y)
sol=qqq$SOLNS
kgrp=c(30:35)
ghat=sol[kgrp, ]
kIQ=c(16,22,28)
topIQ=sol[kIQ, ]

anim=c(1:24)
sire=c(rep(0,9),
   1,0,1,2,0,0,3,4,5,0,1,10,3,4,5)
dam=c(rep(0,9),
  0,6,7,0,9,7,0,0,6,8,0,9,12,15,18)

bi=c(1,1,1,1,1,1,1,1,1,.75,.75,.5,.75,.75,
  .75,.75,.75,.5,.75,.75,.5,.5,.5,.5)
anwr=c(10:24)
AI=AINV(sire,dam,bi)
Xyr=desgn(year,0)
Xg=desgn(gend,0)
h2=0.2
vratio = (1 - h2)/h2
X=cbind(Xyr,Xg)
Z=desgn(anwr,0)
N=length(y)
RI=id(N)
GI=AI*vratio
qqq=MME(X,Z,GI,RI,y)
sol2=qqq$SOLNS
kIQ=c(16,22,28)
topIQ2=sol2[kIQ, ]