source("AllFunctions.R") K_NB <- 32 K_B <- 4 K <- K_NB + K_B what <- c(1,1,1,2,1,3,3,2,1,3,2,3,3,3,3,3,1,3,2,2,3,2,3,3,2,3,3,3,2,3,1,1,0,0,0,0) #type of sampling used: 1 is capture only, 2 is both, # 3 is resighting only, 0 indicates an occasion during the breeding season effort <- read.csv("effort.csv", header=TRUE) effortc <- scale(as.vector(effort$catch))[,1] # capture effort effortr <- scale(as.vector(effort$resight))[,1] # resight effort age <- scale(1:K_NB)[,1] # age as a covariate UCH_M <- as.matrix(read.csv("PMdata.csv", header=F)) # data on already marked birds colnames(UCH_M) <- NULL UCH_U <- as.matrix(read.csv("PUdata.csv", header=F)) # data on previously unmarked birds colnames(UCH_U) <- NULL # create capture history matrices x.matM <- UCH_M[,1:K] y.vectM <- UCH_M[,K+1] DDM <- sum(y.vectM) ##409 marked birds minLogNM <- log(DDM) x.matU <- UCH_U[,1:K] y.vectU <- UCH_U[,K+1] DDU <- sum(y.vectU) ##232 unmarked birds minLogNU <- log(DDU) # Fit all models mentioned in the model selection process, extract estimates and asymptotic CIs of estimates #--------------------------------------------------------------------------------------------------------------------------------------------# #1. beta(t)phi(t)p(eiis)s(eiis)psi(lb)xi(lb)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(K_NB-1, 0.01, 0.99)), rep(0, 16)) npar <- length(pars.start) betat.phit.peiis.seiis.psilb.xilb.etalb.out <- fit.betat.phit.peiis.seiis.psilb.xilb.etalb(pars.start) VC <- diag(solve(betat.phit.peiis.seiis.psilb.xilb.etalb.out$hessian)) #phi expo(betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+2):(2*K_NB)] - 1.96*sqrt(VC[(K_NB+2):(2*K_NB)])) expo(betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+2):(2*K_NB)]) expo(betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+2):(2*K_NB)] + 1.96*sqrt(VC[(K_NB+2):(2*K_NB)])) #p, p1, s, s1 betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+1):(2*K_NB+8)] - 1.96*sqrt(VC[(2*K_NB+1):(2*K_NB+8)]) betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+1):(2*K_NB+8)] betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+1):(2*K_NB+8)] + 1.96*sqrt(VC[(2*K_NB+1):(2*K_NB+8)]) #psi betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+9):(2*K_NB+10)] - 1.96*sqrt(VC[(2*K_NB+9):(2*K_NB+10)]) betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+9):(2*K_NB+10)] betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+9):(2*K_NB+10)] + 1.96*sqrt(VC[(2*K_NB+9):(2*K_NB+10)]) #xi betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+11):(2*K_NB+12)] - 1.96*sqrt(VC[(2*K_NB+11):(2*K_NB+12)]) betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+11):(2*K_NB+12)] betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+11):(2*K_NB+12)] + 1.96*sqrt(VC[(2*K_NB+11):(2*K_NB+12)]) #pb expo(betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+13)] - 1.96*sqrt(VC[(2*K_NB+13)])) expo(betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+13)]) expo(betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+13)] + 1.96*sqrt(VC[(2*K_NB+13)])) #sb expo(betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+14)] - 1.96*sqrt(VC[(2*K_NB+14)])) expo(betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+14)]) expo(betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+14)] + 1.96*sqrt(VC[(2*K_NB+14)])) #eta betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+15):(2*K_NB+16)] - 1.96*sqrt(VC[(2*K_NB+15):(2*K_NB+16)]) betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+15):(2*K_NB+16)] betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+15):(2*K_NB+16)] + 1.96*sqrt(VC[(2*K_NB+15):(2*K_NB+16)]) betat.phit.peiis.seiis.psilb.xilb.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #2. beta(t)phi(c)p(eiis)s(eiis)psi(lb)xi(lb)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 16)) npar <- length(pars.start) betat.phic.peiis.seiis.psilb.xilb.etalb.out <- fit.betat.phic.peiis.seiis.psilb.xilb.etalb(pars.start) VC <- diag(solve(betat.phic.peiis.seiis.psilb.xilb.etalb.out$hessian)) #phi expo(betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+2)]) expo(betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+3):(K_NB+10)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+10)]) betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+3):(K_NB+10)] betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+3):(K_NB+10)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+10)]) #psi betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+11):(K_NB+12)] - 1.96*sqrt(VC[(K_NB+11):(K_NB+12)]) betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+11):(K_NB+12)] betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+11):(K_NB+12)] + 1.96*sqrt(VC[(K_NB+11):(K_NB+12)]) #xi betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+13):(K_NB+14)] - 1.96*sqrt(VC[(K_NB+13):(K_NB+14)]) betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+13):(K_NB+14)] betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+13):(K_NB+14)] + 1.96*sqrt(VC[(K_NB+13):(K_NB+14)]) #pb expo(betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+15)] - 1.96*sqrt(VC[(K_NB+15)])) expo(betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+15)]) expo(betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+15)] + 1.96*sqrt(VC[(K_NB+15)])) #sb expo(betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+16)] - 1.96*sqrt(VC[(K_NB+16)])) expo(betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+16)]) expo(betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+16)] + 1.96*sqrt(VC[(K_NB+16)])) #eta betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+17):(K_NB+18)] - 1.96*sqrt(VC[(K_NB+17):(K_NB+18)]) betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+17):(K_NB+18)] betat.phic.peiis.seiis.psilb.xilb.etalb.out$pars[(K_NB+17):(K_NB+18)] + 1.96*sqrt(VC[(K_NB+17):(K_NB+18)]) betat.phic.peiis.seiis.psilb.xilb.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #3. beta(t)phi(t)p(e)s(eiis)psi(lb)xi(lb)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(K_NB-1, 0.01, 0.99)), rep(0, 14)) npar <- length(pars.start) betat.phit.pe.seiis.psilb.xilb.etalb.out <- fit.betat.phit.pe.seiis.psilb.xilb.etalb(pars.start) VC <- diag(solve(betat.phit.pe.seiis.psilb.xilb.etalb.out$hessian)) #phi expo(betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+2):(2*K_NB)] - 1.96*sqrt(VC[(K_NB+2):(2*K_NB)])) expo(betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+2):(2*K_NB)]) expo(betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+2):(2*K_NB)] + 1.96*sqrt(VC[(K_NB+2):(2*K_NB)])) #p, p1, s, s1 betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+1):(2*K_NB+6)] - 1.96*sqrt(VC[(2*K_NB+1):(2*K_NB+6)]) betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+1):(2*K_NB+6)] betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+1):(2*K_NB+6)] + 1.96*sqrt(VC[(2*K_NB+1):(2*K_NB+6)]) #psi betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+7):(2*K_NB+8)] - 1.96*sqrt(VC[(2*K_NB+7):(2*K_NB+8)]) betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+7):(2*K_NB+8)] betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+7):(2*K_NB+8)] + 1.96*sqrt(VC[(2*K_NB+7):(2*K_NB+8)]) #xi betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+9):(2*K_NB+10)] - 1.96*sqrt(VC[(2*K_NB+9):(2*K_NB+10)]) betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+9):(2*K_NB+10)] betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+9):(2*K_NB+10)] + 1.96*sqrt(VC[(2*K_NB+9):(2*K_NB+10)]) #pb expo(betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+11)] - 1.96*sqrt(VC[(2*K_NB+11)])) expo(betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+11)]) expo(betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+11)] + 1.96*sqrt(VC[(2*K_NB+11)])) #sb expo(betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+12)] - 1.96*sqrt(VC[(2*K_NB+12)])) expo(betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+12)]) expo(betat.phit.pe.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+12)] + 1.96*sqrt(VC[(2*K_NB+12)])) #eta betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+13):(2*K_NB+14)] - 1.96*sqrt(VC[(2*K_NB+13):(2*K_NB+14)]) betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+13):(2*K_NB+14)] betat.phit.peiis.seiis.psilb.xilb.etalb.out$pars[(2*K_NB+13):(2*K_NB+14)] + 1.96*sqrt(VC[(2*K_NB+13):(2*K_NB+14)]) betat.phit.pe.seiis.psilb.xilb.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #4. beta(t)phi(t)p(eiis)s(e)psi(lb)xi(lb)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(K_NB-1, 0.01, 0.99)), rep(0, 14)) npar <- length(pars.start) betat.phit.peiis.se.psilb.xilb.etalb.out <- fit.betat.phit.peiis.se.psilb.xilb.etalb(pars.start) VC <- diag(solve(betat.phit.peiis.se.psilb.xilb.etalb.out$hessian)) #phi expo(betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+2):(2*K_NB)] - 1.96*sqrt(VC[(K_NB+2):(2*K_NB)])) expo(betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+2):(2*K_NB)]) expo(betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+2):(2*K_NB)] + 1.96*sqrt(VC[(K_NB+2):(2*K_NB)])) #p, p1, s, s1 betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+1):(2*K_NB+6)] - 1.96*sqrt(VC[(2*K_NB+1):(2*K_NB+6)]) betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+1):(2*K_NB+6)] betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+1):(2*K_NB+6)] + 1.96*sqrt(VC[(2*K_NB+1):(2*K_NB+6)]) #psi betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+7):(2*K_NB+8)] - 1.96*sqrt(VC[(2*K_NB+7):(2*K_NB+8)]) betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+7):(2*K_NB+8)] betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+7):(2*K_NB+8)] + 1.96*sqrt(VC[(2*K_NB+7):(2*K_NB+8)]) #xi betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+9):(2*K_NB+10)] - 1.96*sqrt(VC[(2*K_NB+9):(2*K_NB+10)]) betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+9):(2*K_NB+10)] betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+9):(2*K_NB+10)] + 1.96*sqrt(VC[(2*K_NB+9):(2*K_NB+10)]) #pb expo(betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+11)] - 1.96*sqrt(VC[(2*K_NB+11)])) expo(betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+11)]) expo(betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+11)] + 1.96*sqrt(VC[(2*K_NB+11)])) #sb expo(betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+12)] - 1.96*sqrt(VC[(2*K_NB+12)])) expo(betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+12)]) expo(betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+12)] + 1.96*sqrt(VC[(2*K_NB+12)])) #eta betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+13):(2*K_NB+14)] - 1.96*sqrt(VC[(2*K_NB+13):(2*K_NB+14)]) betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+13):(2*K_NB+14)] betat.phit.peiis.se.psilb.xilb.etalb.out$pars[(2*K_NB+13):(2*K_NB+14)] + 1.96*sqrt(VC[(2*K_NB+13):(2*K_NB+14)]) betat.phit.peiis.se.psilb.xilb.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #5. beta(t)phi(t)p(eiis)s(eiis)psi(c)xi(lb)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(K_NB-1, 0.01, 0.99)), rep(0, 15)) npar <- length(pars.start) betat.phit.peiis.seiis.psic.xilb.etalb.out <- fit.betat.phit.peiis.seiis.psic.xilb.etalb(pars.start) VC <- diag(solve(betat.phit.peiis.seiis.psic.xilb.etalb.out$hessian)) #phi expo(betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+2):(2*K_NB)] - 1.96*sqrt(VC[(K_NB+2):(2*K_NB)])) expo(betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+2):(2*K_NB)]) expo(betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+2):(2*K_NB)] + 1.96*sqrt(VC[(K_NB+2):(2*K_NB)])) #p, p1, s, s1 betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+1):(2*K_NB+8)] - 1.96*sqrt(VC[(2*K_NB+1):(2*K_NB+8)]) betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+1):(2*K_NB+8)] betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+1):(2*K_NB+8)] + 1.96*sqrt(VC[(2*K_NB+1):(2*K_NB+8)]) #psi expo(betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+9)] - 1.96*sqrt(VC[(2*K_NB+9)])) expo(betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+9)]) expo(betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+9)] + 1.96*sqrt(VC[(2*K_NB+9)])) #xi betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+10):(2*K_NB+11)] - 1.96*sqrt(VC[(2*K_NB+10):(2*K_NB+11)]) betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+10):(2*K_NB+11)] betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+10):(2*K_NB+11)] + 1.96*sqrt(VC[(2*K_NB+10):(2*K_NB+11)]) #pb expo(betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+12)] - 1.96*sqrt(VC[(2*K_NB+12)])) expo(betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+12)]) expo(betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+12)] + 1.96*sqrt(VC[(2*K_NB+12)])) #sb expo(betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+13)] - 1.96*sqrt(VC[(2*K_NB+13)])) expo(betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+13)]) expo(betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+13)] + 1.96*sqrt(VC[(2*K_NB+13)])) #eta betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+14):(2*K_NB+15)] - 1.96*sqrt(VC[(2*K_NB+14):(2*K_NB+15)]) betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+14):(2*K_NB+15)] betat.phit.peiis.seiis.psic.xilb.etalb.out$pars[(2*K_NB+14):(2*K_NB+15)] + 1.96*sqrt(VC[(2*K_NB+14):(2*K_NB+15)]) betat.phit.peiis.seiis.psic.xilb.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #6. beta(t)phi(t)p(eiis)s(eiis)psi(lb)xi(c)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(K_NB-1, 0.01, 0.99)), rep(0, 15)) npar <- length(pars.start) betat.phit.peiis.seiis.psilb.xic.etalb.out <- fit.betat.phit.peiis.seiis.psilb.xic.etalb(pars.start) VC <- diag(solve(betat.phit.peiis.seiis.psilb.xic.etalb.out$hessian)) #phi expo(betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+2):(2*K_NB)] - 1.96*sqrt(VC[(K_NB+2):(2*K_NB)])) expo(betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+2):(2*K_NB)]) expo(betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+2):(2*K_NB)] + 1.96*sqrt(VC[(K_NB+2):(2*K_NB)])) #p, p1, s, s1 betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+1):(2*K_NB+8)] - 1.96*sqrt(VC[(2*K_NB+1):(2*K_NB+8)]) betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+1):(2*K_NB+8)] betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+1):(2*K_NB+8)] + 1.96*sqrt(VC[(2*K_NB+1):(2*K_NB+8)]) #psi betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+9):(2*K_NB+10)] - 1.96*sqrt(VC[(2*K_NB+9):(2*K_NB+10)]) betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+9):(2*K_NB+10)] betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+9):(2*K_NB+10)] + 1.96*sqrt(VC[(2*K_NB+9):(2*K_NB+10)]) #xi expo(betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+11)] - 1.96*sqrt(VC[(2*K_NB+11)])) expo(betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+11)]) expo(betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+11)] + 1.96*sqrt(VC[(2*K_NB+11)])) #pb expo(betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+12)] - 1.96*sqrt(VC[(2*K_NB+12)])) expo(betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+12)]) expo(betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+12)] + 1.96*sqrt(VC[(2*K_NB+12)])) #sb expo(betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+13)] - 1.96*sqrt(VC[(2*K_NB+13)])) expo(betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+13)]) expo(betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+13)] + 1.96*sqrt(VC[(2*K_NB+13)])) #eta betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+14):(2*K_NB+15)] - 1.96*sqrt(VC[(2*K_NB+14):(2*K_NB+15)]) betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+14):(2*K_NB+15)] betat.phit.peiis.seiis.psilb.xic.etalb.out$pars[(2*K_NB+14):(2*K_NB+15)] + 1.96*sqrt(VC[(2*K_NB+14):(2*K_NB+15)]) betat.phit.peiis.seiis.psilb.xic.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #7. beta(t)phi(t)p(eiis)s(eiis)psi(lb)xi(lb)eta(c) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(K_NB-1, 0.01, 0.99)), rep(0, 15)) npar <- length(pars.start) betat.phit.peiis.seiis.psilb.xilb.etac.out <- fit.betat.phit.peiis.seiis.psilb.xilb.etac(pars.start) VC <- diag(solve(betat.phit.peiis.seiis.psilb.xilb.etac.out$hessian)) #phi expo(betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+2):(2*K_NB)] - 1.96*sqrt(VC[(K_NB+2):(2*K_NB)])) expo(betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+2):(2*K_NB)]) expo(betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+2):(2*K_NB)] + 1.96*sqrt(VC[(K_NB+2):(2*K_NB)])) #p, p1, s, s1 betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+1):(2*K_NB+8)] - 1.96*sqrt(VC[(2*K_NB+1):(2*K_NB+8)]) betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+1):(2*K_NB+8)] betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+1):(2*K_NB+8)] + 1.96*sqrt(VC[(2*K_NB+1):(2*K_NB+8)]) #psi betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+9):(2*K_NB+10)] - 1.96*sqrt(VC[(2*K_NB+9):(2*K_NB+10)]) betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+9):(2*K_NB+10)] betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+9):(2*K_NB+10)] + 1.96*sqrt(VC[(2*K_NB+9):(2*K_NB+10)]) #xi betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+11):(2*K_NB+12)] - 1.96*sqrt(VC[(2*K_NB+11):(2*K_NB+12)]) betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+11):(2*K_NB+12)] betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+11):(2*K_NB+12)] + 1.96*sqrt(VC[(2*K_NB+11):(2*K_NB+12)]) #pb expo(betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+13)] - 1.96*sqrt(VC[(2*K_NB+13)])) expo(betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+13)]) expo(betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+13)] + 1.96*sqrt(VC[(2*K_NB+13)])) #sb expo(betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+14)] - 1.96*sqrt(VC[(2*K_NB+14)])) expo(betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+14)]) expo(betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+14)] + 1.96*sqrt(VC[(2*K_NB+14)])) #eta expo(betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+15)] - 1.96*sqrt(VC[(2*K_NB+15)])) expo(betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+15)]) expo(betat.phit.peiis.seiis.psilb.xilb.etac.out$pars[(2*K_NB+15)] + 1.96*sqrt(VC[(2*K_NB+15)])) betat.phit.peiis.seiis.psilb.xilb.etac.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #8. beta(t)phi(c)p(e)s(eiis)psi(lb)xi(lb)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 14)) npar <- length(pars.start) betat.phic.pe.seiis.psilb.xilb.etalb.out <- fit.betat.phic.pe.seiis.psilb.xilb.etalb(pars.start) VC <- diag(solve(betat.phic.pe.seiis.psilb.xilb.etalb.out$hessian)) #N_M exp(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(1)] - 1.96*sqrt(VC[1])) exp(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(1)]) exp(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(1)] + 1.96*sqrt(VC[1])) #N_U exp(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(2)] - 1.96*sqrt(VC[2])) exp(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(2)]) exp(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(2)] + 1.96*sqrt(VC[2])) #phi expo(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+2)]) expo(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+3):(K_NB+8)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+3):(K_NB+8)] betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+3):(K_NB+8)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) #psi betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+9):(K_NB+10)] - 1.96*sqrt(VC[(K_NB+9):(K_NB+10)]) betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+9):(K_NB+10)] betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+9):(K_NB+10)] + 1.96*sqrt(VC[(K_NB+9):(K_NB+10)]) #xi betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+11):(K_NB+12)] - 1.96*sqrt(VC[(K_NB+11):(K_NB+12)]) betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+11):(K_NB+12)] betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+11):(K_NB+12)] + 1.96*sqrt(VC[(K_NB+11):(K_NB+12)]) #pb expo(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+13)] - 1.96*sqrt(VC[(K_NB+13)])) expo(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+13)]) expo(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+13)] + 1.96*sqrt(VC[(K_NB+13)])) #sb expo(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+14)] - 1.96*sqrt(VC[(K_NB+14)])) expo(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+14)]) expo(betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+14)] + 1.96*sqrt(VC[(K_NB+14)])) #eta betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+15):(K_NB+16)] - 1.96*sqrt(VC[(K_NB+15):(K_NB+16)]) betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+15):(K_NB+16)] betat.phic.pe.seiis.psilb.xilb.etalb.out$pars[(K_NB+15):(K_NB+16)] + 1.96*sqrt(VC[(K_NB+15):(K_NB+16)]) betat.phic.pe.seiis.psilb.xilb.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #9. beta(t)phi(c)p(eiis)s(e)psi(lb)xi(lb)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 14)) npar <- length(pars.start) betat.phic.peiis.se.psilb.xilb.etalb.out <- fit.betat.phic.peiis.se.psilb.xilb.etalb(pars.start) VC <- diag(solve(betat.phic.peiis.se.psilb.xilb.etalb.out$hessian)) #phi expo(betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+2)]) expo(betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+3):(K_NB+8)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+3):(K_NB+8)] betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+3):(K_NB+8)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) #psi betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+9):(K_NB+10)] - 1.96*sqrt(VC[(K_NB+9):(K_NB+10)]) betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+9):(K_NB+10)] betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+9):(K_NB+10)] + 1.96*sqrt(VC[(K_NB+9):(K_NB+10)]) #xi betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+11):(K_NB+12)] - 1.96*sqrt(VC[(K_NB+11):(K_NB+12)]) betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+11):(K_NB+12)] betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+11):(K_NB+12)] + 1.96*sqrt(VC[(K_NB+11):(K_NB+12)]) #pb expo(betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+13)] - 1.96*sqrt(VC[(K_NB+13)])) expo(betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+13)]) expo(betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+13)] + 1.96*sqrt(VC[(K_NB+13)])) #sb expo(betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+14)] - 1.96*sqrt(VC[(K_NB+14)])) expo(betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+14)]) expo(betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+14)] + 1.96*sqrt(VC[(K_NB+14)])) #eta betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+15):(K_NB+16)] - 1.96*sqrt(VC[(K_NB+15):(K_NB+16)]) betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+15):(K_NB+16)] betat.phic.peiis.se.psilb.xilb.etalb.out$pars[(K_NB+15):(K_NB+16)] + 1.96*sqrt(VC[(K_NB+15):(K_NB+16)]) betat.phic.peiis.se.psilb.xilb.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #10. beta(t)phi(c)p(eiis)s(eiis)psi(c)xi(lb)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 15)) npar <- length(pars.start) betat.phic.peiis.seiis.psic.xilb.etalb.out <- fit.betat.phic.peiis.seiis.psic.xilb.etalb(pars.start) VC <- diag(solve(betat.phic.peiis.seiis.psic.xilb.etalb.out$hessian)) #phi expo(betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+2)]) expo(betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+3):(K_NB+10)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+10)]) betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+3):(K_NB+10)] betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+3):(K_NB+10)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+10)]) #psi betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+11)] - 1.96*sqrt(VC[(K_NB+11)]) betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+11)] betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+11)] + 1.96*sqrt(VC[(K_NB+11)]) #xi betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+12):(K_NB+13)] - 1.96*sqrt(VC[(K_NB+12):(K_NB+13)]) betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+12):(K_NB+13)] betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+12):(K_NB+13)] + 1.96*sqrt(VC[(K_NB+12):(K_NB+13)]) #pb expo(betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+14)] - 1.96*sqrt(VC[(K_NB+14)])) expo(betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+14)]) expo(betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+14)] + 1.96*sqrt(VC[(K_NB+14)])) #sb expo(betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+15)] - 1.96*sqrt(VC[(K_NB+15)])) expo(betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+15)]) expo(betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+15)] + 1.96*sqrt(VC[(K_NB+15)])) #eta betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+16):(K_NB+17)] - 1.96*sqrt(VC[(K_NB+16):(K_NB+17)]) betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+16):(K_NB+17)] betat.phic.peiis.seiis.psic.xilb.etalb.out$pars[(K_NB+16):(K_NB+17)] + 1.96*sqrt(VC[(K_NB+16):(K_NB+17)]) betat.phic.peiis.seiis.psic.xilb.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #11. beta(t)phi(c)p(eiis)s(eiis)psi(lb)xi(c)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 15)) npar <- length(pars.start) betat.phic.peiis.seiis.psilb.xic.etalb.out <- fit.betat.phic.peiis.seiis.psilb.xic.etalb(pars.start) VC <- diag(solve(betat.phic.peiis.seiis.psilb.xic.etalb.out$hessian)) #phi expo(betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+2)]) expo(betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+3):(K_NB+10)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+10)]) betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+3):(K_NB+10)] betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+3):(K_NB+10)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+10)]) #psi betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+11):(K_NB+12)] - 1.96*sqrt(VC[(K_NB+11):(K_NB+12)]) betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+11):(K_NB+12)] betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+11):(K_NB+12)] + 1.96*sqrt(VC[(K_NB+11):(K_NB+12)]) #xi betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+13)] - 1.96*sqrt(VC[(K_NB+13)]) betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+13)] betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+13)] + 1.96*sqrt(VC[(K_NB+13)]) #pb expo(betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+14)] - 1.96*sqrt(VC[(K_NB+14)])) expo(betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+14)]) expo(betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+14)] + 1.96*sqrt(VC[(K_NB+14)])) #sb expo(betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+15)] - 1.96*sqrt(VC[(K_NB+15)])) expo(betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+15)]) expo(betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+15)] + 1.96*sqrt(VC[(K_NB+15)])) #eta betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+16):(K_NB+17)] - 1.96*sqrt(VC[(K_NB+16):(K_NB+17)]) betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+16):(K_NB+17)] betat.phic.peiis.seiis.psilb.xic.etalb.out$pars[(K_NB+16):(K_NB+17)] + 1.96*sqrt(VC[(K_NB+16):(K_NB+17)]) betat.phic.peiis.seiis.psilb.xic.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #12. beta(t)phi(c)p(eiis)s(eiis)psi(lb)xi(lb)eta(c) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 15)) npar <- length(pars.start) betat.phic.peiis.seiis.psilb.xilb.etac.out <- fit.betat.phic.peiis.seiis.psilb.xilb.etac(pars.start) VC <- diag(solve(betat.phic.peiis.seiis.psilb.xilb.etac.out$hessian)) #phi expo(betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+2)]) expo(betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+3):(K_NB+10)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+10)]) betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+3):(K_NB+10)] betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+3):(K_NB+10)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+10)]) #psi betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+11):(K_NB+12)] - 1.96*sqrt(VC[(K_NB+11):(K_NB+12)]) betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+11):(K_NB+12)] betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+11):(K_NB+12)] + 1.96*sqrt(VC[(K_NB+11):(K_NB+12)]) #xi betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+13):(K_NB+14)] - 1.96*sqrt(VC[(K_NB+13):(K_NB+14)]) betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+13):(K_NB+14)] betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+13):(K_NB+14)] + 1.96*sqrt(VC[(K_NB+13):(K_NB+14)]) #pb expo(betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+15)] - 1.96*sqrt(VC[(K_NB+15)])) expo(betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+15)]) expo(betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+15)] + 1.96*sqrt(VC[(K_NB+15)])) #sb expo(betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+16)] - 1.96*sqrt(VC[(K_NB+16)])) expo(betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+16)]) expo(betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+16)] + 1.96*sqrt(VC[(K_NB+16)])) #eta betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+17)] - 1.96*sqrt(VC[(K_NB+17)]) betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+17)] betat.phic.peiis.seiis.psilb.xilb.etac.out$pars[(K_NB+17)] + 1.96*sqrt(VC[(K_NB+17)]) betat.phic.peiis.seiis.psilb.xilb.etac.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #13. beta(t)phi(c)p(e)s(e)psi(lb)xi(lb)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 12)) npar <- length(pars.start) betat.phic.pe.se.psilb.xilb.etalb.out <- fit.betat.phic.pe.se.psilb.xilb.etalb(pars.start) VC <- diag(solve(betat.phic.pe.se.psilb.xilb.etalb.out$hessian)) #phi expo(betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+2)]) expo(betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+3):(K_NB+6)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+6)]) betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+3):(K_NB+6)] betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+3):(K_NB+6)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+6)]) #psi betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+7):(K_NB+8)] - 1.96*sqrt(VC[(K_NB+7):(K_NB+8)]) betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+7):(K_NB+8)] betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+7):(K_NB+8)] + 1.96*sqrt(VC[(K_NB+7):(K_NB+8)]) #xi betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+9):(K_NB+10)] - 1.96*sqrt(VC[(K_NB+9):(K_NB+10)]) betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+9):(K_NB+10)] betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+9):(K_NB+10)] + 1.96*sqrt(VC[(K_NB+9):(K_NB+10)]) #pb expo(betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+11)] - 1.96*sqrt(VC[(K_NB+11)])) expo(betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+11)]) expo(betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+11)] + 1.96*sqrt(VC[(K_NB+11)])) #sb expo(betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+12)] - 1.96*sqrt(VC[(K_NB+12)])) expo(betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+12)]) expo(betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+12)] + 1.96*sqrt(VC[(K_NB+12)])) #eta betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+13):(K_NB+14)] - 1.96*sqrt(VC[(K_NB+13):(K_NB+14)]) betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+13):(K_NB+14)] betat.phic.pe.se.psilb.xilb.etalb.out$pars[(K_NB+13):(K_NB+14)] + 1.96*sqrt(VC[(K_NB+13):(K_NB+14)]) betat.phic.pe.se.psilb.xilb.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #14. beta(t)phi(c)p(e)s(eiis)psi(c)xi(lb)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 13)) npar <- length(pars.start) betat.phic.pe.seiis.psic.xilb.etalb.out <- fit.betat.phic.pe.seiis.psic.xilb.etalb(pars.start) VC <- diag(solve(betat.phic.pe.seiis.psic.xilb.etalb.out$hessian)) #N_M exp(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(1)] - 1.96*sqrt(VC[1])) exp(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(1)]) exp(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(1)] + 1.96*sqrt(VC[1])) #N_U exp(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(2)] - 1.96*sqrt(VC[2])) exp(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(2)]) exp(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(2)] + 1.96*sqrt(VC[2])) #phi expo(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+2)]) expo(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+3):(K_NB+8)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+3):(K_NB+8)] betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+3):(K_NB+8)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) #psi betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+9)] - 1.96*sqrt(VC[(K_NB+9)]) betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+9)] betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+9)] + 1.96*sqrt(VC[(K_NB+9)]) #xi betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+10):(K_NB+11)] - 1.96*sqrt(VC[(K_NB+10):(K_NB+11)]) betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+10):(K_NB+11)] betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+10):(K_NB+11)] + 1.96*sqrt(VC[(K_NB+10):(K_NB+11)]) #pb expo(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+12)] - 1.96*sqrt(VC[(K_NB+12)])) expo(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+12)]) expo(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+12)] + 1.96*sqrt(VC[(K_NB+12)])) #sb expo(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+13)] - 1.96*sqrt(VC[(K_NB+13)])) expo(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+13)]) expo(betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+13)] + 1.96*sqrt(VC[(K_NB+13)])) #eta betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+14):(K_NB+15)] - 1.96*sqrt(VC[(K_NB+14):(K_NB+15)]) betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+14):(K_NB+15)] betat.phic.pe.seiis.psic.xilb.etalb.out$pars[(K_NB+14):(K_NB+15)] + 1.96*sqrt(VC[(K_NB+14):(K_NB+15)]) betat.phic.pe.seiis.psic.xilb.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #15. beta(t)phi(c)p(e)s(eiis)psi(lb)xi(c)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 13)) npar <- length(pars.start) betat.phic.pe.seiis.psilb.xic.etalb.out <- fit.betat.phic.pe.seiis.psilb.xic.etalb(pars.start) VC <- diag(solve(betat.phic.pe.seiis.psilb.xic.etalb.out$hessian)) #N_M exp(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(1)] - 1.96*sqrt(VC[1])) exp(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(1)]) exp(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(1)] + 1.96*sqrt(VC[1])) #N_U exp(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(2)] - 1.96*sqrt(VC[2])) exp(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(2)]) exp(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(2)] + 1.96*sqrt(VC[2])) #phi expo(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+2)]) expo(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+3):(K_NB+8)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+3):(K_NB+8)] betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+3):(K_NB+8)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) #psi betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+9):(K_NB+10)] - 1.96*sqrt(VC[(K_NB+9):(K_NB+10)]) betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+9):(K_NB+10)] betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+9):(K_NB+10)] + 1.96*sqrt(VC[(K_NB+9):(K_NB+10)]) #xi betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+11)] - 1.96*sqrt(VC[(K_NB+11)]) betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+11)] betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+11)] + 1.96*sqrt(VC[(K_NB+11)]) #pb expo(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+12)] - 1.96*sqrt(VC[(K_NB+12)])) expo(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+12)]) expo(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+12)] + 1.96*sqrt(VC[(K_NB+12)])) #sb expo(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+13)] - 1.96*sqrt(VC[(K_NB+13)])) expo(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+13)]) expo(betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+13)] + 1.96*sqrt(VC[(K_NB+13)])) #eta betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+14):(K_NB+15)] - 1.96*sqrt(VC[(K_NB+14):(K_NB+15)]) betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+14):(K_NB+15)] betat.phic.pe.seiis.psilb.xic.etalb.out$pars[(K_NB+14):(K_NB+15)] + 1.96*sqrt(VC[(K_NB+14):(K_NB+15)]) betat.phic.pe.seiis.psilb.xic.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #16. beta(t)phi(c)p(e)s(eiis)psi(lb)xi(lb)eta(c) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 13)) npar <- length(pars.start) betat.phic.pe.seiis.psilb.xilb.etac.out <- fit.betat.phic.pe.seiis.psilb.xilb.etac(pars.start) VC <- diag(solve(betat.phic.pe.seiis.psilb.xilb.etac.out$hessian)) #N_M exp(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(1)] - 1.96*sqrt(VC[1])) exp(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(1)]) exp(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(1)] + 1.96*sqrt(VC[1])) #N_U exp(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(2)] - 1.96*sqrt(VC[2])) exp(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(2)]) exp(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(2)] + 1.96*sqrt(VC[2])) #phi expo(betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+2)]) expo(betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+3):(K_NB+8)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+3):(K_NB+8)] betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+3):(K_NB+8)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) #psi betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+9):(K_NB+10)] - 1.96*sqrt(VC[(K_NB+9):(K_NB+10)]) betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+9):(K_NB+10)] betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+9):(K_NB+10)] + 1.96*sqrt(VC[(K_NB+9):(K_NB+10)]) #xi betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+11):(K_NB+12)] - 1.96*sqrt(VC[(K_NB+11):(K_NB+12)]) betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+11):(K_NB+12)] betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+11):(K_NB+12)] + 1.96*sqrt(VC[(K_NB+11):(K_NB+12)]) #pb expo(betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+13)] - 1.96*sqrt(VC[(K_NB+13)])) expo(betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+13)]) expo(betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+13)] + 1.96*sqrt(VC[(K_NB+13)])) #sb expo(betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+14)] - 1.96*sqrt(VC[(K_NB+14)])) expo(betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+14)]) expo(betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+14)] + 1.96*sqrt(VC[(K_NB+14)])) #eta betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+15)] - 1.96*sqrt(VC[(K_NB+15)]) betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+15)] betat.phic.pe.seiis.psilb.xilb.etac.out$pars[(K_NB+15)] + 1.96*sqrt(VC[(K_NB+15)]) betat.phic.pe.seiis.psilb.xilb.etac.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #16. beta(t)phi(c)p(e)s(e)psi(lb)xi(c)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 11)) npar <- length(pars.start) betat.phic.pe.se.psilb.xic.etalb.out <- fit.betat.phic.pe.se.psilb.xic.etalb(pars.start) VC <- diag(solve(betat.phic.pe.seiis.psilb.xic.etalb.out$hessian)) #phi expo(betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+2)]) expo(betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+3):(K_NB+6)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+6)]) betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+3):(K_NB+6)] betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+3):(K_NB+6)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+6)]) #psi betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+7):(K_NB+8)] - 1.96*sqrt(VC[(K_NB+7):(K_NB+8)]) betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+7):(K_NB+8)] betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+7):(K_NB+8)] + 1.96*sqrt(VC[(K_NB+7):(K_NB+8)]) #xi betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+9)] - 1.96*sqrt(VC[(K_NB+9)]) betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+9)] betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+9)] + 1.96*sqrt(VC[(K_NB+9)]) #pb expo(betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+10)] - 1.96*sqrt(VC[(K_NB+10)])) expo(betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+10)]) expo(betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+10)] + 1.96*sqrt(VC[(K_NB+10)])) #sb expo(betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+11)] - 1.96*sqrt(VC[(K_NB+11)])) expo(betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+11)]) expo(betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+11)] + 1.96*sqrt(VC[(K_NB+11)])) #eta betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+12):(K_NB+13)] - 1.96*sqrt(VC[(K_NB+12):(K_NB+13)]) betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+12):(K_NB+13)] betat.phic.pe.se.psilb.xic.etalb.out$pars[(K_NB+12):(K_NB+13)] + 1.96*sqrt(VC[(K_NB+12):(K_NB+13)]) betat.phic.pe.se.psilb.xic.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #18. beta(t)phi(c)p(e)s(eiis)psi(c)xi(c)eta(lb) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 12)) npar <- length(pars.start) betat.phic.pe.seiis.psic.xic.etalb.out <- fit.betat.phic.pe.seiis.psic.xic.etalb(pars.start) VC <- diag(solve(betat.phic.pe.seiis.psic.xic.etalb.out$hessian)) #phi expo(betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+2)]) expo(betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+3):(K_NB+8)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+3):(K_NB+8)] betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+3):(K_NB+8)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) #psi betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+9)] - 1.96*sqrt(VC[(K_NB+9)]) betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+9)] betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+9)] + 1.96*sqrt(VC[(K_NB+9)]) #xi betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+10)] - 1.96*sqrt(VC[(K_NB+10)]) betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+10)] betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+10)] + 1.96*sqrt(VC[(K_NB+10)]) #pb expo(betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+11)] - 1.96*sqrt(VC[(K_NB+11)])) expo(betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+11)]) expo(betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+11)] + 1.96*sqrt(VC[(K_NB+11)])) #sb expo(betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+12)] - 1.96*sqrt(VC[(K_NB+12)])) expo(betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+12)]) expo(betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+12)] + 1.96*sqrt(VC[(K_NB+12)])) #eta betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+13):(K_NB+14)] - 1.96*sqrt(VC[(K_NB+13):(K_NB+14)]) betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+13):(K_NB+14)] betat.phic.pe.seiis.psic.xic.etalb.out$pars[(K_NB+13):(K_NB+14)] + 1.96*sqrt(VC[(K_NB+13):(K_NB+14)]) betat.phic.pe.seiis.psic.xic.etalb.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #19. beta(t)phi(c)p(e)s(eiis)psi(lb)xi(c)eta(c) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 12)) npar <- length(pars.start) betat.phic.pe.seiis.psilb.xic.etac.out <- fit.betat.phic.pe.seiis.psilb.xic.etac(pars.start) VC <- diag(solve(betat.phic.pe.seiis.psilb.xic.etac.out$hessian)) #N exp(betat.phic.pe.seiis.psilb.xic.etac.out$pars[1] - 1.96*sqrt(VC[(1)])) exp(betat.phic.pe.seiis.psilb.xic.etac.out$pars[1]) exp(betat.phic.pe.seiis.psilb.xic.etac.out$pars[1] + 1.96*sqrt(VC[(1)])) exp(betat.phic.pe.seiis.psilb.xic.etac.out$pars[2] - 1.96*sqrt(VC[(2)])) exp(betat.phic.pe.seiis.psilb.xic.etac.out$pars[2]) exp(betat.phic.pe.seiis.psilb.xic.etac.out$pars[2] + 1.96*sqrt(VC[(2)])) #phi expo(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+2)]) expo(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+3):(K_NB+8)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+3):(K_NB+8)] betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+3):(K_NB+8)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) #psi betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+9):(K_NB+10)] - 1.96*sqrt(VC[(K_NB+9):(K_NB+10)]) betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+9):(K_NB+10)] betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+9):(K_NB+10)] + 1.96*sqrt(VC[(K_NB+9):(K_NB+10)]) #xi betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+11)] - 1.96*sqrt(VC[(K_NB+11)]) betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+11)] betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+11)] + 1.96*sqrt(VC[(K_NB+11)]) #pb expo(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+12)] - 1.96*sqrt(VC[(K_NB+12)])) expo(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+12)]) expo(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+12)] + 1.96*sqrt(VC[(K_NB+12)])) #sb expo(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+13)] - 1.96*sqrt(VC[(K_NB+13)])) expo(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+13)]) expo(betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+13)] + 1.96*sqrt(VC[(K_NB+13)])) #eta betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+14)] - 1.96*sqrt(VC[(K_NB+14)]) betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+14)] betat.phic.pe.seiis.psilb.xic.etac.out$pars[(K_NB+14)] + 1.96*sqrt(VC[(K_NB+14)]) betat.phic.pe.seiis.psilb.xic.etac.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #20. beta(t)phi(c)p(e)s(e)psi(lb)xi(c)eta(c) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 10)) npar <- length(pars.start) betat.phic.pe.se.psilb.xic.etac.out <- fit.betat.phic.pe.se.psilb.xic.etac(pars.start) VC <- diag(solve(betat.phic.pe.se.psilb.xic.etac.out$hessian)) #N exp(betat.phic.pe.se.psilb.xic.etac.out$pars[1] - 1.96*sqrt(VC[(1)])) exp(betat.phic.pe.se.psilb.xic.etac.out$pars[1]) exp(betat.phic.pe.se.psilb.xic.etac.out$pars[1] + 1.96*sqrt(VC[(1)])) exp(betat.phic.pe.se.psilb.xic.etac.out$pars[2] - 1.96*sqrt(VC[(2)])) exp(betat.phic.pe.se.psilb.xic.etac.out$pars[2]) exp(betat.phic.pe.se.psilb.xic.etac.out$pars[2] + 1.96*sqrt(VC[(2)])) #phi expo(betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+2)]) expo(betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+3):(K_NB+6)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+6)]) betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+3):(K_NB+6)] betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+3):(K_NB+6)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+6)]) #psi betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+7):(K_NB+8)] - 1.96*sqrt(VC[(K_NB+7):(K_NB+8)]) betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+7):(K_NB+8)] betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+7):(K_NB+8)] + 1.96*sqrt(VC[(K_NB+7):(K_NB+8)]) #xi betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+9)] - 1.96*sqrt(VC[(K_NB+9)]) betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+9)] betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+9)] + 1.96*sqrt(VC[(K_NB+9)]) #pb expo(betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+10)] - 1.96*sqrt(VC[(K_NB+10)])) expo(betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+10)]) expo(betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+10)] + 1.96*sqrt(VC[(K_NB+10)])) #sb expo(betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+11)] - 1.96*sqrt(VC[(K_NB+11)])) expo(betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+11)]) expo(betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+11)] + 1.96*sqrt(VC[(K_NB+11)])) #eta betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+12)] - 1.96*sqrt(VC[(K_NB+12)]) betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+12)] betat.phic.pe.se.psilb.xic.etac.out$pars[(K_NB+12)] + 1.96*sqrt(VC[(K_NB+12)]) betat.phic.pe.se.psilb.xic.etac.out$AIC #--------------------------------------------------------------------------------------------------------------------------------------------# #21. beta(t)phi(c)p(e)s(eiis)psi(c)xi(c)eta(c) pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(1, 0.01, 0.99)), rep(0, 11)) npar <- length(pars.start) betat.phic.pe.seiis.psic.xic.etac.out <- fit.betat.phic.pe.seiis.psic.xic.etac(pars.start) VC <- diag(solve(betat.phic.pe.seiis.psic.xic.etac.out$hessian)) #N exp(betat.phic.pe.seiis.psic.xic.etac.out$pars[1] - 1.96*sqrt(VC[(1)])) exp(betat.phic.pe.seiis.psic.xic.etac.out$pars[1]) exp(betat.phic.pe.seiis.psic.xic.etac.out$pars[1] + 1.96*sqrt(VC[(1)])) exp(betat.phic.pe.seiis.psic.xic.etac.out$pars[2] - 1.96*sqrt(VC[(2)])) exp(betat.phic.pe.seiis.psic.xic.etac.out$pars[2]) exp(betat.phic.pe.seiis.psic.xic.etac.out$pars[2] + 1.96*sqrt(VC[(2)])) #phi expo(betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+2)] - 1.96*sqrt(VC[(K_NB+2)])) expo(betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+2)]) expo(betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+2)] + 1.96*sqrt(VC[(K_NB+2)])) #p, p1, s, s1 betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+3):(K_NB+8)] - 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+3):(K_NB+8)] betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+3):(K_NB+8)] + 1.96*sqrt(VC[(K_NB+3):(K_NB+8)]) #psi betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+9)] - 1.96*sqrt(VC[(K_NB+9)]) betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+9)] betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+9)] + 1.96*sqrt(VC[(K_NB+9)]) #xi betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+10)] - 1.96*sqrt(VC[(K_NB+10)]) betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+10)] betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+10)] + 1.96*sqrt(VC[(K_NB+10)]) #pb expo(betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+11)] - 1.96*sqrt(VC[(K_NB+11)])) expo(betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+11)]) expo(betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+11)] + 1.96*sqrt(VC[(K_NB+11)])) #sb expo(betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+12)] - 1.96*sqrt(VC[(K_NB+12)])) expo(betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+12)]) expo(betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+12)] + 1.96*sqrt(VC[(K_NB+12)])) #eta betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+13)] - 1.96*sqrt(VC[(K_NB+13)]) betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+13)] betat.phic.pe.seiis.psic.xic.etac.out$pars[(K_NB+13)] + 1.96*sqrt(VC[(K_NB+13)]) betat.phic.pe.seiis.psic.xic.etac.out$AIC