Initial commit of occuRNMulti function

6 files changed, 419 insertions, 1 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 87d8c4f..260434d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -46,6 +46,7 @@ Collate: 'classes.R' 'unmarkedEstimate.R' 'mapInfo.R' 'unmarkedFrame.R'
     'unmarkedCrossVal.R' 'piFun.R' 'vif.R' 'makePiFun.R' 'posteriorSamples.R'
     'nmixTTD.R'
     'gdistremoval.R'
+    'occuRNMulti.R'
     'plotEffects.R'
     'mixedModelTools.R'
     'power.R' 
diff --git a/NAMESPACE b/NAMESPACE
index 3234918..61da06b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -24,7 +24,7 @@ importFrom(pbapply, pbsapply, pblapply)
 export(occu, occuFP, occuRN, pcount, pcountOpen, multinomPois, distsamp,
        colext, gmultmix, gdistsamp, gpcount, occuPEN, occuPEN_CV, occuMulti,
        occuMS, computeMPLElambda, pcount.spHDS, occuTTD, distsampOpen,
-       multmixOpen, nmixTTD, gdistremoval, goccu)
+       multmixOpen, nmixTTD, gdistremoval, goccu, occuRNMulti)
 export(removalPiFun, doublePiFun)
 export(makeRemPiFun, makeCrPiFun, makeCrPiFunMb, makeCrPiFunMh)
 
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 79ab967..8b206e5 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -69,6 +69,10 @@ 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_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)
 }
diff --git a/R/occuRNMulti.R b/R/occuRNMulti.R
new file mode 100644
index 0000000..5653f28
--- /dev/null
+++ b/R/occuRNMulti.R
@@ -0,0 +1,173 @@
+setClass("unmarkedFrameOccuRNMulti",
+         contains="unmarkedFrameOccuMulti")
+
+setClass("unmarkedFitOccuRNMulti",
+         representation(
+            detformulas = "list",
+            stateformulas = "list"),
+         contains = "unmarkedFit")
+
+occuRNMulti <- function(detformulas, stateformulas, data, 
+                        K = rep(10, length(stateformulas)),
+                        se = TRUE, threads=1){
+ 
+  stopifnot(all(names(data@ylist) == names(stateformulas)))
+
+  S <- length(stateformulas)
+  dep <- create_dep_matrix(stateformulas)
+  
+  data <- as(data, "unmarkedFrameOccuRNMulti")
+  gd <- getDesign(data, detformulas, stateformulas)
+
+  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)
+  det_rng <- min(gd$det_ind):max(gd$det_ind)
+  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,
+                    gd$ylist, gd$state_dm, gd$det_dm, dep, K, Kmin, threads)
+  }
+
+  fm <- optim(starts, nll_C, method="BFGS", hessian=TRUE)
+  covMat <- invertHessian(fm, nP, se)
+
+  fmAIC <- 2 * fm$value + 2 * nP
+  #----------------------------------------------------------------------------
+
+  #Format output---------------------------------------------------------------
+  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))
+
+  umfit <- new("unmarkedFitOccuRNMulti", fitType = "occuRNMulti", call = match.call(),
+                detformulas = detformulas, stateformulas = stateformulas,
+                formula = ~1, data = data,
+                #sitesRemoved = designMats$removed.sites,
+                estimates = estimateList, AIC = fmAIC, opt = fm,
+                negLogLike = fm$value, nllFun = nll_C)
+
+  return(umfit)
+}
+
+setMethod("getDesign", "unmarkedFrameOccuRNMulti",
+  function(umf, detformulas, stateformulas){
+  M <- nrow(umf@ylist[[1]])
+  J <- ncol(umf@ylist[[1]])
+  
+  ylist <- umf@ylist
+
+  # Site covs
+  sc <- umf@siteCovs
+  if(is.null(sc)) sc <- data.frame(dummy_=rep(1, M))
+
+  mm <- lapply(stateformulas, function(s){
+    if(!is.list(s)){
+      return(model.matrix(s, sc))
+    }
+    lapply(s, function(x) model.matrix(x, sc))
+  })
+
+  # Obs covs
+  sc_long_ind <- rep(1:M, each=J)
+  sc_long <- sc[sc_long_ind,]
+  oc <- cbind(sc_long, umf@obsCovs)
+
+  vv <- lapply(detformulas, function(x){
+    model.matrix(x, oc)
+  })
+
+  # Indices
+  state_n <- unlist(lapply(mm, function(x){
+    if(is.list(x)){
+      lapply(x, function(y) ncol(y))
+    } else {
+      ncol(x)
+    }
+  }))
+  state_end <- cumsum(state_n)
+  state_start <- state_end - state_n + 1
+  state_ind <- cbind(state_start, state_end)
+
+  state_prefix <- rep(rownames(state_ind), state_n)
+  state_prefix <- gsub(".", ":", state_prefix, fixed=TRUE)
+  state_prefix <- paste0("[",state_prefix,"]")
+  
+  state_par <- unlist(lapply(mm, function(x){
+    if(is.matrix(x)) return(colnames(x))
+    lapply(x, function(y){
+      return(colnames(y))
+    })
+  }))
+  stopifnot(length(state_prefix) == length(state_par))
+  state_par <- paste(state_prefix, state_par)
+  
+  det_n <- sapply(vv, ncol)
+  det_end <- cumsum(det_n) + max(state_end)
+  det_start <- det_end - det_n + 1
+  det_ind <- cbind(det_end, det_start)
+
+  det_rng <- min(det_ind):max(det_ind)
+
+  det_prefix <- paste0("[",rep(rownames(det_ind), det_n),"]")
+  det_par <- unlist(lapply(vv, colnames))
+  stopifnot(length(det_prefix) == length(det_par))
+  det_par <- paste(det_prefix, det_par)
+  par_names <- c(state_par, det_par)
+
+  list(par_names = par_names, state_ind = state_ind, det_ind = det_ind,
+              ylist = ylist, state_dm = mm, det_dm = vv)
+})
+
+create_dep_matrix <- function(stateformulas){
+  S <- length(stateformulas)
+  dep <- matrix(0, S, S)
+  rownames(dep) <- colnames(dep) <- names(stateformulas)
+  allowed_species <- character(0)
+
+  # Species 1
+  stopifnot(is.call(stateformulas[[1]]))
+  allowed_species <- c(allowed_species, names(stateformulas)[1])
+
+# Species 2
+  if(is.list(stateformulas[[2]])){
+    dep_sp <- names(stateformulas[[2]])[2:length(stateformulas[[2]])]
+    stopifnot(all(dep_sp %in% allowed_species))
+    for (i in dep_sp){
+      dep[names(stateformulas)[2], dep_sp] <- 1
+    }
+  }
+  allowed_species <- c(allowed_species, names(stateformulas)[2])
+
+  # Species 3
+  if(S > 2){
+    if(is.list(stateformulas[[3]])){
+      dep_sp <- names(stateformulas[[3]])[2:length(stateformulas[[3]])]
+      stopifnot(all(dep_sp %in% allowed_species))
+      for (i in dep_sp){
+        dep[names(stateformulas)[3], dep_sp] <- 1
+      }
+    }
+    allowed_species <- c(allowed_species, names(stateformulas)[3])
+  }
+
+  dep
+}
+
+
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 6d89809..5d2eae3 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -453,6 +453,27 @@ BEGIN_RCPP
     return rcpp_result_gen;
 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) {
+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::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< Rcpp::List >::type ylist(ylistSEXP);
+    Rcpp::traits::input_parameter< Rcpp::List >::type Xlist(XlistSEXP);
+    Rcpp::traits::input_parameter< Rcpp::List >::type Vlist(VlistSEXP);
+    Rcpp::traits::input_parameter< arma::imat >::type dep(depSEXP);
+    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));
+    return rcpp_result_gen;
+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) {
@@ -568,6 +589,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_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
new file mode 100644
index 0000000..d405926
--- /dev/null
+++ b/src/nll_occuRNMulti.cpp
@@ -0,0 +1,218 @@
+#include <RcppArmadillo.h>
+#include "utils.h"
+
+#ifdef _OPENMP
+  #include <omp.h>
+#endif
+
+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);
+  double p;
+  double g = 0;
+  for (int j = 0; j<J; j++){
+    p = 1 - pow(1 - r(j), k);
+    g += R::dbinom(y(j), 1, p, 1);
+  }
+  return f * exp(g);
+}
+
+//Calculate abundance likelihood for subordinate species
+double lik_subord_abun(int Kmin, int Kmax, double lam, int J, arma::vec r,
+    arma::rowvec y){
+
+  double lik = 0.0;
+
+  for (int k = Kmin; k < (Kmax+1); k++){
+    lik += calc_fg(k, lam, J, r, y);
+  }
+
+  return lik;
+}
+
+// [[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){
+
+  #ifdef _OPENMP
+    omp_set_num_threads(threads);
+  #endif
+
+  double nll = 0;
+
+  mat y1 = as<mat>(ylist[0]);
+  unsigned M = y1.n_rows;
+  unsigned J = y1.n_cols;
+
+  //Parameters for dominant species
+  mat X1 = as<mat>(Xlist[0]);
+  vec lam1 = exp(X1 * beta.subvec(state_ind(0, 0), state_ind(0, 1)));
+
+  mat V1 = as<mat>(Vlist[0]);
+  vec r1 = inv_logit(V1 * beta.subvec(det_ind(0, 0), det_ind(0,1)));
+
+  //Parameters for 2nd species
+  mat y2 = as<mat>(ylist[1]);
+
+  int sind_row = 1;
+
+  vec lam2;
+  vec lp2, lp2_sp1;
+  if(dep(1,0)){
+    List Xlist_sp2 = Xlist[1];
+    mat X2 = as<mat>(Xlist_sp2[0]);
+    lp2 = X2 * beta.subvec(state_ind(1,0), state_ind(1,1));
+    mat X2_sp1 = as<mat>(Xlist_sp2[1]);
+    lp2_sp1 = X2_sp1 * beta.subvec(state_ind(2,0), state_ind(2,1));
+    sind_row += 2;
+  } else {
+    mat X2 = as<mat>(Xlist[1]);
+    lam2 = exp(X2 * beta.subvec(state_ind(1,0), state_ind(1,1)));
+    sind_row += 1;
+  }
+
+  mat V2 = as<mat>(Vlist[1]);
+  vec r2 = inv_logit(V2 * beta.subvec(det_ind(1, 0), det_ind(1,1)));
+
+  //Parameters for 3rd species (if needed)
+  mat y3;
+  vec r3;
+  vec lam3, lp3, lp3_sp1, lp3_sp2;
+
+  //Very janky...
+  if(S > 2){
+    y3 = as<mat>(ylist[2]);
+
+    if(dep(2,0) | dep(2,1)){
+      List Xlist_sp3 = Xlist[2];
+      mat X3 = as<mat>(Xlist_sp3[0]);
+      lp3 = X3 * beta.subvec(state_ind(sind_row, 0), state_ind(sind_row, 1));
+      sind_row += 1;
+
+      int mat_idx = 1;
+      if(dep(2,0)){
+        mat X3_sp1 = as<mat>(Xlist_sp3[mat_idx]);
+        mat_idx += 1;
+        lp3_sp1 = X3_sp1 * beta.subvec(state_ind(sind_row, 0), state_ind(sind_row, 1));
+        sind_row += 1;
+      }
+      if(dep(2,1)){
+        mat X3_sp2 = as<mat>(Xlist_sp3[mat_idx]);
+        lp3_sp2 = X3_sp2 * beta.subvec(state_ind(sind_row, 0), state_ind(sind_row, 1));
+      }
+    } else {
+      mat X3 = as<mat>(Xlist[2]);
+      lam3 = exp(X3 * beta.subvec(state_ind(sind_row, 0), state_ind(sind_row, 1)));
+    }
+
+    mat V3 = as<mat>(Vlist[2]);
+    r3 = inv_logit(V3 * beta.subvec(det_ind(2, 0), det_ind(2,1)));
+  }
+
+  #pragma omp parallel for reduction(-: nll) if(threads > 1)
+  for (int m = 0; m<M; m++){
+    double fg1, fg2, lik_sp2, lik_sp3;
+
+    //Indices for detection probability
+    int r_st = m * J;
+    int r_end = r_st + J - 1;
+    vec r1_sub = r1.subvec(r_st, r_end);
+    vec r2_sub = r2.subvec(r_st, r_end);
+
+    vec r3_sub;
+    if(S==3){
+      r3_sub = r3.subvec(r_st, r_end);
+    }
+
+    double lik_m = 0.0;
+
+    double lam2_m, lam3_m;
+
+    for (int k1 = Kmin(m, 0); k1<(K(0)+1); k1++){
+
+      //Dominant species
+      fg1 = calc_fg(k1, lam1(m), J, r1_sub, y1.row(m));
+
+      //-----------------------------------------------------------------------
+      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));
+
+        lik_m += fg1 * lik_sp2;
+
+      //-----------------------------------------------------------------------
+      } else if((S == 3) & (dep(1,0) & dep(2,1))){ // sp1 --> sp2 --> sp3
+
+        lik_sp2 = 0;
+
+        lam2_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));
+
+          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));
+
+          lik_sp2 += fg2 * lik_sp3;
+        }
+
+        lik_m += fg1 * lik_sp2;
+
+      //-----------------------------------------------------------------------
+      } else if((S == 3) & (dep(2,0) & dep(2,1))){ // sp1 --> sp3 <-- sp2
+
+        lik_sp2 = 0;
+
+        for (int k2 = Kmin(m, 1); k2<(K(1)+1); k2++){
+
+          fg2 = calc_fg(k2, lam2(m), J, r2_sub, y2.row(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));
+
+          lik_sp2 += fg2 * lik_sp3;
+        }
+
+        lik_m += fg1 * lik_sp2;
+
+      //-----------------------------------------------------------------------
+      } else if((S == 3) & (dep(1,0) & dep(2,0))){ // sp2 <-- sp1 --> sp3
+
+        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));
+
+        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));
+
+        lik_m += fg1 * lik_sp2 * lik_sp3;
+      }
+    }
+
+
+    nll -= log(lik_m);
+  }
+  return nll;
+}