R/fun4Kin.R In optiSel: Optimum Contribution Selection and Population Genetics

```### This function defines the linear component (a) and the constant term (d).
### The quadratic term (Q) remains unchanged.

"fun4Kin" <- function(obj, ageClass, phen, N){
if(obj\$breed=="across breeds"){
obj\$d <- 0
obj\$a <- 0*obj\$a
return(obj)
}
if(is.na(N)){stop("Parameter N is NA.\n")}

ageClass  <- ageClass[ageClass\$Breed==obj\$breed,]
classes   <- ageClass\$Class
r0        <- ageClass\$rcont0
r1        <- ageClass\$rcont1
diffF     <- obj\$mkin\$diffF
NIndiv0   <- pmax(N*r0,1)
NIndiv1   <- pmax(N*r1,1)
nIndiv    <- pmax(ageClass\$n,1)
ra        <- 1 - sum(r1)

Delta00  <- ra^2*(1-sum((phen\$c0*diag(obj\$Q))[phen\$Breed==obj\$breed]))/(2*ra*N)

if(ra>0.9999){
obj\$a    <- setNames(rep(0,nrow(phen)), phen\$Indiv)
obj\$d    <- Delta00
}else{
NPO       <- pPOpop(ageClass, Breed=obj\$breed, N0=ra*N)
nPO       <- NPOdata(phen, Classes=classes, symmetric=TRUE, quiet=TRUE)
NPO       <- NPO[classes,classes, drop=FALSE]
diag(NPO) <- 1/NIndiv0
diag(nPO) <- 1/nIndiv

u <- ageClass[phen\$Class,"rcont1"]
u <- u*(1/pmax(ageClass[phen\$Class,"n"],1) - 1/pmax(N*ageClass[phen\$Class,"rcont0"],1))
u <- u*(obj\$mkin\$selfkin[phen\$Indiv]-obj\$mkin\$kinwac[phen\$Indiv])
u[is.na(u)|phen\$Breed!=obj\$breed] <- 0

Delta01  <- 2*sum(phen\$c1*u)
Delta11  <- c(r1%*%(diag(1/NIndiv1-1/NIndiv0, nrow=length(NIndiv1))*diffF - (nPO-NPO)*diffF)%*%r1)
obj\$a    <- setNames(-2*c(u), phen\$Indiv)
obj\$d    <- Delta00 + Delta01 + Delta11
}

#cat(paste0("Delta00 =", Delta00, "\n"))
#cat(paste0("Delta01 =", Delta01, "\n"))
#cat(paste0("Delta11 =", Delta11, "\n\n"))

return(obj)
}
```

Try the optiSel package in your browser

Any scripts or data that you put into this service are public.

optiSel documentation built on Jan. 16, 2021, 5:30 p.m.