reorganized pre2-codes, added intercept only treatments for nhmms

helske · Nov 8, 2024 · 7a17606 · 7a17606
1 parent 5a56bff
commit 7a17606
Show file tree

Hide file tree

Showing 64 changed files with 2,994 additions and 3,075 deletions.
diff --git a/R/RcppExports.R b/R/RcppExports.R
diff --git a/R/check_build_arguments.R b/R/check_build_arguments.R
@@ -286,3 +286,30 @@
   data <- data[order(data[[id]], data[[time]]), ]
   fill_time(data, id, time)
 }
+#' Checks that the design matrix is of full rank
+#' @noRd
+.check_identifiability <- function(X, type) {
+  n <- nrow(X)
+  nc <- ncol(X)
+
+  # Check if matrix is full rank
+  if (!identical(qr(X)$rank, min(n, nc))) {
+    # Check for zero-only columns
+    zero_col <- apply(X, 2L, function(x) all(x == 0))
+
+    if (any(zero_col)) {
+      k <- sum(zero_col)
+      stop_(
+        "{cli::qty(k)} Predictor{?s} {.var {colnames(X)[zero_col]}}
+         {cli::qty(k)} contain{?s/} only zeros in the design matrix for {.var {type}_formula}."
+      )
+    } else {
+      stop_(
+        "Perfect collinearity found between predictor variables of
+        {.var {type}_formula}."
+      )
+    }
+    return(FALSE)
+  }
+  TRUE
+}
diff --git a/R/create_base_nhmm.R b/R/create_base_nhmm.R
@@ -67,12 +67,12 @@ create_base_nhmm <- function(observations, data, time, id, n_states,
     }
   }
 
-  icp_only_i <- intercept_only(initial_formula)
-  icp_only_s <- intercept_only(transition_formula)
-  icp_only_o <- intercept_only(emission_formula)
-  icp_only_d <- ifelse(mixture, intercept_only(cluster_formula), TRUE)
+  icpt_only_i <- intercept_only(initial_formula)
+  icpt_only_s <- intercept_only(transition_formula)
+  icpt_only_o <- intercept_only(emission_formula)
+  icpt_only_d <- ifelse(mixture, intercept_only(cluster_formula), TRUE)
   y_in_data <- checkmate::test_character(observations)
-  if (!icp_only_i || !icp_only_s || !icp_only_o || !icp_only_d || y_in_data) {
+  if (!icpt_only_i || !icpt_only_s || !icpt_only_o || !icpt_only_d || y_in_data) {
     data <- .check_data(data, time, id)
     ids <- unique(data[[id]])
     times <- unique(data[[time]])
@@ -174,7 +174,7 @@ create_base_nhmm <- function(observations, data, time, id, n_states,
       np_B = B$n_pars,
       np_omega = omega$n_pars,
       multichannel = ifelse(n_channels > 1, "multichannel_", ""),
-      intercept_only = icp_only_i && icp_only_s && icp_only_o && icp_only_d
+      intercept_only = icpt_only_i && icpt_only_s && icpt_only_o && icpt_only_d
     )
   )
 }
diff --git a/R/fit_mnhmm.R b/R/fit_mnhmm.R
@@ -47,10 +47,12 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, lambda, method,
   n_s <- attr(model, "np_A")
   n_o <- attr(model, "np_B")
   n_d <- attr(model, "np_omega")
-  iv_pi <- attr(model$X_pi, "iv")
+  icpt_only_omega <- attr(model$X_omega, "icpt_only")
+  icpt_only_pi <- attr(model$X_pi, "icpt_only")
+  icpt_only_A <- attr(model$X_A, "icpt_only")
+  icpt_only_B <- attr(model$X_B, "icpt_only")
   iv_A <- attr(model$X_A, "iv")
   iv_B <- attr(model$X_B, "iv")
-  iv_omega <- attr(model$X_omega, "iv")
   tv_A <- attr(model$X_A, "tv")
   tv_B <- attr(model$X_B, "tv")
   X_pi <- model$X_pi
@@ -117,7 +119,8 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, lambda, method,
         )
         out <- log_objective_mnhmm_singlechannel(
           eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B,
-          obs, iv_omega, iv_pi, iv_A, iv_B, tv_A, tv_B, Ti
+          obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, 
+          iv_A, iv_B, tv_A, tv_B
         )
         list(
           objective = - (out$loglik - 0.5 * lambda * sum(pars^2)) / n_obs, 
@@ -136,7 +139,8 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, lambda, method,
         )
         out <- forward_mnhmm_singlechannel(
           eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, 
-          obs, iv_omega, iv_pi, iv_A, iv_B, tv_A, tv_B, Ti
+          obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, 
+          iv_A, iv_B, tv_A, tv_B
         )
 
         - (sum(apply(out[, T_, ], 2, logSumExp)) - 0.5 * lambda * sum(pars^2)) / n_obs
@@ -161,7 +165,8 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, lambda, method,
         )
         out <- log_objective_mnhmm_multichannel(
           eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B,
-          obs, iv_omega, iv_pi, iv_A, iv_B, tv_A, tv_B, Ti
+          obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, 
+          iv_A, iv_B, tv_A, tv_B
         )
         list(
           objective = - (out$loglik - 0.5 * lambda * sum(pars^2)) / n_obs, 
@@ -184,11 +189,10 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, lambda, method,
           pars[n_i + n_s + n_o + seq_len(n_d)], D, K_omega
         )
         out <- forward_mnhmm_multichannel(
-          eta_pi, X_pi,
-          eta_A, X_A,
-          eta_B, X_B,
-          eta_omega, X_omega,
-          obs, M)
+          eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, 
+          obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, 
+          iv_A, iv_B, tv_A, tv_B
+        )
 
         - (sum(apply(out[, T_, ], 2, logSumExp)) - 0.5 * lambda * sum(pars^2)) / n_obs
       }

diff --git a/R/fit_nhmm.R b/R/fit_nhmm.R
@@ -46,7 +46,9 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method,
   n_i <- attr(model, "np_pi")
   n_s <- attr(model, "np_A")
   n_o <- attr(model, "np_B")
-  iv_pi <- attr(model$X_pi, "iv")
+  icpt_only_pi <- attr(model$X_pi, "icpt_only")
+  icpt_only_A <- attr(model$X_A, "icpt_only")
+  icpt_only_B <- attr(model$X_B, "icpt_only")
   iv_A <- attr(model$X_A, "iv")
   iv_B <- attr(model$X_B, "iv")
   tv_A <- attr(model$X_A, "tv")
@@ -93,8 +95,8 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method,
             pars[n_i + n_s + seq_len(n_o)], S, M, K_B
           )
           out <- log_objective_nhmm_singlechannel(
-            eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs,
-            iv_pi, iv_A, iv_B, tv_A, tv_B, Ti
+            eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti,
+            icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B
           )
           list(
             objective = - (out$loglik - 0.5 * lambda * sum(pars^2)) / n_obs, 
@@ -109,7 +111,8 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method,
             pars[n_i + n_s + seq_len(n_o)], S, M, K_B
           )
           out <- forward_nhmm_singlechannel(
-            eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs
+            eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti,
+            icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B
           )
           - (sum(apply(out[, T_, ], 2, logSumExp)) - 0.5 * lambda * sum(pars^2)) / n_obs
         }
@@ -123,8 +126,8 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method,
             pars[n_i + n_s + seq_len(n_o)], S, M, K_B
           )
           out <- log_objective_nhmm_multichannel(
-            eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs,
-            iv_pi, iv_A, iv_B, tv_A, tv_B, Ti
+            eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti,
+            icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B
           )
           list(
             objective = - (out$loglik - 0.5 * lambda * sum(pars^2)) / n_obs, 
@@ -139,7 +142,8 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method,
             pars[n_i + n_s + seq_len(n_o)], S, M, K_B
           )
           out <- forward_nhmm_multichannel(
-            eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs
+            eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti,
+            icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B
           )
           - (sum(apply(out[, T_, ], 2, logSumExp)) - 0.5 * lambda * sum(pars^2)) / n_obs
         }
@@ -217,8 +221,8 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method,
         )
         if (C == 1) {
           EM_LBFGS_nhmm_singlechannel(
-            init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B,
-            obs, iv_pi, iv_A, iv_B, tv_A, tv_B, Ti,
+            init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B, obs,
+            Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B,
             n_obs, control_restart$maxeval,
             control_restart$ftol_abs, control_restart$ftol_rel,
             control_restart$xtol_abs, control_restart$xtol_rel, 
@@ -228,8 +232,8 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method,
             control_mstep$print_level, lambda)
         } else {
           EM_LBFGS_nhmm_multichannel(
-            init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B,
-            obs, iv_pi, iv_A, iv_B, tv_A, tv_B, Ti,
+            init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B, obs,
+            Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B,
             n_obs, control_restart$maxeval,
             control_restart$ftol_abs, control_restart$ftol_rel,
             control_restart$xtol_abs, control_restart$xtol_rel, 
@@ -256,8 +260,8 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method,
     }
     if (C == 1) {
       out <- EM_LBFGS_nhmm_singlechannel(
-        init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B,
-        obs, iv_pi, iv_A, iv_B, tv_A, tv_B, Ti,
+        init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B, obs,
+        Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B,
         n_obs, control$maxeval,
         control$ftol_abs, control$ftol_rel,
         control$xtol_abs, control$xtol_rel, 
@@ -267,8 +271,8 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method,
         control_mstep$print_level, lambda)
     } else {
       out <- EM_LBFGS_nhmm_multichannel(
-        init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B,
-        obs, iv_pi, iv_A, iv_B, tv_A, tv_B, Ti,
+        init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B, obs,
+        Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B,
         n_obs, control$maxeval,
         control$ftol_abs, control$ftol_rel,
         control$xtol_abs, control$xtol_rel,