dyn.load("WithBreeding.dll") #this only works in Windows R 32 bit #------------------------------------------------------------------------------ # logit logit <- function(inval) { low.val <- 0.00001 high.val <- 0.99999 inval[invalhigh.val] <- high.val # Output vector: log(inval/(1-inval)) } #------------------------------------------------------------------------------ # Inverse of logit expo <- function(inval) { 1/(1+exp(-inval)) } #------------------------------------------------------------------------------ #proportions to logits plg <- function(invec) { # How many elements? nel <- length(invec) # Case of 2 els: if (nel==2) outvec <- invec[1] # Case of more than 2 els: if (nel>2) { # Cumulative sum: cs <- cumsum(invec)+0.000000001 # Stop NaNs # Output vector: outvec <- invec[1:(nel-1)]/(1-c(0,cs[1:(nel-2)])) } # return logit: logit(outvec) } #------------------------------------------------------------------------------ # logit to proportions lgp <- function(invec) { # How many elements? nel <- length(invec) # Case of 1 el: if (nel==1) p.v <- expo(invec) # Case of more than 1 el: if (nel>1) { # Get inverse-logit vector: thisvec <- expo(invec) # Convert to proportions vector: need cum product of 1-thisvec cp <- cumprod(1-thisvec) # Return proportions vector: p.v <- thisvec*c(1,cp[1:(nel-1)]) } c(p.v,1-sum(p.v)) } #------------------------------------------------------------------------------ LL.WithBreeding.c <- function(N.valM, N.valU, beta.v, phi.m, p.m, p1.m, s.m, s1.m, psi.v, xi.v, pb.v, sb.v, c.v) { .C("WithBreedingLL", length(y.vectM), length(y.vectU), as.integer(y.vectM), as.integer(y.vectU), as.integer(K_NB), as.integer(K_B), as.integer(K), as.integer(x.matM), as.integer(x.matU), as.double(N.valM), as.double(N.valU), as.double(beta.v), as.double(phi.m), as.double(p.m), as.double(p1.m), as.double(s.m), as.double(s1.m), as.double(psi.v), as.double(xi.v), as.double(pb.v), as.double(sb.v), as.double(c.v), LLM = double (1), LLU = double (1), LLcombined = double (1), resultcode = integer(1) )$LLcombined } # ------------------- Functions to fit all models ---------------------------- # #------------------------------------------------------------------------------------------------# # 1. beta(t)phi(t)p(eiis)s(eiis)psi(lb)xi(lb)eta(lb) mLL.betat.phit.peiis.seiis.psilb.xilb.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2):(2*K_NB)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[2*K_NB + 1] + invect[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[2*K_NB + 3] + invect[2*K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[2*K_NB + 5] + invect[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[2*K_NB + 7] + invect[2*K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[2*K_NB + 9] + invect[2*K_NB + 10]*age) this.xi <- expo(invect[2*K_NB + 11] + invect[2*K_NB + 12]*age) this.pb <- expo(invect[2*K_NB + 13]) this.sb <- expo(invect[2*K_NB + 14]) this.c <- expo(invect[2*K_NB + 15] + invect[2*K_NB + 16]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phit.peiis.seiis.psilb.xilb.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phit.peiis.seiis.psilb.xilb.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2):(2*K_NB)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[2*K_NB + 1] + pars[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[2*K_NB + 3] + pars[2*K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[2*K_NB + 5] + pars[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[2*K_NB + 7] + pars[2*K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[2*K_NB + 9] + pars[2*K_NB + 10]*age) xi.est <- expo(pars[2*K_NB + 11] + pars[2*K_NB + 12]*age) pb.est <- expo(pars[2*K_NB + 13]) sb.est <- expo(pars[2*K_NB + 14]) c.est <- expo(pars[2*K_NB + 15] + pars[2*K_NB + 16]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 2. beta(t)phi(c)p(eiis)s(eiis)psi(lb)xi(lb)eta(lb) mLL.betat.phic.peiis.seiis.psilb.xilb.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 9] + invect[K_NB + 10]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[K_NB + 11] + invect[K_NB + 12]*age) this.xi <- expo(invect[K_NB + 13] + invect[K_NB + 14]*age) this.pb <- expo(invect[K_NB + 15]) this.sb <- expo(invect[K_NB + 16]) this.c <- expo(invect[K_NB + 17] + invect[K_NB + 18]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.peiis.seiis.psilb.xilb.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.peiis.seiis.psilb.xilb.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 9] + pars[K_NB + 10]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[K_NB + 11] + pars[K_NB + 12]*age) xi.est <- expo(pars[K_NB + 13] + pars[K_NB + 14]*age) pb.est <- expo(pars[K_NB + 15]) sb.est <- expo(pars[K_NB + 16]) c.est <- expo(pars[K_NB + 17] + pars[K_NB + 18]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 3. beta(t)phi(t)p(e)s(eiis)psi(lb)xi(lb)eta(lb) mLL.betat.phit.pe.seiis.psilb.xilb.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2):(2*K_NB)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[2*K_NB + 1] + invect[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[2*K_NB + 1] + invect[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[2*K_NB + 3] + invect[2*K_NB + 4]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[2*K_NB + 5] + invect[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[2*K_NB + 7] + invect[2*K_NB + 8]*age) this.xi <- expo(invect[2*K_NB + 9] + invect[2*K_NB + 10]*age) this.pb <- expo(invect[2*K_NB + 11]) this.sb <- expo(invect[2*K_NB + 12]) this.c <- expo(invect[2*K_NB + 13] + invect[2*K_NB + 14]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phit.pe.seiis.psilb.xilb.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phit.pe.seiis.psilb.xilb.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2):(2*K_NB)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[2*K_NB + 1] + pars[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[2*K_NB + 1] + pars[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[2*K_NB + 3] + pars[2*K_NB + 4]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[2*K_NB + 5] + pars[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[2*K_NB + 7] + pars[2*K_NB + 8]*age) xi.est <- expo(pars[2*K_NB + 9] + pars[2*K_NB + 10]*age) pb.est <- expo(pars[2*K_NB + 11]) sb.est <- expo(pars[2*K_NB + 12]) c.est <- expo(pars[2*K_NB + 13] + pars[2*K_NB + 14]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 4. beta(t)phi(t)p(eiis)s(e)psi(lb)xi(lb)eta(lb) mLL.betat.phit.peiis.se.psilb.xilb.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2):(2*K_NB)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[2*K_NB + 1] + invect[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[2*K_NB + 3] + invect[2*K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[2*K_NB + 5] + invect[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[2*K_NB + 5] + invect[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[2*K_NB + 7] + invect[2*K_NB + 8]*age) this.xi <- expo(invect[2*K_NB + 9] + invect[2*K_NB + 10]*age) this.pb <- expo(invect[2*K_NB + 11]) this.sb <- expo(invect[2*K_NB + 12]) this.c <- expo(invect[2*K_NB + 13] + invect[2*K_NB + 14]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phit.peiis.se.psilb.xilb.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phit.peiis.se.psilb.xilb.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2):(2*K_NB)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[2*K_NB + 1] + pars[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[2*K_NB + 3] + pars[2*K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[2*K_NB + 5] + pars[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[2*K_NB + 5] + pars[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[2*K_NB + 7] + pars[2*K_NB + 8]*age) xi.est <- expo(pars[2*K_NB + 9] + pars[2*K_NB + 10]*age) pb.est <- expo(pars[2*K_NB + 11]) sb.est <- expo(pars[2*K_NB + 12]) c.est <- expo(pars[2*K_NB + 13] + pars[2*K_NB + 14]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 5. beta(t)phi(t)p(eiis)s(eiis)psi(c)xi(lb)eta(lb) mLL.betat.phit.peiis.seiis.psic.xilb.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2):(2*K_NB)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[2*K_NB + 1] + invect[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[2*K_NB + 3] + invect[2*K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[2*K_NB + 5] + invect[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[2*K_NB + 7] + invect[2*K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- rep(expo(invect[2*K_NB + 9]), K_NB) this.xi <- expo(invect[2*K_NB + 10] + invect[2*K_NB + 11]*age) this.pb <- expo(invect[2*K_NB + 12]) this.sb <- expo(invect[2*K_NB + 13]) this.c <- expo(invect[2*K_NB + 14] + invect[2*K_NB + 15]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phit.peiis.seiis.psic.xilb.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phit.peiis.seiis.psic.xilb.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2):(2*K_NB)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[2*K_NB + 1] + pars[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[2*K_NB + 3] + pars[2*K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[2*K_NB + 5] + pars[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[2*K_NB + 7] + pars[2*K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- rep(expo(pars[2*K_NB + 9]), K_NB) xi.est <- expo(pars[2*K_NB + 10] + pars[2*K_NB + 11]*age) pb.est <- expo(pars[2*K_NB + 12]) sb.est <- expo(pars[2*K_NB + 13]) c.est <- expo(pars[2*K_NB + 14] + pars[2*K_NB + 15]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 6. beta(t)phi(t)p(eiis)s(eiis)psi(lb)xi(c)eta(lb) mLL.betat.phit.peiis.seiis.psilb.xic.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2):(2*K_NB)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[2*K_NB + 1] + invect[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[2*K_NB + 3] + invect[2*K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[2*K_NB + 5] + invect[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[2*K_NB + 7] + invect[2*K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[2*K_NB + 9] + invect[2*K_NB + 10]*age) this.xi <- rep(expo(invect[2*K_NB + 11]), K_NB) this.pb <- expo(invect[2*K_NB + 12]) this.sb <- expo(invect[2*K_NB + 13]) this.c <- expo(invect[2*K_NB + 14] + invect[2*K_NB + 15]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phit.peiis.seiis.psilb.xic.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phit.peiis.seiis.psilb.xic.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2):(2*K_NB)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[2*K_NB + 1] + pars[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[2*K_NB + 3] + pars[2*K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[2*K_NB + 5] + pars[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[2*K_NB + 7] + pars[2*K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[2*K_NB + 9] + pars[2*K_NB + 10]*age) xi.est <- rep(expo(pars[2*K_NB + 11]), K_NB) pb.est <- expo(pars[2*K_NB + 12]) sb.est <- expo(pars[2*K_NB + 13]) c.est <- expo(pars[2*K_NB + 14] + pars[2*K_NB + 15]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 7. beta(t)phi(t)p(eiis)s(eiis)psi(lb)ksi(lb)eta(c) mLL.betat.phit.peiis.seiis.psilb.xilb.etac <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2):(2*K_NB)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[2*K_NB + 1] + invect[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[2*K_NB + 3] + invect[2*K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[2*K_NB + 5] + invect[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[2*K_NB + 7] + invect[2*K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[2*K_NB + 9] + invect[2*K_NB + 10]*age) this.xi <- expo(invect[2*K_NB + 11] + invect[2*K_NB + 12]*age) this.pb <- expo(invect[2*K_NB + 13]) this.sb <- expo(invect[2*K_NB + 14]) this.c <- rep(expo(invect[2*K_NB + 15]), K_NB) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phit.peiis.seiis.psilb.xilb.etac <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phit.peiis.seiis.psilb.xilb.etac, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2):(2*K_NB)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[2*K_NB + 1] + pars[2*K_NB + 2]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[2*K_NB + 3] + pars[2*K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[2*K_NB + 5] + pars[2*K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[2*K_NB + 7] + pars[2*K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[2*K_NB + 9] + pars[2*K_NB + 10]*age) xi.est <- expo(pars[2*K_NB + 11] + pars[2*K_NB + 12]*age) pb.est <- expo(pars[2*K_NB + 13]) sb.est <- expo(pars[2*K_NB + 14]) c.est <- rep(expo(pars[2*K_NB + 15]), K_NB) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 8. beta(t)phi(c)p(e)s(eiis)psi(lb)xi(lb)eta(lb) mLL.betat.phic.pe.seiis.psilb.xilb.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[K_NB + 9] + invect[K_NB + 10]*age) this.xi <- expo(invect[K_NB + 11] + invect[K_NB + 12]*age) this.pb <- expo(invect[K_NB + 13]) this.sb <- expo(invect[K_NB + 14]) this.c <- expo(invect[K_NB + 15] + invect[K_NB + 16]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.pe.seiis.psilb.xilb.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.pe.seiis.psilb.xilb.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[K_NB + 9] + pars[K_NB + 10]*age) xi.est <- expo(pars[K_NB + 11] + pars[K_NB + 12]*age) pb.est <- expo(pars[K_NB + 13]) sb.est <- expo(pars[K_NB + 14]) c.est <- expo(pars[K_NB + 15] + pars[K_NB + 16]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 9. beta(t)phi(c)p(eiis)s(e)psi(lb)xi(lb)eta(lb) mLL.betat.phic.peiis.se.psilb.xilb.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[K_NB + 9] + invect[K_NB + 10]*age) this.xi <- expo(invect[K_NB + 11] + invect[K_NB + 12]*age) this.pb <- expo(invect[K_NB + 13]) this.sb <- expo(invect[K_NB + 14]) this.c <- expo(invect[K_NB + 15] + invect[K_NB + 16]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.peiis.se.psilb.xilb.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.peiis.se.psilb.xilb.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[K_NB + 9] + pars[K_NB + 10]*age) xi.est <- expo(pars[K_NB + 11] + pars[K_NB + 12]*age) pb.est <- expo(pars[K_NB + 13]) sb.est <- expo(pars[K_NB + 14]) c.est <- expo(pars[K_NB + 15] + pars[K_NB + 16]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 10. beta(t)phi(c)p(eiis)s(eiis)psi(c)xi(lb)eta(lb) mLL.betat.phic.peiis.seiis.psic.xilb.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 9] + invect[K_NB + 10]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- rep(expo(invect[K_NB + 11]), K_NB) this.xi <- expo(invect[K_NB + 12] + invect[K_NB + 13]*age) this.pb <- expo(invect[K_NB + 14]) this.sb <- expo(invect[K_NB + 15]) this.c <- expo(invect[K_NB + 16] + invect[K_NB + 17]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.peiis.seiis.psic.xilb.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.peiis.seiis.psic.xilb.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 9] + pars[K_NB + 10]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- rep(expo(pars[K_NB + 11]), K_NB) xi.est <- expo(pars[K_NB + 12] + pars[K_NB + 13]*age) pb.est <- expo(pars[K_NB + 14]) sb.est <- expo(pars[K_NB + 15]) c.est <- expo(pars[K_NB + 16] + pars[K_NB + 17]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 11. beta(t)phi(c)p(eiis)s(eiis)psi(lb)xi(c)eta(lb) mLL.betat.phic.peiis.seiis.psilb.xic.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 9] + invect[K_NB + 10]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[K_NB + 11] + invect[K_NB + 12]*age) this.xi <- rep(expo(invect[K_NB + 13]), K_NB) this.pb <- expo(invect[K_NB + 14]) this.sb <- expo(invect[K_NB + 15]) this.c <- expo(invect[K_NB + 16] + invect[K_NB + 17]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.peiis.seiis.psilb.xic.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.peiis.seiis.psilb.xic.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 9] + pars[K_NB + 10]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[K_NB + 11] + pars[K_NB + 12]*age) xi.est <- rep(expo(pars[K_NB + 13]), K_NB) pb.est <- expo(pars[K_NB + 14]) sb.est <- expo(pars[K_NB + 15]) c.est <- expo(pars[K_NB + 16] + pars[K_NB + 17]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 12. beta(t)phi(c)p(eiis)s(eiis)psi(lb)xi(lb)eta(c) mLL.betat.phic.peiis.seiis.psilb.xilb.etac <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 9] + invect[K_NB + 10]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[K_NB + 11] + invect[K_NB + 12]*age) this.xi <- expo(invect[K_NB + 13] + invect[K_NB + 14]*age) this.pb <- expo(invect[K_NB + 15]) this.sb <- expo(invect[K_NB + 16]) this.c <- rep(expo(invect[K_NB + 17]), K_NB) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.peiis.seiis.psilb.xilb.etac <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.peiis.seiis.psilb.xilb.etac, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 9] + pars[K_NB + 10]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[K_NB + 11] + pars[K_NB + 12]*age) xi.est <- expo(pars[K_NB + 13] + pars[K_NB + 14]*age) pb.est <- expo(pars[K_NB + 15]) sb.est <- expo(pars[K_NB + 16]) c.est <- rep(expo(pars[K_NB + 17]), K_NB) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 13. beta(t)phi(c)p(e)s(e)psi(lb)xi(lb)eta(lb) mLL.betat.phic.pe.se.psilb.xilb.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[K_NB + 7] + invect[K_NB + 8]*age) this.xi <- expo(invect[K_NB + 9] + invect[K_NB + 10]*age) this.pb <- expo(invect[K_NB + 11]) this.sb <- expo(invect[K_NB + 12]) this.c <- expo(invect[K_NB + 13] + invect[K_NB + 14]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.pe.se.psilb.xilb.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.pe.se.psilb.xilb.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[K_NB + 7] + pars[K_NB + 8]*age) xi.est <- expo(pars[K_NB + 9] + pars[K_NB + 10]*age) pb.est <- expo(pars[K_NB + 11]) sb.est <- expo(pars[K_NB + 12]) c.est <- expo(pars[K_NB + 13] + pars[K_NB + 14]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 14. beta(t)phi(c)p(e)s(eiis)psi(c)xi(lb)eta(lb) mLL.betat.phic.pe.seiis.psic.xilb.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- rep(expo(invect[K_NB + 9]), K_NB) this.xi <- expo(invect[K_NB + 10] + invect[K_NB + 11]*age) this.pb <- expo(invect[K_NB + 12]) this.sb <- expo(invect[K_NB + 13]) this.c <- expo(invect[K_NB + 14] + invect[K_NB + 15]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.pe.seiis.psic.xilb.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.pe.seiis.psic.xilb.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- rep(expo(pars[K_NB + 9]), K_NB) xi.est <- expo(pars[K_NB + 10] + pars[K_NB + 11]*age) pb.est <- expo(pars[K_NB + 12]) sb.est <- expo(pars[K_NB + 13]) c.est <- expo(pars[K_NB + 14] + pars[K_NB + 15]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 15. beta(t)phi(c)p(e)s(eiis)psi(lb)xi(c)eta(lb) mLL.betat.phic.pe.seiis.psilb.xic.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[K_NB + 9] + invect[K_NB + 10]*age) this.xi <- rep(expo(invect[K_NB + 11]), K_NB) this.pb <- expo(invect[K_NB + 12]) this.sb <- expo(invect[K_NB + 13]) this.c <- expo(invect[K_NB + 14] + invect[K_NB + 15]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.pe.seiis.psilb.xic.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.pe.seiis.psilb.xic.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[K_NB + 9] + pars[K_NB + 10]*age) xi.est <- rep(expo(pars[K_NB + 11]), K_NB) pb.est <- expo(pars[K_NB + 12]) sb.est <- expo(pars[K_NB + 13]) c.est <- expo(pars[K_NB + 14] + pars[K_NB + 15]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 16. beta(t)phi(c)p(e)s(eiis)psi(lb)xi(lb)eta(c) mLL.betat.phic.pe.seiis.psilb.xilb.etac <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[K_NB + 9] + invect[K_NB + 10]*age) this.xi <- expo(invect[K_NB + 11] + invect[K_NB + 12]*age) this.pb <- expo(invect[K_NB + 13]) this.sb <- expo(invect[K_NB + 14]) this.c <- rep(expo(invect[K_NB + 15]), K_NB) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.pe.seiis.psilb.xilb.etac <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.pe.seiis.psilb.xilb.etac, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[K_NB + 9] + pars[K_NB + 10]*age) xi.est <- expo(pars[K_NB + 11] + pars[K_NB + 12]*age) pb.est <- expo(pars[K_NB + 13]) sb.est <- expo(pars[K_NB + 14]) c.est <- rep(expo(pars[K_NB + 15]), K_NB) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 17. beta(t)phi(c)p(e)s(e)psi(lb)xi(c)eta(lb) mLL.betat.phic.pe.se.psilb.xic.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[K_NB + 7] + invect[K_NB + 8]*age) this.xi <- rep(expo(invect[K_NB + 9]), K_NB) this.pb <- expo(invect[K_NB + 10]) this.sb <- expo(invect[K_NB + 11]) this.c <- expo(invect[K_NB + 12] + invect[K_NB + 13]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.pe.se.psilb.xic.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.pe.se.psilb.xic.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 100000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[K_NB + 7] + pars[K_NB + 8]*age) xi.est <- rep(expo(pars[K_NB + 9]), K_NB) pb.est <- expo(pars[K_NB + 10]) sb.est <- expo(pars[K_NB + 11]) c.est <- expo(pars[K_NB + 12] + pars[K_NB + 13]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 18. beta(t)phi(c)p(e)s(eiis)psi(c)xi(c)eta(lb) mLL.betat.phic.pe.seiis.psic.xic.etalb <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- rep(expo(invect[K_NB + 9]), K_NB) this.xi <- rep(expo(invect[K_NB + 10]), K_NB) this.pb <- expo(invect[K_NB + 11]) this.sb <- expo(invect[K_NB + 12]) this.c <- expo(invect[K_NB + 13] + invect[K_NB + 14]*age) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.pe.seiis.psic.xic.etalb <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.pe.seiis.psic.xic.etalb, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 10000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- rep(expo(pars[K_NB + 9]), K_NB) xi.est <- rep(expo(pars[K_NB + 10]), K_NB) pb.est <- expo(pars[K_NB + 11]) sb.est <- expo(pars[K_NB + 12]) c.est <- expo(pars[K_NB + 13] + pars[K_NB + 14]*age) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 19. beta(t)phi(c)p(e)s(eiis)psi(lb)xi(c)eta(c) mLL.betat.phic.pe.seiis.psilb.xic.etac <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[K_NB + 9] + invect[K_NB + 10]*age) this.xi <- rep(expo(invect[K_NB + 11]), K_NB) this.pb <- expo(invect[K_NB + 12]) this.sb <- expo(invect[K_NB + 13]) this.c <- rep(expo(invect[K_NB + 14]), K_NB) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.pe.seiis.psilb.xic.etac <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.pe.seiis.psilb.xic.etac, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 10000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[K_NB + 9] + pars[K_NB + 10]*age) xi.est <- rep(expo(pars[K_NB + 11]), K_NB) pb.est <- expo(pars[K_NB + 12]) sb.est <- expo(pars[K_NB + 13]) c.est <- rep(expo(pars[K_NB + 14]), K_NB) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 20. beta(t)phi(c)p(e)s(e)psi(lb)xi(c)eta(c) mLL.betat.phic.pe.se.psilb.xic.etac <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- expo(invect[K_NB + 7] + invect[K_NB + 8]*age) this.xi <- rep(expo(invect[K_NB + 9]), K_NB) this.pb <- expo(invect[K_NB + 10]) this.sb <- expo(invect[K_NB + 11]) this.c <- rep(expo(invect[K_NB + 12]), K_NB) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.pe.se.psilb.xic.etac <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.pe.se.psilb.xic.etac, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 10000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- expo(pars[K_NB + 7] + pars[K_NB + 8]*age) xi.est <- rep(expo(pars[K_NB + 9]), K_NB) pb.est <- expo(pars[K_NB + 10]) sb.est <- expo(pars[K_NB + 11]) c.est <- rep(expo(pars[K_NB + 12]), K_NB) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #------------------------------------------------------------------------------------------------# # 21. beta(t)phi(c)p(e)s(eiis)psi(c)xi(c)eta(c) mLL.betat.phic.pe.seiis.psic.xic.etac <- function(invect) { this.NM <- exp(invect[1]) this.NU <- exp(invect[2]) this.beta <- c(lgp(invect[3:(K_NB + 1)])) this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE); this.phi[lower.tri(this.phi)] <- 0 this.p <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p[ ,what[1:K_NB] == 2] <- 0 this.p[lower.tri(this.p)] <- 0 this.p1 <- matrix(expo(invect[K_NB + 3] + invect[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.p1[ ,what[1:K_NB] == 2] <- 0 this.p1[lower.tri(this.p1)] <- 0 this.s <- matrix(expo(invect[K_NB + 5] + invect[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s[ ,what[1:K_NB] == 1] <- 0 this.s[lower.tri(this.s)] <- 0 this.s1 <- matrix(expo(invect[K_NB + 7] + invect[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE); this.s1[ ,what[1:K_NB] == 1] <- 0 this.s1[lower.tri(this.s1)] <- 0 this.psi <- rep(expo(invect[K_NB + 9]), K_NB) this.xi <- rep(expo(invect[K_NB + 10]), K_NB) this.pb <- expo(invect[K_NB + 11]) this.sb <- expo(invect[K_NB + 12]) this.c <- rep(expo(invect[K_NB + 13]), K_NB) LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1, this.psi, this.xi, this.pb, this.sb, this.c) } fit.betat.phic.pe.seiis.psic.xic.etac <- function(pars.start){ this.fit <- try(optim(par = pars.start, fn = mLL.betat.phic.pe.seiis.psic.xic.etac, lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)), method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 10000)), TRUE) minvalue <- this.fit$value hess <- this.fit$hessian npar <- length(this.fit$par) AIC <- 2*minvalue + 2*npar pars <- this.fit$par NM.est <- exp(pars[1]) NU.est <- exp(pars[2]) beta.est <- c(lgp(pars[3:(K_NB + 1)])) phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE) phi.est[lower.tri(phi.est)] <- 0 p.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p.est[ ,what[1:K_NB] == 2] <- 0 p.est[lower.tri(p.est)] <- 0 p1.est <- matrix(expo(pars[K_NB + 3] + pars[K_NB + 4]*effortc), nrow = K_NB, ncol = K_NB, byrow = TRUE) p1.est[ ,what[1:K_NB] == 2] <- 0 p1.est[lower.tri(p1.est)] <- 0 s.est <- matrix(expo(pars[K_NB + 5] + pars[K_NB + 6]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s.est[ ,what[1:K_NB] == 1] <- 0 s.est[lower.tri(s.est)] <- 0 s1.est <- matrix(expo(pars[K_NB + 7] + pars[K_NB + 8]*effortr), nrow = K_NB, ncol = K_NB, byrow = TRUE) s1.est[ ,what[1:K_NB] == 1] <- 0 s1.est[lower.tri(s1.est)] <- 0 psi.est <- rep(expo(pars[K_NB + 9]), K_NB) xi.est <- rep(expo(pars[K_NB + 10]), K_NB) pb.est <- expo(pars[K_NB + 11]) sb.est <- expo(pars[K_NB + 12]) c.est <- rep(expo(pars[K_NB + 13]), K_NB) list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar, "NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3), "phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3), "s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3), "xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3), "c.est" = round(c.est, 3), "pars" = pars) } #---------------------------- data simulation function ---------------------------------------------------------# ######################################## ALL<-function(x,y){ all(x==y)} ######################################### APPLY<-function(x,y){ apply(y,1,ALL,x=x)} ######################################### Length<-function(x){ length(x[x==T])} ####################################### PASTE<-function(x){paste(x,collapse="")} ######################################## data.simulation <- function(K_NB,N_M,N_U,betta,phi,p,p1,s,s1,psi,xi,pb,sb,c,what) ## { Ind <- round((N_M+N_U)*betta) #numbers entering at each occasion K <- K_NB + 4 b <- rep(1:K,Ind) #times of entry of each individual cumInd <- cumsum(Ind) cumInd <- c(0,cumInd) Alive <- matrix(0,nrow=max(cumInd),ncol=K) #matrix of life histories CH <- matrix(0,nrow=max(cumInd),ncol=K) #matrix of capture histories RH <- matrix(0,nrow=max(cumInd),ncol=K) #matrix of resighting histories marked <- sample(c(rep(1, N_M), rep(0, N_U)), N_M+N_U, replace = F) #indicator vector of being previously marked breed <- rep(0, max(cumInd)) #indicator vector of breeding success <- rep(0, max(cumInd)) chicks <- rep(0, max(cumInd)) for(j in 1:K){ if((cumInd[(j)]!=cumInd[j+1])==T){ Alive[((cumInd[(j)]+1):cumInd[j+1]),j] <- 1 } } for(i in 1:max(cumInd)){ for(j in 1:(K-1)){ if(Alive[i,j]==1) (Alive[i,(j+1)] <- rbinom(1,1,phi[(b[i]),j])); } } for(i in 1:max(cumInd)){ for(j in 1:K_NB){ if(Alive[i,j]==1) { if(any((CH[i,1:(j-1)]==1) | (CH[i,1:(j-1)]==3))) CH[i,j] <- rbinom(1,1, p1[b[i],j]) else CH[i,j] <- rbinom(1,1,p[b[i],j]); if(any((RH[i,1:(j-1)]==1))) RH[i,j] <- rbinom(1,1, s1[b[i],j]) else RH[i,j] <- rbinom(1,1,s[b[i],j]); if(((sum(CH[i,1:j])>0)==TRUE)|(marked[i]==1)) if(RH[i,j]==1){if(CH[i,j]==1) CH[i,j] <- 3 else CH[i,j] <- 2}; } } } for(i in 1:max(cumInd)){ if(Alive[i,(K_NB)]==1) { breed[i] <- rbinom(1,1,psi[(b[i])]); #did individual i breed? if((breed[i]==1)) success[i] <- rbinom(1,1,xi[(b[i])]); #did individual i breed successfully? if (success[i]==1) chicks[i] <- rbinom(1,1, c[(b[i])]); #did individual i have chicks? if((breed[i]==1)&((sum(CH[i,1:(K_NB)])>0)==T)|(marked[i]==1)) ##for birds (previously marked or marked this season) that breed { CH[i,(K_NB+1)] <- 2*rbinom(1,1,sb); #if it did breed, and it was (previously marked or only marked this season) then it will be detected if (success[i]==1) CH[i,(K_NB+2)]<- rbinom(1,1, pb); if (any(CH[i,(K_NB+1):(K_NB+2)]>0)) CH[i,(K_NB+4)] <- success[i] else CH[i,(K_NB+4)] <- 0; ##this is if we observe nesting for birds we know are there. 4 is we did, 14 is if we didn't if (any(CH[i,(K_NB+1):(K_NB+2)]>0)) CH[i,(K_NB+3)] <- chicks[i] else CH[i,(K_NB+3)] <- 0; }; if((breed[i]==1)&((sum(CH[i,1:(K_NB)])>0)==F)&(marked[i]==0)) ##for birds (unmarked and undetected till K_NB+1) that breed { CH[i,(K_NB+1)]<- 0 #if it did breed and it wasn't (previously marked or marked this season) then it will not be detected if (success[i]==1) CH[i,(K_NB+2)]<- rbinom(1,1, pb); if (any(CH[i,(K_NB+1):(K_NB+2)]>0)) CH[i,(K_NB+4)] <- success[i] else CH[i,(K_NB+4)] <- 0; ##this is if we observe nesting for birds we know are there. 4 is we did, 14 is if we didn't if (any(CH[i,(K_NB+1):(K_NB+2)]>0)) CH[i,(K_NB+3)] <- chicks[i] else CH[i,(K_NB+3)] <- 0; }; if(breed[i]==0) CH[i,(K_NB+1):(K_NB+4)] <- 0; #if it didn't breed then it will not be detected during breeding season } } CH1 <- cbind(CH, breed, success, chicks) CH_M <- CH[marked==1,][(apply(CH[marked==1,][,(1:K)],1,sum)>0),] CH_U <- CH[marked==0,][(apply(CH[marked==0,][,(1:K)],1,sum)>0),] UCH_M <- matrix(unique(CH_M)[apply(unique(CH_M),1,ALL,y=0)==FALSE],ncol=K,byrow=FALSE) y_M <- matrix(CH_M[apply(CH_M,1,ALL,y=0)==FALSE],ncol=K,byrow=FALSE) n_M <- apply(apply(UCH_M,1,APPLY,y_M),2,Length) UCH_U <- matrix(unique(CH_U)[apply(unique(CH_U),1,ALL,y=0)==FALSE],ncol=K,byrow=FALSE) y_U <- matrix(CH_U[apply(CH_U,1,ALL,y=0)==FALSE],ncol=K,byrow=FALSE) n_U <- apply(apply(UCH_U,1,APPLY,y_U),2,Length) UCH_M <- cbind(UCH_M,n_M) UCH_U <- cbind(UCH_U,n_U) NM <- sum(marked) temp <- apply(Alive[marked==0,what[1:K_NB]==1],1,sum) NU <- length(temp[temp>0]) list("UCH_M"=UCH_M,"UCH_U"=UCH_U,"NM"=NM, "NU"=NU) }