Fix incorrect likelihood calculation in goccu when primary periods we…

…re missing

R/goccu.R

@@ -82,11 +82,23 @@ goccu <- function(psiformula, phiformula, pformula, data,
     }
   }
 
+  # Determine which configurations of available states should not be
+  # included in likelihood because relevant primary periods were missing
+  alpha_drop <- matrix(NA, M, nrow(alpha_potential))
+  for (i in 1:M){
+    dropped <- rep(0, nrow(alpha_potential))
+    for (j in 1:nrow(alpha_potential)){
+      check_drop <- alpha_potential[j,] * missing_session[i,]
+      if(sum(check_drop) > 0) dropped[j] <- 1
+    }
+    alpha_drop[i,] <- dropped
+  }
+
   # Bundle data for TMB
   dataList <- list(y=y, T=T, link=ifelse(linkPsi=='cloglog', 1, 0), 
                    Xpsi=Xpsi, Xphi=Xphi, Xp=Xp,
                    n_possible=n_possible,
-                   alpha_potential=alpha_potential,
+                   alpha_potential=alpha_potential, alpha_drop = alpha_drop,
                    known_present=known_present, known_available=known_available, 
                    missing_session=missing_session)
 

src/TMB/tmb_goccu.hpp

@@ -19,6 +19,7 @@ Type tmb_goccu(objective_function<Type>* obj) {
 
   DATA_INTEGER(n_possible);
   DATA_MATRIX(alpha_potential);
+  DATA_MATRIX(alpha_drop);
   DATA_VECTOR(known_present);
   DATA_MATRIX(known_available);
   DATA_MATRIX(missing_session);
@@ -93,6 +94,12 @@ Type tmb_goccu(objective_function<Type>* obj) {
       exp_poss_lp = 0.0;
 
       for (int k=0; k<n_possible; k++){
+
+        //If this configuration of available states is not possible
+        //because some of the included primary periods were dropped,
+        //then don't add anything to the likelihood
+        if(alpha_drop(i, k) == 1) continue;
+
         poss_lp = log(psi(i));
 
         for (int t=0; t<T; t++){