R/02_trans_prob.R

In vsritter/sim.etm: What the Package Does (One Line, Title Case)

#' Title
#'
#' @param t
#' @param gen
#'
#' @return
#' @export
#'
#' @examples
mcgen <- function(t, gen) {
  idx <- which(!is.na(gen), arr.ind = T)

  a <- gen[subset(idx, idx[,3] == 1)]
  b <- gen[subset(idx, idx[,3] == 2)]

  cumhaz <- function(t) a*t^b
  Ht <- Vectorize(cumhaz)(t)

  qt <- array(0, dim = c(dim(gen)[-3], length(t)-1))
  for (i in 1:length(a)) {
    qt[idx[i,1], idx[i,2], ] <- diff(Ht[i,])
  }

  qt
}

illness_death_prob <- function(t, gen) {
  idx <- which(!is.na(gen), arr.ind = T)

  a <- gen[subset(idx, idx[,3] == 1)]
  b <- gen[subset(idx, idx[,3] == 2)]

  cumhaz <- function(t) a*t^b
  Ht <- Vectorize(cumhaz)(t)

  p11 <- exp(-(Ht[1,] + Ht[2,]))
  p22 <- exp(-Ht[3,])

  a12 = a[1]*b[1]
  a13 = a[2]*b[2]
  a23 = a[3]*b[3]
  b12 = b[1]-1
  b13 = b[2]-1
  b23 = b[3]-1

  s=0
  g <- function(u) {
    exp(-((a12/(b12+1)*u^(b12+1))+(a13/(b13+1)*u^(b13+1)))
        +((a12/(b12+1)*s^(b12+1))+(a13/(b13+1)*s^(b13+1))))*
      a12*u^b12*exp(-(a23/(b23+1)*ti^(b23+1))+
                      (a23/(b23+1)*u^(b23+1)))
  }

  p12=rep(0,length(t))
  e=1
  for(ti in t)
  {
    if(ti<=s)
    {p12[e]=0}
    else
    {p12[e]=integrate(g,s,ti)$value}
    e=e+1
  }

  pt <- array(0, dim = c(dim(gen)[-3], length(t)))
  pt[1, 1, ] <- p11
  pt[1, 2, ] <- p12
  pt[1, 3, ] <- 1 - p11 - p12
  pt[2, 2, ] <- p22
  pt[2, 3, ] <- 1 - p22

  pt
}