diff --git a/DESCRIPTION b/DESCRIPTION index ec24d2e7..1c5bcf76 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: unmarked -Version: 1.5.0.9006 -Date: 2025-08-01 +Version: 1.5.0.9008 +Date: 2025-08-26 Type: Package Title: Models for Data from Unmarked Animals Authors@R: c( diff --git a/R/IDS.R b/R/IDS.R index 094770de..3ff994cc 100644 --- a/R/IDS.R +++ b/R/IDS.R @@ -102,20 +102,20 @@ IDS <- function(lambdaformula = ~1, } ds_hds <- get_ds_info(dataDS@dist.breaks) Xavail_ds <- matrix(0,0,0) - if(has_avail) Xavail_ds <- getDesign(dataDS, form_avail)$X + if(has_avail) Xavail_ds <- getDesign(dataDS, form_avail)$X_state if(is.null(durationDS)) durationDS <- rep(0,0) - gd_pc <- list(y=matrix(0,0,0), X=matrix(0,0,0), V=matrix(0,0,0)) + gd_pc <- list(y=matrix(0,0,0), X_state=matrix(0,0,0), X_det=matrix(0,0,0)) ds_pc <- list(total_area=0, db=c(0,0), a=0, w=0, u=0) Xavail_pc <- matrix(0,0,0) if(is.null(durationPC)) durationPC <- rep(0,0) if(!is.null(dataPC)){ gd_pc <- getDesign(dataPC, form_pc) ds_pc <- get_ds_info(c(0, maxDistPC)) - if(has_avail) Xavail_pc <- getDesign(dataPC, form_avail)$X + if(has_avail) Xavail_pc <- getDesign(dataPC, form_avail)$X_state } - gd_oc <- list(y=matrix(0,0,0), X=matrix(0,0,0), V=matrix(0,0,0)) + gd_oc <- list(y=matrix(0,0,0), X_state=matrix(0,0,0), X_det=matrix(0,0,0)) ds_oc <- list(total_area=0, db=c(0,0), a=0, w=0, u=0) Kmin_oc <- rep(0,0) Xavail_oc <- matrix(0,0,0) @@ -124,7 +124,7 @@ IDS <- function(lambdaformula = ~1, gd_oc <- getDesign(dataOC, form_oc) ds_oc <- get_ds_info(c(0, maxDistOC)) Kmin_oc <- apply(gd_oc$y, 1, max, na.rm=T) - if(has_avail) Xavail_oc <- getDesign(dataOC, form_avail)$X + if(has_avail) Xavail_oc <- getDesign(dataOC, form_avail)$X_state } # Density conversion and unequal area correction @@ -143,13 +143,13 @@ IDS <- function(lambdaformula = ~1, # Parameter stuff------------------------------------------------------------ pind_mat <- matrix(0, nrow=8, ncol=2) - pind_mat[1,] <- c(1, ncol(gd_hds$X)) - pind_mat[2,] <- max(pind_mat) + c(1, ncol(gd_hds$V)) + pind_mat[1,] <- c(1, ncol(gd_hds$X_state)) + pind_mat[2,] <- max(pind_mat) + c(1, ncol(gd_hds$X_det)) if(!is.null(detformulaPC) & !is.null(dataPC)){ - pind_mat[3,] <- max(pind_mat) + c(1, ncol(gd_pc$V)) + pind_mat[3,] <- max(pind_mat) + c(1, ncol(gd_pc$X_det)) } if(!is.null(detformulaOC) & !is.null(dataOC)){ - pind_mat[4,] <- max(pind_mat) + c(1, ncol(gd_oc$V)) + pind_mat[4,] <- max(pind_mat) + c(1, ncol(gd_oc$X_det)) } if(has_avail){ pind_mat[5,] <- max(pind_mat) + c(1, ncol(Xavail_ds)) @@ -167,8 +167,8 @@ IDS <- function(lambdaformula = ~1, if(is.null(starts)){ lam_init <- log(mean(apply(dataDS@y, 1, sum, na.rm=TRUE)) / lam_adjust[1]) - params_tmb <- list(beta_lam = c(lam_init, rep(0, ncol(gd_hds$X)-1)), - beta_hds = c(log(median(dataDS@dist.breaks)),rep(0, ncol(gd_hds$V)-1)), + params_tmb <- list(beta_lam = c(lam_init, rep(0, ncol(gd_hds$X_state)-1)), + beta_hds = c(log(median(dataDS@dist.breaks)),rep(0, ncol(gd_hds$X_det)-1)), beta_pc = rep(0,0), beta_oc = rep(0,0), beta_avail = rep(0,0), @@ -179,11 +179,11 @@ IDS <- function(lambdaformula = ~1, if(keyfun == "hazard") params_tmb$beta_schds <- 0 if(!is.null(detformulaPC) & !is.null(dataPC)){ - params_tmb$beta_pc <- c(log(maxDistPC/2), rep(0, ncol(gd_pc$V)-1)) + params_tmb$beta_pc <- c(log(maxDistPC/2), rep(0, ncol(gd_pc$X_det)-1)) if(keyfun == "hazard") params_tmb$beta_scpc <- 0 } if(!is.null(detformulaOC) & !is.null(dataOC)){ - params_tmb$beta_oc <- c(log(maxDistOC/2), rep(0, ncol(gd_oc$V)-1)) + params_tmb$beta_oc <- c(log(maxDistOC/2), rep(0, ncol(gd_oc$X_det)-1)) if(keyfun == "hazard") params_tmb$beta_scoc <- 0 } if(has_avail){ @@ -224,16 +224,16 @@ IDS <- function(lambdaformula = ~1, pind = pind_mat, lam_adjust = lam_adjust, # HDS data - y_hds = gd_hds$y, X_hds = gd_hds$X, V_hds = gd_hds$V, key_hds = keyidx, + y_hds = gd_hds$y, X_hds = gd_hds$X_state, V_hds = gd_hds$X_det, key_hds = keyidx, db_hds = dataDS@dist.breaks, a_hds = ds_hds$a, w_hds = ds_hds$w, u_hds = ds_hds$u, # PC data - y_pc = gd_pc$y, X_pc = gd_pc$X, V_pc = gd_pc$V, key_pc = keyidx, + y_pc = gd_pc$y, X_pc = gd_pc$X_state, V_pc = gd_pc$X_det, key_pc = keyidx, db_pc = c(0, maxDistPC), a_pc = ds_pc$a, w_pc = ds_pc$w, u_pc = ds_pc$u, # occ data - y_oc = gd_oc$y, X_oc = gd_oc$X, V_oc = gd_oc$V, key_oc = keyidx, + y_oc = gd_oc$y, X_oc = gd_oc$X_state, V_oc = gd_oc$X_det, key_oc = keyidx, db_oc = c(0, maxDistOC), a_oc = ds_oc$a, w_oc = ds_oc$w, u_oc = ds_oc$u, K_oc = K, Kmin_oc = Kmin_oc, @@ -251,37 +251,37 @@ IDS <- function(lambdaformula = ~1, sdr <- TMB::sdreport(tmb_obj) - lam_coef <- get_coef_info(sdr, "lam", colnames(gd_hds$X), + lam_coef <- get_coef_info(sdr, "lam", colnames(gd_hds$X_state), pind_mat[1,1]:pind_mat[1,2]) lam_est <- unmarkedEstimate(name="Density", short.name="lam", - estimates = lam_coef$ests, covMat = lam_coef$cov, fixed=1:ncol(gd_hds$X), + estimates = lam_coef$ests, covMat = lam_coef$cov, fixed=1:ncol(gd_hds$X_state), invlink = "exp", invlinkGrad = "exp") - dist_coef <- get_coef_info(sdr, "hds", colnames(gd_hds$V), + dist_coef <- get_coef_info(sdr, "hds", colnames(gd_hds$X_det), pind_mat[2,1]:pind_mat[2,2]) dist_est <- unmarkedEstimate(name="Distance sampling detection", short.name="ds", - estimates = dist_coef$ests, covMat = dist_coef$cov, fixed=1:ncol(gd_hds$V), + estimates = dist_coef$ests, covMat = dist_coef$cov, fixed=1:ncol(gd_hds$X_det), invlink = "exp", invlinkGrad = "exp") est_list <- list(lam=lam_est, ds=dist_est) if(!is.null(detformulaPC) & !is.null(dataPC)){ - pc_coef <- get_coef_info(sdr, "pc", colnames(gd_pc$V), + pc_coef <- get_coef_info(sdr, "pc", colnames(gd_pc$X_det), pind_mat[3,1]:pind_mat[3,2]) pc_est <- unmarkedEstimate(name="Point count detection", short.name="pc", - estimates = pc_coef$ests, covMat = pc_coef$cov, fixed=1:ncol(gd_pc$V), + estimates = pc_coef$ests, covMat = pc_coef$cov, fixed=1:ncol(gd_pc$X_det), invlink = "exp", invlinkGrad = "exp") est_list <- c(est_list, list(pc=pc_est)) } if(!is.null(detformulaOC) & !is.null(dataOC)){ - oc_coef <- get_coef_info(sdr, "oc", colnames(gd_oc$V), + oc_coef <- get_coef_info(sdr, "oc", colnames(gd_oc$X_det), pind_mat[4,1]:pind_mat[4,2]) oc_est <- unmarkedEstimate(name="Presence/absence detection", short.name="oc", - estimates = oc_coef$ests, covMat = oc_coef$cov, fixed=1:ncol(gd_oc$V), + estimates = oc_coef$ests, covMat = oc_coef$cov, fixed=1:ncol(gd_oc$X_det), invlink = "exp", invlinkGrad = "exp") est_list <- c(est_list, list(oc=oc_est)) } diff --git a/R/RcppExports.R b/R/RcppExports.R index 70b57ac2..a0eeeb63 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -13,36 +13,36 @@ nll_distsamp <- function(y, lam, sig, scale, a, u, w, db, keyfun, survey) { .Call(`_unmarked_nll_distsamp`, y, lam, sig, scale, a, u, w, db, keyfun, survey) } -nll_distsampOpen <- function(y, yt, Xlam, Xgam, Xom, Xsig, Xiota, beta, bi, Xlam_offset, Xgam_offset, Xom_offset, Xsig_offset, Xiota_offset, ytna, lk, mixture, first, last, first1, M, T, delta, dynamics, survey, fix, go_dims, immigration, I, I1, Ib, Ip, a, u, w, db, keyfun, lfac_k, kmyt, lfac_kmyt, fin, A) { - .Call(`_unmarked_nll_distsampOpen`, y, yt, Xlam, Xgam, Xom, Xsig, Xiota, beta, bi, Xlam_offset, Xgam_offset, Xom_offset, Xsig_offset, Xiota_offset, ytna, lk, mixture, first, last, first1, M, T, delta, dynamics, survey, fix, go_dims, immigration, I, I1, Ib, Ip, a, u, w, db, keyfun, lfac_k, kmyt, lfac_kmyt, fin, A) +nll_distsampOpen <- function(y, yt, X_lambda, X_gamma, X_omega, X_det, X_iota, beta, bi, offset_lambda, offset_gamma, offset_omega, offset_det, offset_iota, ytna, lk, mixture, first, last, first1, M, T, delta, dynamics, survey, fix, go_dims, immigration, I, I1, Ib, Ip, a, u, w, db, keyfun, lfac_k, kmyt, lfac_kmyt, fin, A) { + .Call(`_unmarked_nll_distsampOpen`, y, yt, X_lambda, X_gamma, X_omega, X_det, X_iota, beta, bi, offset_lambda, offset_gamma, offset_omega, offset_det, offset_iota, ytna, lk, mixture, first, last, first1, M, T, delta, dynamics, survey, fix, go_dims, immigration, I, I1, Ib, Ip, a, u, w, db, keyfun, lfac_k, kmyt, lfac_kmyt, fin, A) } -nll_gdistremoval <- function(beta, n_param, yDistance, yRemoval, ysum, mixture, keyfun, Xlam, A, Xphi, Xrem, Xdist, db, a, u, w, pl, K, Kmin, threads) { - .Call(`_unmarked_nll_gdistremoval`, beta, n_param, yDistance, yRemoval, ysum, mixture, keyfun, Xlam, A, Xphi, Xrem, Xdist, db, a, u, w, pl, K, Kmin, threads) +nll_gdistremoval <- function(beta, n_param, yDistance, yRemoval, ysum, mixture, keyfun, X_lambda, A, X_phi, X_rem, X_dist, db, a, u, w, pl, K, Kmin, threads) { + .Call(`_unmarked_nll_gdistremoval`, beta, n_param, yDistance, yRemoval, ysum, mixture, keyfun, X_lambda, A, X_phi, X_rem, X_dist, db, a, u, w, pl, K, Kmin, threads) } -nll_gdistsamp <- function(beta, n_param, y, mixture, keyfun, survey, Xlam, Xlam_offset, A, Xphi, Xphi_offset, Xdet, Xdet_offset, db, a, u, w, k, lfac_k, lfac_kmyt, kmyt, Kmin, threads) { - .Call(`_unmarked_nll_gdistsamp`, beta, n_param, y, mixture, keyfun, survey, Xlam, Xlam_offset, A, Xphi, Xphi_offset, Xdet, Xdet_offset, db, a, u, w, k, lfac_k, lfac_kmyt, kmyt, Kmin, threads) +nll_gdistsamp <- function(beta, n_param, y, mixture, keyfun, survey, X_lambda, offset_lambda, A, X_phi, offset_phi, X_det, offset_det, db, a, u, w, k, lfac_k, lfac_kmyt, kmyt, Kmin, threads) { + .Call(`_unmarked_nll_gdistsamp`, beta, n_param, y, mixture, keyfun, survey, X_lambda, offset_lambda, A, X_phi, offset_phi, X_det, offset_det, db, a, u, w, k, lfac_k, lfac_kmyt, kmyt, Kmin, threads) } -nll_gmultmix <- function(beta, n_param, y, mixture, pi_fun, Xlam, Xlam_offset, Xphi, Xphi_offset, Xdet, Xdet_offset, k, lfac_k, lfac_kmyt, kmyt, Kmin, threads) { - .Call(`_unmarked_nll_gmultmix`, beta, n_param, y, mixture, pi_fun, Xlam, Xlam_offset, Xphi, Xphi_offset, Xdet, Xdet_offset, k, lfac_k, lfac_kmyt, kmyt, Kmin, threads) +nll_gmultmix <- function(beta, n_param, y, mixture, pi_fun, X_lambda, offset_lambda, X_phi, offset_phi, X_det, offset_det, k, lfac_k, lfac_kmyt, kmyt, Kmin, threads) { + .Call(`_unmarked_nll_gmultmix`, beta, n_param, y, mixture, pi_fun, X_lambda, offset_lambda, X_phi, offset_phi, X_det, offset_det, k, lfac_k, lfac_kmyt, kmyt, Kmin, threads) } -nll_gpcount <- function(ym, Xlam, Xphi, Xp, beta_lam, beta_phi, beta_p, log_alpha, Xlam_offset, Xphi_offset, Xp_offset, M, mixture, T, threads) { - .Call(`_unmarked_nll_gpcount`, ym, Xlam, Xphi, Xp, beta_lam, beta_phi, beta_p, log_alpha, Xlam_offset, Xphi_offset, Xp_offset, M, mixture, T, threads) +nll_gpcount <- function(ym, X_lambda, X_phi, X_det, beta_lambda, beta_phi, beta_det, log_alpha, offset_lambda, offset_phi, offset_det, M, mixture, T, threads) { + .Call(`_unmarked_nll_gpcount`, ym, X_lambda, X_phi, X_det, beta_lambda, beta_phi, beta_det, log_alpha, offset_lambda, offset_phi, offset_det, M, mixture, T, threads) } nll_multinomPois <- function(beta, pi_fun, Xlam, Xlam_offset, Xdet, Xdet_offset, y, navec, nP, nAP) { .Call(`_unmarked_nll_multinomPois`, beta, pi_fun, Xlam, Xlam_offset, Xdet, Xdet_offset, y, navec, nP, nAP) } -nll_multmixOpen <- function(y, yt, Xlam, Xgam, Xom, Xp, Xiota, beta, bi, Xlam_offset, Xgam_offset, Xom_offset, Xp_offset, Xiota_offset, ytna, yna, lk, mixture, first, last, first1, M, T, J, R, delta, dynamics, fix, go_dims, immigration, I, I1, Ib, Ip, pi_fun, lfac_k, kmyt, lfac_kmyt, fin) { - .Call(`_unmarked_nll_multmixOpen`, y, yt, Xlam, Xgam, Xom, Xp, Xiota, beta, bi, Xlam_offset, Xgam_offset, Xom_offset, Xp_offset, Xiota_offset, ytna, yna, lk, mixture, first, last, first1, M, T, J, R, delta, dynamics, fix, go_dims, immigration, I, I1, Ib, Ip, pi_fun, lfac_k, kmyt, lfac_kmyt, fin) +nll_multmixOpen <- function(y, yt, X_lambda, X_gamma, X_omega, X_det, X_iota, beta, bi, offset_lambda, offset_gamma, offset_omega, offset_det, offset_iota, ytna, yna, lk, mixture, first, last, first1, M, T, J, R, delta, dynamics, fix, go_dims, immigration, I, I1, Ib, Ip, pi_fun, lfac_k, kmyt, lfac_kmyt, fin) { + .Call(`_unmarked_nll_multmixOpen`, y, yt, X_lambda, X_gamma, X_omega, X_det, X_iota, beta, bi, offset_lambda, offset_gamma, offset_omega, offset_det, offset_iota, ytna, yna, lk, mixture, first, last, first1, M, T, J, R, delta, dynamics, fix, go_dims, immigration, I, I1, Ib, Ip, pi_fun, lfac_k, kmyt, lfac_kmyt, fin) } -nll_nmixTTD <- function(beta, y, delta, W, V, pinds, mixture, tdist, N, J, K, naflag, threads) { - .Call(`_unmarked_nll_nmixTTD`, beta, y, delta, W, V, pinds, mixture, tdist, N, J, K, naflag, threads) +nll_nmixTTD <- function(beta, y, delta, X_state, X_det, pinds, mixture, tdist, N, J, K, naflag, threads) { + .Call(`_unmarked_nll_nmixTTD`, beta, y, delta, X_state, X_det, pinds, mixture, tdist, N, J, K, naflag, threads) } nll_occu <- function(y, X, V, beta_psi, beta_p, nd, knownOcc, navec, X_offset, V_offset, link_psi) { @@ -73,15 +73,15 @@ nll_occuRN <- function(beta, n_param, y, X, V, X_offset, V_offset, K, Kmin, thre .Call(`_unmarked_nll_occuRN`, beta, n_param, y, X, V, X_offset, V_offset, K, Kmin, threads) } -nll_occuTTD <- function(beta, y, delta, W, V, Xgam, Xeps, pind, dind, cind, eind, lpsi, tdist, N, T, J, naflag) { - .Call(`_unmarked_nll_occuTTD`, beta, y, delta, W, V, Xgam, Xeps, pind, dind, cind, eind, lpsi, tdist, N, T, J, naflag) +nll_occuTTD <- function(beta, y, delta, X_psi, X_det, X_col, X_ext, pind, dind, cind, eind, lpsi, tdist, N, T, J, naflag) { + .Call(`_unmarked_nll_occuTTD`, beta, y, delta, X_psi, X_det, X_col, X_ext, pind, dind, cind, eind, lpsi, tdist, N, T, J, naflag) } nll_pcount <- function(beta, n_param, y, X, V, X_offset, V_offset, K, Kmin, mixture, threads) { .Call(`_unmarked_nll_pcount`, beta, n_param, y, X, V, X_offset, V_offset, K, Kmin, mixture, threads) } -nll_pcountOpen <- function(ym, Xlam, Xgam, Xom, Xp, Xiota, beta_lam, beta_gam, beta_om, beta_p, beta_iota, log_alpha, Xlam_offset, Xgam_offset, Xom_offset, Xp_offset, Xiota_offset, ytna, ynam, lk, mixture, first, last, M, J, T, delta, dynamics, fix, go_dims, immigration, I, I1, Ib, Ip) { - .Call(`_unmarked_nll_pcountOpen`, ym, Xlam, Xgam, Xom, Xp, Xiota, beta_lam, beta_gam, beta_om, beta_p, beta_iota, log_alpha, Xlam_offset, Xgam_offset, Xom_offset, Xp_offset, Xiota_offset, ytna, ynam, lk, mixture, first, last, M, J, T, delta, dynamics, fix, go_dims, immigration, I, I1, Ib, Ip) +nll_pcountOpen <- function(ym, X_lambda, X_gamma, X_omega, X_det, X_iota, beta_lam, beta_gam, beta_om, beta_p, beta_iota, log_alpha, offset_lambda, offset_gamma, offset_omega, offset_det, offset_iota, ytna, ynam, lk, mixture, first, last, M, J, T, delta, dynamics, fix, go_dims, immigration, I, I1, Ib, Ip) { + .Call(`_unmarked_nll_pcountOpen`, ym, X_lambda, X_gamma, X_omega, X_det, X_iota, beta_lam, beta_gam, beta_om, beta_p, beta_iota, log_alpha, offset_lambda, offset_gamma, offset_omega, offset_det, offset_iota, ytna, ynam, lk, mixture, first, last, M, J, T, delta, dynamics, fix, go_dims, immigration, I, I1, Ib, Ip) } diff --git a/R/colext.R b/R/colext.R index e1fcc384..4560154e 100644 --- a/R/colext.R +++ b/R/colext.R @@ -8,23 +8,20 @@ colext <- function(psiformula = ~ 1, gammaformula = ~ 1, formula <- list(psiformula = psiformula, gammaformula = gammaformula, epsilonformula = epsilonformula, pformula = pformula) check_no_support(formula) - designMats <- getDesign(data, formula = as.formula(paste(unlist(formula), collapse=" "))) - X_psi <- designMats$W - X_gam <- designMats$X.gam - X_eps <- designMats$X.eps - X_det <- designMats$V - y <- designMats$y + dm <- getDesign(data, formula = as.formula(paste(unlist(formula), collapse=" "))) + X_psi <- dm$X_state; X_col <- dm$X_col; X_ext <- dm$X_ext + y <- dm$y M <- nrow(y) T <- data@numPrimary J <- ncol(y) / T psiParms <- colnames(X_psi) - gamParms <- colnames(X_gam) - epsParms <- colnames(X_eps) - detParms <- colnames(X_det) + colParms <- colnames(X_col) + extParms <- colnames(X_ext) + detParms <- colnames(dm$X_det) ## remove final year from transition prob design matrices - X_gam <- as.matrix(X_gam[-seq(T,M*T,by=T),]) - X_eps <- as.matrix(X_eps[-seq(T,M*T,by=T),]) + X_col <- as.matrix(X_col[-seq(T,M*T,by=T),]) + X_ext <- as.matrix(X_ext[-seq(T,M*T,by=T),]) # Determine which periods were sampled at each site site_sampled <- matrix(1, M, T) @@ -55,18 +52,18 @@ colext <- function(psiformula = ~ 1, gammaformula = ~ 1, # Parameter indices pind_mat <- matrix(0, 4, 2) pind_mat[1,] <- c(1, length(psiParms)) - pind_mat[2,] <- max(pind_mat[1,]) + c(1, length(gamParms)) - pind_mat[3,] <- max(pind_mat[2,]) + c(1, length(epsParms)) + pind_mat[2,] <- max(pind_mat[1,]) + c(1, length(colParms)) + pind_mat[3,] <- max(pind_mat[2,]) + c(1, length(extParms)) pind_mat[4,] <- max(pind_mat[3,]) + c(1, length(detParms)) - tmb_dat <- list(y = as.vector(t(y)), X_psi = X_psi, X_gam = X_gam, - X_eps = X_eps, X_det = X_det, + tmb_dat <- list(y = as.vector(t(y)), X_psi = X_psi, X_col = X_col, + X_ext = X_ext, X_det = dm$X_det, M = M, T = T, J = J, site_sampled = site_sampled, nd = no_detects) tmb_pars <- list(beta_psi = rep(0, length(psiParms)), - beta_gam = rep(0, length(gamParms)), - beta_eps = rep(0, length(epsParms)), + beta_col = rep(0, length(colParms)), + beta_ext = rep(0, length(extParms)), beta_det = rep(0, length(detParms))) tmb_obj <- TMB::MakeADFun(data = c(model = "tmb_colext", tmb_dat), @@ -81,8 +78,8 @@ colext <- function(psiformula = ~ 1, gammaformula = ~ 1, sdr <- TMB::sdreport(tmb_obj) psi_coef <- get_coef_info(sdr, "psi", psiParms, pind_mat[1,1]:pind_mat[1,2]) - gam_coef <- get_coef_info(sdr, "gam", gamParms, pind_mat[2,1]:pind_mat[2,2]) - eps_coef <- get_coef_info(sdr, "eps", epsParms, pind_mat[3,1]:pind_mat[3,2]) + col_coef <- get_coef_info(sdr, "col", colParms, pind_mat[2,1]:pind_mat[2,2]) + ext_coef <- get_coef_info(sdr, "ext", extParms, pind_mat[3,1]:pind_mat[3,2]) det_coef <- get_coef_info(sdr, "det", detParms, pind_mat[4,1]:pind_mat[4,2]) psi <- unmarkedEstimate(name = "Initial", short.name = "psi", @@ -92,14 +89,14 @@ colext <- function(psiformula = ~ 1, gammaformula = ~ 1, invlinkGrad = "logistic.grad") col <- unmarkedEstimate(name = "Colonization", short.name = "col", - estimates = gam_coef$ests, - covMat = gam_coef$cov, + estimates = col_coef$ests, + covMat = col_coef$cov, invlink = "logistic", invlinkGrad = "logistic.grad") ext <- unmarkedEstimate(name = "Extinction", short.name = "ext", - estimates = eps_coef$ests, - covMat = eps_coef$cov, + estimates = ext_coef$ests, + covMat = ext_coef$cov, invlink = "logistic", invlinkGrad = "logistic.grad") @@ -110,14 +107,14 @@ colext <- function(psiformula = ~ 1, gammaformula = ~ 1, invlinkGrad = "logistic.grad") estimateList <- unmarkedEstimateList(list(psi = psi, col = col, - ext = ext, det=det)) + ext = ext, det = det)) psis <- plogis(X_psi %*% psi_coef$ests) # Compute projected estimates phis <- array(NA,c(2,2,T-1,M)) - phis[,1,,] <- plogis(X_gam %x% c(-1,1) %*% gam_coef$ests) - phis[,2,,] <- plogis(X_eps %x% c(-1,1) %*% -eps_coef$ests) + phis[,1,,] <- plogis(X_col %x% c(-1,1) %*% col_coef$ests) + phis[,2,,] <- plogis(X_ext %x% c(-1,1) %*% -ext_coef$ests) projected <- array(NA, c(2, T, M)) projected[1,1,] <- 1 - psis @@ -133,9 +130,9 @@ colext <- function(psiformula = ~ 1, gammaformula = ~ 1, # Compute smoothed estimates smoothed <- calculate_smooth(y = y, psi = psis, - gam = plogis(X_gam %*% gam_coef$ests), - eps = plogis(X_eps %*% eps_coef$ests), - p = plogis(X_det %*% det_coef$ests), + col = plogis(X_col %*% col_coef$ests), + ext = plogis(X_ext %*% ext_coef$ests), + p = plogis(dm$X_det %*% det_coef$ests), M = M, T = T, J = J) smoothed.mean <- apply(smoothed, 1:2, mean) rownames(smoothed.mean) <- c("unoccupied","occupied") @@ -148,7 +145,7 @@ colext <- function(psiformula = ~ 1, gammaformula = ~ 1, gamformula = gammaformula, epsformula = epsilonformula, detformula = pformula, - data = data, sitesRemoved = designMats$removed.sites, + data = data, sitesRemoved = dm$removed.sites, estimates = estimateList, AIC = fmAIC, opt = opt, negLogLike = opt$value, nllFun = tmb_obj$fn, @@ -161,14 +158,14 @@ colext <- function(psiformula = ~ 1, gammaformula = ~ 1, # Based on Weir, Fiske, Royle 2009 "TRENDS IN ANURAN OCCUPANCY" # Appendix 1 -calculate_smooth <- function(y, psi, gam, eps, p, M, T, J){ +calculate_smooth <- function(y, psi, col, ext, p, M, T, J){ smoothed <- array(NA, c(2, T, M)) # Turn parameters into matrices p <- matrix(p, M, T*J, byrow=TRUE) - gam <- matrix(gam, M, (T-1), byrow=TRUE) - eps <- matrix(eps, M, (T-1), byrow = TRUE) + col <- matrix(col, M, (T-1), byrow=TRUE) + ext <- matrix(ext, M, (T-1), byrow = TRUE) tind <- rep(1:T, each = J) @@ -202,16 +199,16 @@ calculate_smooth <- function(y, psi, gam, eps, p, M, T, J){ psub <- p[i, tind == t] if(all(is.na(ysub))){ - alpha1[i,t] <- alpha0[i,t-1] * gam[i,t-1] + alpha1[i,t-1] * (1 - eps[i,t-1]) - alpha0[i,t] <- alpha0[i,t-1] * (1-gam[i,t-1]) + alpha1[i,t-1] * eps[i,t-1] + alpha1[i,t] <- alpha0[i,t-1] * col[i,t-1] + alpha1[i,t-1] * (1 - ext[i,t-1]) + alpha0[i,t] <- alpha0[i,t-1] * (1-col[i,t-1]) + alpha1[i,t-1] * ext[i,t-1] } else { # Case when z = 1 cp <- prod(na.omit(dbinom(ysub, 1, psub))) - alpha1[i,t] <- (alpha0[i,t-1] * gam[i,t-1] + alpha1[i,t-1] * (1 - eps[i,t-1])) * cp + alpha1[i,t] <- (alpha0[i,t-1] * col[i,t-1] + alpha1[i,t-1] * (1 - ext[i,t-1])) * cp # Case when z = 0 cp <- prod(na.omit(dbinom(ysub, 1, 0))) - alpha0[i,t] <- (alpha0[i,t-1] * (1-gam[i,t-1]) + alpha1[i,t-1] * eps[i,t-1]) * cp + alpha0[i,t] <- (alpha0[i,t-1] * (1-col[i,t-1]) + alpha1[i,t-1] * ext[i,t-1]) * cp } } @@ -228,16 +225,16 @@ calculate_smooth <- function(y, psi, gam, eps, p, M, T, J){ psub <- p[i, tind == t+1] if(all(is.na(ysub))){ - beta1[i, t] <- eps[i,t] * beta0[i, t+1] + (1-eps[i,t]) * beta1[i, t+1] - beta0[i, t] <- (1-gam[i,t]) * beta0[i, t+1] + gam[i,t] * beta1[i, t+1] + beta1[i, t] <- ext[i,t] * beta0[i, t+1] + (1-ext[i,t]) * beta1[i, t+1] + beta0[i, t] <- (1-col[i,t]) * beta0[i, t+1] + col[i,t] * beta1[i, t+1] } else { cp1 <- prod(na.omit(dbinom(ysub, 1, psub))) cp0 <- prod(na.omit(dbinom(ysub, 1, 0))) # Case when z = 1 - beta1[i, t] <- eps[i,t] * cp0 * beta0[i, t+1] + (1-eps[i,t]) * cp1 * beta1[i, t+1] + beta1[i, t] <- ext[i,t] * cp0 * beta0[i, t+1] + (1-ext[i,t]) * cp1 * beta1[i, t+1] # Case when z = 0 - beta0[i, t] <- (1-gam[i,t]) * cp0 * beta0[i, t+1] + gam[i,t] * cp1 * beta1[i, t+1] + beta0[i, t] <- (1-col[i,t]) * cp0 * beta0[i, t+1] + col[i,t] * cp1 * beta1[i, t+1] } } diff --git a/R/distsamp.R b/R/distsamp.R index 0d3bf49d..6de81329 100644 --- a/R/distsamp.R +++ b/R/distsamp.R @@ -16,14 +16,8 @@ distsamp <- function(formula, data, if(missing(starts)) starts <- NULL #Generate design matrix - designMats <- getDesign(data, formula) - X <- designMats$X; V <- designMats$V; y <- designMats$y - X.offset <- designMats$X.offset; V.offset <- designMats$V.offset - if(is.null(X.offset)) - X.offset <- rep(0, nrow(X)) - if(is.null(V.offset)) - V.offset <- rep(0, nrow(V)) - + dm <- getDesign(data, formula) + y <- dm$y M <- nrow(y) J <- ncol(y) @@ -63,8 +57,8 @@ distsamp <- function(formula, data, kmsq = A <- A) # Set up parameters - lamParms <- colnames(X) - detParms <- colnames(V) + lamParms <- colnames(dm$X_state) + detParms <- colnames(dm$X_det) scaleParms <- character(0) nAP <- length(lamParms) nDP <- length(detParms) @@ -98,8 +92,8 @@ distsamp <- function(formula, data, switch(keyfun, halfnorm = { nll <- function(param) { - sigma <- drop(exp(V %*% param[(nAP+1):nP] + V.offset)) - lambda <- drop(exp(X %*% param[1:nAP] + X.offset)) + sigma <- drop(exp(dm$X_det %*% param[(nAP+1):nP] + dm$offset_det)) + lambda <- drop(exp(dm$X_state %*% param[1:nAP] + dm$offset_state)) if(identical(output, "density")) lambda <- lambda * A for(i in 1:M) { @@ -129,8 +123,8 @@ distsamp <- function(formula, data, }}, exp = { nll <- function(param) { - rate <- drop(exp(V %*% param[(nAP+1):nP] + V.offset)) - lambda <- drop(exp(X %*% param[1:nAP] + X.offset)) + rate <- drop(exp(dm$X_det %*% param[(nAP+1):nP] + dm$offset_det)) + lambda <- drop(exp(dm$X_state %*% param[1:nAP] + dm$offset_state)) if(identical(output, "density")) lambda <- lambda * A for(i in 1:M) { @@ -165,9 +159,9 @@ distsamp <- function(formula, data, }}, hazard = { nll <- function(param) { - shape <- drop(exp(V %*% param[(nAP+1):(nP-1)] + V.offset)) + shape <- drop(exp(dm$X_det %*% param[(nAP+1):(nP-1)] + dm$offset_det)) scale <- drop(exp(param[nP])) - lambda <- drop(exp(X %*% param[1:nAP] + X.offset)) + lambda <- drop(exp(dm$X_state %*% param[1:nAP] + dm$offset_state)) if(identical(output, "density")) lambda <- lambda * A for(i in 1:M) { @@ -203,7 +197,7 @@ distsamp <- function(formula, data, }}, uniform = { nll <- function(param) { - lambda <- drop(exp(X %*% param + X.offset)) + lambda <- drop(exp(dm$X_state %*% param + dm$offset_state)) if(identical(output, "density")) lambda <- lambda * A ll <- dpois(y, lambda * u, log=TRUE) @@ -220,10 +214,10 @@ distsamp <- function(formula, data, beta.sig <- param[(nAP+1):nP] scale <- -99.0 } - lambda <- drop(exp(X %*% beta.lam + X.offset)) + lambda <- drop(exp(dm$X_state %*% beta.lam + dm$offset_state)) if(identical(output, "density")) lambda <- lambda * A - sigma <- drop(exp(V %*% beta.sig + V.offset)) + sigma <- drop(exp(dm$X_det %*% beta.sig + dm$offset_det)) nll_distsamp( y, lambda, sigma, scale, a, u, w, db, @@ -261,14 +255,14 @@ distsamp <- function(formula, data, if(output == "abund") A <- rep(1, length(A)) forms <- split_formula(formula) inps <- get_ranef_inputs(forms, list(det=siteCovs(data), state=siteCovs(data)), - list(V, X), designMats[c("Z_det","Z_state")]) + list(dm$X_det, dm$X_state), dm[c("Z_det","Z_state")]) keyfun_type <- switch(keyfun, uniform={0}, halfnorm={1}, exp={2}, hazard={3}) survey_type <- switch(survey, line={0}, point={1}) tmb_dat <- c(list(y=y, survey_type=survey_type, keyfun_type=keyfun_type, - A=A, db=db, a=a, w=w, u=u, offset_state=X.offset, - offset_det=V.offset), inps$data) + A=A, db=db, a=a, w=w, u=u, offset_state=dm$offset_state, + offset_det=dm$offset_det), inps$data) tmb_param <- c(inps$pars, list(beta_scale=rep(0,0))) @@ -324,7 +318,7 @@ distsamp <- function(formula, data, dsfit <- new("unmarkedFitDS", fitType = "distsamp", call = match.call(), opt = fm, formula = formula, data = data, keyfun=keyfun, - sitesRemoved = designMats$removed.sites, unitsOut=unitsOut, + sitesRemoved = dm$removed.sites, unitsOut=unitsOut, estimates = estimateList, AIC = fmAIC, negLogLike = fm$value, nllFun = nll, output=output, TMB=tmb_mod) return(dsfit) diff --git a/R/distsampOpen.R b/R/distsampOpen.R index 551035b7..d253d810 100644 --- a/R/distsampOpen.R +++ b/R/distsampOpen.R @@ -52,24 +52,11 @@ distsampOpen <- function(lambdaformula, gammaformula, omegaformula, pformula, D <- getDesign(data, formula) y <- D$y - Xlam <- D$Xlam - Xgam <- D$Xgam - Xom <- D$Xom - Xsig <- D$Xp - Xiota<- D$Xiota - - delta <- D$delta; go.dims <- D$go.dims - deltamax <- max(delta, na.rm=TRUE) + deltamax <- max(D$delta, na.rm=TRUE) M <- nrow(y) T <- data@numPrimary J <- ncol(getY(data)) / T - Xlam.offset <- D$Xlam.offset - Xgam.offset <- D$Xgam.offset - Xom.offset <- D$Xom.offset - Xsig.offset <- D$Xp.offset - Xiota.offset<- D$Xiota.offset - y <- array(y, c(M, J, T)) yt <- apply(y, c(1,3), function(x) { if(all(is.na(x))) @@ -124,21 +111,21 @@ distsampOpen <- function(lambdaformula, gammaformula, omegaformula, pformula, A <- rep(1, M) } - lamParms <- colnames(Xlam) - gamParms <- colnames(Xgam) - omParms <- colnames(Xom) - nAP <- ncol(Xlam) - nGP <- ncol(Xgam) - nOP <- ncol(Xom) + lamParms <- colnames(D$X_lambda) + gamParms <- colnames(D$X_gamma) + omParms <- colnames(D$X_omega) + nAP <- ncol(D$X_lambda) + nGP <- ncol(D$X_gamma) + nOP <- ncol(D$X_omega) #No parameters if uniform key function - nDP <- ifelse(keyfun == "uniform", 0, ncol(Xsig)) + nDP <- ifelse(keyfun == "uniform", 0, ncol(D$X_det)) detParms <- character(0) - if(keyfun != "uniform") detParms <- colnames(Xsig) + if(keyfun != "uniform") detParms <- colnames(D$X_det) - nIP <- ifelse(immigration, ncol(Xiota), 0) + nIP <- ifelse(immigration, ncol(D$X_iota), 0) iotaParms <- character(0) - if(immigration) iotaParms <- colnames(Xiota) + if(immigration) iotaParms <- colnames(D$X_iota) if(identical(fix, "gamma")) { if(!identical(dynamics, "constant")) @@ -215,12 +202,12 @@ distsampOpen <- function(lambdaformula, gammaformula, omegaformula, pformula, nll <- function(parms) { nll_distsampOpen( yperm, yt, - Xlam, Xgam, Xom, Xsig, Xiota, + D$X_lambda, D$X_gamma, D$X_omega, D$X_det, D$X_iota, parms, beta_ind - 1, - Xlam.offset, Xgam.offset, Xom.offset, Xsig.offset, Xiota.offset, + D$offset_lambda, D$offset_gamma, D$offset_omega, D$offset_det, D$offset_iota, ytna, lk, mixture, first - 1, last - 1, first1 - 1, M, T, - delta, dynamics, survey, fix, go.dims, immigration, + D$delta, dynamics, survey, fix, D$go.dims, immigration, I, I1, lik_trans$Ib, lik_trans$Ip, a, u, w, db, keyfun, lfac.k, kmyt, lfac.kmyt, fin, A diff --git a/R/gdistremoval.R b/R/gdistremoval.R index e26cca6a..3bb91a8d 100644 --- a/R/gdistremoval.R +++ b/R/gdistremoval.R @@ -83,19 +83,19 @@ gdistremoval <- function(lambdaformula=~1, phiformula=~1, removalformula=~1, Kmin = apply(ysum, 1, max, na.rm=T) # Parameters----------------------------------------------------------------- - n_param <- c(ncol(gd$Xlam), ifelse(mixture=="P",0,1), - ifelse(T>1,ncol(gd$Xphi),0), - ifelse(keyfun=="uniform", 0, ncol(gd$Xdist)), + n_param <- c(ncol(gd$X_lambda), ifelse(mixture=="P",0,1), + ifelse(T>1,ncol(gd$X_phi),0), + ifelse(keyfun=="uniform", 0, ncol(gd$X_dist)), ifelse(keyfun=="hazard",1,0), - ncol(gd$Xrem)) + ncol(gd$X_rem)) nP <- sum(n_param) - pnames <- colnames(gd$Xlam) + pnames <- colnames(gd$X_lambda) if(mixture!="P") pnames <- c(pnames, "alpha") - if(data@numPrimary > 1) pnames <- c(pnames, colnames(gd$Xphi)) - if(keyfun!="uniform") pnames <- c(pnames, colnames(gd$Xdist)) + if(data@numPrimary > 1) pnames <- c(pnames, colnames(gd$X_phi)) + if(keyfun!="uniform") pnames <- c(pnames, colnames(gd$X_dist)) if(keyfun=="hazard") pnames <- c(pnames, "scale") - pnames <- c(pnames, colnames(gd$Xrem)) + pnames <- c(pnames, colnames(gd$X_rem)) lam_ind <- 1:n_param[1] a_ind <- n_param[1]+1 @@ -134,7 +134,7 @@ gdistremoval <- function(lambdaformula=~1, phiformula=~1, removalformula=~1, nll <- function(param){ nll_gdistremoval(param, n_param, gd$yDist, gd$yRem, ysum, mixture_code, keyfun, - gd$Xlam, A, gd$Xphi, gd$Xrem, gd$Xdist, db, a, t(u), w, pl, + gd$X_lambda, A, gd$X_phi, gd$X_rem, gd$X_dist, db, a, t(u), w, pl, K, Kmin, threads=threads) } @@ -175,8 +175,8 @@ gdistremoval <- function(lambdaformula=~1, phiformula=~1, removalformula=~1, dlist <- list(lambda=siteCovs(data), phi=yearlySiteCovs(data), dist=siteCovs(data), rem=obsCovs(data)) inps <- get_ranef_inputs(formlist, dlist, - gd[c("Xlam","Xphi","Xdist","Xrem")], - gd[c("Zlam","Zphi","Zdist","Zrem")]) + gd[c("X_lambda","X_phi","X_dist","X_rem")], + gd[c("Z_lambda","Z_phi","Z_dist","Z_rem")]) keyfun_type <- switch(keyfun, uniform={0}, halfnorm={1}, exp={2}, hazard={3}) diff --git a/R/gdistsamp.R b/R/gdistsamp.R index d239f2f9..9f7ca562 100644 --- a/R/gdistsamp.R +++ b/R/gdistsamp.R @@ -26,19 +26,11 @@ formlist <- list(lambdaformula = lambdaformula, phiformula = phiformula, check_no_support(formlist) form <- as.formula(paste(unlist(formlist), collapse=" ")) D <- getDesign(data, formula = form) +X_lambda <- D$X_state +offset_lambda <- D$offset_state -Xlam <- D$Xlam -Xphi <- D$Xphi -Xdet <- D$Xdet y <- D$y # MxJT -Xlam.offset <- D$Xlam.offset -Xphi.offset <- D$Xphi.offset -Xdet.offset <- D$Xdet.offset -if(is.null(Xlam.offset)) Xlam.offset <- rep(0, nrow(Xlam)) -if(is.null(Xphi.offset)) Xphi.offset <- rep(0, nrow(Xphi)) -if(is.null(Xdet.offset)) Xdet.offset <- rep(0, nrow(Xdet)) - M <- nrow(y) T <- data@numPrimary R <- ncol(y) @@ -90,22 +82,22 @@ switch(unitsOut, kmsq = A <- A) -lamPars <- colnames(Xlam) +lamPars <- colnames(X_lambda) if(T==1) { phiPars <- character(0) nPP <- 0 } else { - phiPars <- colnames(Xphi) - nPP <- ncol(Xphi) + phiPars <- colnames(D$X_phi) + nPP <- ncol(D$X_phi) } if(identical(keyfun, "uniform")) { nDP <- 0 detPars <- character(0) } else { - nDP <- ncol(Xdet) - detPars <- colnames(Xdet) + nDP <- ncol(D$X_det) + detPars <- colnames(D$X_det) } if(identical(keyfun, "hazard")) { nSP <- 1 @@ -126,7 +118,7 @@ if(identical(mixture, "NB")) { nbPar <- character(0) } -nLP <- ncol(Xlam) +nLP <- ncol(X_lambda) nP <- nLP + nPP + nDP + nSP + nOP cp <- array(as.numeric(NA), c(M, T, J+1)) @@ -150,24 +142,24 @@ for(i in 1:M) { switch(keyfun, halfnorm = { - altdetParms <- paste("sigma", colnames(Xdet), sep="") + altdetParms <- paste("sigma", colnames(D$X_det), sep="") if(missing(starts)) { starts <- rep(0, nP) starts[nLP+nPP+1] <- log(max(db)) } nll_R <- function(pars) { - lambda <- exp(Xlam %*% pars[1:nLP] + Xlam.offset) + lambda <- exp(X_lambda %*% pars[1:nLP] + offset_lambda) if(identical(output, "density")) lambda <- lambda * A if(T==1) phi <- matrix(1, M, T) else { - phi <- plogis(Xphi %*% pars[(nLP+1):(nLP+nPP)] + Xphi.offset) + phi <- plogis(D$X_phi %*% pars[(nLP+1):(nLP+nPP)] + D$offset_phi) phi <- matrix(phi, M, T, byrow=TRUE) } - sigma <- exp(Xdet %*% pars[(nLP+nPP+1):(nLP+nPP+nDP)]+Xdet.offset) + sigma <- exp(D$X_det %*% pars[(nLP+nPP+1):(nLP+nPP+nDP)]+D$offset_det) sigma <- matrix(sigma, M, T, byrow=TRUE) switch(mixture, @@ -220,23 +212,23 @@ halfnorm = { } }, exp = { - altdetParms <- paste("rate", colnames(Xdet), sep="") + altdetParms <- paste("rate", colnames(D$X_det), sep="") if(missing(starts)) { starts <- rep(0, nP) starts[nLP+nPP+1] <- log(max(db)) } nll_R <- function(pars) { - lambda <- exp(Xlam %*% pars[1:nLP] + Xlam.offset) + lambda <- exp(X_lambda %*% pars[1:nLP] + offset_lambda) if(identical(output, "density")) lambda <- lambda * A if(T==1) phi <- matrix(1, M, T) else { - phi <- plogis(Xphi %*% pars[(nLP+1):(nLP+nPP)] + Xphi.offset) + phi <- plogis(D$X_phi %*% pars[(nLP+1):(nLP+nPP)] + D$offset_phi) phi <- matrix(phi, M, T, byrow=TRUE) } - rate <- exp(Xdet %*% pars[(nLP+nPP+1):(nLP+nPP+nDP)] + Xdet.offset) + rate <- exp(D$X_det %*% pars[(nLP+nPP+1):(nLP+nPP+nDP)] + D$offset_det) rate <- matrix(rate, M, T, byrow=TRUE) switch(mixture, @@ -290,22 +282,22 @@ exp = { } }, hazard = { - altdetParms <- paste("shape", colnames(Xdet), sep="") + altdetParms <- paste("shape", colnames(D$X_det), sep="") if(missing(starts)) { starts <- rep(0, nP) } nll_R <- function(pars) { - lambda <- exp(Xlam %*% pars[1:nLP] + Xlam.offset) + lambda <- exp(X_lambda %*% pars[1:nLP] + offset_lambda) if(identical(output, "density")) lambda <- lambda * A if(T==1) phi <- matrix(1, M, T) else { - phi <- plogis(Xphi %*% pars[(nLP+1):(nLP+nPP)] + Xphi.offset) + phi <- plogis(D$X_phi %*% pars[(nLP+1):(nLP+nPP)] + D$offset_phi) phi <- matrix(phi, M, T, byrow=TRUE) } - shape <- exp(Xdet %*% pars[(nLP+nPP+1):(nLP+nPP+nDP)]+Xdet.offset) + shape <- exp(D$X_det %*% pars[(nLP+nPP+1):(nLP+nPP+nDP)]+D$offset_det) shape <- matrix(shape, M, T, byrow=TRUE) scale <- exp(pars[nLP+nPP+nDP+1]) @@ -364,13 +356,13 @@ uniform = { starts <- rep(0, nP) } nll_R <- function(pars) { - lambda <- exp(Xlam %*% pars[1:nLP] + Xlam.offset) + lambda <- exp(X_lambda %*% pars[1:nLP] + offset_lambda) if(identical(output, "density")) lambda <- lambda * A if(T==1) phi <- matrix(1, M, T) else { - phi <- plogis(Xphi %*% pars[(nLP+1):(nLP+nPP)] + Xphi.offset) + phi <- plogis(D$X_phi %*% pars[(nLP+1):(nLP+nPP)] + D$offset_phi) phi <- matrix(phi, M, T, byrow=TRUE) } p <- 1 @@ -418,7 +410,7 @@ if(engine =="C"){ nll <- function(params){ nll_gdistsamp(params, n_param, y_long, mixture_code, keyfun, survey, - Xlam, Xlam.offset, A, Xphi, Xphi.offset, Xdet, Xdet.offset, + X_lambda, offset_lambda, A, D$X_phi, D$offset_phi, D$X_det, D$offset_det, db, a, t(u), w, k, lfac.k, lfac.kmytC, kmytC, Kmin, threads) } diff --git a/R/getDesign.R b/R/getDesign.R index 2252c485..0d875312 100644 --- a/R/getDesign.R +++ b/R/getDesign.R @@ -61,24 +61,24 @@ setMethod("getDesign", "unmarkedFrame", } list(y = y, - X = X_state, X.offset = offset_state, Z_state = Z_state, - V = X_det, V.offset = offset_det, Z_det = Z_det, + X_state = X_state, offset_state = offset_state, Z_state = Z_state, + X_det = X_det, offset_det = offset_det, Z_det = Z_det, removed.sites = which(drop_sites)) }) # UnmarkedMultFrame -# used by colext and base class for G3 and DailMadsen +# used by colext, occuTTD, nmixTTD and base class for G3 and DailMadsen setMethod("getDesign", "unmarkedMultFrame", function(umf, formula, na.rm = TRUE){ - psiformula <- as.formula(formula[[2]][[2]][[2]]) + stateformula <- as.formula(formula[[2]][[2]][[2]]) gamformula <- as.formula(as.call(list(as.name("~"), formula[[2]][[2]][[3]]))) epsformula <- as.formula(as.call(list(as.name("~"), formula[[2]][[3]]))) detformula <- as.formula(as.call(list(as.name("~"), formula[[3]]))) # Process state and detection with generic umf method - comb_form <- list(as.name("~"), detformula, psiformula[[2]]) + comb_form <- list(as.name("~"), detformula, stateformula[[2]]) comb_form <- as.formula(as.call(comb_form)) out <- methods::callNextMethod(umf, formula = comb_form, na.rm = FALSE) @@ -123,16 +123,16 @@ setMethod("getDesign", "unmarkedMultFrame", warning("Site(s) ", paste(which(drop_sites), collapse = ","), " dropped due to missing values", call.=FALSE) y <- y[!drop_sites,,drop=FALSE] - out$X <- out$X[!drop_sites,,drop=FALSE] + out$X_state <- out$X_state[!drop_sites,,drop=FALSE] drop_sites_per <- rep(drop_sites, each = T) X_col <- X_col[!drop_sites_per,,drop=FALSE] X_ext <- X_ext[!drop_sites_per,,drop=FALSE] drop_sites_obs <- rep(drop_sites, each = R) - out$V <- out$V[!drop_sites_obs,,drop=FALSE] + out$X_det <- out$X_det[!drop_sites_obs,,drop=FALSE] } # Combine outputs - list(y = y, W = out$X, X.gam = X_col, X.eps = X_ext, V = out$V, + list(y = y, X_state = out$X_state, X_col = X_col, X_ext = X_ext, X_det = out$X_det, removed.sites = which(drop_sites)) }) @@ -141,12 +141,12 @@ setMethod("getDesign", "unmarkedMultFrame", setMethod("getDesign", "unmarkedFrameG3", function(umf, formula, na.rm = TRUE){ - lamformula <- as.formula(formula[[2]][[2]]) + stateformula <- as.formula(formula[[2]][[2]]) phiformula <- as.formula(as.call(list(as.name("~"), formula[[2]][[3]]))) detformula <- as.formula(as.call(list(as.name("~"), formula[[3]]))) # Process state and detection with generic umf method - comb_form <- list(as.name("~"), detformula, lamformula[[2]]) + comb_form <- list(as.name("~"), detformula, stateformula[[2]]) comb_form <- as.formula(as.call(comb_form)) # Have to use getMethod because this inherits from unmarkedMultFrame getDesign_generic <- methods::getMethod("getDesign", "unmarkedFrame") @@ -179,20 +179,20 @@ setMethod("getDesign", "unmarkedFrameG3", warning("Site(s) ", paste(which(drop_sites), collapse = ","), " dropped due to missing values", call.=FALSE) y <- y[!drop_sites,,drop=FALSE] - out$X <- out$X[!drop_sites,,drop=FALSE] - out$X.offset <- out$X.offset[!drop_sites] + out$X_state <- out$X_state[!drop_sites,,drop=FALSE] + out$offset_state <- out$offset_state[!drop_sites] drop_sites_per <- rep(drop_sites, each = T) X_phi <- X_phi[!drop_sites_per,,drop=FALSE] offset_phi <- offset_phi[!drop_sites_per] drop_sites_obs <- rep(drop_sites, each = R) - out$V <- out$V[!drop_sites_obs,,drop=FALSE] - out$V.offset <- out$V.offset[!drop_sites_obs] + out$X_det <- out$X_det[!drop_sites_obs,,drop=FALSE] + out$offset_det <- out$offset_det[!drop_sites_obs] } # Combine outputs - test = list(y = y, Xlam = out$X, Xlam.offset = out$X.offset, - Xphi = X_phi, Xphi.offset = offset_phi, - Xdet = out$V, Xdet.offset = out$V.offset, + list(y = y, X_state = out$X_state, offset_state = out$offset_state, + X_phi = X_phi, offset_phi = offset_phi, + X_det = out$X_det, offset_det = out$offset_det, removed.sites = which(drop_sites)) }) @@ -359,11 +359,11 @@ setMethod("getDesign", "unmarkedFrameDailMadsen", } } - list(y = y, Xlam = X_lambda, Xlam.offset = offset_lambda, - Xgam = X_gamma, Xgam.offset = offset_gamma, - Xom = X_omega, Xom.offset = offset_omega, - Xiota = X_iota, Xiota.offset = offset_iota, - Xp = X_det, Xp.offset = offset_det, + list(y = y, X_lambda = X_lambda, offset_lambda = offset_lambda, + X_gamma = X_gamma, offset_gamma = offset_gamma, + X_omega = X_omega, offset_omega = offset_omega, + X_iota = X_iota, offset_iota = offset_iota, + X_det = X_det, offset_det = offset_det, delta = delta, removed.sites = which(drop_sites), go.dims = go.dims) }) @@ -427,56 +427,56 @@ setMethod("getDesign", "unmarkedFrameGDR", if(return.frames) return(list(sc=sc, ysc=ysc, oc=oc)) lam_fixed <- reformulas::nobars(formula$lambdaformula) - Xlam <- model.matrix(lam_fixed, + X_lambda <- model.matrix(lam_fixed, model.frame(lam_fixed, sc, na.action=NULL)) phi_fixed <- reformulas::nobars(formula$phiformula) - Xphi <- model.matrix(phi_fixed, + X_phi<- model.matrix(phi_fixed, model.frame(phi_fixed, ysc, na.action=NULL)) dist_fixed <- reformulas::nobars(formula$distanceformula) - Xdist <- model.matrix(dist_fixed, + X_dist <- model.matrix(dist_fixed, model.frame(dist_fixed, ysc, na.action=NULL)) rem_fixed <- reformulas::nobars(formula$removalformula) - Xrem <- model.matrix(rem_fixed, + X_rem <- model.matrix(rem_fixed, model.frame(rem_fixed, oc, na.action=NULL)) - Zlam <- get_Z(formula$lambdaformula, sc) - Zphi <- get_Z(formula$phiformula, ysc) - Zdist <- get_Z(formula$distanceformula, ysc) - Zrem <- get_Z(formula$removalformula, oc) + Z_lambda <- get_Z(formula$lambdaformula, sc) + Z_phi <- get_Z(formula$phiformula, ysc) + Z_dist <- get_Z(formula$distanceformula, ysc) + Z_rem <- get_Z(formula$removalformula, oc) # Check if there are missing yearlySiteCovs ydist_mat <- apply(matrix(yDist, nrow=M*T, byrow=TRUE), 1, function(x) any(is.na(x))) yrem_mat <- apply(matrix(yRem, nrow=M*T, byrow=TRUE), 1, function(x) any(is.na(x))) ok_missing_phi_covs <- ydist_mat | yrem_mat - missing_phi_covs <- apply(Xphi, 1, function(x) any(is.na(x))) + missing_phi_covs <- apply(X_phi, 1, function(x) any(is.na(x))) if(!all(which(missing_phi_covs) %in% which(ok_missing_phi_covs))){ stop("Missing yearlySiteCovs values for some observations that are not missing", call.=FALSE) } # Check if there are missing dist covs - missing_dist_covs <- apply(Xdist, 1, function(x) any(is.na(x))) + missing_dist_covs <- apply(X_dist, 1, function(x) any(is.na(x))) ok_missing_dist_covs <- ydist_mat if(!all(which(missing_dist_covs) %in% which(ok_missing_dist_covs))){ stop("Missing yearlySiteCovs values for some distance observations that are not missing", call.=FALSE) } # Check if there are missing rem covs - missing_obs_covs <- apply(Xrem, 1, function(x) any(is.na(x))) + missing_obs_covs <- apply(X_rem, 1, function(x) any(is.na(x))) missing_obs_covs <- apply(matrix(missing_obs_covs, nrow=M*T, byrow=TRUE), 1, function(x) any(x)) ok_missing_obs_covs <- yrem_mat if(!all(which(missing_obs_covs) %in% which(ok_missing_obs_covs))){ stop("Missing obsCovs values for some removal observations that are not missing", call.=FALSE) } - if(any(is.na(Xlam))){ + if(any(is.na(X_lambda))){ stop("gdistremoval does not currently handle missing values in siteCovs", call.=FALSE) } - list(yDist=yDist, yRem=yRem, Xlam=Xlam, Xphi=Xphi, Xdist=Xdist, Xrem=Xrem, - Zlam=Zlam, Zphi=Zphi, Zdist=Zdist, Zrem=Zrem) + list(yDist=yDist, yRem=yRem, X_lambda=X_lambda, X_phi=X_phi, X_dist=X_dist, X_rem=X_rem, + Z_lambda=Z_lambda, Z_phi=Z_phi, Z_dist=Z_dist, Z_rem=Z_rem) }) @@ -521,11 +521,11 @@ setMethod("getDesign", "unmarkedFrameOccuFP", " dropped due to missing values", call.=FALSE) y <- y[!drop_sites,,drop=FALSE] - out$X <- out$X[!drop_sites,,drop=FALSE] - out$X.offset <- out$X.offset[!drop_sites] + out$X_state <- out$X_state[!drop_sites,,drop=FALSE] + out$offset_state <- out$offset_state[!drop_sites] drop_sites_obs <- rep(drop_sites, each = J) - out$V <- out$V[!drop_sites_obs,,drop=FALSE] - out$V.offset <- out$V.offset[!drop_sites_obs] + out$X_det <- out$X_det[!drop_sites_obs,,drop=FALSE] + out$offset_det <- out$offset_det[!drop_sites_obs] X_fp <- X_fp[!drop_sites_obs,,drop=FALSE] offset_fp <- offset_fp[!drop_sites_obs] @@ -533,8 +533,8 @@ setMethod("getDesign", "unmarkedFrameOccuFP", offset_b <- offset_b[!drop_sites_obs] } - # Combine outputs - new_out <- list(U = X_fp, U.offset = offset_fp, W = X_b, W.offset = offset_b) + # Combine outputs (U=fp, W=b) + new_out <- list(X_fp = X_fp, offset_fp = offset_fp, X_b = X_b, offset_b = offset_b) out <- c(out, new_out) out$y <- y out$removed.sites <- which(drop_sites) diff --git a/R/gmultmix.R b/R/gmultmix.R index a809a4cd..29406d7e 100644 --- a/R/gmultmix.R +++ b/R/gmultmix.R @@ -16,19 +16,10 @@ formlist <- list(lambdaformula = lambdaformula, phiformula = phiformula, check_no_support(formlist) form <- as.formula(paste(unlist(formlist), collapse=" ")) D <- getDesign(data, formula = form) - -Xlam <- D$Xlam -Xphi <- D$Xphi -Xdet <- D$Xdet +X_lambda <- D$X_state +offset_lambda <- D$offset_state y <- D$y # MxJT -Xlam.offset <- D$Xlam.offset -Xphi.offset <- D$Xphi.offset -Xdet.offset <- D$Xdet.offset -if(is.null(Xlam.offset)) Xlam.offset <- rep(0, nrow(Xlam)) -if(is.null(Xphi.offset)) Xphi.offset <- rep(0, nrow(Xphi)) -if(is.null(Xdet.offset)) Xdet.offset <- rep(0, nrow(Xdet)) - K <- check_K_multinomial(K, K_adjust = 100, y, data@numPrimary) k <- 0:K lk <- length(k) @@ -48,17 +39,17 @@ yt <- apply(y, 1:2, function(x) { piFun <- data@piFun -lamPars <- colnames(Xlam) -detPars <- colnames(Xdet) -nLP <- ncol(Xlam) +lamPars <- colnames(X_lambda) +detPars <- colnames(D$X_det) +nLP <- ncol(X_lambda) if(T==1) { nPP <- 0 phiPars <- character(0) } else if(T>1) { - nPP <- ncol(Xphi) - phiPars <- colnames(Xphi) + nPP <- ncol(D$X_phi) + phiPars <- colnames(D$X_phi) } -nDP <- ncol(Xdet) +nDP <- ncol(D$X_det) nP <- nLP + nPP + nDP + (mixture%in%c('NB','ZIP')) if(!missing(starts) && length(starts) != nP) stop(paste("The number of starting values should be", nP)) @@ -81,12 +72,12 @@ for(i in 1:M) { } nll_R <- function(pars) { - lambda <- exp(Xlam %*% pars[1:nLP] + Xlam.offset) + lambda <- exp(X_lambda %*% pars[1:nLP] + offset_lambda) if(T==1) phi <- 1 else if(T>1) - phi <- drop(plogis(Xphi %*% pars[(nLP+1):(nLP+nPP)] + Xphi.offset)) - p <- plogis(Xdet %*% pars[(nLP+nPP+1):(nLP+nPP+nDP)] + Xdet.offset) + phi <- drop(plogis(D$X_phi %*% pars[(nLP+1):(nLP+nPP)] + D$offset_phi)) + p <- plogis(D$X_det %*% pars[(nLP+nPP+1):(nLP+nPP+nDP)] + D$offset_det) phi.mat <- matrix(phi, M, T, byrow=TRUE) phi <- as.numeric(phi.mat) @@ -143,8 +134,8 @@ if(engine=="R"){ Kmin <- apply(yt, 1, max, na.rm=TRUE) nll <- function(params) { - nll_gmultmix(params, n_param, y_long, mixture_code, piFun, Xlam, Xlam.offset, - Xphi, Xphi.offset, Xdet, Xdet.offset, k, lfac.k, lfac.kmytC, + nll_gmultmix(params, n_param, y_long, mixture_code, piFun, X_lambda, offset_lambda, + D$X_phi, D$offset_phi, D$X_det, D$offset_det, k, lfac.k, lfac.kmytC, kmytC, Kmin, threads) } diff --git a/R/goccu.R b/R/goccu.R index 9db8eead..75c922de 100644 --- a/R/goccu.R +++ b/R/goccu.R @@ -31,11 +31,8 @@ goccu <- function(psiformula, phiformula, pformula, data, # handle offsets gd <- getDesign(data, formula = formula) - + X_psi <- gd$X_state y <- gd$y - Xpsi <- gd$Xlam # The G3 getDesign method identifies state as lam - Xphi <- gd$Xphi - Xp <- gd$Xdet M <- nrow(y) T <- data@numPrimary @@ -81,7 +78,7 @@ goccu <- function(psiformula, phiformula, pformula, data, # Bundle data for TMB dataList <- list(y=y, T=T, link=ifelse(linkPsi=='cloglog', 1, 0), - Xpsi=Xpsi, Xphi=Xphi, Xp=Xp, + X_psi=X_psi, X_phi=gd$X_phi, X_det=gd$X_det, n_possible=n_possible, alpha_potential=alpha_potential, alpha_drop = alpha_drop, known_present=known_present, known_available=known_available, @@ -89,13 +86,13 @@ goccu <- function(psiformula, phiformula, pformula, data, # Provide dimensions and starting values for parameters # This part should change to be more like occu() if we add random effects - psi_ind <- 1:ncol(Xpsi) - phi_ind <- 1:ncol(Xphi) + max(psi_ind) - p_ind <- 1:ncol(Xp) + max(phi_ind) - nP <- max(p_ind) + psi_ind <- 1:ncol(X_psi) + phi_ind <- 1:ncol(gd$X_phi) + max(psi_ind) + det_ind <- 1:ncol(gd$X_det) + max(phi_ind) + nP <- max(det_ind) params <- list(beta_psi = rep(0, length(psi_ind)), beta_phi = rep(0, length(phi_ind)), - beta_p = rep(0, length(p_ind))) + beta_det = rep(0, length(det_ind))) # Create TMB object tmb_mod <- TMB::MakeADFun(data = c(model = "tmb_goccu", dataList), @@ -109,14 +106,14 @@ goccu <- function(psiformula, phiformula, pformula, data, covMat <- invertHessian(opt, nP, se) ests <- opt$par - names(ests) <- c(colnames(Xpsi), colnames(Xphi), colnames(Xp)) + names(ests) <- c(colnames(X_psi), colnames(gd$X_phi), colnames(gd$X_det)) fmAIC <- 2 * opt$value + 2 * nP psi_est <- unmarkedEstimate(name = "Occupancy", short.name = "psi", estimates = ests[psi_ind], covMat = covMat[psi_ind, psi_ind, drop=FALSE], - fixed = 1:ncol(Xpsi), + fixed = 1:ncol(X_psi), invlink = psiInvLink, invlinkGrad = psiLinkGrad, randomVarInfo=list() @@ -125,23 +122,22 @@ goccu <- function(psiformula, phiformula, pformula, data, phi_est <- unmarkedEstimate(name = "Availability", short.name = "phi", estimates = ests[phi_ind], covMat = covMat[phi_ind, phi_ind, drop=FALSE], - fixed = 1:ncol(Xphi), + fixed = 1:ncol(gd$X_phi), invlink = "logistic", invlinkGrad = "logistic.grad", randomVarInfo=list() ) - p_est <- unmarkedEstimate(name = "Detection", short.name = "p", - estimates = ests[p_ind], - covMat = covMat[p_ind, p_ind, drop=FALSE], - fixed = 1:ncol(Xp), + det_est <- unmarkedEstimate(name = "Detection", short.name = "p", + estimates = ests[det_ind], + covMat = covMat[det_ind, det_ind, drop=FALSE], + fixed = 1:ncol(gd$X_det), invlink = "logistic", invlinkGrad = "logistic.grad", randomVarInfo=list() ) - estimate_list <- unmarkedEstimateList(list(psi=psi_est, phi=phi_est, - det=p_est)) + estimate_list <- unmarkedEstimateList(list(psi=psi_est, phi=phi_est, det=det_est)) # Create unmarkedFit object-------------------------------------------------- umfit <- new("unmarkedFitGOccu", fitType = "goccu", call = match.call(), @@ -160,8 +156,7 @@ goccu <- function(psiformula, phiformula, pformula, data, setMethod("predict_inputs_from_umf", "unmarkedFitGOccu", function(object, type, newdata, na.rm, re.form=NA){ designMats <- getDesign(newdata, object@formula, na.rm=na.rm) - X_idx <- switch(type, psi="Xlam", phi="Xphi", det="Xdet") # Xlam = state design matrix - off_idx <- paste0(X_idx, ".offset") + X_idx <- switch(type, psi="X_state", phi="X_phi", det="X_det") list(X=designMats[[X_idx]], offset=NULL) }) @@ -190,12 +185,12 @@ setMethod("fitted_internal", "unmarkedFitGOccu", function(object){ M <- numSites(object@data) JT <- obsNum(object@data) gd <- getDesign(object@data, object@formula, na.rm=FALSE) - Xpsi <- gd$Xlam + X_psi <- gd$X_state - psi <- drop(plogis(Xpsi %*% coef(object, "psi"))) + psi <- drop(plogis(X_psi %*% coef(object, "psi"))) psi <- matrix(psi, nrow=M, ncol=JT) - phi <- drop(plogis(gd$Xphi %*% coef(object, "phi"))) + phi <- drop(plogis(gd$X_phi %*% coef(object, "phi"))) phi <- rep(phi, each = JT / object@data@numPrimary) phi <- matrix(phi, nrow=M, ncol=JT, byrow=TRUE) @@ -215,10 +210,10 @@ setMethod("ranef_internal", "unmarkedFitGOccu", function(object, ...){ gd <- getDesign(object@data, object@formula, na.rm=FALSE) y_array <- array(t(gd$y), c(J, T, M)) - Xpsi <- gd$Xlam + X_psi <- gd$X_state - psi <- drop(plogis(Xpsi %*% coef(object, "psi"))) - phi <- drop(plogis(gd$Xphi %*% coef(object, "phi"))) + psi <- drop(plogis(X_psi %*% coef(object, "psi"))) + phi <- drop(plogis(gd$X_phi %*% coef(object, "phi"))) phi <- matrix(phi, nrow=M, ncol=T, byrow=TRUE) p <- getP(object) p_array <- array(t(p), c(J, T, M)) diff --git a/R/gpcount.R b/R/gpcount.R index af4dade2..d0b97483 100644 --- a/R/gpcount.R +++ b/R/gpcount.R @@ -14,19 +14,10 @@ formlist <- list(lambdaformula = lambdaformula, phiformula = phiformula, check_no_support(formlist) form <- as.formula(paste(unlist(formlist), collapse=" ")) D <- getDesign(data, formula = form) - -Xlam <- D$Xlam -Xphi <- D$Xphi -Xdet <- D$Xdet +X_lambda <- D$X_state +offset_lambda <- D$offset_state ym <- D$y # MxJT -Xlam.offset <- D$Xlam.offset -Xphi.offset <- D$Xphi.offset -Xdet.offset <- D$Xdet.offset -if(is.null(Xlam.offset)) Xlam.offset <- rep(0, nrow(Xlam)) -if(is.null(Xphi.offset)) Xphi.offset <- rep(0, nrow(Xphi)) -if(is.null(Xdet.offset)) Xdet.offset <- rep(0, nrow(Xdet)) - if(missing(K) || is.null(K)) { K <- max(ym, na.rm=TRUE) + 100 warning("K was not specified, so was set to max(y)+100 =", K) @@ -41,12 +32,12 @@ J <- numY(data) / T y <- array(ym, c(I, J, T)) -lamPars <- colnames(Xlam) -detPars <- colnames(Xdet) -nLP <- ncol(Xlam) -nPP <- ncol(Xphi) -phiPars <- colnames(Xphi) -nDP <- ncol(Xdet) +lamPars <- colnames(X_lambda) +detPars <- colnames(D$X_det) +nLP <- ncol(X_lambda) +nPP <- ncol(D$X_phi) +phiPars <- colnames(D$X_phi) +nDP <- ncol(D$X_det) nP <- nLP + nPP + nDP + (mixture%in%c('NB','ZIP')) if(!missing(starts) && length(starts) != nP) stop("There should be", nP, "starting values, not", length(starts)) @@ -54,10 +45,10 @@ if(!missing(starts) && length(starts) != nP) if(identical(engine, "R")) { # Minus negative log-likelihood nll <- function(pars) { - lam <- exp(Xlam %*% pars[1:nLP] + Xlam.offset) - phi <- plogis(Xphi %*% pars[(nLP+1):(nLP+nPP)] + Xphi.offset) + lam <- exp(X_lambda %*% pars[1:nLP] + offset_lambda) + phi <- plogis(D$X_phi %*% pars[(nLP+1):(nLP+nPP)] + D$offset_phi) phi <- matrix(phi, I, T, byrow=TRUE) - p <- plogis(Xdet %*% pars[(nLP+nPP+1):(nLP+nPP+nDP)] + Xdet.offset) + p <- plogis(D$X_det %*% pars[(nLP+nPP+1):(nLP+nPP+nDP)] + D$offset_det) p <- matrix(p, I, byrow=TRUE) p <- array(p, c(I, J, T)) # byrow? L <- rep(NA, I) @@ -99,14 +90,14 @@ nll <- function(pars) { } else if(identical(engine, "C")) { nll <- function(pars) { - beta.lam <- pars[1:nLP] - beta.phi <- pars[(nLP+1):(nLP+nPP)] - beta.p <- pars[(nLP+nPP+1):(nLP+nPP+nDP)] + beta_lambda <- pars[1:nLP] + beta_phi <- pars[(nLP+1):(nLP+nPP)] + beta_det <- pars[(nLP+nPP+1):(nLP+nPP+nDP)] log.alpha <- 1 if(mixture %in% c("NB", "ZIP")) log.alpha <- pars[nP] - nll_gpcount(ym, Xlam, Xphi, Xdet, beta.lam, beta.phi, beta.p, - log.alpha, Xlam.offset, Xphi.offset, Xdet.offset, + nll_gpcount(ym, X_lambda, D$X_phi, D$X_det, beta_lambda, beta_phi, beta_det, + log.alpha, offset_lambda, D$offset_phi, D$offset_det, as.integer(K), mixture, T, threads) } } diff --git a/R/multinomPois.R b/R/multinomPois.R index 076410f3..7c144df3 100644 --- a/R/multinomPois.R +++ b/R/multinomPois.R @@ -8,24 +8,18 @@ multinomPois <- function(formula, data, starts, method = "BFGS", stop("Data is not a data frame or unmarkedFrame.") engine <- match.arg(engine, c("C", "R", "TMB")) if(any(sapply(split_formula(formula), has_random))) engine <- "TMB" - designMats <- getDesign(data, formula) - X <- designMats$X; V <- designMats$V; y <- designMats$y - X.offset <- designMats$X.offset; V.offset <- designMats$V.offset - if (is.null(X.offset)) { - X.offset <- rep(0, nrow(X)) - } - if (is.null(V.offset)) { - V.offset <- rep(0, nrow(V)) - } + dm <- getDesign(data, formula) + y <- dm$y + J <- ncol(y) R <- obsNum(data) M <- nrow(y) piFun <- data@piFun - lamParms <- colnames(X) - detParms <- colnames(V) - nDP <- ncol(V) - nAP <- ncol(X) + lamParms <- colnames(dm$X_state) + detParms <- colnames(dm$X_det) + nDP <- ncol(dm$X_det) + nAP <- ncol(dm$X_state) lamIdx <- 1:nAP pIdx <- (nAP+1):(nAP+nDP) nP <- nDP + nAP @@ -36,8 +30,8 @@ multinomPois <- function(formula, data, starts, method = "BFGS", navec <- is.na(yvec) nll_R <- function(parms) { - lambda <- exp(X %*% parms[1 : nAP] + X.offset) - p <- plogis(V %*% parms[(nAP + 1) : nP] + V.offset) + lambda <- exp(dm$X_state %*% parms[1 : nAP] + dm$offset_state) + p <- plogis(dm$X_det %*% parms[(nAP + 1) : nP] + dm$offset_det) p.matrix <- matrix(p, M, R, byrow = TRUE) pi <- do.call(piFun, list(p = p.matrix)) logLikeSite <- dpois(y, matrix(lambda, M, J) * pi, log = TRUE) @@ -48,7 +42,7 @@ multinomPois <- function(formula, data, starts, method = "BFGS", nll_C <- function(params) { nll_multinomPois( params,piFun, - X, X.offset, V, V.offset, + dm$X_state, dm$offset_state, dm$X_det, dm$offset_det, yC, navecC, nP,nAP ) } @@ -86,15 +80,15 @@ multinomPois <- function(formula, data, starts, method = "BFGS", forms <- split_formula(formula) obs_all <- add_covariates(obsCovs(data), siteCovs(data), numSites(data)*obsNum(data)) inps <- get_ranef_inputs(forms, list(det=obs_all, state=siteCovs(data)), - list(V, X), designMats[c("Z_det","Z_state")]) + list(dm$X_det, dm$X_state), dm[c("Z_det","Z_state")]) if(!piFun%in%c('doublePiFun','removalPiFun','depDoublePiFun')){ stop("Custom pi functions are not supported by TMB engine.") } pifun_type <- switch(piFun, removalPiFun={0}, doublePiFun={1}, depDoublePiFun={2}) - tmb_dat <- c(list(y=y, pifun_type=pifun_type, offset_state=X.offset, - offset_det=V.offset), inps$data) + tmb_dat <- c(list(y=y, pifun_type=pifun_type, offset_state=dm$offset_state, + offset_det=dm$offset_det), inps$data) # Fit model in TMB if(missing(starts)) starts <- NULL @@ -130,7 +124,7 @@ multinomPois <- function(formula, data, starts, method = "BFGS", umfit <- new("unmarkedFitMPois", fitType = "multinomPois", call = match.call(), formula = formula, data = data, - estimates = estimateList, sitesRemoved = designMats$removed.sites, + estimates = estimateList, sitesRemoved = dm$removed.sites, AIC = fmAIC, opt = fm, negLogLike = fm$value, nllFun = nll, TMB=tmb_mod) return(umfit) diff --git a/R/multmixOpen.R b/R/multmixOpen.R index 56bec1a8..fca1448c 100644 --- a/R/multmixOpen.R +++ b/R/multmixOpen.R @@ -52,17 +52,6 @@ multmixOpen <- function(lambdaformula, gammaformula, omegaformula, pformula, last <- apply(!ytna, 1, function(x) max(which(x))) first1 <- which(first==1)[1] - Xlam.offset <- D$Xlam.offset - Xgam.offset <- D$Xgam.offset - Xom.offset <- D$Xom.offset - Xp.offset <- D$Xp.offset - Xiota.offset <- D$Xiota.offset - if(is.null(Xlam.offset)) Xlam.offset <- rep(0, M) - if(is.null(Xgam.offset)) Xgam.offset <- rep(0, M*(T-1)) - if(is.null(Xom.offset)) Xom.offset <- rep(0, M*(T-1)) - if(is.null(Xp.offset)) Xp.offset <- rep(0, M*T*R) - if(is.null(Xiota.offset)) Xiota.offset <- rep(0, M*(T-1)) - #K stuff K <- check_K_multinomial(K, K_adjust = 20, D$y, T) k <- 0:K @@ -82,18 +71,18 @@ multmixOpen <- function(lambdaformula, gammaformula, omegaformula, pformula, } } - lamParms <- colnames(D$Xlam) - gamParms <- colnames(D$Xgam) - omParms <- colnames(D$Xom) - detParms <- colnames(D$Xp) - nAP <- ncol(D$Xlam) - nGP <- ncol(D$Xgam) - nOP <- ncol(D$Xom) - nDP <- ncol(D$Xp) + lamParms <- colnames(D$X_lambda) + gamParms <- colnames(D$X_gamma) + omParms <- colnames(D$X_omega) + detParms <- colnames(D$X_det) + nAP <- ncol(D$X_lambda) + nGP <- ncol(D$X_gamma) + nOP <- ncol(D$X_omega) + nDP <- ncol(D$X_det) - nIP <- ifelse(immigration, ncol(D$Xiota), 0) + nIP <- ifelse(immigration, ncol(D$X_iota), 0) iotaParms <- character(0) - if(immigration) iotaParms <- colnames(D$Xiota) + if(immigration) iotaParms <- colnames(D$X_iota) if(identical(fix, "gamma")) { if(!identical(dynamics, "constant")) @@ -168,9 +157,9 @@ multmixOpen <- function(lambdaformula, gammaformula, omegaformula, pformula, nll <- function(parms) { nll_multmixOpen( yperm, yt, - D$Xlam, D$Xgam, D$Xom, D$Xp, D$Xiota, + D$X_lambda, D$X_gamma, D$X_omega, D$X_det, D$X_iota, parms, beta_ind - 1, - Xlam.offset, Xgam.offset, Xom.offset, Xp.offset, Xiota.offset, + D$offset_lambda, D$offset_gamma, D$offset_omega, D$offset_det, D$offset_iota, ytna, yna, lk, mixture, first - 1, last - 1, first1 - 1, M, T, J, R, D$delta, dynamics, fix, D$go.dims, immigration, diff --git a/R/nmixTTD.R b/R/nmixTTD.R index 6d2fcf7e..2192e8fc 100644 --- a/R/nmixTTD.R +++ b/R/nmixTTD.R @@ -24,11 +24,9 @@ nmixTTD <- function(stateformula=~1, detformula=~1, data, K=100, formula <- as.formula(paste(unlist(formula),collapse=" ")) #Process input data---------------------------------------------------------- - designMats <- getDesign(data, formula) - #V = detection; W = abundance - V <- designMats$V; W <- designMats$W - y <- designMats$y - removed <- designMats$removed.sites + dm <- getDesign(data, formula) + y <- dm$y + removed <- dm$removed.sites N <- nrow(y) R <- ncol(y) @@ -44,8 +42,8 @@ nmixTTD <- function(stateformula=~1, detformula=~1, data, K=100, delta <- as.numeric(yvec0) knownOccLog <- knownOccLog[-removed] + if(length(dm$removed.sites)>0) knownOccLog <- knownOccLog[-dm$removed.sites] # Set up parameter names and indices----------------------------------------- - occParms <- colnames(X) - detParms <- colnames(V) - nDP <- ncol(V) - nOP <- ncol(X) + occParms <- colnames(dm$X_state) + detParms <- colnames(dm$X_det) + nDP <- ncol(dm$X_det) + nOP <- ncol(dm$X_state) nP <- nDP + nOP psiIdx <- 1:nOP pIdx <- (nOP+1):nP @@ -53,8 +49,8 @@ occu <- function(formula, data, knownOcc = numeric(0), beta.psi <- params[1:nOP] beta.p <- params[(nOP+1):nP] nll_occu( - yvec, X, V, beta.psi, beta.p, nd, knownOccLog, navec, - X.offset, V.offset, linkPsi + yvec, dm$X_state, dm$X_det, beta.psi, beta.p, nd, knownOccLog, navec, + dm$offset_state, dm$offset_det, linkPsi ) } } else if (identical(engine, "R")){ @@ -63,9 +59,9 @@ occu <- function(formula, data, knownOcc = numeric(0), M <- nrow(y) nll <- function(params) { - psi <- psiLinkFunc(X %*% params[1 : nOP] + X.offset) + psi <- psiLinkFunc(dm$X_state %*% params[1 : nOP] + dm$offset_state) psi[knownOccLog] <- 1 - pvec <- plogis(V %*% params[(nOP + 1) : nP] + V.offset) + pvec <- plogis(dm$X_det %*% params[(nOP + 1) : nP] + dm$offset_det) cp <- (pvec^yvec) * ((1 - pvec)^(1 - yvec)) cp[navec] <- 1 # so that NA's don't modify likelihood cpmat <- matrix(cp, M, J, byrow = TRUE) # @@ -102,10 +98,10 @@ occu <- function(formula, data, knownOcc = numeric(0), forms <- split_formula(formula) obs_all <- add_covariates(obsCovs(data), siteCovs(data), length(getY(data))) inps <- get_ranef_inputs(forms, list(det=obs_all, state=siteCovs(data)), - list(V, X), designMats[c("Z_det","Z_state")]) + list(dm$X_det, dm$X_state), dm[c("Z_det","Z_state")]) tmb_dat <- c(list(y=y, no_detect=nd, link=ifelse(linkPsi=="cloglog",1,0), - offset_state=X.offset, offset_det=V.offset), inps$data) + offset_state=dm$offset_state, offset_det=dm$offset_det), inps$data) # Fit model with TMB if(missing(starts)) starts <- NULL @@ -152,7 +148,7 @@ occu <- function(formula, data, knownOcc = numeric(0), # Create unmarkedFit object-------------------------------------------------- umfit <- new("unmarkedFitOccu", fitType = "occu", call = match.call(), formula = formula, data = data, - sitesRemoved = designMats$removed.sites, + sitesRemoved = dm$removed.sites, estimates = estimateList, AIC = fmAIC, opt = fm, negLogLike = fm$value, nllFun = nll, knownOcc = knownOccLog, TMB=tmb_mod) diff --git a/R/occuCOP.R b/R/occuCOP.R index 288ad631..0b223dbb 100644 --- a/R/occuCOP.R +++ b/R/occuCOP.R @@ -96,8 +96,8 @@ setMethod("predict_inputs_from_umf", "unmarkedFitOccuCOP", designMats = getDesign(umf = newdata, formula = comb_form, na.rm = na.rm) - if (type == "psi") list_els <- c("X", "Z_state") - if (type == "lambda") list_els <- c("V", "Z_det") + if (type == "psi") list_els <- c("X_state", "Z_state") + if (type == "lambda") list_els <- c("X_det", "Z_det") X <- designMats[[list_els[1]]] if (is.null(re.form)) X <- cbind(X, designMats[[list_els[2]]]) return(list(X = X, offset = NULL)) @@ -139,10 +139,10 @@ setMethod("fitted_internal", "unmarkedFitOccuCOP", function(object) { comb_form <- list(as.name("~"), fl$lambdaformula, fl$psiformula[[2]]) des <- getDesign(data, comb_form, na.rm = FALSE) estim_psi = as.numeric(do.call(object["psi"]@invlink, - list(as.matrix(des$X %*% coef(object, 'psi'))))) + list(as.matrix(des$X_state %*% coef(object, 'psi'))))) estim_lambda = do.call(object["lambda"]@invlink, list(matrix( - as.numeric(des$V %*% coef(object, 'lambda')), + as.numeric(des$X_det %*% coef(object, 'lambda')), nrow = M, ncol = J, byrow = T))) return(estim_psi * estim_lambda) }) @@ -461,10 +461,10 @@ occuCOP <- function(data, } # Xpsi is the fixed effects design matrix for occupancy - Xpsi <- designMats$X + Xpsi <- designMats$X_state # Xlambda is the fixed effects design matrix for detection rate - Xlambda <- designMats$V + Xlambda <- designMats$X_det # Zpsi is the random effects design matrix for occupancy Zpsi <- designMats$Z_state diff --git a/R/occuFP.R b/R/occuFP.R index 14711b5a..da9f83e4 100644 --- a/R/occuFP.R +++ b/R/occuFP.R @@ -13,44 +13,25 @@ occuFP <- function(detformula = ~ 1,FPformula = ~ 1,Bformula = ~ 1,stateformula if(!is(data, "unmarkedFrameOccuFP")) stop("Data is not an unmarkedFrameOccuFP object.") type <- data@type - if(sum(type[2:3])==0) stop("Only type 1 data. No data types with false positives. Use occu instead.") - designMats <- getDesign(data, detformula,FPformula,Bformula,stateformula) - X <- designMats$X; V <- designMats$V; U <- designMats$U; W <- designMats$W; - y <- designMats$y - + dm <- getDesign(data, detformula,FPformula,Bformula,stateformula) + y <- dm$y if(any(type[1:2]>0)) if(any(y[,1:sum(type[1:2])]>1,na.rm = TRUE)) stop("Values of y for type 1 and type 2 data must be 0 or 1.") if(type[3]>0) if(any(y[,1:sum(type[3])]>2,na.rm = TRUE)) stop("Values of y for type 3 data must be 0, 1, or 2.") - removed <- designMats$removed.sites - X.offset <- designMats$X.offset; V.offset <- designMats$V.offset; U.offset <- designMats$U.offset; W.offset <- designMats$W.offset - if(is.null(X.offset)) { - X.offset <- rep(0, nrow(X)) - } - if(is.null(V.offset)) { - V.offset <- rep(0, nrow(V)) - } - if(is.null(U.offset)) { - U.offset <- rep(0, nrow(U)) - } - if(is.null(W.offset)) { - W.offset <- rep(0, nrow(W)) - } - J <- ncol(y) M <- nrow(y) - - occParms <- colnames(X) - detParms <- colnames(V) - FPParms <- colnames(U) - if(type[3]!=0) BParms <- colnames(W) else BParms = NULL - nDP <- ncol(V) - nFP <- ncol(U) - nBP <- ifelse(type[3]!=0,ncol(W),0) - nOP <- ncol(X) + occParms <- colnames(dm$X_state) + detParms <- colnames(dm$X_det) + FPParms <- colnames(dm$X_fp) + if(type[3]!=0) BParms <- colnames(dm$X_b) else BParms = NULL + nDP <- ncol(dm$X_det) + nFP <- ncol(dm$X_fp) + nBP <- ifelse(type[3]!=0,ncol(dm$X_b),0) + nOP <- ncol(dm$X_state) nP <- nDP + nOP + nFP + nBP if(!missing(starts) && length(starts) != nP) @@ -63,12 +44,12 @@ occuFP <- function(detformula = ~ 1,FPformula = ~ 1,Bformula = ~ 1,stateformula nll <- function(params) { - psi <- plogis(X %*% params[1 : nOP] + X.offset) - pvec <- plogis(V %*% params[(nOP + 1) : (nOP + nDP)] + V.offset) - fvec <- plogis(U %*% params[(nOP + nDP + 1) : (nOP + nDP + nFP)] + U.offset) + psi <- plogis(dm$X_state %*% params[1 : nOP] + dm$offset_state) + pvec <- plogis(dm$X_det %*% params[(nOP + 1) : (nOP + nDP)] + dm$offset_det) + fvec <- plogis(dm$X_fp %*% params[(nOP + nDP + 1) : (nOP + nDP + nFP)] + dm$offset_fp) if (type[1]!=0) fvec[rep(c(rep(TRUE,type[1]),rep(FALSE,sum(type[2:3]))),M)] = 0 if (type[3]!=0){ - bvec <- plogis(W %*% params[(nOP + nDP + nFP + 1) : nP] + W.offset) + bvec <- plogis(dm$X_b %*% params[(nOP + nDP + nFP + 1) : nP] + dm$offset_b) if (type[1]!=0|type[2]!=0) bvec[rep(c(rep(TRUE,sum(type[1:2])),rep(FALSE,type[3])),M)] = 0} if (type[3]==0){ bvec <- matrix(0,M*J,1) @@ -130,7 +111,7 @@ occuFP <- function(detformula = ~ 1,FPformula = ~ 1,Bformula = ~ 1,stateformula umfit <- new("unmarkedFitOccuFP", fitType = "occuFP", call = match.call(), detformula = detformula,FPformula = FPformula,Bformula = Bformula, stateformula = stateformula, formula = ~1, type = type, data = data, - sitesRemoved = designMats$removed.sites, + sitesRemoved = dm$removed.sites, estimates = estimateList, AIC = fmAIC, opt = fm, negLogLike = fm$value, nllFun = nll) diff --git a/R/occuPEN.R b/R/occuPEN.R index 55c75637..ce2b5202 100644 --- a/R/occuPEN.R +++ b/R/occuPEN.R @@ -3,7 +3,7 @@ computeMPLElambda = function(formula, data, knownOcc = numeric(0), starts, metho engine <- match.arg(engine, c("C", "R")) designMats <- getDesign(data, formula) - X <- designMats$X; V <- designMats$V; y <- designMats$y + X_state <- designMats$X_state; X_det <- designMats$X_det; y <- designMats$y removed <- designMats$removed.sites y <- truncateToBinary(y) @@ -13,18 +13,18 @@ computeMPLElambda = function(formula, data, knownOcc = numeric(0), starts, metho if(length(removed)>0) knownOccLog <- knownOccLog[-removed] - nDP <- ncol(V) - nOP <- ncol(X) + nDP <- ncol(X_det) + nOP <- ncol(X_state) nP <- nDP + nOP if(!missing(starts) && length(starts) != nP) stop(paste("The number of starting values should be", nP)) if(missing(starts)) starts <- rep(0, nP) - LRparams = glm.fit(x=X,y=apply(y,1,max),family=binomial(),intercept=F,start=starts[1:nOP]) + LRparams = glm.fit(x=X_state,y=apply(y,1,max),family=binomial(),intercept=F,start=starts[1:nOP]) naiveOcc = mean(LRparams$fitted.values) occuOutMLE = occu(formula,data,knownOcc = knownOcc, starts = starts, method = "BFGS", engine = engine, se = TRUE) - meanDet = mean((1+exp(-occuOutMLE[2]@estimates%*%t(V)))^-1) + meanDet = mean((1+exp(-occuOutMLE[2]@estimates%*%t(X_det)))^-1) MPLElambda = sqrt(sum(diag(occuOutMLE[2]@covMat)))*(1-(1-meanDet)^(dim(y)[2]))*(1-naiveOcc) # what if there are different numbers of visits to different sites? return(MPLElambda) } @@ -49,7 +49,7 @@ occuPEN_CV <- function(formula, data, knownOcc = numeric(0), starts, engine <- match.arg(engine, c("C", "R")) designMats <- getDesign(data, formula) - X <- designMats$X; V <- designMats$V; y <- designMats$y + X_state <- designMats$X_state; X_det <- designMats$X_det; y <- designMats$y y <- truncateToBinary(y) J <- ncol(y) M <- nrow(y) @@ -78,15 +78,15 @@ occuPEN_CV <- function(formula, data, knownOcc = numeric(0), starts, beta.psi <- params[1:nOP] beta.p <- params[(nOP+1):nP] nll_occu( - yvec, X, V, beta.psi, beta.p, nd, knownOccLog, navec, - X.offset, V.offset, "logit" + yvec, X_state, X_det, beta.psi, beta.p, nd, knownOccLog, navec, + offset_state, offset_det, "logit" ) } } else { nll <- function(params) { # penalize this function - psi <- plogis(X %*% params[1 : nOP] + X.offset) + psi <- plogis(X_state %*% params[1 : nOP] + offset_state) psi[knownOccLog] <- 1 - pvec <- plogis(V %*% params[(nOP + 1) : nP] + V.offset) + pvec <- plogis(X_det %*% params[(nOP + 1) : nP] + offset_det) cp <- (pvec^yvec) * ((1 - pvec)^(1 - yvec)) cp[navec] <- 1 # so that NA's don't modify likelihood cpmat <- matrix(cp, M, J, byrow = TRUE) # @@ -103,15 +103,9 @@ occuPEN_CV <- function(formula, data, knownOcc = numeric(0), starts, occuTest = data[which(foldAssignments==fold),] # test on this fold designMats <- getDesign(occuTest, formula) - X <- designMats$X; V <- designMats$V; y <- designMats$y + X_state <- designMats$X_state; X_det <- designMats$X_det; y <- designMats$y removed <- designMats$removed.sites - X.offset <- designMats$X.offset; V.offset <- designMats$V.offset - if(is.null(X.offset)) { - X.offset <- rep(0, nrow(X)) - } - if(is.null(V.offset)) { - V.offset <- rep(0, nrow(V)) - } + offset_state <- designMats$offset_state; offset_det <- designMats$offset_det y <- truncateToBinary(y) J <- ncol(y) @@ -123,10 +117,10 @@ occuPEN_CV <- function(formula, data, knownOcc = numeric(0), starts, if(length(removed)>0) knownOccLog <- knownOccLog[-removed] - occParms <- colnames(X) - detParms <- colnames(V) - nDP <- ncol(V) - nOP <- ncol(X) + occParms <- colnames(X_state) + detParms <- colnames(X_det) + nDP <- ncol(X_det) + nOP <- ncol(X_state) nP <- nDP + nOP if(!missing(starts) && length(starts) != nP) @@ -186,20 +180,14 @@ occuPEN <- function(formula, data, knownOcc = numeric(0), starts, engine <- match.arg(engine, c("C", "R")) designMats <- getDesign(data, formula) - X <- designMats$X; V <- designMats$V; y <- designMats$y + X_state <- designMats$X_state; X_det <- designMats$X_det; y <- designMats$y - if (ncol(X)==1 & pen.type=="MPLE") stop("MPLE requires occupancy covariates.") + if (ncol(X_state)==1 & pen.type=="MPLE") stop("MPLE requires occupancy covariates.") - if (ncol(X)==1 & ncol(V)==1 & pen.type=="Ridge") stop("Ridge requires covariates.") + if (ncol(X_state)==1 & ncol(X_det)==1 & pen.type=="Ridge") stop("Ridge requires covariates.") removed <- designMats$removed.sites - X.offset <- designMats$X.offset; V.offset <- designMats$V.offset - if(is.null(X.offset)) { - X.offset <- rep(0, nrow(X)) - } - if(is.null(V.offset)) { - V.offset <- rep(0, nrow(V)) - } + offset_state <- designMats$offset_state; offset_det <- designMats$offset_det y <- truncateToBinary(y) J <- ncol(y) @@ -211,10 +199,10 @@ occuPEN <- function(formula, data, knownOcc = numeric(0), starts, if(length(removed)>0) knownOccLog <- knownOccLog[-removed] - occParms <- colnames(X) - detParms <- colnames(V) - nDP <- ncol(V) - nOP <- ncol(X) + occParms <- colnames(X_state) + detParms <- colnames(X_det) + nDP <- ncol(X_det) + nOP <- ncol(X_state) nP <- nDP + nOP if(!missing(starts) && length(starts) != nP) @@ -226,11 +214,11 @@ occuPEN <- function(formula, data, knownOcc = numeric(0), starts, nd <- ifelse(rowSums(y,na.rm=TRUE) == 0, 1, 0) # no det at site i ## need to add offsets !!!!!!!!!!!!!! - ## and fix bug causing crash when NAs are in V + ## and fix bug causing crash when NAs are in X_det ## compute logistic regression MPLE targets and lambda: if (pen.type=="MPLE") { - LRparams = glm.fit(x=X,y=apply(y,1,max),family=binomial(),intercept=F,start=starts[1:nOP]) + LRparams = glm.fit(x=X_state,y=apply(y,1,max),family=binomial(),intercept=F,start=starts[1:nOP]) MPLElambda = computeMPLElambda(formula, data, knownOcc = numeric(0), starts, method = "BFGS", engine = engine) if (MPLElambda != lambda) warning("Supplied lambda does not match the computed value. Proceeding with the supplied lambda.") } @@ -252,14 +240,14 @@ occuPEN <- function(formula, data, knownOcc = numeric(0), starts, } else { stop("pen.type not found") } - nll_occuPEN(yvec, X, V, beta.psi, beta.p, nd, knownOccLog, navec, - X.offset, V.offset, penalty) + nll_occuPEN(yvec, X_state, X_det, beta.psi, beta.p, nd, knownOccLog, navec, + offset_state, offset_det, penalty) } } else { nll <- function(params) { # penalize this function - psi <- plogis(X %*% params[1 : nOP] + X.offset) + psi <- plogis(X_state %*% params[1 : nOP] + offset_state) psi[knownOccLog] <- 1 - pvec <- plogis(V %*% params[(nOP + 1) : nP] + V.offset) + pvec <- plogis(X_det %*% params[(nOP + 1) : nP] + offset_det) cp <- (pvec^yvec) * ((1 - pvec)^(1 - yvec)) cp[navec] <- 1 # so that NA's don't modify likelihood cpmat <- matrix(cp, M, J, byrow = TRUE) # diff --git a/R/occuRN.R b/R/occuRN.R index 62b09249..71f2649d 100644 --- a/R/occuRN.R +++ b/R/occuRN.R @@ -4,29 +4,23 @@ occuRN <- function(formula, data, K = 25, starts, method = "BFGS", se = TRUE, engine=c("C","R"), threads = 1, ...) { - check_no_support(split_formula(formula)) + check_no_support(split_formula(formula)) - if(!is(data, "unmarkedFrameOccu")) - stop("Data is not an unmarkedFrameOccu object.") - - engine <- match.arg(engine, c("C", "R")) - designMats <- getDesign(data, formula) - X <- designMats$X; V <- designMats$V; y <- designMats$y - X.offset <- designMats$X.offset; V.offset <- designMats$V.offset - if(is.null(X.offset)) - X.offset <- rep(0, nrow(X)) - if (is.null(V.offset)) - V.offset <- rep(0, nrow(V)) + if(!is(data, "unmarkedFrameOccu")) + stop("Data is not an unmarkedFrameOccu object.") + engine <- match.arg(engine, c("C", "R")) + dm <- getDesign(data, formula) + y <- dm$y y <- truncateToBinary(y) J <- ncol(y) M <- nrow(y) - occParms <- colnames(X) - detParms <- colnames(V) - nDP <- ncol(V) - nOP <- ncol(X) + occParms <- colnames(dm$X_state) + detParms <- colnames(dm$X_det) + nDP <- ncol(dm$X_det) + nOP <- ncol(dm$X_state) nP <- nDP + nOP if(!missing(starts) && length(starts) != nP) @@ -40,7 +34,7 @@ occuRN <- function(formula, data, K = 25, starts, method = "BFGS", { ## compute individual level detection probabilities - r.ij <- matrix(plogis(V %*% parms[(nOP + 1) : nP] + V.offset), M, J, + r.ij <- matrix(plogis(dm$X_det %*% parms[(nOP + 1) : nP] + dm$offset_det), M, J, byrow = TRUE) ## compute list of detection probabilities along N @@ -59,7 +53,7 @@ occuRN <- function(formula, data, K = 25, starts, method = "BFGS", cp.in <- sapply(cp.ij.list, rowProds) ## compute P(N = n | lambda_i) along i - lambda.i <- exp(X %*% parms[1 : nOP] + X.offset) + lambda.i <- exp(dm$X_state %*% parms[1 : nOP] + dm$offset_state) lambda.in <- sapply(n, function(x) dpois(x, lambda.i)) ## integrate over P(y_i | N = n) * P(N = n | lambda_i) wrt n @@ -74,7 +68,7 @@ occuRN <- function(formula, data, K = 25, starts, method = "BFGS", n_param <- c(nOP, nDP) Kmin <- apply(y, 1, function(x) max(x, na.rm=TRUE)) nll <- function(params){ - nll_occuRN(params, n_param, y, X, V, X.offset, V.offset, + nll_occuRN(params, n_param, y, dm$X_state, dm$X_det, dm$offset_state, dm$offset_det, K, Kmin, threads) } } @@ -102,7 +96,7 @@ occuRN <- function(formula, data, K = 25, starts, method = "BFGS", umfit <- new("unmarkedFitOccuRN", fitType = "occuRN", call = match.call(), formula = formula, data = data, - sitesRemoved = designMats$removed.sites, estimates = estimateList, + sitesRemoved = dm$removed.sites, estimates = estimateList, AIC = fmAIC, opt = fm, negLogLike = fm$value, nllFun = nll, K = K) return(umfit) diff --git a/R/occuTTD.R b/R/occuTTD.R index 57514506..fde1b346 100644 --- a/R/occuTTD.R +++ b/R/occuTTD.R @@ -31,12 +31,10 @@ occuTTD <- function(psiformula=~1, gammaformula=~1, epsilonformula=~1, } #Process input data---------------------------------------------------------- - designMats <- getDesign(data, formula) - #V = detection; W = occupancy - V <- designMats$V; W <- designMats$W - X.gam <- designMats$X.gam; X.eps <- designMats$X.eps - y <- designMats$y - removed <- designMats$removed.sites + dm <- getDesign(data, formula) + X_psi <- dm$X_state + X_col <- dm$X_col; X_ext <- dm$X_ext + y <- dm$y N <- nrow(y) R <- ncol(y) @@ -47,13 +45,13 @@ occuTTD <- function(psiformula=~1, gammaformula=~1, epsilonformula=~1, yvec <- as.numeric(t(y)) naflag <- as.numeric(is.na(yvec)) surveyLength <- data@surveyLength - if(length(removed>0)) surveyLength <- surveyLength[-removed,] + if(length(dm$removed.sites>0)) surveyLength <- surveyLength[-dm$removed.sites,] ymax <- as.numeric(t(surveyLength)) delta <- as.numeric(yvec1){ #Remove final year from col/ext design matrices - X.gam <- as.matrix(X.gam[-seq(T,N*T,by=T),,drop=FALSE]) - X.eps <- as.matrix(X.eps[-seq(T,N*T,by=T),,drop=FALSE]) + X_col <- as.matrix(X_col[-seq(T,N*T,by=T),,drop=FALSE]) + X_ext <- as.matrix(X_ext[-seq(T,N*T,by=T),,drop=FALSE]) - gamParms <- colnames(X.gam); nGP <- ncol(X.gam) - epsParms <- colnames(X.eps); nEP <- ncol(X.eps) + gamParms <- colnames(X_col); nGP <- ncol(X_col) + epsParms <- colnames(X_ext); nEP <- ncol(X_ext) col_inds <- (nOP+1):(nOP+nGP) ext_inds <- (nOP+nGP+1):(nOP+nGP+nEP) } @@ -80,9 +78,9 @@ occuTTD <- function(psiformula=~1, gammaformula=~1, epsilonformula=~1, nll_R <- function(params){ #Get occupancy and detection parameters - psi <- linkFunc(W %*% params[psi_inds]) + psi <- linkFunc(X_psi %*% params[psi_inds]) psi <- cbind(1-psi, psi) - lam <- exp(V %*% params[det_inds]) + lam <- exp(dm$X_det %*% params[det_inds]) #Simplified version of Garrard et al. 2013 eqn 5 #Extended to Weibull @@ -101,8 +99,8 @@ occuTTD <- function(psiformula=~1, gammaformula=~1, epsilonformula=~1, #If dynamic, get col/ext probs and transition prob matrix if(T>1){ - col <- plogis(X.gam %*% params[col_inds]) - ext <- plogis(X.eps %*% params[ext_inds]) + col <- plogis(X_col %*% params[col_inds]) + ext <- plogis(X_ext %*% params[ext_inds]) phi <- cbind(1-col, col, ext, 1-ext) } @@ -136,7 +134,7 @@ occuTTD <- function(psiformula=~1, gammaformula=~1, epsilonformula=~1, nll_C <- function(params){ nll_occuTTD( - params, yvec, delta, W, V, X.gam, X.eps, + params, yvec, delta, X_psi, dm$X_det, X_col, X_ext, range(psi_inds)-1, range(det_inds)-1, range(col_inds)-1, range(ext_inds)-1, linkPsi, ttdDist, N, T, J, naflag @@ -206,7 +204,7 @@ occuTTD <- function(psiformula=~1, gammaformula=~1, epsilonformula=~1, gamformula = gammaformula, epsformula = epsilonformula, detformula = detformula, - data = data, sitesRemoved = removed, + data = data, sitesRemoved = dm$removed.sites, estimates = estimateList, AIC = fmAIC, opt = fm, negLogLike = fm$value, nllFun = nll) diff --git a/R/pcount.R b/R/pcount.R index bbce47e0..955903b1 100644 --- a/R/pcount.R +++ b/R/pcount.R @@ -19,16 +19,14 @@ pcount <- function(formula, data, K, mixture = c("P", "NB", "ZIP"), starts, stop("ZIP mixture not available for R engine") # Generate design matrices------------------------------------------------- - designMats <- getDesign(data, formula) - X <- designMats$X; V <- designMats$V; y <- designMats$y - X.offset <- designMats$X.offset; V.offset <- designMats$V.offset - Z_state <- designMats$Z_state; Z_det <- designMats$Z_det + dm <- getDesign(data, formula) + y <- dm$y # Set up parameter names and indices--------------------------------------- - lamParms <- colnames(X) - detParms <- colnames(V) - nDP <- ncol(V) - nAP <- ncol(X) + lamParms <- colnames(dm$X_state) + detParms <- colnames(dm$X_det) + nDP <- ncol(dm$X_det) + nAP <- ncol(dm$X_state) lamIdx <- 1:nAP pIdx <- (nAP+1):(nAP+nDP) n_param <- c(nAP, nDP, ifelse(mixture != "P", 1, 0)) @@ -59,8 +57,8 @@ pcount <- function(formula, data, K, mixture = c("P", "NB", "ZIP"), starts, ijk <- expand.grid(k = 0:K, j = 1:J, i = 1:M) ijk.to.ikj <- with(ijk, order(i, k, j)) nll <- function(parms) { - theta.i <- exp(X %*% parms[lamIdx] + X.offset) - p.ij <- plogis(V %*% parms[pIdx] + V.offset) + theta.i <- exp(dm$X_state %*% parms[lamIdx] + dm$offset_state) + p.ij <- plogis(dm$X_det %*% parms[pIdx] + dm$offset_det) theta.ik <- rep(theta.i, each = K + 1) p.ijk <- rep(p.ij, each = K + 1) @@ -83,8 +81,8 @@ pcount <- function(formula, data, K, mixture = c("P", "NB", "ZIP"), starts, } } else if(identical(engine, "C")) { nll <- function(parms) { - nll_pcount(parms, n_param, y, X, V, X.offset, V.offset, K, Kmin, - mixture_code, threads) + nll_pcount(parms, n_param, y, dm$X_state, dm$X_det, dm$offset_state, dm$offset_det, + K, Kmin, mixture_code, threads) } } @@ -120,10 +118,10 @@ pcount <- function(formula, data, K, mixture = c("P", "NB", "ZIP"), starts, forms <- split_formula(formula) obs_all <- add_covariates(obsCovs(data), siteCovs(data), length(getY(data))) inps <- get_ranef_inputs(forms, list(det=obs_all, state=siteCovs(data)), - list(V, X), designMats[c("Z_det","Z_state")]) + list(dm$X_det, dm$X_state), dm[c("Z_det","Z_state")]) tmb_dat <- c(list(y=y, K=K, Kmin=Kmin, mixture=mixture_code, - offset_state=X.offset, offset_det=V.offset), inps$data) + offset_state=dm$offset_state, offset_det=dm$offset_det), inps$data) tmb_param <- c(inps$pars, list(beta_scale=rep(0,0))) if(mixture_code > 1) tmb_param$beta_scale <- rep(0,1) @@ -184,7 +182,7 @@ pcount <- function(formula, data, K, mixture = c("P", "NB", "ZIP"), starts, # Create unmarkedFit object------------------------------------------------ umfit <- new("unmarkedFitPCount", fitType="pcount", call=match.call(), formula = formula, data = data, - sitesRemoved = designMats$removed.sites, + sitesRemoved = dm$removed.sites, estimates = estimateList, AIC = fmAIC, opt = fm, negLogLike = fm$value, nllFun = nll, K = K, mixture = mixture, TMB=tmb_mod) diff --git a/R/pcount.spHDS.R b/R/pcount.spHDS.R index 093a1ee7..489e492f 100644 --- a/R/pcount.spHDS.R +++ b/R/pcount.spHDS.R @@ -20,25 +20,16 @@ function (formula, data, K, mixture = c("P", "NB", "ZIP"), starts, #engine <- match.arg(engine, c("C", "R")) if (identical(mixture, "ZIP") & identical(engine, "R")) stop("ZIP mixture not available for engine='R'") - designMats <- getDesign(data, formula) + dm <- getDesign(data, formula) + y <- dm$y - X <- designMats$X - V <- designMats$V - y <- designMats$y - X.offset <- designMats$X.offset - V.offset <- designMats$V.offset - - if (is.null(X.offset)) - X.offset <- rep(0, nrow(X)) - if (is.null(V.offset)) - V.offset <- rep(0, nrow(V)) NAmat <- is.na(y) J <- ncol(y) M <- nrow(y) - lamParms <- colnames(X) - detParms <- colnames(V) - nDP <- ncol(V) - nAP <- ncol(X) + lamParms <- colnames(dm$X_state) + detParms <- colnames(dm$X_det) + nDP <- ncol(dm$X_det) + nAP <- ncol(dm$X_state) if (missing(K)) { K <- max(y, na.rm = TRUE) + 100 warning("K was not specified and was set to ", K, ".") @@ -68,9 +59,9 @@ function (formula, data, K, mixture = c("P", "NB", "ZIP"), starts, tmp.parms[1]<- -1*exp(tmp.parms[1]) parms[(nAP + 1):(nAP + nDP)]<- tmp.parms - theta.i <- exp(X %*% parms[1:nAP] + X.offset) - p.ij <- 2*plogis(V %*% (parms[(nAP + 1):(nAP + nDP)]) + - V.offset) + theta.i <- exp(dm$X_state %*% parms[1:nAP] + dm$offset_state) + p.ij <- 2*plogis(dm$X_det %*% (parms[(nAP + 1):(nAP + nDP)]) + + dm$offset_det) theta.ik <- rep(theta.i, each = K + 1) p.ijk <- rep(p.ij, each = K + 1) bin.ijk <- dbinom(y.ijk, k.ijk, p.ijk) @@ -120,7 +111,7 @@ function (formula, data, K, mixture = c("P", "NB", "ZIP"), starts, nP]), invlink = "logistic", invlinkGrad = "logistic.grad") } umfit <- new("unmarkedFitPCount", fitType = "pcount", call = match.call(), - formula = formula, data = data, sitesRemoved = designMats$removed.sites, + formula = formula, data = data, sitesRemoved = dm$removed.sites, estimates = estimateList, AIC = fmAIC, opt = fm, negLogLike = fm$value, nllFun = nll, K = K, mixture = mixture) return(umfit) diff --git a/R/pcountOpen.R b/R/pcountOpen.R index 9156aada..2e22860a 100644 --- a/R/pcountOpen.R +++ b/R/pcountOpen.R @@ -22,29 +22,12 @@ formula <- as.formula(paste(unlist(formlist), collapse=" ")) D <- getDesign(data, formula) y <- D$y -Xlam <- D$Xlam -Xgam <- D$Xgam -Xom <- D$Xom -Xp <- D$Xp -Xiota <- D$Xiota - delta <- D$delta; go.dims <- D$go.dims deltamax <- max(delta, na.rm=TRUE) M <- nrow(y) T <- data@numPrimary J <- ncol(getY(data)) / T -Xlam.offset <- D$Xlam.offset -Xgam.offset <- D$Xgam.offset -Xom.offset <- D$Xom.offset -Xp.offset <- D$Xp.offset -Xiota.offset <- D$Xiota.offset -if(is.null(Xlam.offset)) Xlam.offset <- rep(0, M) -if(is.null(Xgam.offset)) Xgam.offset <- rep(0, M*(T-1)) -if(is.null(Xom.offset)) Xom.offset <- rep(0, M*(T-1)) -if(is.null(Xp.offset)) Xp.offset <- rep(0, M*T*J) -if(is.null(Xiota.offset)) Xiota.offset <- rep(0, M*(T-1)) - yna <- is.na(y) yna[] <- as.integer(yna) y <- array(y, c(M, J, T)) @@ -64,18 +47,18 @@ if(K <= max(y, na.rm = TRUE)) k <- 0:K lk <- length(k) -lamParms <- colnames(Xlam) -gamParms <- colnames(Xgam) -omParms <- colnames(Xom) -detParms <- colnames(Xp) -nAP <- ncol(Xlam) -nGP <- ncol(Xgam) -nOP <- ncol(Xom) -nDP <- ncol(Xp) +lamParms <- colnames(D$X_lambda) +gamParms <- colnames(D$X_gamma) +omParms <- colnames(D$X_omega) +detParms <- colnames(D$X_det) +nAP <- ncol(D$X_lambda) +nGP <- ncol(D$X_gamma) +nOP <- ncol(D$X_omega) +nDP <- ncol(D$X_det) if(immigration) { - iotaParms <- colnames(Xiota) - nIP <- ncol(Xiota) + iotaParms <- colnames(D$X_iota) + nIP <- ncol(D$X_iota) } else { nIP <- 0 iotaParms <- character(0) @@ -150,9 +133,9 @@ nll <- function(parms) { log.alpha <- parms[nP] nll_pcountOpen( ym, - Xlam, Xgam, Xom, Xp, Xiota, + D$X_lambda, D$X_gamma, D$X_omega, D$X_det, D$X_iota, beta.lam, beta.gam, beta.om, beta.p, beta.iota, log.alpha, - Xlam.offset, Xgam.offset, Xom.offset, Xp.offset, Xiota.offset, + D$offset_lambda, D$offset_gamma, D$offset_omega, D$offset_det, D$offset_iota, ytna, yna, lk, mixture, first, last, M, J, T, delta, dynamics, fix, go.dims, immigration, diff --git a/R/predict.R b/R/predict.R index 86aeac1b..0f08e73c 100644 --- a/R/predict.R +++ b/R/predict.R @@ -130,10 +130,10 @@ check_type <- function(mod, type){ setMethod("predict_inputs_from_umf", "unmarkedFit", function(object, type, newdata, na.rm, re.form){ designMats <- getDesign(newdata, object@formula, na.rm = na.rm) - if(type == "state") list_els <- c("X","Z_state","X.offset") - if(type == "det") list_els <- c("V","Z_det","V.offset") + if(type == "state") list_els <- c("X_state","Z_state","offset_state") + if(type == "det") list_els <- c("X_det","Z_det","offset_det") if(type == "scale"){ # no covariates - n <- nrow(designMats$V) + n <- nrow(designMats$X_det) return(list(X = matrix(1, nrow=n, ncol=1), offset = rep(0, n))) } @@ -323,7 +323,7 @@ setMethod("predict_by_chunk", "unmarkedFitPCount", setMethod("predict_inputs_from_umf", "unmarkedFitColExt", function(object, type, newdata, na.rm, re.form=NA){ designMats <- getDesign(newdata, object@formula, na.rm = na.rm) - list_el <- switch(type, psi="W", col="X.gam", ext="X.eps", det="V") + list_el <- switch(type, psi="X_state", col="X_col", ext="X_ext", det="X_det") # colext doesn't support offsets list(X=designMats[[list_el]], offset=NULL) }) @@ -347,8 +347,8 @@ setMethod("predict_inputs_from_umf", "unmarkedFitOccuFP", function(object, type, newdata, na.rm, re.form=NA){ designMats <- getDesign(newdata, object@detformula, object@FPformula, object@Bformula, object@stateformula, na.rm=na.rm) - X_idx <- switch(type, state="X", det="V", fp="U", b="W") - off_idx <- paste0(X_idx, ".offset") + X_idx <- switch(type, state="X_state", det="X_det", fp="X_fp", b="X_b") + off_idx <- paste0("offset_", type) list(X=designMats[[X_idx]], offset=designMats[[off_idx]]) }) @@ -387,9 +387,9 @@ setMethod("check_predict_arguments", "unmarkedFitDailMadsen", setMethod("predict_inputs_from_umf", "unmarkedFitDailMadsen", function(object, type, newdata, na.rm, re.form=NA){ designMats <- getDesign(newdata, object@formula, na.rm=na.rm) - X_idx <- switch(type, lambda="Xlam", gamma="Xgam", omega="Xom", - iota="Xiota", det="Xp") - off_idx <- paste0(X_idx, ".offset") + X_idx <- switch(type, lambda="X_lambda", gamma="X_gamma", omega="X_omega", + iota="X_iota", det="X_det") + off_idx <- paste0("offset_", type) list(X=designMats[[X_idx]], offset=designMats[[off_idx]]) }) @@ -444,8 +444,8 @@ setMethod("predict_by_chunk", "unmarkedFitDailMadsen", setMethod("predict_inputs_from_umf", "unmarkedFitGMM", function(object, type, newdata, na.rm, re.form=NA){ designMats <- getDesign(newdata, object@formula, na.rm=na.rm) - X_idx <- switch(type, lambda="Xlam", phi="Xphi", det="Xdet") - off_idx <- paste0(X_idx, ".offset") + X_idx <- switch(type, lambda="X_state", phi="X_phi", det="X_det") + off_idx <- switch(type, lambda="offset_state", phi="offset_phi", det="offset_det") list(X=designMats[[X_idx]], offset=designMats[[off_idx]]) }) @@ -469,7 +469,7 @@ setMethod("get_orig_data", "unmarkedFitGMM", function(object, type, ...){ setMethod("predict_inputs_from_umf", "unmarkedFitOccuTTD", function(object, type, newdata, na.rm, re.form=NA){ designMats <- getDesign(newdata, object@formula, na.rm = na.rm) - list_el <- switch(type, psi="W", col="X.gam", ext="X.eps", det="V") + list_el <- switch(type, psi="X_state", col="X_col", ext="X_ext", det="X_det") list(X=designMats[[list_el]], offset=NULL) }) @@ -491,7 +491,7 @@ setMethod("get_orig_data", "unmarkedFitOccuTTD", function(object, type, ...){ setMethod("predict_inputs_from_umf", "unmarkedFitNmixTTD", function(object, type, newdata, na.rm, re.form=NA){ designMats <- getDesign(newdata, object@formula, na.rm = na.rm) - list_el <- switch(type, state="W", det="V") + list_el <- switch(type, state="X_state", det="X_det") list(X=designMats[[list_el]], offset=NULL) }) @@ -511,10 +511,10 @@ setMethod("get_orig_data", "unmarkedFitNmixTTD", function(object, type, ...){ setMethod("predict_inputs_from_umf", "unmarkedFitGDR", function(object, type, newdata, na.rm, re.form=NA){ designMats <- getDesign(newdata, object@formlist) - if(type == "lambda") list_els <- c("Xlam","Zlam") - if(type == "phi") list_els <- c("Xphi","Zphi") - if(type == "dist") list_els <- c("Xdist","Zdist") - if(type == "rem") list_els <- c("Xrem", "Zrem") + if(type == "lambda") list_els <- c("X_lambda","Z_lambda") + if(type == "phi") list_els <- c("X_phi","Z_phi") + if(type == "dist") list_els <- c("X_dist","Z_dist") + if(type == "rem") list_els <- c("X_rem", "Z_rem") X <- designMats[[list_els[1]]] if(is.null(re.form)) X <- cbind(X, designMats[[list_els[2]]]) list(X=X, offset=NULL) diff --git a/R/unmarkedFit.R b/R/unmarkedFit.R index f09594fb..555133a3 100644 --- a/R/unmarkedFit.R +++ b/R/unmarkedFit.R @@ -473,14 +473,10 @@ setMethod("getFP", "unmarkedFitOccuFP", function(object, na.rm = TRUE) designMats <- getDesign(umf, detformula,FPformula,Bformula,stateformula, na.rm = na.rm) type = object@type y <- designMats$y - U <- designMats$U - U.offset <- designMats$U.offset - if (is.null(U.offset)) - U.offset <- rep(0, nrow(U)) M <- nrow(y) J <- ncol(y) fpars <- coef(object, type = "fp") - f <- plogis(U %*% fpars + U.offset) + f <- plogis(designMats$X_fp %*% fpars + designMats$offset_fp) f <- matrix(f, M, J, byrow = TRUE) if (type[1]!=0){ f[,1:type[1]] = 0 @@ -498,16 +494,12 @@ setMethod("getB", "unmarkedFitOccuFP", function(object, na.rm = TRUE) umf <- object@data designMats <- getDesign(umf, detformula,FPformula,Bformula,stateformula, na.rm = na.rm) y <- designMats$y - W <- designMats$W - W.offset <- designMats$W.offset - if (is.null(W.offset)) - W.offset <- rep(0, nrow(W)) M <- nrow(y) J <- ncol(y) type = object@type if (type[3]!=0){ bpars <- coef(object, type = "b") - b <- plogis(W %*% bpars + W.offset) + b <- plogis(designMats$X_b %*% bpars + designMats$offset_b) b <- matrix(b, M, J, byrow = TRUE) } if (type[3]==0){ diff --git a/R/zzz.R b/R/zzz.R deleted file mode 100644 index 063c71ce..00000000 --- a/R/zzz.R +++ /dev/null @@ -1 +0,0 @@ -.shiny_env <- new.env(parent=emptyenv()) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index dcf8ea60..f31ff0db 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -59,25 +59,25 @@ BEGIN_RCPP END_RCPP } // nll_distsampOpen -double nll_distsampOpen(arma::ucube y, arma::imat yt, arma::mat Xlam, arma::mat Xgam, arma::mat Xom, arma::mat Xsig, arma::mat Xiota, arma::vec beta, arma::umat bi, arma::colvec Xlam_offset, arma::colvec Xgam_offset, arma::colvec Xom_offset, arma::colvec Xsig_offset, arma::colvec Xiota_offset, arma::imat ytna, int lk, std::string mixture, Rcpp::IntegerVector first, Rcpp::IntegerVector last, int first1, int M, int T, arma::imat delta, std::string dynamics, std::string survey, std::string fix, std::string go_dims, bool immigration, arma::imat I, arma::imat I1, Rcpp::List Ib, Rcpp::List Ip, arma::mat a, arma::mat u, arma::vec w, arma::vec db, std::string keyfun, arma::vec lfac_k, arma::cube kmyt, arma::cube lfac_kmyt, arma::icube fin, arma::vec A); -RcppExport SEXP _unmarked_nll_distsampOpen(SEXP ySEXP, SEXP ytSEXP, SEXP XlamSEXP, SEXP XgamSEXP, SEXP XomSEXP, SEXP XsigSEXP, SEXP XiotaSEXP, SEXP betaSEXP, SEXP biSEXP, SEXP Xlam_offsetSEXP, SEXP Xgam_offsetSEXP, SEXP Xom_offsetSEXP, SEXP Xsig_offsetSEXP, SEXP Xiota_offsetSEXP, SEXP ytnaSEXP, SEXP lkSEXP, SEXP mixtureSEXP, SEXP firstSEXP, SEXP lastSEXP, SEXP first1SEXP, SEXP MSEXP, SEXP TSEXP, SEXP deltaSEXP, SEXP dynamicsSEXP, SEXP surveySEXP, SEXP fixSEXP, SEXP go_dimsSEXP, SEXP immigrationSEXP, SEXP ISEXP, SEXP I1SEXP, SEXP IbSEXP, SEXP IpSEXP, SEXP aSEXP, SEXP uSEXP, SEXP wSEXP, SEXP dbSEXP, SEXP keyfunSEXP, SEXP lfac_kSEXP, SEXP kmytSEXP, SEXP lfac_kmytSEXP, SEXP finSEXP, SEXP ASEXP) { +double nll_distsampOpen(arma::ucube y, arma::imat yt, arma::mat X_lambda, arma::mat X_gamma, arma::mat X_omega, arma::mat X_det, arma::mat X_iota, arma::vec beta, arma::umat bi, arma::colvec offset_lambda, arma::colvec offset_gamma, arma::colvec offset_omega, arma::colvec offset_det, arma::colvec offset_iota, arma::imat ytna, int lk, std::string mixture, Rcpp::IntegerVector first, Rcpp::IntegerVector last, int first1, int M, int T, arma::imat delta, std::string dynamics, std::string survey, std::string fix, std::string go_dims, bool immigration, arma::imat I, arma::imat I1, Rcpp::List Ib, Rcpp::List Ip, arma::mat a, arma::mat u, arma::vec w, arma::vec db, std::string keyfun, arma::vec lfac_k, arma::cube kmyt, arma::cube lfac_kmyt, arma::icube fin, arma::vec A); +RcppExport SEXP _unmarked_nll_distsampOpen(SEXP ySEXP, SEXP ytSEXP, SEXP X_lambdaSEXP, SEXP X_gammaSEXP, SEXP X_omegaSEXP, SEXP X_detSEXP, SEXP X_iotaSEXP, SEXP betaSEXP, SEXP biSEXP, SEXP offset_lambdaSEXP, SEXP offset_gammaSEXP, SEXP offset_omegaSEXP, SEXP offset_detSEXP, SEXP offset_iotaSEXP, SEXP ytnaSEXP, SEXP lkSEXP, SEXP mixtureSEXP, SEXP firstSEXP, SEXP lastSEXP, SEXP first1SEXP, SEXP MSEXP, SEXP TSEXP, SEXP deltaSEXP, SEXP dynamicsSEXP, SEXP surveySEXP, SEXP fixSEXP, SEXP go_dimsSEXP, SEXP immigrationSEXP, SEXP ISEXP, SEXP I1SEXP, SEXP IbSEXP, SEXP IpSEXP, SEXP aSEXP, SEXP uSEXP, SEXP wSEXP, SEXP dbSEXP, SEXP keyfunSEXP, SEXP lfac_kSEXP, SEXP kmytSEXP, SEXP lfac_kmytSEXP, SEXP finSEXP, SEXP ASEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< arma::ucube >::type y(ySEXP); Rcpp::traits::input_parameter< arma::imat >::type yt(ytSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xlam(XlamSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xgam(XgamSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xom(XomSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xsig(XsigSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xiota(XiotaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_lambda(X_lambdaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_gamma(X_gammaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_omega(X_omegaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_det(X_detSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_iota(X_iotaSEXP); Rcpp::traits::input_parameter< arma::vec >::type beta(betaSEXP); Rcpp::traits::input_parameter< arma::umat >::type bi(biSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xlam_offset(Xlam_offsetSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xgam_offset(Xgam_offsetSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xom_offset(Xom_offsetSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xsig_offset(Xsig_offsetSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xiota_offset(Xiota_offsetSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_lambda(offset_lambdaSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_gamma(offset_gammaSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_omega(offset_omegaSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_det(offset_detSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_iota(offset_iotaSEXP); Rcpp::traits::input_parameter< arma::imat >::type ytna(ytnaSEXP); Rcpp::traits::input_parameter< int >::type lk(lkSEXP); Rcpp::traits::input_parameter< std::string >::type mixture(mixtureSEXP); @@ -106,13 +106,13 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::cube >::type lfac_kmyt(lfac_kmytSEXP); Rcpp::traits::input_parameter< arma::icube >::type fin(finSEXP); Rcpp::traits::input_parameter< arma::vec >::type A(ASEXP); - rcpp_result_gen = Rcpp::wrap(nll_distsampOpen(y, yt, Xlam, Xgam, Xom, Xsig, Xiota, beta, bi, Xlam_offset, Xgam_offset, Xom_offset, Xsig_offset, Xiota_offset, ytna, lk, mixture, first, last, first1, M, T, delta, dynamics, survey, fix, go_dims, immigration, I, I1, Ib, Ip, a, u, w, db, keyfun, lfac_k, kmyt, lfac_kmyt, fin, A)); + rcpp_result_gen = Rcpp::wrap(nll_distsampOpen(y, yt, X_lambda, X_gamma, X_omega, X_det, X_iota, beta, bi, offset_lambda, offset_gamma, offset_omega, offset_det, offset_iota, ytna, lk, mixture, first, last, first1, M, T, delta, dynamics, survey, fix, go_dims, immigration, I, I1, Ib, Ip, a, u, w, db, keyfun, lfac_k, kmyt, lfac_kmyt, fin, A)); return rcpp_result_gen; END_RCPP } // nll_gdistremoval -double nll_gdistremoval(arma::vec beta, arma::uvec n_param, arma::vec yDistance, arma::vec yRemoval, arma::mat ysum, int mixture, std::string keyfun, arma::mat Xlam, arma::vec A, arma::mat Xphi, arma::mat Xrem, arma::mat Xdist, arma::vec db, arma::mat a, arma::mat u, arma::vec w, arma::uvec pl, int K, arma::uvec Kmin, int threads); -RcppExport SEXP _unmarked_nll_gdistremoval(SEXP betaSEXP, SEXP n_paramSEXP, SEXP yDistanceSEXP, SEXP yRemovalSEXP, SEXP ysumSEXP, SEXP mixtureSEXP, SEXP keyfunSEXP, SEXP XlamSEXP, SEXP ASEXP, SEXP XphiSEXP, SEXP XremSEXP, SEXP XdistSEXP, SEXP dbSEXP, SEXP aSEXP, SEXP uSEXP, SEXP wSEXP, SEXP plSEXP, SEXP KSEXP, SEXP KminSEXP, SEXP threadsSEXP) { +double nll_gdistremoval(arma::vec beta, arma::uvec n_param, arma::vec yDistance, arma::vec yRemoval, arma::mat ysum, int mixture, std::string keyfun, arma::mat X_lambda, arma::vec A, arma::mat X_phi, arma::mat X_rem, arma::mat X_dist, arma::vec db, arma::mat a, arma::mat u, arma::vec w, arma::uvec pl, int K, arma::uvec Kmin, int threads); +RcppExport SEXP _unmarked_nll_gdistremoval(SEXP betaSEXP, SEXP n_paramSEXP, SEXP yDistanceSEXP, SEXP yRemovalSEXP, SEXP ysumSEXP, SEXP mixtureSEXP, SEXP keyfunSEXP, SEXP X_lambdaSEXP, SEXP ASEXP, SEXP X_phiSEXP, SEXP X_remSEXP, SEXP X_distSEXP, SEXP dbSEXP, SEXP aSEXP, SEXP uSEXP, SEXP wSEXP, SEXP plSEXP, SEXP KSEXP, SEXP KminSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -123,11 +123,11 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::mat >::type ysum(ysumSEXP); Rcpp::traits::input_parameter< int >::type mixture(mixtureSEXP); Rcpp::traits::input_parameter< std::string >::type keyfun(keyfunSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xlam(XlamSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_lambda(X_lambdaSEXP); Rcpp::traits::input_parameter< arma::vec >::type A(ASEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xphi(XphiSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xrem(XremSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xdist(XdistSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_phi(X_phiSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_rem(X_remSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_dist(X_distSEXP); Rcpp::traits::input_parameter< arma::vec >::type db(dbSEXP); Rcpp::traits::input_parameter< arma::mat >::type a(aSEXP); Rcpp::traits::input_parameter< arma::mat >::type u(uSEXP); @@ -136,13 +136,13 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type K(KSEXP); Rcpp::traits::input_parameter< arma::uvec >::type Kmin(KminSEXP); Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); - rcpp_result_gen = Rcpp::wrap(nll_gdistremoval(beta, n_param, yDistance, yRemoval, ysum, mixture, keyfun, Xlam, A, Xphi, Xrem, Xdist, db, a, u, w, pl, K, Kmin, threads)); + rcpp_result_gen = Rcpp::wrap(nll_gdistremoval(beta, n_param, yDistance, yRemoval, ysum, mixture, keyfun, X_lambda, A, X_phi, X_rem, X_dist, db, a, u, w, pl, K, Kmin, threads)); return rcpp_result_gen; END_RCPP } // nll_gdistsamp -double nll_gdistsamp(arma::vec beta, arma::uvec n_param, arma::vec y, int mixture, std::string keyfun, std::string survey, arma::mat Xlam, arma::vec Xlam_offset, arma::vec A, arma::mat Xphi, arma::vec Xphi_offset, arma::mat Xdet, arma::vec Xdet_offset, arma::vec db, arma::mat a, arma::mat u, arma::vec w, arma::vec k, arma::vec lfac_k, arma::vec lfac_kmyt, arma::vec kmyt, arma::uvec Kmin, int threads); -RcppExport SEXP _unmarked_nll_gdistsamp(SEXP betaSEXP, SEXP n_paramSEXP, SEXP ySEXP, SEXP mixtureSEXP, SEXP keyfunSEXP, SEXP surveySEXP, SEXP XlamSEXP, SEXP Xlam_offsetSEXP, SEXP ASEXP, SEXP XphiSEXP, SEXP Xphi_offsetSEXP, SEXP XdetSEXP, SEXP Xdet_offsetSEXP, SEXP dbSEXP, SEXP aSEXP, SEXP uSEXP, SEXP wSEXP, SEXP kSEXP, SEXP lfac_kSEXP, SEXP lfac_kmytSEXP, SEXP kmytSEXP, SEXP KminSEXP, SEXP threadsSEXP) { +double nll_gdistsamp(arma::vec beta, arma::uvec n_param, arma::vec y, int mixture, std::string keyfun, std::string survey, arma::mat X_lambda, arma::vec offset_lambda, arma::vec A, arma::mat X_phi, arma::vec offset_phi, arma::mat X_det, arma::vec offset_det, arma::vec db, arma::mat a, arma::mat u, arma::vec w, arma::vec k, arma::vec lfac_k, arma::vec lfac_kmyt, arma::vec kmyt, arma::uvec Kmin, int threads); +RcppExport SEXP _unmarked_nll_gdistsamp(SEXP betaSEXP, SEXP n_paramSEXP, SEXP ySEXP, SEXP mixtureSEXP, SEXP keyfunSEXP, SEXP surveySEXP, SEXP X_lambdaSEXP, SEXP offset_lambdaSEXP, SEXP ASEXP, SEXP X_phiSEXP, SEXP offset_phiSEXP, SEXP X_detSEXP, SEXP offset_detSEXP, SEXP dbSEXP, SEXP aSEXP, SEXP uSEXP, SEXP wSEXP, SEXP kSEXP, SEXP lfac_kSEXP, SEXP lfac_kmytSEXP, SEXP kmytSEXP, SEXP KminSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -152,13 +152,13 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type mixture(mixtureSEXP); Rcpp::traits::input_parameter< std::string >::type keyfun(keyfunSEXP); Rcpp::traits::input_parameter< std::string >::type survey(surveySEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xlam(XlamSEXP); - Rcpp::traits::input_parameter< arma::vec >::type Xlam_offset(Xlam_offsetSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_lambda(X_lambdaSEXP); + Rcpp::traits::input_parameter< arma::vec >::type offset_lambda(offset_lambdaSEXP); Rcpp::traits::input_parameter< arma::vec >::type A(ASEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xphi(XphiSEXP); - Rcpp::traits::input_parameter< arma::vec >::type Xphi_offset(Xphi_offsetSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xdet(XdetSEXP); - Rcpp::traits::input_parameter< arma::vec >::type Xdet_offset(Xdet_offsetSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_phi(X_phiSEXP); + Rcpp::traits::input_parameter< arma::vec >::type offset_phi(offset_phiSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_det(X_detSEXP); + Rcpp::traits::input_parameter< arma::vec >::type offset_det(offset_detSEXP); Rcpp::traits::input_parameter< arma::vec >::type db(dbSEXP); Rcpp::traits::input_parameter< arma::mat >::type a(aSEXP); Rcpp::traits::input_parameter< arma::mat >::type u(uSEXP); @@ -169,13 +169,13 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::vec >::type kmyt(kmytSEXP); Rcpp::traits::input_parameter< arma::uvec >::type Kmin(KminSEXP); Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); - rcpp_result_gen = Rcpp::wrap(nll_gdistsamp(beta, n_param, y, mixture, keyfun, survey, Xlam, Xlam_offset, A, Xphi, Xphi_offset, Xdet, Xdet_offset, db, a, u, w, k, lfac_k, lfac_kmyt, kmyt, Kmin, threads)); + rcpp_result_gen = Rcpp::wrap(nll_gdistsamp(beta, n_param, y, mixture, keyfun, survey, X_lambda, offset_lambda, A, X_phi, offset_phi, X_det, offset_det, db, a, u, w, k, lfac_k, lfac_kmyt, kmyt, Kmin, threads)); return rcpp_result_gen; END_RCPP } // nll_gmultmix -double nll_gmultmix(arma::vec beta, arma::uvec n_param, arma::vec y, int mixture, std::string pi_fun, arma::mat Xlam, arma::vec Xlam_offset, arma::mat Xphi, arma::vec Xphi_offset, arma::mat Xdet, arma::vec Xdet_offset, arma::vec k, arma::vec lfac_k, arma::vec lfac_kmyt, arma::vec kmyt, arma::vec Kmin, int threads); -RcppExport SEXP _unmarked_nll_gmultmix(SEXP betaSEXP, SEXP n_paramSEXP, SEXP ySEXP, SEXP mixtureSEXP, SEXP pi_funSEXP, SEXP XlamSEXP, SEXP Xlam_offsetSEXP, SEXP XphiSEXP, SEXP Xphi_offsetSEXP, SEXP XdetSEXP, SEXP Xdet_offsetSEXP, SEXP kSEXP, SEXP lfac_kSEXP, SEXP lfac_kmytSEXP, SEXP kmytSEXP, SEXP KminSEXP, SEXP threadsSEXP) { +double nll_gmultmix(arma::vec beta, arma::uvec n_param, arma::vec y, int mixture, std::string pi_fun, arma::mat X_lambda, arma::vec offset_lambda, arma::mat X_phi, arma::vec offset_phi, arma::mat X_det, arma::vec offset_det, arma::vec k, arma::vec lfac_k, arma::vec lfac_kmyt, arma::vec kmyt, arma::vec Kmin, int threads); +RcppExport SEXP _unmarked_nll_gmultmix(SEXP betaSEXP, SEXP n_paramSEXP, SEXP ySEXP, SEXP mixtureSEXP, SEXP pi_funSEXP, SEXP X_lambdaSEXP, SEXP offset_lambdaSEXP, SEXP X_phiSEXP, SEXP offset_phiSEXP, SEXP X_detSEXP, SEXP offset_detSEXP, SEXP kSEXP, SEXP lfac_kSEXP, SEXP lfac_kmytSEXP, SEXP kmytSEXP, SEXP KminSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -184,44 +184,44 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::vec >::type y(ySEXP); Rcpp::traits::input_parameter< int >::type mixture(mixtureSEXP); Rcpp::traits::input_parameter< std::string >::type pi_fun(pi_funSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xlam(XlamSEXP); - Rcpp::traits::input_parameter< arma::vec >::type Xlam_offset(Xlam_offsetSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xphi(XphiSEXP); - Rcpp::traits::input_parameter< arma::vec >::type Xphi_offset(Xphi_offsetSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xdet(XdetSEXP); - Rcpp::traits::input_parameter< arma::vec >::type Xdet_offset(Xdet_offsetSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_lambda(X_lambdaSEXP); + Rcpp::traits::input_parameter< arma::vec >::type offset_lambda(offset_lambdaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_phi(X_phiSEXP); + Rcpp::traits::input_parameter< arma::vec >::type offset_phi(offset_phiSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_det(X_detSEXP); + Rcpp::traits::input_parameter< arma::vec >::type offset_det(offset_detSEXP); Rcpp::traits::input_parameter< arma::vec >::type k(kSEXP); Rcpp::traits::input_parameter< arma::vec >::type lfac_k(lfac_kSEXP); Rcpp::traits::input_parameter< arma::vec >::type lfac_kmyt(lfac_kmytSEXP); Rcpp::traits::input_parameter< arma::vec >::type kmyt(kmytSEXP); Rcpp::traits::input_parameter< arma::vec >::type Kmin(KminSEXP); Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); - rcpp_result_gen = Rcpp::wrap(nll_gmultmix(beta, n_param, y, mixture, pi_fun, Xlam, Xlam_offset, Xphi, Xphi_offset, Xdet, Xdet_offset, k, lfac_k, lfac_kmyt, kmyt, Kmin, threads)); + rcpp_result_gen = Rcpp::wrap(nll_gmultmix(beta, n_param, y, mixture, pi_fun, X_lambda, offset_lambda, X_phi, offset_phi, X_det, offset_det, k, lfac_k, lfac_kmyt, kmyt, Kmin, threads)); return rcpp_result_gen; END_RCPP } // nll_gpcount -double nll_gpcount(arma::mat ym, arma::mat Xlam, arma::mat Xphi, arma::mat Xp, arma::vec beta_lam, arma::vec beta_phi, arma::vec beta_p, double log_alpha, arma::vec Xlam_offset, arma::vec Xphi_offset, arma::vec Xp_offset, int M, std::string mixture, int T, int threads); -RcppExport SEXP _unmarked_nll_gpcount(SEXP ymSEXP, SEXP XlamSEXP, SEXP XphiSEXP, SEXP XpSEXP, SEXP beta_lamSEXP, SEXP beta_phiSEXP, SEXP beta_pSEXP, SEXP log_alphaSEXP, SEXP Xlam_offsetSEXP, SEXP Xphi_offsetSEXP, SEXP Xp_offsetSEXP, SEXP MSEXP, SEXP mixtureSEXP, SEXP TSEXP, SEXP threadsSEXP) { +double nll_gpcount(arma::mat ym, arma::mat X_lambda, arma::mat X_phi, arma::mat X_det, arma::vec beta_lambda, arma::vec beta_phi, arma::vec beta_det, double log_alpha, arma::vec offset_lambda, arma::vec offset_phi, arma::vec offset_det, int M, std::string mixture, int T, int threads); +RcppExport SEXP _unmarked_nll_gpcount(SEXP ymSEXP, SEXP X_lambdaSEXP, SEXP X_phiSEXP, SEXP X_detSEXP, SEXP beta_lambdaSEXP, SEXP beta_phiSEXP, SEXP beta_detSEXP, SEXP log_alphaSEXP, SEXP offset_lambdaSEXP, SEXP offset_phiSEXP, SEXP offset_detSEXP, SEXP MSEXP, SEXP mixtureSEXP, SEXP TSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< arma::mat >::type ym(ymSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xlam(XlamSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xphi(XphiSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xp(XpSEXP); - Rcpp::traits::input_parameter< arma::vec >::type beta_lam(beta_lamSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_lambda(X_lambdaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_phi(X_phiSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_det(X_detSEXP); + Rcpp::traits::input_parameter< arma::vec >::type beta_lambda(beta_lambdaSEXP); Rcpp::traits::input_parameter< arma::vec >::type beta_phi(beta_phiSEXP); - Rcpp::traits::input_parameter< arma::vec >::type beta_p(beta_pSEXP); + Rcpp::traits::input_parameter< arma::vec >::type beta_det(beta_detSEXP); Rcpp::traits::input_parameter< double >::type log_alpha(log_alphaSEXP); - Rcpp::traits::input_parameter< arma::vec >::type Xlam_offset(Xlam_offsetSEXP); - Rcpp::traits::input_parameter< arma::vec >::type Xphi_offset(Xphi_offsetSEXP); - Rcpp::traits::input_parameter< arma::vec >::type Xp_offset(Xp_offsetSEXP); + Rcpp::traits::input_parameter< arma::vec >::type offset_lambda(offset_lambdaSEXP); + Rcpp::traits::input_parameter< arma::vec >::type offset_phi(offset_phiSEXP); + Rcpp::traits::input_parameter< arma::vec >::type offset_det(offset_detSEXP); Rcpp::traits::input_parameter< int >::type M(MSEXP); Rcpp::traits::input_parameter< std::string >::type mixture(mixtureSEXP); Rcpp::traits::input_parameter< int >::type T(TSEXP); Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); - rcpp_result_gen = Rcpp::wrap(nll_gpcount(ym, Xlam, Xphi, Xp, beta_lam, beta_phi, beta_p, log_alpha, Xlam_offset, Xphi_offset, Xp_offset, M, mixture, T, threads)); + rcpp_result_gen = Rcpp::wrap(nll_gpcount(ym, X_lambda, X_phi, X_det, beta_lambda, beta_phi, beta_det, log_alpha, offset_lambda, offset_phi, offset_det, M, mixture, T, threads)); return rcpp_result_gen; END_RCPP } @@ -246,25 +246,25 @@ BEGIN_RCPP END_RCPP } // nll_multmixOpen -double nll_multmixOpen(arma::ucube y, arma::imat yt, arma::mat Xlam, arma::mat Xgam, arma::mat Xom, arma::mat Xp, arma::mat Xiota, arma::vec beta, arma::umat bi, arma::colvec Xlam_offset, arma::colvec Xgam_offset, arma::colvec Xom_offset, arma::colvec Xp_offset, arma::colvec Xiota_offset, arma::imat ytna, arma::ucube yna, int lk, std::string mixture, Rcpp::IntegerVector first, Rcpp::IntegerVector last, int first1, int M, int T, int J, int R, arma::imat delta, std::string dynamics, std::string fix, std::string go_dims, bool immigration, arma::imat I, arma::imat I1, Rcpp::List Ib, Rcpp::List Ip, std::string pi_fun, arma::vec lfac_k, arma::cube kmyt, arma::cube lfac_kmyt, arma::icube fin); -RcppExport SEXP _unmarked_nll_multmixOpen(SEXP ySEXP, SEXP ytSEXP, SEXP XlamSEXP, SEXP XgamSEXP, SEXP XomSEXP, SEXP XpSEXP, SEXP XiotaSEXP, SEXP betaSEXP, SEXP biSEXP, SEXP Xlam_offsetSEXP, SEXP Xgam_offsetSEXP, SEXP Xom_offsetSEXP, SEXP Xp_offsetSEXP, SEXP Xiota_offsetSEXP, SEXP ytnaSEXP, SEXP ynaSEXP, SEXP lkSEXP, SEXP mixtureSEXP, SEXP firstSEXP, SEXP lastSEXP, SEXP first1SEXP, SEXP MSEXP, SEXP TSEXP, SEXP JSEXP, SEXP RSEXP, SEXP deltaSEXP, SEXP dynamicsSEXP, SEXP fixSEXP, SEXP go_dimsSEXP, SEXP immigrationSEXP, SEXP ISEXP, SEXP I1SEXP, SEXP IbSEXP, SEXP IpSEXP, SEXP pi_funSEXP, SEXP lfac_kSEXP, SEXP kmytSEXP, SEXP lfac_kmytSEXP, SEXP finSEXP) { +double nll_multmixOpen(arma::ucube y, arma::imat yt, arma::mat X_lambda, arma::mat X_gamma, arma::mat X_omega, arma::mat X_det, arma::mat X_iota, arma::vec beta, arma::umat bi, arma::colvec offset_lambda, arma::colvec offset_gamma, arma::colvec offset_omega, arma::colvec offset_det, arma::colvec offset_iota, arma::imat ytna, arma::ucube yna, int lk, std::string mixture, Rcpp::IntegerVector first, Rcpp::IntegerVector last, int first1, int M, int T, int J, int R, arma::imat delta, std::string dynamics, std::string fix, std::string go_dims, bool immigration, arma::imat I, arma::imat I1, Rcpp::List Ib, Rcpp::List Ip, std::string pi_fun, arma::vec lfac_k, arma::cube kmyt, arma::cube lfac_kmyt, arma::icube fin); +RcppExport SEXP _unmarked_nll_multmixOpen(SEXP ySEXP, SEXP ytSEXP, SEXP X_lambdaSEXP, SEXP X_gammaSEXP, SEXP X_omegaSEXP, SEXP X_detSEXP, SEXP X_iotaSEXP, SEXP betaSEXP, SEXP biSEXP, SEXP offset_lambdaSEXP, SEXP offset_gammaSEXP, SEXP offset_omegaSEXP, SEXP offset_detSEXP, SEXP offset_iotaSEXP, SEXP ytnaSEXP, SEXP ynaSEXP, SEXP lkSEXP, SEXP mixtureSEXP, SEXP firstSEXP, SEXP lastSEXP, SEXP first1SEXP, SEXP MSEXP, SEXP TSEXP, SEXP JSEXP, SEXP RSEXP, SEXP deltaSEXP, SEXP dynamicsSEXP, SEXP fixSEXP, SEXP go_dimsSEXP, SEXP immigrationSEXP, SEXP ISEXP, SEXP I1SEXP, SEXP IbSEXP, SEXP IpSEXP, SEXP pi_funSEXP, SEXP lfac_kSEXP, SEXP kmytSEXP, SEXP lfac_kmytSEXP, SEXP finSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< arma::ucube >::type y(ySEXP); Rcpp::traits::input_parameter< arma::imat >::type yt(ytSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xlam(XlamSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xgam(XgamSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xom(XomSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xp(XpSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xiota(XiotaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_lambda(X_lambdaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_gamma(X_gammaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_omega(X_omegaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_det(X_detSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_iota(X_iotaSEXP); Rcpp::traits::input_parameter< arma::vec >::type beta(betaSEXP); Rcpp::traits::input_parameter< arma::umat >::type bi(biSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xlam_offset(Xlam_offsetSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xgam_offset(Xgam_offsetSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xom_offset(Xom_offsetSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xp_offset(Xp_offsetSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xiota_offset(Xiota_offsetSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_lambda(offset_lambdaSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_gamma(offset_gammaSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_omega(offset_omegaSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_det(offset_detSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_iota(offset_iotaSEXP); Rcpp::traits::input_parameter< arma::imat >::type ytna(ytnaSEXP); Rcpp::traits::input_parameter< arma::ucube >::type yna(ynaSEXP); Rcpp::traits::input_parameter< int >::type lk(lkSEXP); @@ -290,21 +290,21 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::cube >::type kmyt(kmytSEXP); Rcpp::traits::input_parameter< arma::cube >::type lfac_kmyt(lfac_kmytSEXP); Rcpp::traits::input_parameter< arma::icube >::type fin(finSEXP); - rcpp_result_gen = Rcpp::wrap(nll_multmixOpen(y, yt, Xlam, Xgam, Xom, Xp, Xiota, beta, bi, Xlam_offset, Xgam_offset, Xom_offset, Xp_offset, Xiota_offset, ytna, yna, lk, mixture, first, last, first1, M, T, J, R, delta, dynamics, fix, go_dims, immigration, I, I1, Ib, Ip, pi_fun, lfac_k, kmyt, lfac_kmyt, fin)); + rcpp_result_gen = Rcpp::wrap(nll_multmixOpen(y, yt, X_lambda, X_gamma, X_omega, X_det, X_iota, beta, bi, offset_lambda, offset_gamma, offset_omega, offset_det, offset_iota, ytna, yna, lk, mixture, first, last, first1, M, T, J, R, delta, dynamics, fix, go_dims, immigration, I, I1, Ib, Ip, pi_fun, lfac_k, kmyt, lfac_kmyt, fin)); return rcpp_result_gen; END_RCPP } // nll_nmixTTD -double nll_nmixTTD(const arma::vec beta, const arma::vec y, const arma::vec delta, const arma::mat W, const arma::mat V, const arma::umat pinds, const std::string mixture, const std::string tdist, int N, int J, int K, const arma::vec naflag, int threads); -RcppExport SEXP _unmarked_nll_nmixTTD(SEXP betaSEXP, SEXP ySEXP, SEXP deltaSEXP, SEXP WSEXP, SEXP VSEXP, SEXP pindsSEXP, SEXP mixtureSEXP, SEXP tdistSEXP, SEXP NSEXP, SEXP JSEXP, SEXP KSEXP, SEXP naflagSEXP, SEXP threadsSEXP) { +double nll_nmixTTD(const arma::vec beta, const arma::vec y, const arma::vec delta, const arma::mat X_state, const arma::mat X_det, const arma::umat pinds, const std::string mixture, const std::string tdist, int N, int J, int K, const arma::vec naflag, int threads); +RcppExport SEXP _unmarked_nll_nmixTTD(SEXP betaSEXP, SEXP ySEXP, SEXP deltaSEXP, SEXP X_stateSEXP, SEXP X_detSEXP, SEXP pindsSEXP, SEXP mixtureSEXP, SEXP tdistSEXP, SEXP NSEXP, SEXP JSEXP, SEXP KSEXP, SEXP naflagSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const arma::vec >::type beta(betaSEXP); Rcpp::traits::input_parameter< const arma::vec >::type y(ySEXP); Rcpp::traits::input_parameter< const arma::vec >::type delta(deltaSEXP); - Rcpp::traits::input_parameter< const arma::mat >::type W(WSEXP); - Rcpp::traits::input_parameter< const arma::mat >::type V(VSEXP); + Rcpp::traits::input_parameter< const arma::mat >::type X_state(X_stateSEXP); + Rcpp::traits::input_parameter< const arma::mat >::type X_det(X_detSEXP); Rcpp::traits::input_parameter< const arma::umat >::type pinds(pindsSEXP); Rcpp::traits::input_parameter< const std::string >::type mixture(mixtureSEXP); Rcpp::traits::input_parameter< const std::string >::type tdist(tdistSEXP); @@ -313,7 +313,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type K(KSEXP); Rcpp::traits::input_parameter< const arma::vec >::type naflag(naflagSEXP); Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); - rcpp_result_gen = Rcpp::wrap(nll_nmixTTD(beta, y, delta, W, V, pinds, mixture, tdist, N, J, K, naflag, threads)); + rcpp_result_gen = Rcpp::wrap(nll_nmixTTD(beta, y, delta, X_state, X_det, pinds, mixture, tdist, N, J, K, naflag, threads)); return rcpp_result_gen; END_RCPP } @@ -471,18 +471,18 @@ BEGIN_RCPP END_RCPP } // nll_occuTTD -double nll_occuTTD(arma::vec beta, arma::vec y, arma::vec delta, arma::mat W, arma::mat V, arma::mat Xgam, arma::mat Xeps, arma::vec pind, arma::vec dind, arma::vec cind, arma::vec eind, std::string lpsi, std::string tdist, int N, int T, int J, arma::vec naflag); -RcppExport SEXP _unmarked_nll_occuTTD(SEXP betaSEXP, SEXP ySEXP, SEXP deltaSEXP, SEXP WSEXP, SEXP VSEXP, SEXP XgamSEXP, SEXP XepsSEXP, SEXP pindSEXP, SEXP dindSEXP, SEXP cindSEXP, SEXP eindSEXP, SEXP lpsiSEXP, SEXP tdistSEXP, SEXP NSEXP, SEXP TSEXP, SEXP JSEXP, SEXP naflagSEXP) { +double nll_occuTTD(arma::vec beta, arma::vec y, arma::vec delta, arma::mat X_psi, arma::mat X_det, arma::mat X_col, arma::mat X_ext, arma::vec pind, arma::vec dind, arma::vec cind, arma::vec eind, std::string lpsi, std::string tdist, int N, int T, int J, arma::vec naflag); +RcppExport SEXP _unmarked_nll_occuTTD(SEXP betaSEXP, SEXP ySEXP, SEXP deltaSEXP, SEXP X_psiSEXP, SEXP X_detSEXP, SEXP X_colSEXP, SEXP X_extSEXP, SEXP pindSEXP, SEXP dindSEXP, SEXP cindSEXP, SEXP eindSEXP, SEXP lpsiSEXP, SEXP tdistSEXP, SEXP NSEXP, SEXP TSEXP, SEXP JSEXP, SEXP naflagSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< arma::vec >::type beta(betaSEXP); Rcpp::traits::input_parameter< arma::vec >::type y(ySEXP); Rcpp::traits::input_parameter< arma::vec >::type delta(deltaSEXP); - Rcpp::traits::input_parameter< arma::mat >::type W(WSEXP); - Rcpp::traits::input_parameter< arma::mat >::type V(VSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xgam(XgamSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xeps(XepsSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_psi(X_psiSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_det(X_detSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_col(X_colSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_ext(X_extSEXP); Rcpp::traits::input_parameter< arma::vec >::type pind(pindSEXP); Rcpp::traits::input_parameter< arma::vec >::type dind(dindSEXP); Rcpp::traits::input_parameter< arma::vec >::type cind(cindSEXP); @@ -493,7 +493,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type T(TSEXP); Rcpp::traits::input_parameter< int >::type J(JSEXP); Rcpp::traits::input_parameter< arma::vec >::type naflag(naflagSEXP); - rcpp_result_gen = Rcpp::wrap(nll_occuTTD(beta, y, delta, W, V, Xgam, Xeps, pind, dind, cind, eind, lpsi, tdist, N, T, J, naflag)); + rcpp_result_gen = Rcpp::wrap(nll_occuTTD(beta, y, delta, X_psi, X_det, X_col, X_ext, pind, dind, cind, eind, lpsi, tdist, N, T, J, naflag)); return rcpp_result_gen; END_RCPP } @@ -519,28 +519,28 @@ BEGIN_RCPP END_RCPP } // nll_pcountOpen -double nll_pcountOpen(arma::imat ym, arma::mat Xlam, arma::mat Xgam, arma::mat Xom, arma::mat Xp, arma::mat Xiota, arma::colvec beta_lam, arma::colvec beta_gam, arma::colvec beta_om, arma::colvec beta_p, arma::colvec beta_iota, double log_alpha, arma::colvec Xlam_offset, arma::colvec Xgam_offset, arma::colvec Xom_offset, arma::colvec Xp_offset, arma::colvec Xiota_offset, arma::imat ytna, arma::imat ynam, int lk, std::string mixture, Rcpp::IntegerVector first, Rcpp::IntegerVector last, int M, int J, int T, arma::imat delta, std::string dynamics, std::string fix, std::string go_dims, bool immigration, arma::imat I, arma::imat I1, Rcpp::List Ib, Rcpp::List Ip); -RcppExport SEXP _unmarked_nll_pcountOpen(SEXP ymSEXP, SEXP XlamSEXP, SEXP XgamSEXP, SEXP XomSEXP, SEXP XpSEXP, SEXP XiotaSEXP, SEXP beta_lamSEXP, SEXP beta_gamSEXP, SEXP beta_omSEXP, SEXP beta_pSEXP, SEXP beta_iotaSEXP, SEXP log_alphaSEXP, SEXP Xlam_offsetSEXP, SEXP Xgam_offsetSEXP, SEXP Xom_offsetSEXP, SEXP Xp_offsetSEXP, SEXP Xiota_offsetSEXP, SEXP ytnaSEXP, SEXP ynamSEXP, SEXP lkSEXP, SEXP mixtureSEXP, SEXP firstSEXP, SEXP lastSEXP, SEXP MSEXP, SEXP JSEXP, SEXP TSEXP, SEXP deltaSEXP, SEXP dynamicsSEXP, SEXP fixSEXP, SEXP go_dimsSEXP, SEXP immigrationSEXP, SEXP ISEXP, SEXP I1SEXP, SEXP IbSEXP, SEXP IpSEXP) { +double nll_pcountOpen(arma::imat ym, arma::mat X_lambda, arma::mat X_gamma, arma::mat X_omega, arma::mat X_det, arma::mat X_iota, arma::colvec beta_lam, arma::colvec beta_gam, arma::colvec beta_om, arma::colvec beta_p, arma::colvec beta_iota, double log_alpha, arma::colvec offset_lambda, arma::colvec offset_gamma, arma::colvec offset_omega, arma::colvec offset_det, arma::colvec offset_iota, arma::imat ytna, arma::imat ynam, int lk, std::string mixture, Rcpp::IntegerVector first, Rcpp::IntegerVector last, int M, int J, int T, arma::imat delta, std::string dynamics, std::string fix, std::string go_dims, bool immigration, arma::imat I, arma::imat I1, Rcpp::List Ib, Rcpp::List Ip); +RcppExport SEXP _unmarked_nll_pcountOpen(SEXP ymSEXP, SEXP X_lambdaSEXP, SEXP X_gammaSEXP, SEXP X_omegaSEXP, SEXP X_detSEXP, SEXP X_iotaSEXP, SEXP beta_lamSEXP, SEXP beta_gamSEXP, SEXP beta_omSEXP, SEXP beta_pSEXP, SEXP beta_iotaSEXP, SEXP log_alphaSEXP, SEXP offset_lambdaSEXP, SEXP offset_gammaSEXP, SEXP offset_omegaSEXP, SEXP offset_detSEXP, SEXP offset_iotaSEXP, SEXP ytnaSEXP, SEXP ynamSEXP, SEXP lkSEXP, SEXP mixtureSEXP, SEXP firstSEXP, SEXP lastSEXP, SEXP MSEXP, SEXP JSEXP, SEXP TSEXP, SEXP deltaSEXP, SEXP dynamicsSEXP, SEXP fixSEXP, SEXP go_dimsSEXP, SEXP immigrationSEXP, SEXP ISEXP, SEXP I1SEXP, SEXP IbSEXP, SEXP IpSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< arma::imat >::type ym(ymSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xlam(XlamSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xgam(XgamSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xom(XomSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xp(XpSEXP); - Rcpp::traits::input_parameter< arma::mat >::type Xiota(XiotaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_lambda(X_lambdaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_gamma(X_gammaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_omega(X_omegaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_det(X_detSEXP); + Rcpp::traits::input_parameter< arma::mat >::type X_iota(X_iotaSEXP); Rcpp::traits::input_parameter< arma::colvec >::type beta_lam(beta_lamSEXP); Rcpp::traits::input_parameter< arma::colvec >::type beta_gam(beta_gamSEXP); Rcpp::traits::input_parameter< arma::colvec >::type beta_om(beta_omSEXP); Rcpp::traits::input_parameter< arma::colvec >::type beta_p(beta_pSEXP); Rcpp::traits::input_parameter< arma::colvec >::type beta_iota(beta_iotaSEXP); Rcpp::traits::input_parameter< double >::type log_alpha(log_alphaSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xlam_offset(Xlam_offsetSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xgam_offset(Xgam_offsetSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xom_offset(Xom_offsetSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xp_offset(Xp_offsetSEXP); - Rcpp::traits::input_parameter< arma::colvec >::type Xiota_offset(Xiota_offsetSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_lambda(offset_lambdaSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_gamma(offset_gammaSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_omega(offset_omegaSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_det(offset_detSEXP); + Rcpp::traits::input_parameter< arma::colvec >::type offset_iota(offset_iotaSEXP); Rcpp::traits::input_parameter< arma::imat >::type ytna(ytnaSEXP); Rcpp::traits::input_parameter< arma::imat >::type ynam(ynamSEXP); Rcpp::traits::input_parameter< int >::type lk(lkSEXP); @@ -559,7 +559,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::imat >::type I1(I1SEXP); Rcpp::traits::input_parameter< Rcpp::List >::type Ib(IbSEXP); Rcpp::traits::input_parameter< Rcpp::List >::type Ip(IpSEXP); - rcpp_result_gen = Rcpp::wrap(nll_pcountOpen(ym, Xlam, Xgam, Xom, Xp, Xiota, beta_lam, beta_gam, beta_om, beta_p, beta_iota, log_alpha, Xlam_offset, Xgam_offset, Xom_offset, Xp_offset, Xiota_offset, ytna, ynam, lk, mixture, first, last, M, J, T, delta, dynamics, fix, go_dims, immigration, I, I1, Ib, Ip)); + rcpp_result_gen = Rcpp::wrap(nll_pcountOpen(ym, X_lambda, X_gamma, X_omega, X_det, X_iota, beta_lam, beta_gam, beta_om, beta_p, beta_iota, log_alpha, offset_lambda, offset_gamma, offset_omega, offset_det, offset_iota, ytna, ynam, lk, mixture, first, last, M, J, T, delta, dynamics, fix, go_dims, immigration, I, I1, Ib, Ip)); return rcpp_result_gen; END_RCPP } diff --git a/src/TMB/tmb_colext.hpp b/src/TMB/tmb_colext.hpp index 666f697e..77d35255 100644 --- a/src/TMB/tmb_colext.hpp +++ b/src/TMB/tmb_colext.hpp @@ -7,8 +7,8 @@ Type tmb_colext(objective_function* obj) { DATA_VECTOR(y); DATA_MATRIX(X_psi); - DATA_MATRIX(X_gam); - DATA_MATRIX(X_eps); + DATA_MATRIX(X_col); + DATA_MATRIX(X_ext); DATA_MATRIX(X_det); DATA_INTEGER(M); DATA_INTEGER(T); @@ -17,8 +17,8 @@ Type tmb_colext(objective_function* obj) { DATA_MATRIX(nd); PARAMETER_VECTOR(beta_psi); - PARAMETER_VECTOR(beta_gam); - PARAMETER_VECTOR(beta_eps); + PARAMETER_VECTOR(beta_col); + PARAMETER_VECTOR(beta_ext); PARAMETER_VECTOR(beta_det); Type loglik = 0.0; @@ -26,11 +26,11 @@ Type tmb_colext(objective_function* obj) { vector psi = X_psi * beta_psi; psi = invlogit(psi); - vector gam = X_gam * beta_gam; - gam = invlogit(gam); + vector col = X_col * beta_col; + col = invlogit(col); - vector eps = X_eps * beta_eps; - eps = invlogit(eps); + vector ext = X_ext * beta_ext; + ext = invlogit(ext); vector p = X_det * beta_det; p = invlogit(p); @@ -54,10 +54,10 @@ Type tmb_colext(objective_function* obj) { for (int t=0; t<(T-1); t++){ phi = phi.setZero(); - phi(0, 0) = 1 - eps(Tidx); - phi(0, 1) = eps(Tidx); - phi(1, 0) = gam(Tidx); - phi(1, 1) = 1 - gam(Tidx); + phi(0, 0) = 1 - ext(Tidx); + phi(0, 1) = ext(Tidx); + phi(1, 0) = col(Tidx); + phi(1, 1) = 1 - col(Tidx); Tidx += 1; if(site_sampled(i, t) == 1){ diff --git a/src/TMB/tmb_goccu.hpp b/src/TMB/tmb_goccu.hpp index 6007dc05..546aaad1 100644 --- a/src/TMB/tmb_goccu.hpp +++ b/src/TMB/tmb_goccu.hpp @@ -13,9 +13,9 @@ Type tmb_goccu(objective_function* obj) { DATA_INTEGER(link); - DATA_MATRIX(Xpsi); - DATA_MATRIX(Xphi); - DATA_MATRIX(Xp); + DATA_MATRIX(X_psi); + DATA_MATRIX(X_phi); + DATA_MATRIX(X_det); DATA_INTEGER(n_possible); DATA_MATRIX(alpha_potential); @@ -26,11 +26,11 @@ Type tmb_goccu(objective_function* obj) { PARAMETER_VECTOR(beta_psi); PARAMETER_VECTOR(beta_phi); - PARAMETER_VECTOR(beta_p); + PARAMETER_VECTOR(beta_det); - vector psi = Xpsi * beta_psi; - vector phi = Xphi * beta_phi; - vector p = Xp * beta_p; + vector psi = X_psi * beta_psi; + vector phi = X_phi * beta_phi; + vector p = X_det * beta_det; if(link == 1){ //psi = cloglog(psi); diff --git a/src/nll_distsampOpen.cpp b/src/nll_distsampOpen.cpp index 0f289af9..4a718bba 100644 --- a/src/nll_distsampOpen.cpp +++ b/src/nll_distsampOpen.cpp @@ -9,10 +9,10 @@ using namespace arma; // [[Rcpp::export]] double nll_distsampOpen(arma::ucube y, arma::imat yt, - arma::mat Xlam, arma::mat Xgam, arma::mat Xom, arma::mat Xsig, arma::mat Xiota, + arma::mat X_lambda, arma::mat X_gamma, arma::mat X_omega, arma::mat X_det, arma::mat X_iota, arma::vec beta, arma::umat bi, - arma::colvec Xlam_offset, arma::colvec Xgam_offset, arma::colvec Xom_offset, - arma::colvec Xsig_offset, arma::colvec Xiota_offset, + arma::colvec offset_lambda, arma::colvec offset_gamma, arma::colvec offset_omega, + arma::colvec offset_det, arma::colvec offset_iota, arma::imat ytna, int lk, std::string mixture, Rcpp::IntegerVector first, Rcpp::IntegerVector last, int first1, int M, int T, arma::imat delta, std::string dynamics, std::string survey, @@ -28,7 +28,7 @@ double nll_distsampOpen(arma::ucube y, arma::imat yt, //Lambda vec beta_lam = beta.subvec(bi(0,0), bi(0,1)); - vec lam = exp(Xlam*beta_lam + Xlam_offset) % A; + vec lam = exp(X_lambda*beta_lam + offset_lambda) % A; //Get 2nd abundance dist parameter set up if necessary double alpha = 0.0; @@ -44,10 +44,10 @@ double nll_distsampOpen(arma::ucube y, arma::imat yt, if((fix != "omega") && (dynamics != "trend")) { vec beta_om = beta.subvec(bi(2,0), bi(2,1)); if((dynamics == "ricker") || (dynamics == "gompertz")){ - omv = exp(Xom*beta_om + Xom_offset); + omv = exp(X_omega*beta_om + offset_omega); } else if((dynamics == "constant") || (dynamics == "autoreg") || (dynamics == "notrend")){ - omv = 1.0/(1.0+exp(-1*(Xom*beta_om + Xom_offset))); + omv = 1.0/(1.0+exp(-1*(X_omega*beta_om + offset_omega))); } } omv.reshape(T-1, M); @@ -60,7 +60,7 @@ double nll_distsampOpen(arma::ucube y, arma::imat yt, gam = (1-om) % lamMat; } else if (fix != "gamma") { vec beta_gam = beta.subvec(bi(1,0), bi(1,1)); - mat gamv = exp(Xgam*beta_gam + Xgam_offset); + mat gamv = exp(X_gamma*beta_gam + offset_gamma); gamv.reshape(T-1, M); gam = trans(gamv); } @@ -69,7 +69,7 @@ double nll_distsampOpen(arma::ucube y, arma::imat yt, mat sig(M, T); if(keyfun != "uniform"){ vec beta_sig = beta.subvec(bi(3,0), bi(3,1)); - mat sigv = exp(Xsig*beta_sig + Xsig_offset); + mat sigv = exp(X_det*beta_sig + offset_det); sigv.reshape(T, M); sig = trans(sigv); } @@ -84,7 +84,7 @@ double nll_distsampOpen(arma::ucube y, arma::imat yt, mat iota = zeros(M,T-1); if(immigration){ vec beta_iota = beta.subvec(bi(4,0), bi(4,1)); - mat iotav = exp(Xiota*beta_iota + Xiota_offset); + mat iotav = exp(X_iota*beta_iota + offset_iota); iotav.reshape(T-1, M); iota = trans(iotav); } diff --git a/src/nll_gdistremoval.cpp b/src/nll_gdistremoval.cpp index 8a04471f..9067bdab 100644 --- a/src/nll_gdistremoval.cpp +++ b/src/nll_gdistremoval.cpp @@ -15,41 +15,41 @@ using namespace arma; // [[Rcpp::export]] double nll_gdistremoval(arma::vec beta, arma::uvec n_param, arma::vec yDistance, arma::vec yRemoval, arma::mat ysum, int mixture, std::string keyfun, - arma::mat Xlam, arma::vec A, arma::mat Xphi, arma::mat Xrem, - arma::mat Xdist, arma::vec db, arma::mat a, arma::mat u, arma::vec w, + arma::mat X_lambda, arma::vec A, arma::mat X_phi, arma::mat X_rem, + arma::mat X_dist, arma::vec db, arma::mat a, arma::mat u, arma::vec w, arma::uvec pl, int K, arma::uvec Kmin, int threads){ #ifdef _OPENMP omp_set_num_threads(threads); #endif - int M = Xlam.n_rows; - int T = Xphi.n_rows / M; + int M = X_lambda.n_rows; + int T = X_phi.n_rows / M; int Rdist = yDistance.size() / M; unsigned Jdist = Rdist / T; int Rrem = yRemoval.size() / M; unsigned Jrem = Rrem / T; //Abundance - const vec lambda = exp(Xlam * beta_sub(beta, n_param, 0)) % A; + const vec lambda = exp(X_lambda * beta_sub(beta, n_param, 0)) % A; double log_alpha = beta_sub(beta, n_param, 1)(0); //length 1 vector //Availability vec phi = ones(M*T); if(T > 1){ - phi = inv_logit(Xphi * beta_sub(beta, n_param, 2)); + phi = inv_logit(X_phi * beta_sub(beta, n_param, 2)); } //Distance sampling detection vec dist_param(M*T); if(keyfun != "uniform"){ - dist_param = exp(Xdist * beta_sub(beta, n_param, 3)); + dist_param = exp(X_dist * beta_sub(beta, n_param, 3)); } double scale = exp(beta_sub(beta, n_param, 4)(0)); //Removal sampling detection vec remP(M * T * Jrem); - remP = inv_logit(Xrem * beta_sub(beta, n_param, 5)); + remP = inv_logit(X_rem * beta_sub(beta, n_param, 5)); double loglik = 0.0; diff --git a/src/nll_gdistsamp.cpp b/src/nll_gdistsamp.cpp index 2ae91331..9604964c 100644 --- a/src/nll_gdistsamp.cpp +++ b/src/nll_gdistsamp.cpp @@ -14,8 +14,8 @@ using namespace arma; // [[Rcpp::export]] double nll_gdistsamp(arma::vec beta, arma::uvec n_param, arma::vec y, int mixture, std::string keyfun, std::string survey, - arma::mat Xlam, arma::vec Xlam_offset, arma::vec A, arma::mat Xphi, - arma::vec Xphi_offset, arma::mat Xdet, arma::vec Xdet_offset, arma::vec db, + arma::mat X_lambda, arma::vec offset_lambda, arma::vec A, arma::mat X_phi, + arma::vec offset_phi, arma::mat X_det, arma::vec offset_det, arma::vec db, arma::mat a, arma::mat u, arma::vec w, arma::vec k, arma::vec lfac_k, arma::vec lfac_kmyt, arma::vec kmyt, arma::uvec Kmin, int threads){ @@ -23,27 +23,27 @@ double nll_gdistsamp(arma::vec beta, arma::uvec n_param, arma::vec y, omp_set_num_threads(threads); #endif - int M = Xlam.n_rows; - int T = Xphi.n_rows / M; + int M = X_lambda.n_rows; + int T = X_phi.n_rows / M; int R = y.size() / M; unsigned J = R / T; int lk = k.size(); int K = lk - 1; //Abundance - const vec lambda = exp(Xlam * beta_sub(beta, n_param, 0) + Xlam_offset) % A; + const vec lambda = exp(X_lambda * beta_sub(beta, n_param, 0) + offset_lambda) % A; double log_alpha = beta_sub(beta, n_param, 4)(0); //length 1 vector //Availability vec phi = ones(M*T); if(T > 1){ - phi = inv_logit(Xphi * beta_sub(beta, n_param, 1) + Xphi_offset); + phi = inv_logit(X_phi * beta_sub(beta, n_param, 1) + offset_phi); } //Detection vec det_param(M*T); if(keyfun != "uniform"){ - det_param = exp(Xdet * beta_sub(beta, n_param, 2) + Xdet_offset); + det_param = exp(X_det * beta_sub(beta, n_param, 2) + offset_det); } double scale = exp(beta_sub(beta, n_param, 3)(0)); diff --git a/src/nll_gmultmix.cpp b/src/nll_gmultmix.cpp index f4e90d55..55fb2e25 100644 --- a/src/nll_gmultmix.cpp +++ b/src/nll_gmultmix.cpp @@ -12,9 +12,9 @@ using namespace arma; // [[Rcpp::export]] double nll_gmultmix(arma::vec beta, arma::uvec n_param, arma::vec y, - int mixture, std::string pi_fun, arma::mat Xlam, - arma::vec Xlam_offset, arma::mat Xphi, arma::vec Xphi_offset, - arma::mat Xdet, arma::vec Xdet_offset, + int mixture, std::string pi_fun, arma::mat X_lambda, + arma::vec offset_lambda, arma::mat X_phi, arma::vec offset_phi, + arma::mat X_det, arma::vec offset_det, arma::vec k, arma::vec lfac_k, arma::vec lfac_kmyt, arma::vec kmyt, arma::vec Kmin, int threads){ @@ -22,22 +22,22 @@ double nll_gmultmix(arma::vec beta, arma::uvec n_param, arma::vec y, omp_set_num_threads(threads); #endif - int M = Xlam.n_rows; - int T = Xphi.n_rows / M; - int J = Xdet.n_rows / (M * T); + int M = X_lambda.n_rows; + int T = X_phi.n_rows / M; + int J = X_det.n_rows / (M * T); int R = y.size() / (M * T); int lk = k.size(); int K = k.size() - 1; - vec lambda = exp(Xlam * beta_sub(beta, n_param, 0) + Xlam_offset); + vec lambda = exp(X_lambda * beta_sub(beta, n_param, 0) + offset_lambda); double log_alpha = beta_sub(beta, n_param, 3)(0); vec phi = ones(M*T); if(T > 1){ - phi = inv_logit(Xphi * beta_sub(beta, n_param, 1) + Xphi_offset); + phi = inv_logit(X_phi * beta_sub(beta, n_param, 1) + offset_phi); } - vec p = inv_logit(Xdet * beta_sub(beta, n_param, 2) + Xdet_offset); + vec p = inv_logit(X_det * beta_sub(beta, n_param, 2) + offset_det); double loglik = 0.0; #pragma omp parallel for reduction(+: loglik) if(threads > 1) diff --git a/src/nll_gpcount.cpp b/src/nll_gpcount.cpp index 62e621f3..474c8d43 100644 --- a/src/nll_gpcount.cpp +++ b/src/nll_gpcount.cpp @@ -9,9 +9,9 @@ using namespace Rcpp ; using namespace arma ; // [[Rcpp::export]] -double nll_gpcount(arma::mat ym, arma::mat Xlam, arma::mat Xphi, arma::mat Xp, - arma::vec beta_lam, arma::vec beta_phi, arma::vec beta_p, double log_alpha, - arma::vec Xlam_offset, arma::vec Xphi_offset, arma::vec Xp_offset, +double nll_gpcount(arma::mat ym, arma::mat X_lambda, arma::mat X_phi, arma::mat X_det, + arma::vec beta_lambda, arma::vec beta_phi, arma::vec beta_det, double log_alpha, + arma::vec offset_lambda, arma::vec offset_phi, arma::vec offset_det, int M, std::string mixture, int T, int threads){ #ifdef _OPENMP @@ -27,10 +27,10 @@ double nll_gpcount(arma::mat ym, arma::mat Xlam, arma::mat Xphi, arma::mat Xp, alpha = 1.0/(1.0+exp(-log_alpha)); int R = ym.n_rows; int J = ym.n_cols / T; - arma::colvec lam = exp(Xlam*beta_lam + Xlam_offset); - arma::colvec logit_phi = Xphi*beta_phi + Xphi_offset; + arma::colvec lam = exp(X_lambda*beta_lambda + offset_lambda); + arma::colvec logit_phi = X_phi*beta_phi + offset_phi; arma::mat phiv = 1.0/(1.0+exp(-logit_phi)); - arma::colvec logit_p = Xp*beta_p + Xp_offset; + arma::colvec logit_p = X_det*beta_det + offset_det; arma::mat pv = 1.0/(1.0+exp(-logit_p)); phiv.reshape(T,R); arma::mat phi = trans(phiv); diff --git a/src/nll_multmixOpen.cpp b/src/nll_multmixOpen.cpp index 6fb755aa..250ae5f8 100644 --- a/src/nll_multmixOpen.cpp +++ b/src/nll_multmixOpen.cpp @@ -10,10 +10,10 @@ using namespace arma; // [[Rcpp::export]] double nll_multmixOpen(arma::ucube y, arma::imat yt, - arma::mat Xlam, arma::mat Xgam, arma::mat Xom, arma::mat Xp, arma::mat Xiota, + arma::mat X_lambda, arma::mat X_gamma, arma::mat X_omega, arma::mat X_det, arma::mat X_iota, arma::vec beta, arma::umat bi, - arma::colvec Xlam_offset, arma::colvec Xgam_offset, arma::colvec Xom_offset, - arma::colvec Xp_offset, arma::colvec Xiota_offset, + arma::colvec offset_lambda, arma::colvec offset_gamma, arma::colvec offset_omega, + arma::colvec offset_det, arma::colvec offset_iota, arma::imat ytna, arma::ucube yna, int lk, std::string mixture, Rcpp::IntegerVector first, Rcpp::IntegerVector last, int first1, int M, int T, int J, int R, arma::imat delta, std::string dynamics, @@ -28,7 +28,7 @@ double nll_multmixOpen(arma::ucube y, arma::imat yt, //Lambda vec beta_lam = beta.subvec(bi(0,0), bi(0,1)); - vec lam = exp(Xlam*beta_lam + Xlam_offset); + vec lam = exp(X_lambda*beta_lam + offset_lambda); //Get 2nd abundance dist parameter set up if necessary double alpha = 0.0; @@ -44,10 +44,10 @@ double nll_multmixOpen(arma::ucube y, arma::imat yt, if((fix != "omega") && (dynamics != "trend")) { vec beta_om = beta.subvec(bi(2,0), bi(2,1)); if((dynamics == "ricker") || (dynamics == "gompertz")){ - omv = exp(Xom*beta_om + Xom_offset); + omv = exp(X_omega*beta_om + offset_omega); } else if((dynamics == "constant") || (dynamics == "autoreg") || (dynamics == "notrend")){ - omv = 1.0/(1.0+exp(-1*(Xom*beta_om + Xom_offset))); + omv = 1.0/(1.0+exp(-1*(X_omega*beta_om + offset_omega))); } } omv.reshape(T-1, M); @@ -60,14 +60,14 @@ double nll_multmixOpen(arma::ucube y, arma::imat yt, gam = (1-om) % lamMat; } else if (fix != "gamma") { vec beta_gam = beta.subvec(bi(1,0), bi(1,1)); - mat gamv = exp(Xgam*beta_gam + Xgam_offset); + mat gamv = exp(X_gamma*beta_gam + offset_gamma); gamv.reshape(T-1, M); gam = trans(gamv); } //Detection prob vec beta_p = beta.subvec(bi(3,0), bi(3,1)); - vec pv = inv_logit(Xp * beta_p + Xp_offset); + vec pv = inv_logit(X_det * beta_p + offset_det); //Transform into cube (J x M x T) cube p((const double*) pv.begin(), R, M, T); @@ -75,7 +75,7 @@ double nll_multmixOpen(arma::ucube y, arma::imat yt, mat iota = zeros(M,T-1); if(immigration){ vec beta_iota = beta.subvec(bi(4,0), bi(4,1)); - mat iotav = exp(Xiota*beta_iota + Xiota_offset); + mat iotav = exp(X_iota*beta_iota + offset_iota); iotav.reshape(T-1, M); iota = trans(iotav); } diff --git a/src/nll_nmixTTD.cpp b/src/nll_nmixTTD.cpp index 50aab532..73f75555 100644 --- a/src/nll_nmixTTD.cpp +++ b/src/nll_nmixTTD.cpp @@ -9,7 +9,7 @@ using namespace arma; // [[Rcpp::export]] double nll_nmixTTD(const arma::vec beta, const arma::vec y, const arma::vec delta, - const arma::mat W, const arma::mat V, const arma::umat pinds, + const arma::mat X_state, const arma::mat X_det, const arma::umat pinds, const std::string mixture, const std::string tdist, int N, int J, int K, const arma::vec naflag, int threads){ @@ -18,10 +18,10 @@ double nll_nmixTTD(const arma::vec beta, const arma::vec y, const arma::vec delt #endif //Get abundance lambda values - const vec lamN = exp(W * beta.subvec(pinds(0,0), pinds(0,1))); + const vec lamN = exp(X_state * beta.subvec(pinds(0,0), pinds(0,1))); //Get detection lambda values - const vec lamP = exp(V * beta.subvec(pinds(1,0), pinds(1,1))); + const vec lamP = exp(X_det * beta.subvec(pinds(1,0), pinds(1,1))); //Get alpha if negative binomial double alpha = 1.0; diff --git a/src/nll_occuTTD.cpp b/src/nll_occuTTD.cpp index 8005f2a2..9bd86baf 100644 --- a/src/nll_occuTTD.cpp +++ b/src/nll_occuTTD.cpp @@ -32,7 +32,7 @@ using namespace arma; // [[Rcpp::export]] double nll_occuTTD( arma::vec beta, arma::vec y, arma::vec delta, - arma::mat W, arma::mat V, arma::mat Xgam, arma::mat Xeps, + arma::mat X_psi, arma::mat X_det, arma::mat X_col, arma::mat X_ext, arma::vec pind, arma::vec dind, arma::vec cind, arma::vec eind, std::string lpsi, std::string tdist, int N, int T, int J, @@ -41,7 +41,7 @@ double nll_occuTTD( arma::vec beta, arma::vec y, arma::vec delta, int ys = y.size(); //Get psi values - colvec raw_psi = W * beta.subvec(pind(0), pind(1)); + colvec raw_psi = X_psi * beta.subvec(pind(0), pind(1)); if(lpsi == "cloglog"){ raw_psi = 1 - exp(-exp(raw_psi)); } else { @@ -50,7 +50,7 @@ double nll_occuTTD( arma::vec beta, arma::vec y, arma::vec delta, const mat psi = join_rows(1-raw_psi, raw_psi); //Get lambda values - const vec lam = exp(V * beta.subvec(dind(0), dind(1))); + const vec lam = exp(X_det * beta.subvec(dind(0), dind(1))); vec e_lamt(ys); if(tdist == "weibull"){ @@ -68,8 +68,8 @@ double nll_occuTTD( arma::vec beta, arma::vec y, arma::vec delta, mat phi_raw(N*(T-1), 4); if(T > 1){ - colvec col = Xgam * beta.subvec(cind(0), cind(1)); - colvec ext = Xeps * beta.subvec(eind(0), eind(1)); + colvec col = X_col * beta.subvec(cind(0), cind(1)); + colvec ext = X_ext * beta.subvec(eind(0), eind(1)); ext = 1 / (1 + exp(-ext)); col = 1 / (1 + exp(-col)); phi_raw.col(0) = 1-col; diff --git a/src/nll_pcountOpen.cpp b/src/nll_pcountOpen.cpp index 235976c2..eeb17476 100644 --- a/src/nll_pcountOpen.cpp +++ b/src/nll_pcountOpen.cpp @@ -8,11 +8,11 @@ using namespace Rcpp ; // [[Rcpp::export]] -double nll_pcountOpen(arma::imat ym, arma::mat Xlam, arma::mat Xgam, arma::mat Xom, - arma::mat Xp, arma::mat Xiota, arma::colvec beta_lam, arma::colvec beta_gam, +double nll_pcountOpen(arma::imat ym, arma::mat X_lambda, arma::mat X_gamma, arma::mat X_omega, + arma::mat X_det, arma::mat X_iota, arma::colvec beta_lam, arma::colvec beta_gam, arma::colvec beta_om, arma::colvec beta_p, arma::colvec beta_iota, - double log_alpha, arma::colvec Xlam_offset, arma::colvec Xgam_offset, - arma::colvec Xom_offset, arma::colvec Xp_offset, arma::colvec Xiota_offset, + double log_alpha, arma::colvec offset_lambda, arma::colvec offset_gamma, + arma::colvec offset_omega, arma::colvec offset_det, arma::colvec offset_iota, arma::imat ytna, arma::imat ynam, int lk, std::string mixture, Rcpp::IntegerVector first, Rcpp::IntegerVector last, int M, int J, int T, arma::imat delta, std::string dynamics, std::string fix, std::string go_dims, @@ -28,13 +28,13 @@ double nll_pcountOpen(arma::imat ym, arma::mat Xlam, arma::mat Xgam, arma::mat X psi = 1.0/(1.0+exp(-log_alpha)); // linear predictors - arma::colvec lam = exp(Xlam*beta_lam + Xlam_offset); + arma::colvec lam = exp(X_lambda*beta_lam + offset_lambda); arma::mat omv = arma::ones(M*(T-1)); if((fix != "omega") && (dynamics != "trend")) { if((dynamics == "ricker") || (dynamics == "gompertz")) - omv = exp(Xom*beta_om + Xom_offset); + omv = exp(X_omega*beta_om + offset_omega); else if((dynamics == "constant") || (dynamics == "autoreg") || (dynamics == "notrend")) - omv = 1.0/(1.0+exp(-1*(Xom*beta_om + Xom_offset))); + omv = 1.0/(1.0+exp(-1*(X_omega*beta_om + offset_omega))); } omv.reshape(T-1, M); arma::mat om = arma::trans(omv); @@ -44,18 +44,18 @@ double nll_pcountOpen(arma::imat ym, arma::mat Xlam, arma::mat Xgam, arma::mat X gam = (1-om) % lamMat; } else { if(fix != "gamma") { - arma::mat gamv = exp(Xgam*beta_gam + Xgam_offset); + arma::mat gamv = exp(X_gamma*beta_gam + offset_gamma); gamv.reshape(T-1, M); gam = arma::trans(gamv); } } - arma::mat pv = 1.0/(1.0+exp(-1*(Xp*beta_p + Xp_offset))); + arma::mat pv = 1.0/(1.0+exp(-1*(X_det*beta_p + offset_det))); pv.reshape(J*T, M); arma::mat pm = trans(pv); //Immigration arma::mat iotav = arma::zeros(M*(T-1)); if(immigration) { - iotav = exp(Xiota*beta_iota + Xiota_offset); + iotav = exp(X_iota*beta_iota + offset_iota); } iotav.reshape(T-1, M); arma::mat iota = arma::trans(iotav);