estimateParams()

NAMESPACE

@@ -1,5 +1,15 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(getInit,ExponentialTransition)
+S3method(getInit,WeibullTransition)
+S3method(getResults,ExponentialTransition)
+S3method(getResults,WeibullTransition)
+S3method(getTarget,ExponentialTransition)
+S3method(getTarget,WeibullTransition)
+S3method(haz,ExponentialTransition)
+S3method(haz,WeibullTransition)
+S3method(survTrans,ExponentialTransition)
+S3method(survTrans,WeibullTransition)
 export(ExpHazOS)
 export(ExpQuantOS)
 export(ExpSurvOS)
@@ -18,26 +28,33 @@ export(avgHRExpOS)
 export(avgHRIntegExpOS)
 export(censoringByNumberEvents)
 export(empSignificant)
+export(estimateParams)
 export(exponential_transition)
 export(getCensoredData)
 export(getClinicalTrials)
 export(getDatasetWideFormat)
 export(getEventsAll)
+export(getInit)
 export(getNumberEvents)
 export(getOneClinicalTrial)
 export(getOneToTwoRows)
 export(getPCWDistr)
 export(getPCWHazard)
+export(getResults)
 export(getSimulatedData)
 export(getSumPCW)
+export(getTarget)
 export(getTimePoint)
 export(getWaitTimeSum)
+export(haz)
 export(integrateVector)
 export(logRankTest)
+export(negLogLik)
 export(passedLogRank)
 export(piecewise_exponential)
 export(prepareData)
 export(pwA)
+export(survTrans)
 export(trackEventsPerTrial)
 export(weibull_transition)
 import(checkmate)

NEWS.md

@@ -3,6 +3,7 @@
 ### New Features
 
 -   `prepareData` allows formatting of trial data for log-likelihood computation.
+-   `estimateParams` estimates parameters of exponential and Weibull transition hazards for a given data set.
 
 ### Bug Fixes
 

R/estimateParams.R

@@ -53,3 +53,307 @@ prepareData <- function(data) {
 
   as.data.frame(dataNew[, -which(names(dataNew) == "time")], row.names = seq_len(nrow(dataNew)))
 }
+
+#' Compute the Negative Log-Likelihood for a Given Data Set and Transition Model
+#'
+#' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr
+#'   see [exponential_transition()] or [weibull_transition()] for details.
+#' @param data (`data.frame`)\cr in the format created by [prepareData()].
+#'
+#' @return The value of the negative log-likelihood.
+#' @export
+#'
+#' @details
+#' Calculates the negative log-likelihood for a given data set and transition model. It uses the hazard
+#' and survival functions specific to the transition model.
+#'
+#' @examples
+#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
+#' simData <- getOneClinicalTrial(
+#'   nPat = c(30), transitionByArm = list(transition),
+#'   dropout = list(rate = 0.8, time = 12),
+#'   accrual = list(param = "time", value = 1)
+#' )
+#' negLogLik(transition, prepareData(simData))
+negLogLik <- function(transition, data) {
+  with(data, -sum(log(haz(transition, exit, trans)^status
+    * survTrans(transition, exit, trans) / survTrans(transition, entry, trans))))
+}
+
+#' Hazard Function for Different Transition Models
+#'
+#' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr
+#'   see [exponential_transition()] or [weibull_transition()] for details.
+#' @param t (`numeric`)\cr time at which hazard is to be computed.
+#' @param trans (`integer`)\cr index specifying the transition type.
+#'
+#' @return The hazard rate for the specified transition and time.
+#' @export
+#'
+#' @examples
+#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
+#' haz(transition, 0.4, 2)
+haz <- function(transition, t, trans) {
+  assert_class(transition, "TransitionParameters")
+  assert_numeric(t, lower = 0)
+  assert_subset(trans, c(1, 2, 3))
+  UseMethod("haz")
+}
+
+#' @describeIn haz for an exponential transition model.
+#' @export
+#'
+#' @examples
+#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
+#' haz(transition, 0.4, 2)
+haz.ExponentialTransition <- function(transition, t, trans) {
+  # params (in this order): h01, h02, h12.
+  params <- unlist(transition$hazards, use.names = FALSE)
+  params[trans]
+}
+
+#' @describeIn haz for the Weibull transition model.
+#' @export
+#'
+#' @examples
+#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
+#' haz(transition, 0.4, 2)
+haz.WeibullTransition <- function(transition, t, trans) {
+  # params (in this order): h01, h02, h12, p01, p02, p12.
+  params <- c(unlist(transition$hazards, use.names = FALSE), unlist(transition$weibull_rates, use.names = FALSE))
+  params[trans] * params[trans + 3] * t^(params[trans + 3] - 1)
+}
+
+#' Survival Function for Different Transition Models
+#'
+#' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr
+#'   see [exponential_transition()] or [weibull_transition()] for details.
+#' @param t (`numeric`)\cr time at which survival probability is to be computed.
+#' @param trans (`integer`)\cr index specifying the transition type.
+#'
+#' @return The survival probability for the specified transition and time.
+#' @export
+#'
+#' @examples
+#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
+#' survTrans(transition, 0.4, 2)
+survTrans <- function(transition, t, trans) {
+  assert_class(transition, "TransitionParameters")
+  assert_numeric(t, lower = 0)
+  assert_subset(trans, c(1, 2, 3))
+  UseMethod("survTrans")
+}
+
+#' @describeIn survTrans for the Exponential Transition Model
+#' @export
+#'
+#' @examples
+#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
+#' survTrans(transition, 0.4, 2)
+survTrans.ExponentialTransition <- function(transition, t, trans) {
+  # params (in this order): h01, h02, h12.
+  params <- unlist(transition$hazards, use.names = FALSE)
+  exp(-params[trans] * t)
+}
+
+#' @describeIn survTrans for the Weibull Transition Model
+#' @export
+#'
+#' @examples
+#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
+#' survTrans(transition, 0.4, 2)
+survTrans.WeibullTransition <- function(transition, t, trans) {
+  # params (in this order): h01, h02, h12, p01, p02, p12.
+  params <- c(unlist(transition$hazards, use.names = FALSE), unlist(transition$weibull_rates, use.names = FALSE))
+  exp(-params[trans] * t^params[trans + 3])
+}
+
+#' Retrieve Initial Parameter Vectors for Likelihood Maximization
+#'
+#' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr containing the initial parameters.
+#'   See [exponential_transition()] or [weibull_transition()] for details.
+#'
+#' @return The numeric vector of initial parameters for likelihood maximization.
+#' @export
+#'
+#' @examples
+#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
+#' getInit(transition)
+getInit <- function(transition) {
+  assert_class(transition, "TransitionParameters")
+  UseMethod("getInit")
+}
+
+#' @describeIn getInit for the Exponential Transition Model
+#' @export
+#'
+#' @examples
+#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
+#' getInit(transition)
+getInit.ExponentialTransition <- function(transition) {
+  unlist(transition$hazards, use.names = FALSE)
+}
+
+#' @describeIn getInit for the Weibull Transition Model
+#' @export
+#'
+#' @examples
+#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
+#' getInit(transition)
+getInit.WeibullTransition <- function(transition) {
+  c(unlist(transition$hazards, use.names = FALSE), unlist(transition$weibull_rates, use.names = FALSE))
+}
+
+#' Generate the Target Function for Optimization
+#'
+#' @param transition (`TransitionParameters`)\cr
+#'   specifying the distribution family. See [exponential_transition()] or [weibull_transition()] for details.
+#'
+#' @return Function that calculates the negative log-likelihood for the given parameters.
+#' @export
+#'
+#' @details
+#' This function creates a target function for optimization, computing the negative log-likelihood for given
+#' parameters, data, and transition model type.
+#'
+#' @examples
+#' transition <- exponential_transition(2, 1.3, 0.8)
+#' simData <- getOneClinicalTrial(
+#'   nPat = c(30), transitionByArm = list(transition),
+#'   dropout = list(rate = 0.8, time = 12),
+#'   accrual = list(param = "time", value = 1)
+#' )
+#' params <- c(1.2, 1.5, 1.6) # For ExponentialTransition
+#' data <- prepareData(simData)
+#' transition <- exponential_transition()
+#' fun <- getTarget(transition)
+#' fun(params, data)
+getTarget <- function(transition) {
+  assert_class(transition, "TransitionParameters")
+  UseMethod("getTarget", transition)
+}
+
+#' @describeIn getTarget for the Exponential Transition Model
+#' @export
+#'
+#' @examples
+#' transition <- exponential_transition(2, 1.3, 0.8)
+#' simData <- getOneClinicalTrial(
+#'   nPat = c(30), transitionByArm = list(transition),
+#'   dropout = list(rate = 0.8, time = 12),
+#'   accrual = list(param = "time", value = 1)
+#' )
+#' params <- c(1.2, 1.5, 1.6)
+#' data <- prepareData(simData)
+#' transition <- exponential_transition()
+#' target <- getTarget(transition)
+#' target(params, data)
+getTarget.ExponentialTransition <- function(transition) {
+  function(params, data) {
+    negLogLik(transition = exponential_transition(h01 = params[1], h02 = params[2], h12 = params[3]), data = data)
+  }
+}
+
+#' @describeIn getTarget for the Weibull Transition Model
+#' @export
+#'
+#' @examples
+#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
+#' simData <- getOneClinicalTrial(
+#'   nPat = c(30), transitionByArm = list(transition),
+#'   dropout = list(rate = 0.8, time = 12),
+#'   accrual = list(param = "time", value = 1)
+#' )
+#' params <- c(1.2, 1.5, 1.6, 0.8, 1.3, 1.1)
+#' data <- prepareData(simData)
+#' transition <- weibull_transition()
+#' target <- getTarget(transition)
+#' target(params, data)
+getTarget.WeibullTransition <- function(transition) {
+  function(params, data) {
+    negLogLik(transition = weibull_transition(
+      h01 = params[1], h02 = params[2], h12 = params[3],
+      p01 = params[4], p02 = params[5], p12 = params[6]
+    ), data = data)
+  }
+}
+
+#' Format Results of Parameter Estimation for Different Transition Models
+#'
+#' @param transition (`TransitionParameters`)\cr
+#'   see [exponential_transition()] or [weibull_transition()] for details.
+#' @param res (`numeric` vector)\cr vector of parameter estimates from the likelihood maximization procedure.
+#'
+#' @return Returns a `TransitionParameters` object with parameter estimates.
+#' @export
+#'
+#' @examples
+#' results <- c(1.2, 1.5, 1.6)
+#' getResults(exponential_transition(), results)
+getResults <- function(transition, res) {
+  UseMethod("getResults")
+}
+
+#' @describeIn getResults for the Exponential Transition Model
+#' @export
+#'
+#' @examples
+#' results <- c(1.2, 1.5, 1.6)
+#' getResults(exponential_transition(), results)
+getResults.ExponentialTransition <- function(transition, res) {
+  exponential_transition(h01 = res[1], h02 = res[2], h12 = res[3])
+}
+
+#' @describeIn getResults for the Weibull Transition Model
+#' @export
+#'
+#' @examples
+#' results <- c(1.2, 1.5, 1.6, 2, 2.5, 1)
+#' getResults(weibull_transition(), results)
+getResults.WeibullTransition <- function(transition, res) {
+  weibull_transition(
+    h01 = res[1], h02 = res[2], h12 = res[3],
+    p01 = res[4], p02 = res[5], p12 = res[6]
+  )
+}
+
+#' Estimate Parameters of the Multistate Model Using Clinical Trial Data
+#'
+#' @param data (`data.frame`)\cr in the format produced by [getOneClinicalTrial()].
+#' @param transition (`TransitionParameters` object)\cr specifying the assumed distribution of transition hazards.
+#'   Initial parameters for optimization can be specified here.
+#'   See [exponential_transition()] or [weibull_transition()] for details.
+#'
+#' @return Returns a `TransitionParameters` object with the estimated parameters.
+#' @export
+#'
+#' @details
+#' This function estimates parameters for transition models using clinical trial data.
+#' The `transition` object can be initialized with starting values for parameter estimation.
+#' It uses [stats::optim()] to optimize the parameters.
+#'
+#' @examples
+#' transition <- exponential_transition(h01 = 2, h02 = 1.4, h12 = 1.6)
+#' simData <- getOneClinicalTrial(
+#'   nPat = c(30), transitionByArm = list(transition),
+#'   dropout = list(rate = 0.3, time = 12),
+#'   accrual = list(param = "time", value = 1)
+#' )
+#' # Initialize transition with desired starting values for optimization:
+#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
+#' estimate <- estimateParams(simData, transition)
+estimateParams <- function(data, transition) {
+  data <- prepareData(data)
+  par <- getInit(transition)
+  target <- getTarget(transition)
+
+  res <- stats::optim(
+    par = par,
+    fn = target,
+    method = "L-BFGS-B",
+    lower = 1e-3,
+    data = data
+  )$par
+
+  getResults(transition, res)
+}