Internal developments with phase-type

.Rbuildignore

@@ -13,3 +13,4 @@ vignettes/examples*
 bin
 ^codecov\.yml$
 ^\.github$
+.Rprofile

.gitignore

@@ -8,3 +8,4 @@
 *~
 vignettes/*_cache*
 vignettes/*_files*
+*.exe

DESCRIPTION

@@ -8,10 +8,11 @@ Description: Bayesian multi-state models for intermittently-observed data.  Mark
 License: GPL (>= 3)
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 Depends: 
     R (>= 4.1.0)
 Imports:
+    abind,
     dplyr,
     tidyr,
     magrittr,

NAMESPACE

@@ -18,6 +18,7 @@ export(qmatrix)
 export(standardise_to)
 export(standardize_to)
 export(totlos)
+importFrom(abind,abind)
 importFrom(cli,cli_abort)
 importFrom(cli,cli_progress_bar)
 importFrom(cli,cli_progress_done)
@@ -50,6 +51,7 @@ importFrom(posterior,as_draws)
 importFrom(posterior,as_draws_matrix)
 importFrom(posterior,draws_of)
 importFrom(posterior,ess_bulk)
+importFrom(posterior,is_rvar)
 importFrom(posterior,merge_chains)
 importFrom(posterior,ndraws)
 importFrom(posterior,rdo)

R/msmbayes-package.R

@@ -8,7 +8,7 @@
 #' 
 #' @name msmbayes-package
 #' @importFrom stats delete.response na.omit reshape setNames terms quantile runif qnorm
-#' @importFrom posterior as_draws as_draws_matrix rhat ess_bulk rvar ndraws rvar_sum "%**%" rdo rvar_sum draws_of merge_chains
+#' @importFrom posterior as_draws as_draws_matrix rhat ess_bulk rvar ndraws rvar_sum "%**%" rdo rvar_sum draws_of merge_chains is_rvar
 #' @importFrom cli cli_abort cli_warn qty cli_progress_bar cli_progress_update cli_progress_done
 #' @importFrom glue glue
 #' @importFrom rlang caller_env .data
@@ -17,6 +17,7 @@
 #' @importFrom tidyr pivot_longer
 #' @importFrom utils head
 #' @importFrom stringr str_match
+#' @importFrom abind abind
 #'
 #' @md
 "_PACKAGE"

R/msmbayes.R

@@ -11,9 +11,9 @@
 #'   where K is the number of states, or a factor with these integers
 #'   as level labels.
 #'
-#' @param time Character string naming the observation time variable in the data
+#' @param time Character string naming the observation time variable in the data.
 #'
-#' @param subject Character string naming the individual ID variable in the data
+#' @param subject Character string naming the individual ID variable in the data.
 #'
 #' @param Q Matrix indicating the transition structure.  A zero entry
 #'   indicates that instantaneous transitions from (row) to (column)
@@ -107,13 +107,11 @@ msmbayes <- function(data, state, time, subject,
                      ...){
   qm <- form_qmodel(Q)
 
-  if (!is.null(nphase)){
-    pm <- form_phasetype(nphase, Q)
+  pm <- form_phasetype(nphase, Q)
+  if (pm$phasetype){
     qm <- phase_expand_qmodel(qm, pm)
     E <- pm$E
-  } else {
-    pm <- list(phasetype=FALSE, pdat=NULL)
-  }
+  } 
   em <- form_emodel(E, pm$Efix)
 
   check_data(data, state, time, subject, qm) 

R/msmbayes_internal.R

@@ -1,5 +1,7 @@
 #' @return List of information about the transition structure
 #'
+#' (TODO document fully)
+#'
 #' @noRd
 form_qmodel <- function(Q,Qfix=NULL){
   check_Q(Q)
@@ -14,14 +16,14 @@ form_qmodel <- function(Q,Qfix=NULL){
   } else {
     qfixrow <- qfixcol <- qfix <- as.array(numeric(0))
   }
+  qrow <- row(Q)[Q>0]
+  qcol <- col(Q)[Q>0]
+  qvec <- Q[cbind(qrow,qcol)]
   list(
-    Q = Q,
-    K = nrow(Q),
+    Q = Q, K = nrow(Q),
+    qvec = qvec, qrow=qrow, qcol=qcol,
     nqpars = length(Q[Q>0]),
-    qrow = row(Q)[Q>0],
-    qcol = col(Q)[Q>0],
-    qfixrow = qfixrow,
-    qfixcol = qfixcol,
+    qfixrow = qfixrow, qfixcol = qfixcol,
     qfix = qfix
   )
 }
@@ -89,3 +91,50 @@ transient_states <- function(qm){
 absorbing_states <- function(qm){
   which(rowSums(qm$Q) == 0)
 }
+
+
+## Utilities for transition intensity matrices without reference to
+## msmbayes
+
+## Mean sojourn time 
+Q_to_mst <- function(Q){
+  diag(Q) <- 0
+  1 / rowSums(Q)
+}
+
+##' Convert phase-type transition intensities to mixture representation
+##'
+##' A higher-level wrapper around `phase_mixture` which does the core calculation
+
+##' Currently unused and untested
+##'
+##' @param Qphase Intensity matrix on phased space
+##'
+##' @param nphase Numeric vector concatenating number of phases per state
+##'
+##' @return List with the components:
+##'
+##' \code{mix}: Mixture probs and mean sojourn times conditional on mixture component, from \code{\link{phase_mixture}}
+##'
+##' \code{mst}: Marginal mean sojourn times
+##'
+##'
+##' @noRd
+Qphase_to_mix <- function(Qphase, nphase){
+  K <- sum(nphase)
+  stopifnot(nrow(Qphase)==K && ncol(Qphase)==K && K>0)
+  qm <- form_qmodel(Qphase)
+  pdat <- form_phasedata(nphase)
+  tdat <- form_phasetrans(qm, pdat)
+  mix <- list()
+  mst <- numeric()
+  for (i in seq_along(nphase)){
+    mix[[i]] <-  cbind(
+      state = i,
+      phase = seq(nphase[i]),
+      phase_mixture(qm$qvec, tdat, i)
+    )
+    mst[[i]] <- rvarn_sum(mix[[i]]$mixprob * mix[[i]]$mst)
+  }
+  list(mix=do.call("rbind",mix), mst=mst)
+}