# We are going to try to predict selection response using 
# the breeders equation

# The data are from a selection experiment by Beldade et al. (2002) Nature 416, 844-847
# with 15 selection lines that are followed for 9/10 generations, including three control lines.
# two traits are under selection, A Anterior relative eyespot size P posterior relative eyespot size

# Load a file which contains "prepared" statistics.

prelinedatn<-read.table("prelinedatn")

# the following function takes the prepared statistics as input, a name of a specific selection line
# and a G matrix
# don't go through this in detail. Just try if you can spot some of the equations from the lectures.

liner<-function(predat,linename,gmat){
# means and variances per line
dat<-subset(predat,subset=line==linename)
maxline<-max(dat$g)-1
# calculations on data within those generations where selection occurred
dat$sA<-dat$As-dat$A;
dat$sP<-dat$Ps-dat$P;
dat$BA<-(dat$VP*(dat$sA)-dat$COV*(dat$sP))/(dat$VA*dat$VP-dat$COV^2);
dat$BA<-dat$BA/2;
dat$BP<-(-dat$COV*(dat$sA)+dat$VA*(dat$sP))/(dat$VA*dat$VP-dat$COV^2);
dat$BP<-dat$BP/2;
dat$RA<-(gmat[1]*(dat$VP*(dat$sA)-dat$COV*(dat$sP))+gmat[2]*(-dat$COV*(dat$sA)+dat$VA*(dat$sP)))/(dat$VA*dat$VP-dat$COV^2);
dat$RA<-dat$RA/2; # response in A
dat$RP<-(gmat[2]*(dat$VP*(dat$sA)-dat$COV*(dat$sP))+gmat[3]*(-dat$COV*(dat$sA)+dat$VA*(dat$sP)))/(dat$VA*dat$VP-dat$COV^2);
dat$RP<-dat$RP/2; # response in P
dat$predA<-dat$A+dat$RA;
dat$predP<-dat$P+dat$RP;
# same calculations, but with changing G matrix acccording predictions from infinitesimal model
# the G matrices before and after selection are not added as extra variables to each generation
# there is only one set of variables G11, G12 and G22 that are added to each generation and that 
# correspond to the G matrix before selection in that generation
dat$G11<-c(rep(gmat[1],length(dat$predP)));
dat$G12<-c(rep(gmat[2],length(dat$predP)));
dat$G22<-c(rep(gmat[3],length(dat$predP)));
for(i in (maxline+1):1){
dat$G11[i]<-dat$G11[i+1]+(1/(dat$COV[i+1]^2-dat$VA[i+1]*dat$VP[i+1])^2*(-2*dat$COV[i+1]^3*dat$G11[i+1]*dat$G12[i+1]+
dat$VP[i+1]*(dat$G12[i+1]*dat$VA[i+1]*(2*dat$COVs[i+1]*dat$G11[i+1]-dat$G12[i+1]*dat$VA[i+1])-
dat$G11[i+1]^2*(dat$VA[i+1]-dat$VAs[i+1])*dat$VP[i+1])+dat$G12[i+1]^2*dat$VA[i+1]^2*dat$VPs[i+1]+
dat$COV[i+1]^2*(2*dat$COVs[i+1]*dat$G11[i+1]*dat$G12[i+1]+dat$G12[i+1]^2*(dat$VA[i+1]+dat$VAs[i+1])+
dat$G11[i+1]^2*(dat$VP[i+1]+dat$VPs[i+1]))-
2*dat$COV[i+1]*(dat$COVs[i+1]*(dat$G12[i+1]^2*dat$VA[i+1]+dat$G11[i+1]^2*dat$VP[i+1])+
dat$G11[i+1]*dat$G12[i+1]*(-dat$VA[i+1]*dat$VP[i+1]+dat$VAs[i+1]*dat$VP[i+1]+dat$VA[i+1]*dat$VPs[i+1]))))/4;
dat$G12[i]<-dat$G12[i+1]+(1/(dat$COV[i+1]^2-dat$VA[i+1]*dat$VP[i+1])^2*(-dat$COV[i+1]^3*(dat$G12[i+1]^2+
dat$G11[i+1]*dat$G22[i+1])+dat$COVs[i+1]*(dat$G12[i+1]^2+dat$G11[i+1]*dat$G22[i+1])*dat$VA[i+1]*dat$VP[i+1]+
dat$G12[i+1]*(dat$G11[i+1]*(-dat$VA[i+1]+dat$VAs[i+1])*dat$VP[i+1]^2+
dat$G22[i+1]*dat$VA[i+1]^2*(-dat$VP[i+1]+dat$VPs[i+1]))-
dat$COV[i+1]*(2*dat$COVs[i+1]*dat$G12[i+1]*(dat$G22[i+1]*dat$VA[i+1]+dat$G11[i+1]*dat$VP[i+1])+(dat$G12[i+1]^2+
dat$G11[i+1]*dat$G22[i+1])*(dat$VAs[i+1]*dat$VP[i+1]+dat$VA[i+1]*(-dat$VP[i+1]+dat$VPs[i+1])))+
dat$COV[i+1]^2*(dat$COVs[i+1]*(dat$G12[i+1]^2+dat$G11[i+1]*dat$G22[i+1])+
dat$G12[i+1]*(dat$G22[i+1]*(dat$VA[i+1]+dat$VAs[i+1])+dat$G11[i+1]*(dat$VP[i+1]+dat$VPs[i+1])))))/4;
dat$G22[i]<-dat$G22[i+1]+(1/(dat$COV[i+1]^2-dat$VA[i+1]*dat$VP[i+1])^2*(-2*dat$COV[i+1]^3*dat$G12[i+1]*dat$G22[i+1]+
dat$VP[i+1]*(dat$G22[i+1]*dat$VA[i+1]*(2*dat$COVs[i+1]*dat$G12[i+1]-dat$G22[i+1]*dat$VA[i+1])-
dat$G12[i+1]^2*(dat$VA[i+1]-dat$VAs[i+1])*dat$VP[i+1])+dat$G22[i+1]^2*dat$VA[i+1]^2*dat$VPs[i+1]+
dat$COV[i+1]^2*(2*dat$COVs[i+1]*dat$G12[i+1]*dat$G22[i+1]+dat$G22[i+1]^2*(dat$VA[i+1]+dat$VAs[i+1])+
dat$G12[i+1]^2*(dat$VP[i+1]+dat$VPs[i+1]))-
2*dat$COV[i+1]*(dat$COVs[i+1]*(dat$G22[i+1]^2*dat$VA[i+1]+dat$G12[i+1]^2*dat$VP[i+1])+
dat$G12[i+1]*dat$G22[i+1]*(-dat$VA[i+1]*dat$VP[i+1]+dat$VAs[i+1]*dat$VP[i+1]+dat$VA[i+1]*dat$VPs[i+1]))))/4
}
dat$RAg<-(dat$G11*(dat$VP*(dat$sA)-dat$COV*(dat$sP))+dat$G12*(-dat$COV*(dat$sA)+dat$VA*(dat$sP)))/(dat$VA*dat$VP-dat$COV^2);
dat$RAg<-dat$RAg/2; # response in A allowing G to change
dat$RPg<-(dat$G12*(dat$VP*(dat$sA)-dat$COV*(dat$sP))+dat$G22*(-dat$COV*(dat$sA)+dat$VA*(dat$sP)))/(dat$VA*dat$VP-dat$COV^2);
dat$RPg<-dat$RPg/2; # response in P allowing G to change
dat$predAg<-dat$A+dat$RAg;
dat$predPg<-dat$P+dat$RPg;
dat$eA<-c(dat$A[-(maxline+2)]-dat$predA[-1],NA); # "error"
dat$eP<-c(dat$P[-(maxline+2)]-dat$predP[-1],NA);
dat$eAg<-c(dat$A[-(maxline+2)]-dat$predAg[-1],NA);
dat$ePg<-c(dat$P[-(maxline+2)]-dat$predPg[-1],NA);
dat$eAl<-c(dat$A[-(maxline+2)]-dat$predAl[-1],NA);
dat$ePl<-c(dat$P[-(maxline+2)]-dat$predPl[-1],NA);
# calculating from base population and selection response only
# cumulative responses
dat$cumRA<-c(NA,rep(0,(maxline+1)))
dat$cumRP<-c(NA,rep(0,(maxline+1)))
for(i in 2:(maxline+2)){
dat$cumRA[i]<-sum(dat$RA[i:(maxline+2)]);
dat$cumRP[i]<-sum(dat$RP[i:(maxline+2)]);};
dat$cumRAg<-c(NA,rep(0,(maxline+1)))
dat$cumRPg<-c(NA,rep(0,(maxline+1)))
for(i in 2:(maxline+2)){
dat$cumRAg[i]<-sum(dat$RAg[i:(maxline+2)]);
dat$cumRPg[i]<-sum(dat$RPg[i:(maxline+2)]);};
dat
}

linenames<-levels(prelinedatn$line)

# ML estimate obtained using another routine

# Genetic Variance in A
GA<-17.3
# Genetic correlation (between -1 and 1) 
GC<-0.85
# Genetic variance in P
GP<-29.9

Gm<-c(GA,GC*sqrt(GA*GP),GP) # G matrix
# genetic variance in A, genetic covariance AP, genetic variance in P



	linedat1<-liner(prelinedatn,linenames[1],Gm)
	linedat2<-liner(prelinedatn,linenames[2],Gm)
	linedat3<-liner(prelinedatn,linenames[3],Gm)
	linedat4<-liner(prelinedatn,linenames[4],Gm)
	linedat5<-liner(prelinedatn,linenames[5],Gm)
	linedat6<-liner(prelinedatn,linenames[6],Gm)
	linedat7<-liner(prelinedatn,linenames[7],Gm)
	linedat8<-liner(prelinedatn,linenames[8],Gm)
	linedat9<-liner(prelinedatn,linenames[9],Gm)
	linedat10<-liner(prelinedatn,linenames[10],Gm)
	linedat11<-liner(prelinedatn,linenames[11],Gm)
	linedat12<-liner(prelinedatn,linenames[12],Gm)
	linedat13<-liner(prelinedatn,linenames[13],Gm)
	linedat14<-liner(prelinedatn,linenames[14],Gm)
	linedat15<-liner(prelinedatn,linenames[15],Gm)
	linedat<-rbind(linedat1,linedat2,linedat3,linedat4,linedat5,linedat6,linedat7,linedat8, linedat9,linedat10,linedat11,linedat12,linedat13,linedat14,linedat15)


# make a graph with the data and the model predictions

# for a model where we assume G fixed and no environmental effects per generation


 plot(linedat$P~linedat$A,type="n",xlim=c(-10,60),ylim=c(10,100))
 for( i in c(1:15)){ # control lines removed
 data<-subset(linedat,subset=linedat$line==linenames[i])
 points(data$P~data$A)
 lines(data$A,data$P,lty=3)
 lines(data$A[data$g==0]+data$cumRA,data$P[data$g==0]+data$cumRP,lty=1)
 } 
 
 
 # try a few other G matrices