# R/MLE.popsize.R In MasterBayes: ML and MCMC Methods for Pedigree Reconstruction and Analysis

```"MLE.popsize"<-function(X.list, USdam=FALSE, USsire=FALSE, nUS=NULL, ped=NULL, shrink=NULL){

if(length(USdam)==1){
if(USdam==TRUE){
USdam<-rep(1, length(X.list\$X))
}else{
}
}else{
if(length(USdam)!=length(X.list\$X)){stop("length of USdam does not equal number of offspring")}
}

if(USsire[1]=="USdam"){
USsiredam<-TRUE
USsire<-USdam
}else{
USsiredam<-FALSE
}

if(length(USsire)==1){
if(USsire==TRUE){
nbetaS<-1
USsire<-rep(1, length(X.list\$X))
}else{
nbetaS<-0
}
}else{
if(length(USsire)!=length(X.list\$X)){stop("length of USsire does not equal number of offspring")}
nbetaS<-length(unique(USsire))
}

if(length(nUS)==0){
}else{
warning("beta is wrong size in popsize.loglik")
stop()
}else{
nUS<-as.matrix(nUS)
}
}

nsire<-unlist(lapply(X.list\$X, function(x){length(x\$sire.id)}))
ndam<-unlist(lapply(X.list\$X, function(x){length(x\$dam.id)}))
nudam<-unlist(lapply(X.list\$X, function(x){length(x\$restdam.id)}))
nusire<-unlist(lapply(X.list\$X, function(x){length(x\$restsire.id)}))

npar<-ndam*nsire

if(is.null(ped)==FALSE){
ped<-NULL
}
}

if(is.null(ped)==FALSE){
if(sum(is.na(ped[,3])==FALSE)==0 & nbetaS>0){
ped<-NULL
}
}

if(is.null(ped)==FALSE){
ped<-ped[match(X.list\$id[as.numeric(names(X.list\$X))], ped[,1]),]
}

X<-lapply(X.list\$X, function(x){list(N=NULL,G=NULL)})

for(i in 1:length(X)){
nsampS<-nsire[i]-(nbetaS>0)

sampD<-c(1:npar[i])
sampS<-c(1:npar[i])

sampD<-c(1:npar[i])[-c((1:nsire[i])+((nudam[i]-1)*nsire[i]))]
}
if(nbetaS>0){
sampS<-c(1:npar[i])[-c(((0:(ndam[i]-1))*nsire[i])+nusire[i])]
}
DandS<-intersect(sampS, sampD)
DnotS<-setdiff(sampD, sampS)
SnotD<-setdiff(sampS, sampD)
X[[i]]\$N<-c(nsampD*nsampS, nsampS,  nsampD, 1)
if(is.null(ped)){
X[[i]]\$G<-c(sum(X.list\$X[[i]]\$G[DandS]), sum(X.list\$X[[i]]\$G[SnotD]), sum(X.list\$X[[i]]\$G[DnotS]),sum(X.list\$X[[i]]\$G[notDS]))
}
}

optim.out <- optim(nUS, popsize.loglik, method = "L-BFGS-B",lower=lower, upper=upper, control = list(fnscale = -1, ndeps=ndpes), hessian = TRUE, X=X, USdam=USdam, USsire=USsire, ped=ped, USsiredam=USsiredam, shrink=shrink)

if(USsiredam){
C[1]<-solve(-1 * optim.out\$hessian)
}else{
}
B.start <- rep(optim.out\$par,2)
}else{
C <- solve(-1 * optim.out\$hessian)
B.start <- optim.out\$par
}
list(nUS=B.start, C=C)
}
```

## Try the MasterBayes package in your browser

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

MasterBayes documentation built on July 30, 2017, 1:01 a.m.