Add occupancy support, model sp1-->sp3<--sp2 doesnt work

4 files changed, 114 insertions, 52 deletions

diff --git a/R/RcppExports.R b/R/RcppExports.R
index 8b206e5..3d86707 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -69,8 +69,8 @@ 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_occuRNMulti <- function(beta, state_ind, det_ind, S, ylist, Xlist, Vlist, dep, K, Kmin, threads) {
-    .Call(`_unmarked_nll_occuRNMulti`, beta, state_ind, det_ind, S, ylist, Xlist, Vlist, dep, K, Kmin, threads)
+nll_occuRNMulti <- function(beta, state_ind, det_ind, S, modOcc, ylist, Xlist, Vlist, dep, K, Kmin, threads) {
+    .Call(`_unmarked_nll_occuRNMulti`, beta, state_ind, det_ind, S, modOcc, ylist, Xlist, Vlist, dep, K, Kmin, threads)
 }
 
 nll_occuTTD <- function(beta, y, delta, W, V, Xgam, Xeps, pind, dind, cind, eind, lpsi, tdist, N, T, J, naflag) {
diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R
index 5653f28..a8278c1 100644
--- a/R/occuRNMulti.R
+++ b/R/occuRNMulti.R
@@ -7,14 +7,31 @@ setClass("unmarkedFitOccuRNMulti",
             stateformulas = "list"),
          contains = "unmarkedFit")
 
-occuRNMulti <- function(detformulas, stateformulas, data, 
+occuRNMulti <- function(detformulas, stateformulas, data, modelOccupancy,
                         K = rep(10, length(stateformulas)),
-                        se = TRUE, threads=1){
+                        starts, method="BFGS", se = TRUE, threads=1, ...){
  
   stopifnot(all(names(data@ylist) == names(stateformulas)))
 
   S <- length(stateformulas)
   dep <- create_dep_matrix(stateformulas)
+
+  if(missing(modelOccupancy)){
+    modelOccupancy <- rep(0, S)
+    names(modelOccupancy) <- names(stateformulas)
+  } else {
+    stopifnot(identical(names(modelOccupancy), names(stateformulas)))
+    if(is.logical(modelOccupancy)){
+      modelOccupancy <- modelOccupancy*1
+    }
+    not_allowed_occ <- colSums(dep)
+    sp_mismatch <- modelOccupancy & not_allowed_occ
+    if(any(sp_mismatch)){
+      stop(paste("Species", 
+                 paste(names(modelOccupancy)[sp_mismatch], collapse=", "),
+                 "cannot use occupancy models."), call.=FALSE)
+    }
+  }
   
   data <- as(data, "unmarkedFrameOccuRNMulti")
   gd <- getDesign(data, detformulas, stateformulas)
@@ -22,17 +39,32 @@ occuRNMulti <- function(detformulas, stateformulas, data,
   Kmin <- sapply(gd$ylist, function(x) apply(x, 1, max))
 
   nP <- max(gd$det_ind)
-  state_rng <- min(gd$state_ind):max(gd$state_ind)
+  
+  has_occ <- sum(modelOccupancy) > 0
+  if(has_occ){
+    sp_occ <- names(modelOccupancy)[modelOccupancy == 1]
+    sp_abun <- names(modelOccupancy)[modelOccupancy == 0]
+    occ_ind <- gd$state_ind[grepl(paste(paste0("^",sp_occ), collapse="|"),
+                                  rownames(gd$state_ind)),,drop=FALSE]
+    occ_rng <- min(occ_ind):max(occ_ind)
+    abun_ind <- gd$state_ind[grepl(paste(paste0("^",sp_abun), collapse="|"),
+                                  rownames(gd$state_ind)),,drop=FALSE]
+    abun_rng <- min(abun_ind):max(abun_ind)
+  } else {
+    state_rng <- min(gd$state_ind):max(gd$state_ind)
+  }
+
   det_rng <- min(gd$det_ind):max(gd$det_ind)
-  starts <- rep(0, nP)
+  
+  if(missing(starts)) starts <- rep(0, nP)
   names(starts) <- gd$par_names
 
   nll_C <- function(pars){
-    nll_occuRNMulti(pars, gd$state_ind-1, gd$det_ind-1, S,
+    nll_occuRNMulti(pars, gd$state_ind-1, gd$det_ind-1, S, modelOccupancy,
                     gd$ylist, gd$state_dm, gd$det_dm, dep, K, Kmin, threads)
   }
 
-  fm <- optim(starts, nll_C, method="BFGS", hessian=TRUE)
+  fm <- optim(starts, nll_C, method=method, hessian=se, ...)
   covMat <- invertHessian(fm, nP, se)
 
   fmAIC <- 2 * fm$value + 2 * nP
@@ -42,19 +74,33 @@ occuRNMulti <- function(detformulas, stateformulas, data,
   ests <- fm$par
   names(ests) <- gd$par_names
 
-  state <- unmarkedEstimate(name = "Abundance", short.name = "lam",
-                            estimates = ests[state_rng],
-                            covMat = as.matrix(covMat[state_rng, state_rng]),
-                            invlink = "exp",
-                            invlinkGrad = "exp")
-
   det <- unmarkedEstimate(name = "Detection", short.name = "p",
                           estimates = ests[det_rng],
                           covMat = as.matrix(covMat[det_rng, det_rng]),
                           invlink = "logistic",
                           invlinkGrad = "logistic.grad")
 
-  estimateList <- unmarkedEstimateList(list(state=state, det=det))
+  if(has_occ){
+    lam <- unmarkedEstimate(name = "Abundance", short.name = "lam",
+                            estimates = ests[abun_rng],
+                            covMat = as.matrix(covMat[abun_rng, abun_rng]),
+                            invlink = "exp",
+                            invlinkGrad = "exp")
+    psi <- unmarkedEstimate(name = "Occupancy", short.name = "psi",
+                          estimates = ests[occ_rng],
+                          covMat = as.matrix(covMat[occ_rng, occ_rng]),
+                          invlink = "logistic",
+                          invlinkGrad = "logistic.grad")
+
+    estimateList <- unmarkedEstimateList(list(lam=lam, psi=psi, det=det))
+  } else {
+    lam <- unmarkedEstimate(name = "Abundance", short.name = "lam",
+                            estimates = ests[state_rng],
+                            covMat = as.matrix(covMat[state_rng, state_rng]),
+                            invlink = "exp",
+                            invlinkGrad = "exp")
+    estimateList <- unmarkedEstimateList(list(lam=lam, det=det))
+  }
 
   umfit <- new("unmarkedFitOccuRNMulti", fitType = "occuRNMulti", call = match.call(),
                 detformulas = detformulas, stateformulas = stateformulas,
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 5d2eae3..ec850b2 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -454,8 +454,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // nll_occuRNMulti
-double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, int S, Rcpp::List ylist, Rcpp::List Xlist, Rcpp::List Vlist, arma::imat dep, arma::ivec K, arma::imat Kmin, int threads);
-RcppExport SEXP _unmarked_nll_occuRNMulti(SEXP betaSEXP, SEXP state_indSEXP, SEXP det_indSEXP, SEXP SSEXP, SEXP ylistSEXP, SEXP XlistSEXP, SEXP VlistSEXP, SEXP depSEXP, SEXP KSEXP, SEXP KminSEXP, SEXP threadsSEXP) {
+double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, int S, arma::ivec modOcc, Rcpp::List ylist, Rcpp::List Xlist, Rcpp::List Vlist, arma::imat dep, arma::ivec K, arma::imat Kmin, int threads);
+RcppExport SEXP _unmarked_nll_occuRNMulti(SEXP betaSEXP, SEXP state_indSEXP, SEXP det_indSEXP, SEXP SSEXP, SEXP modOccSEXP, SEXP ylistSEXP, SEXP XlistSEXP, SEXP VlistSEXP, SEXP depSEXP, SEXP KSEXP, SEXP KminSEXP, SEXP threadsSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -463,6 +463,7 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< arma::mat >::type state_ind(state_indSEXP);
     Rcpp::traits::input_parameter< arma::mat >::type det_ind(det_indSEXP);
     Rcpp::traits::input_parameter< int >::type S(SSEXP);
+    Rcpp::traits::input_parameter< arma::ivec >::type modOcc(modOccSEXP);
     Rcpp::traits::input_parameter< Rcpp::List >::type ylist(ylistSEXP);
     Rcpp::traits::input_parameter< Rcpp::List >::type Xlist(XlistSEXP);
     Rcpp::traits::input_parameter< Rcpp::List >::type Vlist(VlistSEXP);
@@ -470,7 +471,7 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< arma::ivec >::type K(KSEXP);
     Rcpp::traits::input_parameter< arma::imat >::type Kmin(KminSEXP);
     Rcpp::traits::input_parameter< int >::type threads(threadsSEXP);
-    rcpp_result_gen = Rcpp::wrap(nll_occuRNMulti(beta, state_ind, det_ind, S, ylist, Xlist, Vlist, dep, K, Kmin, threads));
+    rcpp_result_gen = Rcpp::wrap(nll_occuRNMulti(beta, state_ind, det_ind, S, modOcc, ylist, Xlist, Vlist, dep, K, Kmin, threads));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -589,7 +590,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_unmarked_nll_occuMulti", (DL_FUNC) &_unmarked_nll_occuMulti, 15},
     {"_unmarked_nll_occuPEN", (DL_FUNC) &_unmarked_nll_occuPEN, 11},
     {"_unmarked_nll_occuRN", (DL_FUNC) &_unmarked_nll_occuRN, 10},
-    {"_unmarked_nll_occuRNMulti", (DL_FUNC) &_unmarked_nll_occuRNMulti, 11},
+    {"_unmarked_nll_occuRNMulti", (DL_FUNC) &_unmarked_nll_occuRNMulti, 12},
     {"_unmarked_nll_occuTTD", (DL_FUNC) &_unmarked_nll_occuTTD, 17},
     {"_unmarked_nll_pcount", (DL_FUNC) &_unmarked_nll_pcount, 11},
     {"_unmarked_nll_pcountOpen", (DL_FUNC) &_unmarked_nll_pcountOpen, 35},
diff --git a/src/nll_occuRNMulti.cpp b/src/nll_occuRNMulti.cpp
index d405926..5402a9a 100644
--- a/src/nll_occuRNMulti.cpp
+++ b/src/nll_occuRNMulti.cpp
@@ -8,10 +8,6 @@
 using namespace Rcpp;
 using namespace arma;
 
-//arma::mat inv_logit( arma::mat inp ){
-//  return(1 / (1 + exp(-1 * inp)));
-//}
-
 //Calculate part of abundance likelihood to integrate over (abun f[g] * detect [g])
 double calc_fg(int k, double lam, int J, arma::vec r, arma::rowvec y){
   double f = R::dpois(k, lam, 0);
@@ -29,7 +25,7 @@ double lik_subord_abun(int Kmin, int Kmax, double lam, int J, arma::vec r,
     arma::rowvec y){
 
   double lik = 0.0;
-
+  lam = exp(lam);
   for (int k = Kmin; k < (Kmax+1); k++){
     lik += calc_fg(k, lam, J, r, y);
   }
@@ -37,10 +33,20 @@ double lik_subord_abun(int Kmin, int Kmax, double lam, int J, arma::vec r,
   return lik;
 }
 
+//Calculate occupancy likelihood for subordinate species
+double lik_subord_occ(int Kmin, double psi, int J, arma::vec p, arma::rowvec y){
+  double g = 1;
+  psi = inv_logit(psi);
+  for (int j = 0; j<J; j++){
+    g *= R::dbinom(y(j), 1, p(j), 0);
+  }
+  return g * psi + (1 - Kmin) * (1 - psi);
+}
+
 // [[Rcpp::export]]
 double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, int S,
-    Rcpp::List ylist, Rcpp::List Xlist, Rcpp::List Vlist, arma::imat dep,
-    arma::ivec K, arma::imat Kmin, int threads){
+    arma::ivec modOcc, Rcpp::List ylist, Rcpp::List Xlist, Rcpp::List Vlist,
+    arma::imat dep, arma::ivec K, arma::imat Kmin, int threads){
 
   #ifdef _OPENMP
     omp_set_num_threads(threads);
@@ -109,6 +115,7 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i
         lp3_sp2 = X3_sp2 * beta.subvec(state_ind(sind_row, 0), state_ind(sind_row, 1));
       }
     } else {
+      // Unused at the moment I think
       mat X3 = as<mat>(Xlist[2]);
       lam3 = exp(X3 * beta.subvec(state_ind(sind_row, 0), state_ind(sind_row, 1)));
     }
@@ -134,7 +141,7 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i
 
     double lik_m = 0.0;
 
-    double lam2_m, lam3_m;
+    double par2_m, par3_m; // state parameters for species 2 and 3
 
     for (int k1 = Kmin(m, 0); k1<(K(0)+1); k1++){
 
@@ -145,11 +152,12 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i
       if((S == 2) & dep(1,0)){  //sp1 --> sp2
 
         lik_sp2 = 0;
-
-        //ALLOW OCCUPANCY HERE
-        lam2_m = exp(lp2(m) + k1 * lp2_sp1(m));
-        lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), lam2_m, J, r2_sub, y2.row(m));
-
+        par2_m = lp2(m) + k1 * lp2_sp1(m);
+        if(modOcc(1)){
+          lik_sp2 = lik_subord_occ(Kmin(m, 1), par2_m, J, r2_sub, y2.row(m));
+        } else {
+          lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), par2_m, J, r2_sub, y2.row(m));
+        }
         lik_m += fg1 * lik_sp2;
 
       //-----------------------------------------------------------------------
@@ -157,17 +165,18 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i
 
         lik_sp2 = 0;
 
-        lam2_m = exp(lp2(m) + k1 * lp2_sp1(m));
+        par2_m = exp(lp2(m) + k1 * lp2_sp1(m));
 
         for (int k2 = Kmin(m, 1); k2<(K(1)+1); k2++){
-          fg2 = calc_fg(k2, lam2_m, J, r2_sub, y2.row(m));
+          fg2 = calc_fg(k2, par2_m, J, r2_sub, y2.row(m));
 
           lik_sp3 = 0;
-
-          //ALLOW OCCUPANCY HERE
-          lam3_m = exp(lp3(m) + k2 * lp3_sp2(m));
-          lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), lam3_m, J, r3_sub, y3.row(m));
-
+          par3_m = lp3(m) + k2 * lp3_sp2(m);
+          if(modOcc(2)){
+            lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m));
+          } else {
+            lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m));
+          }
           lik_sp2 += fg2 * lik_sp3;
         }
 
@@ -182,12 +191,13 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i
 
           fg2 = calc_fg(k2, lam2(m), J, r2_sub, y2.row(m));
 
+          par3_m = lp3(m) + k1 * lp3_sp1(m) + k2 * lp3_sp2(m);
           lik_sp3 = 0;
-
-          //ALLOW OCCUPANCY HERE
-          lam3_m = exp(lp3(m) + k1 * lp3_sp1(m) + k2 * lp3_sp2(m));
-          lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), lam3_m, J, r3_sub, y3.row(m));
-
+          if(modOcc(1)){
+            lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m));
+          } else {
+            lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m));
+          }
           lik_sp2 += fg2 * lik_sp3;
         }
 
@@ -196,22 +206,27 @@ double nll_occuRNMulti(arma::vec beta, arma::mat state_ind, arma::mat det_ind, i
       //-----------------------------------------------------------------------
       } else if((S == 3) & (dep(1,0) & dep(2,0))){ // sp2 <-- sp1 --> sp3
 
+
+        par2_m = lp2(m) + k1 * lp2_sp1(m);
         lik_sp2 = 0;
-        //ALLOW OCCUPANCY HERE
-        lam2_m = exp(lp2(m) + k1 * lp2_sp1(m));
-        lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), lam2_m, J, r2_sub, y2.row(m));
+        if(modOcc(1)){
+          lik_sp2 = lik_subord_occ(Kmin(m, 1), par2_m, J, r2_sub, y2.row(m));
+        } else {
+          lik_sp2 = lik_subord_abun(Kmin(m, 1), K(1), par2_m, J, r2_sub, y2.row(m));
+        }
 
+        par3_m = lp3(m) + k1 * lp3_sp1(m);
         lik_sp3 = 0;
-
-        //ALLOW OCCUPANCY HERE
-        lam3_m = exp(lp3(m) + k1 * lp3_sp1(m));
-        lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), lam3_m, J, r3_sub, y3.row(m));
+        if(modOcc(2)){
+          lik_sp3 = lik_subord_occ(Kmin(m, 2), par3_m, J, r3_sub, y3.row(m));
+        } else {
+          lik_sp3 = lik_subord_abun(Kmin(m, 2), K(2), par3_m, J, r3_sub, y3.row(m));
+        }
 
         lik_m += fg1 * lik_sp2 * lik_sp3;
       }
     }
 
-
     nll -= log(lik_m);
   }
   return nll;