NAs in occuCOP: fix the C++ likelihood function

1 files changed, 10 insertions, 5 deletions

diff --git a/src/nll_occuCOP.cpp b/src/nll_occuCOP.cpp
index abae16c..5ce374e 100644
--- a/src/nll_occuCOP.cpp
+++ b/src/nll_occuCOP.cpp
@@ -19,24 +19,29 @@ double nll_occuCOP(arma::icolvec y, arma::icolvec L,
   //Calculate lambda back-transformed from log
   arma::colvec lambda = exp(Xlambda*beta_lambda);
   
-
   double ll=0.0;
   int k=0; // counter
   // for each site i in 1:M
   for(int i=0; i<M; i++) {
     double iLambdaL=0.0; // init sum(lambda_ij * L_ij)
     double iN=0.0; // init sum(y) = total count of detec at site i
+    int NbRemoved=0; // init count of the removed observations at site i
     for(int j=0; j<J; j++) {
       if(!removed(k)) {
+        // If the observation is not removed from the analysis
+        // (removed if there is a NA in y, L or in the relevant covariates for this site and obs)
         iLambdaL += lambda(k)*L(k);
         iN += y(k);
+        NbRemoved += 1;
       }
       k++;
     }
-    if(iN>0) {
-      ll += log(psi(i) * pow(iLambdaL, iN) / tgamma(iN + 1) * exp(-iLambdaL));
-    } else {
-      ll += log(psi(i) * exp(-iLambdaL) + 1-psi(i));
+    if (!NbRemoved < J) {
+      if(iN>0) {
+        ll += log(psi(i) * pow(iLambdaL, iN) / tgamma(iN + 1) * exp(-iLambdaL));
+      } else {
+        ll += log(psi(i) * exp(-iLambdaL) + 1-psi(i));
+      }
     }
   }
   return -ll;