#Additional file 3: R-code for SE(rpc)
#setwd("M:/wrk/programs/Rprograms/rpc/prediction")
rm(list=ls(all=TRUE))

##################################################################################################################################
# Prediction of SE(rpc_hat) according to Equation 1 of 
# "The standard error of the genetic correlation: how much data do we need to estimate a purebred-crossbred genetic correlation"
# by Piter Bijma and John W. M. Bastiaansen; published in GSE in 2014/15.
##################################################################################################################################
res_sd <- data.frame() 

####INPUT VALUES FOR SCHEMES TO BE EVALUATED####
vpp=1					#purebred phenotypic variance (Note: does not affect SE(rpc))
vpc=1					#crossbred phenotypic variance (Note: does not affect SE(rpc))
nsvals=c(100)			#numbers of sires
ndpvals=c(2)			#numbers of dams of sire line per sire
ndcvals=ndpvals
#ndcvals=c(500)			#numbers of dams of other line per sire
nopvals=c(50)			#number of purebred offspring per dam
nocvals=c(50)			#number of crossbred offspring per dam
rpcvals=c(0.8)          	#rpc.
h2pvals=c(0.3)		#heritability of purebred trait
h2cvals=c(0.3)
#h2cvals=c(0.3)		#heritability of crossbred trait
c2pvals=c(0.0,0.05,0.1,0.15,0.2,0.25)			#common litter effect (c^2) for purebred trait
c2cvals=c(0.0,0.05,0.1,0.15,0.2,0.25)		#common litter effect (c^2) for crossbred trait

###LOOPING over alternative schemes###
alt=0
for (i in 1:length(nsvals)) {
  ns=nsvals[i]
  for (j in 1:length(ndpvals)) {
    ndp=ndpvals[j]
    for (jj in 1:length(ndcvals)) {
      ndc=ndcvals[jj]
      for (k in 1:length(nopvals)) {
        nop=nopvals[k]
        #for (kk in 1:length(nocvals)) {
          noc=nop#nocvals[kk]
          for (l in 1:length(rpcvals)) {
            rpc=rpcvals[l]
            for (m in 1:length(h2pvals)) {
              h2p=h2pvals[m]
              vap=h2p*vpp
                #for (mm in 1:length(h2cvals)) {
                h2c=h2p #h2cvals[mm]
                vac=h2c*vpc
                capc=rpc*sqrt(vap*vac)
                for (n in 1:length(c2pvals)) {
                  c2p=c2pvals[n]
                  vcp=c2p*vpp
                  vep=(1-h2p-c2p)*vpp
                  for (nn in 1:length(c2cvals)) {
                    c2c=c2cvals[nn]
                    vcc=c2c*vpc
                    vec=(1-h2c-c2c)*vpc
                    alt=alt+1

			  #predicted SE(rpc) from Equation 1
                    r2p=0.25*vap/(0.25*vap+(0.25*vap+vcp)/ndp+(0.5*vap+vep)/(nop*ndp))
                    r2c=0.25*vac/(0.25*vac+(0.25*vac+vcc)/ndc+(0.5*vac+vec)/(noc*ndc))
                    sdrpc=sqrt( 1/(ns-1) * ( 1/(r2p*r2c) + rpc^2*(1+0.5*(1/r2p^2+1/r2c^2)-2*(1/r2p+1/r2c)) + rpc^4) )

			  #store predictions
                    res_sd <- rbind(res_sd,data.frame(alt=alt,ns=ns,ndp=ndp,ndc=ndc,nop=nop,noc=noc,
                              rpc=rpc,h2p=h2p,h2c=h2c,c2p=c2p,c2c=c2c,sdrpc=round(sdrpc,4) ))

         }
        }
       }
      }
     }
    } 
   }
  }
# }
#}

directory=getwd()
write("OUTPUT: Predicted SD of ESTIMATED VARIANCE COMPONENTS","Psdvce.txt",append=FALSE)
write(c(directory,date()),"Psdvce.txt",append=TRUE)
write.table(res_sd,"Psdvce.txt",row.names=FALSE,quote=FALSE,append=TRUE)
print (paste("File Psdvce.txt written to directory ",directory))

##############END OF PROGRAM##################################################################