R/inverse_mc_sampler.R

In csetraynor/mctte: Compartmental models applied to time-to-event.

#' Simulation of a multi-state model
#' General function to generate a multi-state model requires a compiled model using mrgsolve 
#' for the Kolmogorov ODE, which is used to sample the state occupancy in the model.
#' @param modeldat mrgsolve model with defined the multi-state model.
#' @param id id parameter
#' @param covs covariates of the model
#' @param nstate number of states
#' @param absor absorving states, if more than one pass as a vector
#' @param qfun function to pass q.
#' @param G the G matrix of nstate x nstate.
#' @param A the alpha matrix of nstate x nstate.
#' @export
mcttesim_weibull <- function(modeldat, id, covs, nstate, absor, qfun, G, A){
  
  newdat = mcttesim1_weibull(modeldat, id, nstate, absor, qfun, G, A)
  newdat <- bind_rows(newdat)
  left_join(newdat, covs, by = "ID")
}

mcttesim1_weibull <- function(simdat, id, nstate, absor, q, G, A){
  Time1 = 0
  currstate = 1
  outdata = NULL
  while( (currstate %!in% absor) & Time1 < Inf ) {
    
    U = runif(1)
    pdat = simdat [ ,3:ncol(simdat)] ## only keep p columns
    pcurr = pdat[ , get_pcurrent(nstate, currstate)] ## add two cols for ID and time
    etime = get_tte(U, pcurr)
    
    if(etime != -99){ ### if not administrative right censoring found
      eventtime = get_time( simdat, etime)
      
      Gt = G[currstate, ]
      At = A[currstate, ]
      probs = seq_along(Gt) %>%
        map_dbl(function(i){
          q(eventtime, Gt[i], At[i])
        })
      newstate <- which(rmultinom(1, 1, probs) == 1)
      Time2 = Time1 + eventtime
      
      outdata <- rbind(
        outdata,
        data.frame(ID = id,
                   Time1 = Time1,
                   Time2 = Time2,
                   State1 = currstate,
                   State2 = newstate)
      )
      Time1 = Time2
      currstate = newstate
    } else { ## for administrative censoring
      Time2 = Time1 + get_max_time(simdat)
      outdata <- rbind(
        outdata,
        data.frame(ID = id,
                   Time1 = Time1,
                   Time2 = Time2,
                   State1 = currstate,
                   State2 = currstate)
      )
      Time1 = Inf
    }
  }
  outdata
}



#' Simulation of a multi-state model
#' General function to generate a multi-state model requires a compiled model using mrgsolve 
#' for the Kolmogorov ODE, which is used to sample the state occupancy in the model.
#' @param modeldat mrgsolve model with defined the multi-state model.
#' @param id id parameter
#' @param covs covariates of the model
#' @param nstate number of states
#' @param absor absorving states, if more than one pass as a vector
#' @param qfun function to pass q.
#' @param BS BS matrix of nstate x nstate.
#' @export
mcttesim_bs <- function(modeldat, id, covs, nstate, absor, qfun, BS){
  
  newdat = mcttesim1_weibull(modeldat, id, nstate, absor, qfun, BS)
  newdat <- bind_rows(newdat)
  left_join(newdat, covs, by = "ID")
}

mcttesim1_bs <- function(simdat, id, nstate, absor, q, G, A){
  Time1 = 0
  currstate = 1
  outdata = NULL
  while( (currstate %!in% absor) & Time1 < Inf ) {
    
    U = runif(1)
    pdat = simdat [ ,3:ncol(simdat)] ## only keep p columns
    pcurr = pdat[ , get_pcurrent(nstate, currstate)] ## add two cols for ID and time
    etime = get_tte(U, pcurr)
    
    if(etime != -99){ ### if not administrative right censoring found
      eventtime = get_time( simdat, etime)
      
      BSt = BS[currstate, ]
      probs = seq_along(Gt) %>%
        map_dbl(function(i){
          q(eventtime, BSt[i])
        })
      newstate <- which(rmultinom(1, 1, probs) == 1)
      Time2 = Time1 + eventtime
      
      outdata <- rbind(
        outdata,
        data.frame(ID = id,
                   Time1 = Time1,
                   Time2 = Time2,
                   State1 = currstate,
                   State2 = newstate)
      )
      Time1 = Time2
      currstate = newstate
    } else { ## for administrative censoring
      Time2 = Time1 + get_max_time(simdat)
      outdata <- rbind(
        outdata,
        data.frame(ID = id,
                   Time1 = Time1,
                   Time2 = Time2,
                   State1 = currstate,
                   State2 = currstate)
      )
      Time1 = Inf
    }
  }
  outdata
}




################ Helpers #############################
get_pcurrent <- function(nstate, currstate){
  diag(matrix(1:(nstate^2), nrow = nstate))[currstate]
}
get_tte <- function(U, pcurr){
  match(1, cumsum( pcurr %>% unlist() < U), nomatch = -99L)
}
get_pmatrix <- function(currpdat, nstate){
  matrix(c(currpdat), ncol = nstate, byrow = TRUE)
}

get_Q <- function(q, nstate){
  matrix(c(q), ncol = nstate, byrow = TRUE)
}

getq_probs <- function(q, currstate){
  q[currstate, ]
}

### Internally normalised to sum 1
get_probs <- function(transprob, currstate){
  shortprob = transprob[-currstate]
  sumprob = sum(shortprob)
  newprob = shortprob / sumprob
  newprob = c(newprob[seq_len(currstate-1)], 0, newprob[(currstate):length(newprob)])
  newprob
}

get_time <- function(simdat, etime){
  unlist(simdat[etime, 2])
}

get_max_time <- function(simdat) {
  unlist( simdat[nrow(simdat), 2] )
}


#' Simulation based on homogenous system
#' The homogenous is the basic multi-state model. 
#' @param modeldat mrgsolve model with defined the multi-state model.
#' @param id id parameter
#' @param covs covariates of the model
#' @param nstate number of states
#' @param absor absorving states, if more than one pass as a vector
#' @export
mcttesim <- function(modeldat, id, covs, nstate, absor, q){

  newdat = mcttesim1(modeldat, id, nstate, absor, q)
  newdat <- bind_rows(newdat)
  left_join(newdat, covs, by = "ID")
}

mcttesim1 <- function(simdat, id, nstate, absor, q){
  Time1 = 0
  currstate = 1
  outdata = NULL
  while( (currstate %!in% absor) & Time1 < Inf ) {

    U = runif(1)
    pdat = simdat [ ,3:ncol(simdat)] ## only keep p columns
    pcurr = pdat[ , get_pcurrent(nstate, currstate)] ## add two cols for ID and time
    etime = get_tte(U, pcurr)

    if(etime != -99){ ### if not administrative right censoring found
      eventtime = get_time( simdat, etime)
      qmat = get_Q(q, nstate)
      probs = unlist(qmat[currstate, ])
      newstate <- which(rmultinom(1, 1, probs) == 1)
      Time2 = Time1 + eventtime

      outdata <- rbind(
        outdata,
        data.frame(ID = id,
                   Time1 = Time1,
                   Time2 = Time2,
                   State1 = currstate,
                   State2 = newstate)
      )
      Time1 = Time2
      currstate = newstate
    } else { ## for administrative censoring
      Time2 = Time1 + get_max_time(simdat)
      outdata <- rbind(
        outdata,
        data.frame(ID = id,
                   Time1 = Time1,
                   Time2 = Time2,
                   State1 = currstate,
                   State2 = currstate)
      )
      Time1 = Inf
    }
  }
  outdata
}

# 
# mcttesim1id <- function(simdat, id, nstate, absor){
#   Time1 = 0
#   currstate = 1
#   outdata = NULL
#   while( (currstate != absor) & Time1 < Inf ) {
#     U = runif(1)
#     pdat = simdat [ ,3:ncol(simdat)] ## only keep p columns
#     pcurr = pdat[ , get_pcurrent(nstate, currstate)] ## add two cols for ID and time
#     etime = get_tte(U, pcurr)
# 
#     if(etime != -99){ ### if not administrative right censoring found
#       currpdat = pdat[etime, ]
#       pmatrix = get_pmatrix(currpdat, nstate)
#       transprob = unlist(pmatrix[currstate, ])
#       probs = get_probs(transprob, currstate)
#       newstate <- which(rmultinom(1, 1, probs) == 1)
#       Time2 = Time1 + get_time( simdat, etime)
# 
#       outdata <- rbind(
#         outdata,
#         data.frame(ID = id,
#                    Time1 = Time1,
#                    Time2 = Time2,
#                    State1 = currstate,
#                    State2 = newstate)
#       )
#       Time1 = Time2
#       currstate = newstate
#     } else { ## for administrative censoring
#       Time2 = Time1 + get_max_time(simdat)
#       outdata <- rbind(
#         outdata,
#         data.frame(ID = id,
#                    Time1 = Time1,
#                    Time2 = Time2,
#                    State1 = currstate,
#                    State2 = currstate)
#       )
#       Time1 = Inf
#     }
#   }
#   outdata
# }
#