fan-hmm

NAMESPACE

@@ -85,6 +85,7 @@ export(build_mhmm)
 export(build_mm)
 export(build_mmm)
 export(cluster_names)
+export(estimate_fanhmm)
 export(estimate_mnhmm)
 export(estimate_nhmm)
 export(fit_model)

R/RcppExports.R

@@ -41,6 +41,10 @@ eta_to_gamma_cube_field <- function(eta) {
     .Call(`_seqHMM_eta_to_gamma_cube_field`, eta)
 }
 
+EM_LBFGS_fanhmm_singlechannel <- function(eta_pi, x, eta_A, X_A, eta_B, X_B, rho_A, W_A, rho_B, W_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound, obs_0) {
+    .Call(`_seqHMM_EM_LBFGS_fanhmm_singlechannel`, eta_pi, x, eta_A, X_A, eta_B, X_B, rho_A, W_A, rho_B, W_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound, obs_0)
+}
+
 fast_quantiles <- function(X, probs) {
     .Call(`_seqHMM_fast_quantiles`, X, probs)
 }

R/ame_obs.R

@@ -102,6 +102,7 @@ ame_obs.nhmm <- function(
   newdata[[variable]][newdata[[time]] >= start_time] <- values[2]
   X2 <- update(model, newdata)[c("X_pi", "X_A", "X_B")]
   C <- model$n_channels
+  start <- which(sort(unique(newdata[[time]])) == start_time)
   if (C == 1L) {
     times <- as.numeric(colnames(model$observations))
     symbol_names <- list(model$symbol_names)
@@ -117,7 +118,7 @@ ame_obs.nhmm <- function(
       attr(X2$X_B, "icpt_only"), attr(X2$X_A, "iv"), 
       attr(X2$X_B, "iv"), attr(X2$X_A, "tv"), attr(X2$X_B, "tv"), 
       X2$X_pi, X2$X_A, X2$X_B,
-      model$boot$gamma_pi, model$boot$gamma_A, model$boot$gamma_B, start_time, 
+      model$boot$gamma_pi, model$boot$gamma_A, model$boot$gamma_B, start, 
       probs, model$boot$idx - 1L
     )
     d <- data.frame(
@@ -145,7 +146,7 @@ ame_obs.nhmm <- function(
       attr(X2$X_B, "icpt_only"), attr(X2$X_A, "iv"), 
       attr(X2$X_B, "iv"), attr(X2$X_A, "tv"), attr(X2$X_B, "tv"), 
       X2$X_pi, X2$X_A, X2$X_B,
-      model$boot$gamma_pi, model$boot$gamma_A, model$boot$gamma_B, start_time, 
+      model$boot$gamma_pi, model$boot$gamma_A, model$boot$gamma_B, start, 
       probs, model$boot$idx - 1L
     )
     d <- data.frame(
@@ -159,6 +160,7 @@ ame_obs.nhmm <- function(
       }
     }
   }
+  colnames(d)[2] <- time
   d[d[[time]] >= start_time, ]
 }
 

R/bootstrap.R

@@ -96,7 +96,7 @@ bootstrap_coefs.nhmm <- function(model, nsim = 1000,
     checkmate::test_int(x = nsim, lower = 0L), 
     "Argument {.arg nsim} must be a single positive integer."
   )
-  init <- model$etas
+  init <- setNames(model$etas, c("eta_pi", "eta_A", "eta_B"))
   gammas_mle <- model$gammas
   lambda <- model$estimation_results$lambda
   bound <- model$estimation_results$bound
@@ -199,7 +199,7 @@ bootstrap_coefs.mnhmm <- function(model, nsim = 1000,
     checkmate::test_int(x = nsim, lower = 0L), 
     "Argument {.arg nsim} must be a single positive integer."
   )
-  init <- model$etas
+  init <- setNames(model$etas, c("eta_pi", "eta_A", "eta_B", "eta_omega"))
   gammas_mle <- model$gammas
   pcp_mle <- posterior_cluster_probabilities(model)
   lambda <- model$estimation_results$lambda

R/build_fanhmm.R

@@ -0,0 +1,88 @@
+#' Build a Feedback-Augmented Non-homogeneous Hidden Markov Model
+#'
+#' @noRd
+build_fanhmm <- function(
+    observations, n_states, initial_formula, 
+    transition_formula, emission_formula, autoregression_formula, 
+    feedback_formula, data, time, id, state_names = NULL) {
+
+  stopifnot(
+    !is.na(autoregression_formula) || !is.na(feedback_formula),
+    "Provide {.arg autoregression_formula} and/or {.arg feedback_formula} for FAN-HMM."
+  )
+  out <- create_base_nhmm(
+    observations, data, time, id, n_states, state_names, channel_names = NULL,
+    initial_formula, transition_formula, emission_formula) 
+  out[c("cluster_names", "n_clusters", "X_omega")] <- NULL
+  setNames(
+    create_initial_values(list(), out$model, 0), c("pi", "A", "B")
+  )
+  stopifnot_(
+    out$n_channels == 1L,
+    "Currently only single-channel responses are supported for FAN-HMM.")
+  if(is.na(autoregression_formula)) {
+    out$model$W_A <- array(
+      0, c(0L, out$model$length_of_sequences, out$model$n_sequences)
+    )
+    out$model$rhos$A <- create_rho_A_inits(
+      NULL, n_states, out$model$n_symbols, 0, 0
+    )
+    np_rho_A <- 0
+  } else {
+    W_A <- model_matrix_autoregression_formula(
+      autoregression_formula, data, 
+      out$model$n_sequences, 
+      out$model$length_of_sequences, n_states,
+      out$model$n_symbols, time, id, 
+      out$model$sequence_lengths
+    )
+    out$model$W_A <- W_A$X
+    out$model$rhos$A <- create_rho_A_inits(
+      NULL, n_states, out$model$n_symbols, nrow(out$model$W_A), 0
+    )
+    out$extras$intercept_only <- FALSE
+    np_rho_A <- W_A$n_pars
+  }
+  if(is.na(feedback_formula)) {
+    out$model$W_B <- array(
+      0, c(0L, out$model$length_of_sequences, out$model$n_sequences)
+    )
+    out$model$rhos$B <- create_rho_B_inits(
+      NULL, n_states, out$model$n_symbols, 0, 0
+    )
+    np_rho_B <- 0
+  } else {
+    W_B <- model_matrix_feedback_formula_formula(
+      feedback_formula_formula, data, 
+      out$model$n_sequences, 
+      out$model$length_of_sequences, n_states,
+      out$model$n_symbols, time, id, 
+      out$model$sequence_lengths
+    )
+    out$model$W_B <- W_B$X
+    out$model$rhos$B <- create_rho_B_inits(
+      NULL, n_states, out$model$n_symbols, nrow(out$model$W_B), 0
+    )
+    out$extras$intercept_only <- FALSE
+    np_rho_B <- W_B$n_pars
+  }
+  structure(
+    c(
+      out$model,
+      list(autoregression_formula = autoregression_formula, 
+           feedback_formula_formula = feedback_formula_formula
+      )
+    ),
+    class = "fanhmm",
+    nobs = attr(out$model$observations, "nobs"),
+    df = out$extras$np_pi + out$extras$np_A + out$extras$np_B + np_rho_A + 
+      np_rho_B,
+    type = paste0(out$extras$multichannel, "fanhmm"),
+    intercept_only = out$extras$intercept_only,
+    np_pi = out$extras$np_pi,
+    np_A = out$extras$np_A,
+    np_B = out$extras$np_B,
+    np_rho_A = np_rho_A,
+    np_rho_B = np_rho_B
+  )
+}

R/build_mnhmm.R

@@ -23,7 +23,9 @@ build_mnhmm <- function(
     observations, data, time, id, n_states, state_names, channel_names, 
     initial_formula, transition_formula, emission_formula, cluster_formula, 
     cluster_names)
-  out$model$etas <- create_initial_values(list(), out$model, 0)
+  out$model$etas <- setNames(
+    create_initial_values(list(), out$model, 0), c("pi", "A", "B", "omega")
+  )
   structure(
     out$model,
     class = "mnhmm",

R/build_nhmm.R

@@ -10,7 +10,9 @@ build_nhmm <- function(
     observations, data, time, id, n_states, state_names, channel_names,
     initial_formula, transition_formula, emission_formula) 
   out[c("cluster_names", "n_clusters", "X_omega")] <- NULL
-  out$model$etas <- create_initial_values(list(), out$model, 0)
+  out$model$etas <- setNames(
+    create_initial_values(list(), out$model, 0), c("pi", "A", "B")
+  )
   structure(
     out$model,
     class = "nhmm",

R/create_base_nhmm.R

@@ -121,7 +121,7 @@ create_base_nhmm <- function(observations, data, time, id, n_states,
   )
   B <- model_matrix_emission_formula(
     emission_formula, data, n_sequences, length_of_sequences, n_states, 
-    n_symbols, n_channels, time, id, sequence_lengths
+    n_symbols, time, id, sequence_lengths
   )
   if (mixture) {
     omega <- model_matrix_cluster_formula(

R/create_initial_values.R

@@ -23,84 +23,84 @@ create_initial_values_ <- function(inits, init_sd, S, M, K_pi, K_A, K_B, D = 1,
                                    K_omega = 0) {
   if(!is.null(inits$initial_probs)) {
     if (D > 1) {
-      pi <- lapply(
+      eta_pi <- lapply(
         seq_len(D), function(i) {
           create_inits_vector(inits$initial_probs[[i]], S, K_pi, init_sd)
         }
       )
     } else {
-      pi <- create_inits_vector(inits$initial_probs, S, K_pi, init_sd)
+      eta_pi <- create_inits_vector(inits$initial_probs, S, K_pi, init_sd)
     }
   } else {
-    pi <- create_eta_pi_inits(
-      inits$pi, S, K_pi, init_sd, D
+    eta_pi <- create_eta_pi_inits(
+      inits$eta_pi, S, K_pi, init_sd, D
     )
   }
 
   if(!is.null(inits$transition_probs)) {
     if (D > 1) {
-      A <- lapply(
+      eta_A <- lapply(
         seq_len(D), function(i) {
           create_inits_matrix(inits$transition_probs[[i]], S, S, K_A, init_sd)
         }
       )
     } else {
-      A <- create_inits_matrix(
+      eta_A <- create_inits_matrix(
         inits$transition_probs, S, S, K_A, init_sd
       )
     }
   } else {
-    A <- create_eta_A_inits(inits$A, S, K_A, init_sd, D)
+    eta_A <- create_eta_A_inits(inits$eta_A, S, K_A, init_sd, D)
   }
 
   if(!is.null(inits$emission_probs)) {
     if (D > 1) {
       if (length(M) > 1) {
-        B <- lapply(
+        eta_B <- lapply(
           seq_len(D), function(i) {
             lapply(seq_len(length(M)), function(j) {
               create_inits_matrix(
                 inits$emission_probs[[i]][[j]], S, M[j], K_B, init_sd)
             })
           })
       } else {
-        B <- lapply(
+        eta_B <- lapply(
           seq_len(D), function(i) {
             create_inits_matrix(inits$emission_probs[[i]], S, M, K_B, init_sd)
           }
         )
       }
     } else {
       if (length(M) > 1) {
-        B <- lapply(seq_len(length(M)), function(j) {
+        eta_B <- lapply(seq_len(length(M)), function(j) {
           create_inits_matrix(
             inits$emission_probs[[j]], S, M[j], K_B, init_sd)
         })
       } else {
-        B <- create_inits_matrix(
+        eta_B <- create_inits_matrix(
           inits$emission_probs, S, M, K_B, init_sd
         )
       }
     }
   } else {
-    B <- create_eta_B_inits(inits$B, S, M, K_B, init_sd, D)
+    eta_B <- create_eta_B_inits(inits$eta_B, S, M, K_B, init_sd, D)
   }
   out <- list(
-    pi = pi,
-    A = A,
-    B = B
+    eta_pi = eta_pi,
+    eta_A = eta_A,
+    eta_B = eta_B
   )
   if (D > 1) {
     if(!is.null(inits$cluster_probs)) {
-      omega <- create_inits_vector(
+      eta_omega <- create_inits_vector(
         inits$cluster_probs, D, K_omega, init_sd
       )
     } else {
-      omega <- create_eta_omega_inits(
-        inits$omega, D, K_omega, init_sd
+      eta_omega <- create_eta_omega_inits(
+        inits$eta_omega, D, K_omega, init_sd
       )
     }
-    out$omega <- omega
+    out$eta_omega <- eta_omega
   }
   out
 }
@@ -302,4 +302,47 @@ create_inits_matrix <- function(x, n, m, K, sd = 0) {
   z
 }
 
+create_rho_A_inits <- function(x, S, M, L, init_sd = 0) {
+  if (is.null(x$rho_A)) {
+    create_rho_A(numeric((S - 1) * L * (M - 1) * S), S, M, L, init_sd)
+  } else {
+    stopifnot_(
+      length(x$rho_A) == (S - 1) * L * (M - 1) * S,
+      paste0(
+        "Number of initial values for {.val rho_A} is not equal to ",
+        "(S - 1) * L * (M - 1) * S = {(S - 1) * L * (M - 1) * S}."
+      )
+    )
+    create_rho_A(x$rho_A, S, M, L, init_sd)
+  }
+}
+
+create_rho_A <- function(x, S, M, L, init_sd = 0) {
+  n <- (S - 1) * L * (M - 1)
+  lapply(seq_len(S), function(i) {
+    array(rnorm(n, x[(i - 1) * n + 1:n], init_sd), c(S - 1, L, M - 1))
+  })
+}
+
+create_rho_B_inits <- function(x, S, M, L, init_sd = 0) {
+  if (is.null(x$rho_B)) {
+    create_rho_B(numeric((M - 1) * L * (M - 1) * S), S, M, L, init_sd)
+  } else {
+    stopifnot_(
+      length(x$rho_B) == (M - 1) * L * (M - 1) * S,
+      paste0(
+        "Number of initial values for {.val rho_B} is not equal to ",
+        "(M - 1) * L * (M - 1) * S = {(M - 1) * L * (M - 1) * S}."
+      )
+    )
+    create_rho_B(x$rho_B, S, M, L, init_sd)
+  }
+}
+
+create_rho_B <- function(x, S, M, L, init_sd = 0) {
+  n <- (M - 1) * L * (M - 1)
+  lapply(seq_len(S), function(i) {
+    array(rnorm(n, x[(i - 1) * n + 1:n], init_sd), c(M - 1, L, M - 1))
+  })
+}