Request from Howard around R-HTA

[RobertASmithBresMed]

[RobertASmithBresMed]

# Test that the ArmaMarkov function matches
# the Rloop function for a simple example.

# remotes::install_github("RobertASmith/darkpeak")
# install.packages("microbenchmark")
library(testthat)
library(darkpeak)
library(microbenchmark)


# function to run time dependent markov model in base R
Rloopfunc <- function(m_TR,
                      a_P) {
  # throughout the number of cycles
  # estimate the Markov trace for cycle the next cycle (t + 1)
  for (t in 1:(nrow(m_TR) - 1)) {
    m_TR[t + 1,] <- m_TR[t,] %*% a_P[, , t]
  }
  # return the trace
  m_TR
}

# function to create example array
# number of health states defined using n_hs
makeTransProbArray <- function(n_hs = 4,
                               n_periods = 1000){
  # need the random number stream to be consistent- fix at seed = 100
  set.seed(100)
  # create an array from slicing this matrix
  array(
    data = runif(n_hs  * n_hs * n_periods),
    dim = c(n_hs, n_hs, n_periods)
  )
}

# function to create example markov trace
# number of health states defined using n_hs
makeMarkovTrace <- function(n_hs = 4,
                            n_periods = 1000) {
  # make empty trace
  m_TR <- matrix(data = 0,
                 nrow = n_periods,
                 ncol = n_hs)
  # initialise trace
  m_TR[1,]  <- c(1, rep(0, n_hs - 1))
  # return trace
  m_TR

}



#==================================================#
# EQUAL
#==================================================#
test_that("ArmaMarkov equals BaseRloop", {
  # test three different symmetric matrix dimensions:
  for(s in c(3,12,21)){

    # sample the number of periods - means we are testing this too.
    n_periods <- sample(x = c(1000, 323, 4780),
                        size = 1,
                        replace = T)

    # in each case test that the two approaches get equal answers
    testthat::expect_identical(

      ArmaTDMarkovLoop(m_TR = makeMarkovTrace(s, n_periods),
                       a_P = makeTransProbArray(s, n_periods)),

      Rloopfunc(m_TR = makeMarkovTrace(s, n_periods),
                a_P = makeTransProbArray(s, n_periods))

    ) # end expect equal
  }
}) # end testthat.



#==================================================#
# FASTER
#==================================================#
test_that("ArmaMarkov is faster than BaseRloop", {
  # test three different symmetric matrix dimensions:
  for(s in c(3,12,21)){

    # sample the number of periods - means we are testing this too.
    n_periods <- sample(x = c(1000, 323, 4780),
                        size = 1,
                        replace = T)

    mb <- microbenchmark::microbenchmark(

      Cpp = ArmaTDMarkovLoop(m_TR = makeMarkovTrace(s, n_periods),
                             a_P = makeTransProbArray(s, n_periods)),

      Rloop = Rloopfunc(m_TR = makeMarkovTrace(s, n_periods),
                        a_P = makeTransProbArray(s, n_periods))
    )


    # in each case test that the ArmaMarkov function is faster
    # Mean
    testthat::expect_false(
      mean(mb[mb$expr == "Cpp" ,"time"]) > mean(mb[mb$expr == "Rloop" ,"time"])
    ) # end expect equal

    # Median
    testthat::expect_false(
      median(mb[mb$expr == "Cpp" ,"time"]) > median(mb[mb$expr == "Rloop" ,"time"])
    ) # end expect equal

  }
}) # end testthat.



#==================================================#
# FOR HOWARD THOM
#==================================================#

s <- 5 ; n_periods = 100000

mb <- microbenchmark::microbenchmark(

  Cpp = ArmaTDMarkovLoop(m_TR = makeMarkovTrace(s, n_periods),
                         a_P = makeTransProbArray(s, n_periods)),

  Rloop = Rloopfunc(m_TR = makeMarkovTrace(s, n_periods),
                    a_P = makeTransProbArray(s, n_periods))
)

autoplot(mb)
print(mb)