#' @include helpers.R
#' @include helpers_jags.R
#' @include Model-validity.R
#' @include ModelParams-class.R
#' @include CrmPackClass-class.R
NULL
# GeneralModel-class ----
#' `GeneralModel`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`GeneralModel`] is a general model class, from which all other specific
#' model-like classes inherit.
#'
#' @note The `datamodel` must obey the convention that the data input is
#' called exactly in the same way as in the corresponding data class.
#' All prior distributions for parameters should be contained in the
#' model function `priormodel`. The background is that this can
#' be used to simulate from the prior distribution, before obtaining any data.
#'
#' @slot datamodel (`function`)\cr a function representing the `JAGS` data model
#' specification.
#' @slot priormodel (`function`)\cr a function representing the `JAGS` prior
#' specification.
#' @slot modelspecs (`function`)\cr a function computing the list of the data
#' model and prior model specifications that are required to be specified
#' completely (e.g. prior parameters, reference dose, etc.), based on the data
#' slots that are required as arguments of this function.
#' Apart of data arguments, this function can be specified with one additional
#' (optional) argument `from_prior` of type `logical` and length one. This
#' `from_prior` flag can be used to differentiate the output of the `modelspecs`,
#' as its value is taken directly from the `from_prior` argument of the `mcmc`
#' method that invokes `modelspecs` function. That is, when `from_prior` is
#' `TRUE`, then only `priormodel` JAGS model is used (`datamodel` is not used)
#' by the `mcmc`, and hence `modelspecs` function should return all the parameters
#' that are required by the `priormodel` only. If the value of `from_prior` is
#' `FALSE`, then both JAGS models `datamodel` and `priormodel` are used in the
#' MCMC sampler, and hence `modelspecs` function should return all the parameters
#' required by both `datamodel` and `priormodel`.
#' @slot init (`function`)\cr a function computing the list of starting values
#' for parameters required to be initialized in the MCMC sampler, based on the
#' data slots that are required as arguments of this function.
#' @slot datanames (`character`)\cr the names of all data slots that are used
#' by `datamodel` JAGS function. No other names should be specified here.
#' @slot datanames_prior (`character`)\cr the names of all data slots that are
#' used by `priormodel` JAGS function. No other names should be specified here.
#' @slot sample (`character`)\cr names of all parameters from which you would
#' like to save the MCMC samples.
#'
#' @seealso [`ModelPseudo`].
#'
#' @aliases GeneralModel
#' @export
#'
.GeneralModel <- setClass(
Class = "GeneralModel",
slots = c(
datamodel = "function",
priormodel = "function",
modelspecs = "function",
init = "function",
datanames = "character",
datanames_prior = "character",
sample = "character"
),
prototype = prototype(
datamodel = I,
priormodel = I,
init = function() {
list()
}
),
contains = "CrmPackClass",
validity = v_general_model
)
## default constructor ----
#' @rdname GeneralModel-class
#' @note Typically, end users will not use the `.DefaultGeneralModel()` function.
#' @export
.DefaultGeneralModel <- function() {
stop(paste0("Class GeneralModel should not be instantiated directly. Please use one of its subclasses instead."))
}
# ModelLogNormal ----
## class ----
#' `ModelLogNormal`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`ModelLogNormal`] is the class for a model with a reference dose and bivariate
#' normal prior on the model parameters `alpha0` and natural logarithm of `alpha1`,
#' i.e.: \deqn{(alpha0, log(alpha1)) ~ Normal(mean, cov),}. Transformations other
#' than `log`, e.g. identity, can be specified too in `priormodel` slot.
#' The parameter `alpha1` has a log-normal distribution by default to ensure
#' positivity of `alpha1` which further guarantees `exp(alpha1) > 1`.
#' The slots of this class contain the mean vector, the covariance and
#' precision matrices of the bivariate normal distribution, as well as the
#' reference dose. Note that the precision matrix is an inverse of the
#' covariance matrix in the `JAGS`.
#' All ("normal") model specific classes inherit from this class.
#'
#' @slot params (`ModelParamsNormal`)\cr bivariate normal prior parameters.
#' @slot ref_dose (`positive_number`)\cr the reference dose.
#'
#' @seealso [`ModelParamsNormal`], [`LogisticNormal`], [`LogisticLogNormal`],
#' [`LogisticLogNormalSub`], [`ProbitLogNormal`], [`ProbitLogNormalRel`].
#'
#' @aliases ModelLogNormal
#' @export
#'
.ModelLogNormal <- setClass(
Class = "ModelLogNormal",
contains = "GeneralModel",
slots = c(
params = "ModelParamsNormal",
ref_dose = "positive_number"
)
)
## constructor ----
#' @rdname ModelLogNormal-class
#'
#' @param mean (`numeric`)\cr the prior mean vector.
#' @param cov (`matrix`)\cr the prior covariance matrix. The precision matrix
#' `prec` is internally calculated as an inverse of `cov`.
#' @param ref_dose (`number`)\cr the reference dose \eqn{x*} (strictly positive
#' number).
#'
#' @export
#'
ModelLogNormal <- function(mean, cov, ref_dose = 1) {
params <- ModelParamsNormal(mean, cov)
.ModelLogNormal(
params = params,
ref_dose = positive_number(ref_dose),
priormodel = function() {
theta ~ dmnorm(mean, prec)
alpha0 <- theta[1]
alpha1 <- exp(theta[2])
},
modelspecs = function(from_prior) {
ms <- list(mean = params@mean, prec = params@prec)
if (!from_prior) {
ms$ref_dose <- ref_dose
}
ms
},
init = function() {
list(theta = c(0, 1))
},
datanames = c("nObs", "y", "x"),
sample = c("alpha0", "alpha1")
)
}
## default constructor ----
#' @rdname ModelLogNormal-class
#' @note Typically, end users will not use the `.DefaultModelLogNormal()` function.
#' @export
.DefaultModelLogNormal <- function() {
ModelLogNormal(mean = c(-0.85, 1), cov = matrix(c(1, -0.5, -0.5, 1), nrow = 2))
}
# LogisticNormal ----
## class ----
#' `LogisticNormal`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`LogisticNormal`] is the class for the usual logistic regression model with
#' a bivariate normal prior on the intercept and slope.
#'
#' @details The covariate is the natural logarithm of the dose \eqn{x} divided by
#' the reference dose \eqn{x*}, i.e.:
#' \deqn{logit[p(x)] = alpha0 + alpha1 * log(x/x*),}
#' where \eqn{p(x)} is the probability of observing a DLT for a given dose \eqn{x}.
#' The prior \deqn{(alpha0, alpha1) ~ Normal(mean, cov).}
#'
#' @seealso [`ModelLogNormal`], [`LogisticLogNormal`], [`LogisticLogNormalSub`],
#' [`ProbitLogNormal`], [`ProbitLogNormalRel`], [`LogisticNormalMixture`].
#'
#' @aliases LogisticNormal
#' @export
#'
.LogisticNormal <- setClass(
Class = "LogisticNormal",
contains = "ModelLogNormal"
)
## constructor ----
#' @rdname LogisticNormal-class
#'
#' @inheritParams ModelLogNormal
#'
#' @export
#' @example examples/Model-class-LogisticNormal.R
#'
LogisticNormal <- function(mean, cov, ref_dose = 1) {
model_ln <- ModelLogNormal(mean = mean, cov = cov, ref_dose = ref_dose)
.LogisticNormal(
model_ln,
datamodel = function() {
for (i in 1:nObs) {
logit(p[i]) <- alpha0 + alpha1 * log(x[i] / ref_dose)
y[i] ~ dbern(p[i])
}
},
priormodel = function() {
theta ~ dmnorm(mean, prec)
alpha0 <- theta[1]
alpha1 <- theta[2]
}
)
}
## default constructor ----
#' @rdname LogisticNormal-class
#' @note Typically, end users will not use the `.DefaultLogisticNormal()` function.
#' @export
.DefaultLogisticNormal <- function() {
LogisticNormal(mean = c(-0.85, 1), cov = matrix(c(1, -0.5, -0.5, 1), nrow = 2))
}
# LogisticLogNormal ----
## class ----
#' `LogisticLogNormal`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`LogisticLogNormal`] is the class for the usual logistic regression model
#' with a bivariate normal prior on the intercept and log slope.
#'
#' @details The covariate is the natural logarithm of the dose \eqn{x} divided by
#' the reference dose \eqn{x*}, i.e.:
#' \deqn{logit[p(x)] = alpha0 + alpha1 * log(x/x*),}
#' where \eqn{p(x)} is the probability of observing a DLT for a given dose \eqn{x}.
#' The prior \deqn{(alpha0, log(alpha1)) ~ Normal(mean, cov).}
#'
#' @seealso [`ModelLogNormal`], [`LogisticNormal`], [`LogisticLogNormalSub`],
#' [`ProbitLogNormal`], [`ProbitLogNormalRel`], [`LogisticLogNormalMixture`],
#' [`DALogisticLogNormal`].
#'
#' @aliases LogisticLogNormal
#' @export
#'
.LogisticLogNormal <- setClass(
Class = "LogisticLogNormal",
contains = "ModelLogNormal"
)
## constructor ----
#' @rdname LogisticLogNormal-class
#'
#' @inheritParams ModelLogNormal
#'
#' @export
#' @example examples/Model-class-LogisticLogNormal.R
#'
LogisticLogNormal <- function(mean, cov, ref_dose = 1) {
model_ln <- ModelLogNormal(mean = mean, cov = cov, ref_dose = ref_dose)
.LogisticLogNormal(
model_ln,
datamodel = function() {
for (i in 1:nObs) {
logit(p[i]) <- alpha0 + alpha1 * log(x[i] / ref_dose)
y[i] ~ dbern(p[i])
}
}
)
}
## default constructor ----
#' @rdname LogisticLogNormal-class
#' @note Typically, end users will not use the `.DefaultLogisticLogNormal()` function.
#' @export
.DefaultLogisticLogNormal <- function() {
LogisticLogNormal(
mean = c(-0.85, 1),
cov = matrix(c(1, -0.5, -0.5, 1), nrow = 2),
ref_dose = 50
)
}
# LogisticLogNormalSub ----
## class ----
#' `LogisticLogNormalSub`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`LogisticLogNormalSub`] is the class for a standard logistic model with
#' bivariate (log) normal prior with subtractive dose standardization.
#'
#' @details The covariate is the dose \eqn{x} minus the reference dose \eqn{x*},
#' i.e.:
#' \deqn{logit[p(x)] = alpha0 + alpha1 * (x - x*),}
#' where \eqn{p(x)} is the probability of observing a DLT for a given dose \eqn{x}.
#' The prior \deqn{(alpha0, log(alpha1)) ~ Normal(mean, cov).}
#'
#' @slot params (`ModelParamsNormal`)\cr bivariate normal prior parameters.
#' @slot ref_dose (`number`)\cr the reference dose \eqn{x*}.
#'
#' @seealso [`LogisticNormal`], [`LogisticLogNormal`], [`ProbitLogNormal`],
#' [`ProbitLogNormalRel`].
#'
#' @aliases LogisticLogNormalSub
#' @export
#'
.LogisticLogNormalSub <- setClass(
Class = "LogisticLogNormalSub",
slots = c(
params = "ModelParamsNormal",
ref_dose = "numeric"
),
contains = "GeneralModel"
)
## constructor ----
#' @rdname LogisticLogNormalSub-class
#'
#' @param mean (`numeric`)\cr the prior mean vector.
#' @param cov (`matrix`)\cr the prior covariance matrix. The precision matrix
#' `prec` is internally calculated as an inverse of `cov`.
#' @param ref_dose (`number`)\cr the reference dose \eqn{x*}.
#'
#' @export
#' @example examples/Model-class-LogisticLogNormalSub.R
#'
LogisticLogNormalSub <- function(mean, cov, ref_dose = 0) {
params <- ModelParamsNormal(mean, cov)
.LogisticLogNormalSub(
params = params,
ref_dose = ref_dose,
datamodel = function() {
for (i in 1:nObs) {
logit(p[i]) <- alpha0 + alpha1 * (x[i] - ref_dose)
y[i] ~ dbern(p[i])
}
},
priormodel = function() {
theta ~ dmnorm(mean, prec)
alpha0 <- theta[1]
alpha1 <- exp(theta[2])
},
modelspecs = function(from_prior) {
ms <- list(mean = params@mean, prec = params@prec)
if (!from_prior) {
ms$ref_dose <- ref_dose
}
ms
},
init = function() {
list(theta = c(0, -20))
},
datanames = c("nObs", "y", "x"),
sample = c("alpha0", "alpha1")
)
}
## default constructor ----
#' @rdname LogisticLogNormalSub-class
#' @note Typically, end-users will not use the `.DefaultLogisticLogNormalSub()` function.
#' @export
.DefaultLogisticLogNormalSub <- function() {
LogisticLogNormalSub(
mean = c(-0.85, 1),
cov = matrix(c(1, -0.5, -0.5, 1), nrow = 2),
ref_dose = 50
)
}
# ProbitLogNormal ----
## class ----
#' `ProbitLogNormal`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`ProbitLogNormal`] is the class for probit regression model with a
#' bivariate normal prior on the intercept and log slope.
#'
#' @details The covariate is the natural logarithm of dose \eqn{x} divided by a
#' reference dose \eqn{x*}, i.e.:
#' \deqn{probit[p(x)] = alpha0 + alpha1 * log(x/x*),}
#' where \eqn{p(x)} is the probability of observing a DLT for a given dose \eqn{x}.
#' The prior \deqn{(alpha0, log(alpha1)) ~ Normal(mean, cov).}
#'
#' @note This model is also used in the [`DualEndpoint`] classes, so this class
#' can be used to check the prior assumptions on the dose-toxicity model, even
#' when sampling from the prior distribution of the dual endpoint model is not
#' possible.
#'
#' @seealso [`ModelLogNormal`], [`LogisticNormal`], [`LogisticLogNormal`],
#' [`LogisticLogNormalSub`], [`ProbitLogNormalRel`].
#'
#' @aliases ProbitLogNormalLogDose
#' @export
#'
.ProbitLogNormal <- setClass(
Class = "ProbitLogNormal",
contains = "ModelLogNormal"
)
## constructor ----
#' @rdname ProbitLogNormal-class
#'
#' @inheritParams ModelLogNormal
#'
#' @export
#' @example examples/Model-class-ProbitLogNormal.R
#'
ProbitLogNormal <- function(mean, cov, ref_dose = 1) {
model_ln <- ModelLogNormal(mean = mean, cov = cov, ref_dose = ref_dose)
.ProbitLogNormal(
model_ln,
datamodel = function() {
for (i in 1:nObs) {
probit(p[i]) <- alpha0 + alpha1 * log(x[i] / ref_dose)
y[i] ~ dbern(p[i])
}
}
)
}
## default constructor ----
#' @rdname ProbitLogNormal-class
#' @note Typically, end users will not use the `.DefaultProbitLogNormal()` function.
#' @export
.DefaultProbitLogNormal <- function() {
ProbitLogNormal(
mean = c(-0.85, 1),
cov = matrix(c(1, -0.5, -0.5, 1), nrow = 2),
ref_dose = 7.2
)
}
# ProbitLogNormalRel ----
## class ----
#' `ProbitLogNormalRel`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`ProbitLogNormalRel`] is the class for probit regression model with a bivariate
#' normal prior on the intercept and log slope.
#'
#' @details The covariate is the dose \eqn{x} divided by a reference dose \eqn{x*},
#' i.e.:
#' \deqn{probit[p(x)] = alpha0 + alpha1 * x/x*,}
#' where \eqn{p(x)} is the probability of observing a DLT for a given dose \eqn{x}.
#' The prior \deqn{(alpha0, log(alpha1)) ~ Normal(mean, cov).}
#'
#' @note This model is also used in the [`DualEndpoint`] classes, so this class
#' can be used to check the prior assumptions on the dose-toxicity model, even
#' when sampling from the prior distribution of the dual endpoint model is not
#' possible.
#'
#' @seealso [`ModelLogNormal`], [`LogisticNormal`], [`LogisticLogNormal`],
#' [`LogisticLogNormalSub`], [`ProbitLogNormal`].
#'
#' @aliases ProbitLogNormalRel
#' @export
#'
.ProbitLogNormalRel <- setClass(
Class = "ProbitLogNormalRel",
contains = "ModelLogNormal"
)
## constructor ----
#' @rdname ProbitLogNormalRel-class
#'
#' @inheritParams ModelLogNormal
#'
#' @export
#' @example examples/Model-class-ProbitLogNormalRel.R
#'
ProbitLogNormalRel <- function(mean, cov, ref_dose = 1) {
model_ln <- ModelLogNormal(mean = mean, cov = cov, ref_dose = ref_dose)
.ProbitLogNormalRel(
model_ln,
datamodel = function() {
for (i in 1:nObs) {
probit(p[i]) <- alpha0 + alpha1 * (x[i] / ref_dose)
y[i] ~ dbern(p[i])
}
}
)
}
## default constructor ----
#' @rdname ProbitLogNormalRel-class
#' @note Typically, end users will not use the `.DefaultProbitLogNormalRel()` function.
#' @export
.DefaultProbitLogNormalRel <- function() {
ProbitLogNormalRel(mean = c(-0.85, 1), cov = matrix(c(1, -0.5, -0.5, 1), nrow = 2))
}
# LogisticLogNormalGrouped ----
## class ----
#' `LogisticLogNormalGrouped`
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' [`LogisticLogNormalGrouped`] is the class for a logistic regression model
#' for both the mono and the combo arms of the simultaneous dose escalation
#' design.
#'
#' @details The continuous covariate is the natural logarithm of the dose \eqn{x} divided by
#' the reference dose \eqn{x*} as in [`LogisticLogNormal`]. In addition,
#' \eqn{I_c} is a binary indicator covariate which is 1 for the combo arm and 0 for the mono arm.
#' The model is then defined as:
#' \deqn{logit[p(x)] = (alpha0 + I_c * delta0) + (alpha1 + I_c * delta1) * log(x / x*),}
#' where \eqn{p(x)} is the probability of observing a DLT for a given dose \eqn{x},
#' and `delta0` and `delta1` are the differences in the combo arm compared to the mono intercept
#' and slope parameters `alpha0` and `alpha1`.
#' The prior is defined as \deqn{(alpha0, log(delta0), log(alpha1), log(delta1)) ~ Normal(mean, cov).}
#'
#' @seealso [`ModelLogNormal`], [`LogisticLogNormal`].
#'
#' @aliases LogisticLogNormalGrouped
#' @export
#'
.LogisticLogNormalGrouped <- setClass(
Class = "LogisticLogNormalGrouped",
contains = "ModelLogNormal"
)
## constructor ----
#' @rdname LogisticLogNormalGrouped-class
#'
#' @inheritParams ModelLogNormal
#'
#' @export
#' @example examples/Model-class-LogisticLogNormalGrouped.R
#'
LogisticLogNormalGrouped <- function(mean, cov, ref_dose = 1) {
params <- ModelParamsNormal(mean, cov)
.LogisticLogNormalGrouped(
params = params,
ref_dose = positive_number(ref_dose),
priormodel = function() {
theta ~ dmnorm(mean, prec)
alpha0 <- theta[1]
delta0 <- exp(theta[2])
alpha1 <- exp(theta[3])
delta1 <- exp(theta[4])
},
datamodel = function() {
for (i in 1:nObs) {
logit(p[i]) <- (alpha0 + is_combo[i] * delta0) +
(alpha1 + is_combo[i] * delta1) * log(x[i] / ref_dose)
y[i] ~ dbern(p[i])
}
},
modelspecs = function(group, from_prior) {
ms <- list(
mean = params@mean,
prec = params@prec
)
if (!from_prior) {
ms$ref_dose <- ref_dose
ms$is_combo <- as.integer(group == "combo")
}
ms
},
init = function() {
list(theta = c(0, 1, 1, 1))
},
datanames = c("nObs", "y", "x"),
sample = c("alpha0", "delta0", "alpha1", "delta1")
)
}
## default constructor ----
#' @rdname LogisticLogNormalGrouped-class
#' @note Typically, end users will not use the `.DefaultLogisticLogNormalGrouped()` function.
#' @export
.DefaultLogisticLogNormalGrouped <- function() {
LogisticLogNormalGrouped(
mean = rep(0, 4),
cov = diag(rep(1, 4)),
)
}
# LogisticKadane ----
## class ----
#' `LogisticKadane`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`LogisticKadane`] is the class for the logistic model in the parametrization
#' of Kadane et al. (1980).
#'
#' @details Let `rho0 = p(xmin)` be the probability of a DLT at the minimum dose
#' `xmin`, and let `gamma` be the dose with target toxicity probability `theta`,
#' i.e. \eqn{p(gamma) = theta}. Then it can easily be shown that the logistic
#' regression model has intercept
#' \deqn{[gamma * logit(rho0) - xmin * logit(theta)] / [gamma - xmin]}
#' and slope
#' \deqn{[logit(theta) - logit(rho0)] / [gamma - xmin].}
#'
#' The priors are \deqn{gamma ~ Unif(xmin, xmax).} and
#' \deqn{rho0 ~ Unif(0, theta).}
#'
#' @note The slots of this class, required for creating the model, are the target
#' toxicity, as well as the minimum and maximum of the dose range. Note that
#' these can be different from the minimum and maximum of the dose grid in the
#' data later on.
#'
#' @slot theta (`proportion`)\cr the target toxicity probability.
#' @slot xmin (`number`)\cr the minimum of the dose range.
#' @slot xmax (`number`)\cr the maximum of the dose range.
#'
#' @seealso [`ModelLogNormal`]
#'
#' @aliases LogisticKadane
#' @export
#'
.LogisticKadane <- setClass(
Class = "LogisticKadane",
contains = "GeneralModel",
slots = c(
theta = "numeric",
xmin = "numeric",
xmax = "numeric"
),
prototype = prototype(
theta = 0.3,
xmin = 0.1,
xmax = 1
),
validity = v_model_logistic_kadane
)
## constructor ----
#' @rdname LogisticKadane-class
#'
#' @param theta (`proportion`)\cr the target toxicity probability.
#' @param xmin (`number`)\cr the minimum of the dose range.
#' @param xmax (`number`)\cr the maximum of the dose range.
#'
#' @export
#' @example examples/Model-class-LogisticKadane.R
#'
LogisticKadane <- function(theta, xmin, xmax) {
.LogisticKadane(
theta = theta,
xmin = xmin,
xmax = xmax,
datamodel = function() {
for (i in 1:nObs) {
logit(p[i]) <- (1 / (gamma - xmin)) *
(gamma * logit(rho0) - xmin * logit(theta) + x[i] * (logit(theta) - logit(rho0)))
y[i] ~ dbern(p[i])
}
},
priormodel = function() {
rho0 ~ dunif(0, theta)
gamma ~ dunif(xmin, xmax)
},
modelspecs = function() {
list(theta = theta, xmin = xmin, xmax = xmax)
},
init = function() {
list(rho0 = theta / 10, gamma = (xmax - xmin) / 2)
},
datanames = c("nObs", "y", "x"),
sample = c("rho0", "gamma")
)
}
## default constructor ----
#' @rdname LogisticKadane-class
#' @note Typically, end-users will not use the `.DefaultLogisticKadane()` function.
#' @export
.DefaultLogisticKadane <- function() {
LogisticKadane(theta = 0.33, xmin = 1, xmax = 200)
}
# LogisticKadaneBetaGamma ----
## class ----
#' `LogisticKadaneBetaGamma`
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' [`LogisticKadaneBetaGamma`] is the class for the logistic model in the parametrization
#' of Kadane et al. (1980), using a beta and a gamma distribution as the model priors.
#'
#' @details Let `rho0 = p(xmin)` be the probability of a DLT at the minimum dose
#' `xmin`, and let `gamma` be the dose with target toxicity probability `theta`,
#' i.e. \eqn{p(gamma) = theta}. Then it can easily be shown that the logistic
#' regression model has intercept
#' \deqn{[gamma * logit(rho0) - xmin * logit(theta)] / [gamma - xmin]}
#' and slope
#' \deqn{[logit(theta) - logit(rho0)] / [gamma - xmin].}
#'
#' The prior for `gamma`, is \deqn{gamma ~ Gamma(shape, rate).}.
#' The prior for `rho0 = p(xmin)`, is \deqn{rho0 ~ Beta(alpha, beta).}
#'
#' @note The slots of this class, required for creating the model, are the same
#' as in the `LogisticKadane` class. In addition, the shape parameters of the
#' Beta prior distribution of `rho0` and the shape and rate parameters of the
#' Gamma prior distribution of `gamma`, are required for creating the prior model.
#'
#' @slot theta (`proportion`)\cr the target toxicity probability.
#' @slot xmin (`number`)\cr the minimum of the dose range.
#' @slot xmax (`number`)\cr the maximum of the dose range.
#' @slot alpha (`number`)\cr the first shape parameter of the Beta prior distribution
#' of `rho0 = p(xmin)` the probability of a DLT at the minimum dose `xmin`.
#' @slot beta (`number`)\cr the second shape parameter of the Beta prior distribution
#' of `rho0 = p(xmin)` the probability of a DLT at the minimum dose `xmin`.
#' @slot shape (`number`)\cr the shape parameter of the Gamma prior distribution
#' of `gamma` the dose with target toxicity probability `theta`.
#' @slot rate (`number`)\cr the rate parameter of the Gamma prior distribution
#' of `gamma` the dose with target toxicity probability `theta`.
#'
#' @seealso [`ModelLogNormal`], [`LogisticKadane`].
#'
#' @aliases LogisticKadaneBetaGamma
#' @export
#'
.LogisticKadaneBetaGamma <- setClass(
Class = "LogisticKadaneBetaGamma",
contains = "LogisticKadane",
slots = c(
alpha = "numeric",
beta = "numeric",
shape = "numeric",
rate = "numeric"
),
prototype = prototype(
theta = 0.3,
xmin = 0.1,
xmax = 1,
alpha = 1,
beta = 0.5,
shape = 1.2,
rate = 2.5
),
validity = v_model_logistic_kadane_beta_gamma
)
## constructor ----
#' @rdname LogisticKadaneBetaGamma-class
#'
#' @inheritParams LogisticKadane
#'
#' @param alpha (`number`)\cr the first shape parameter of the Beta prior distribution
#' `rho0 = p(xmin)` the probability of a DLT at the minimum dose `xmin`.
#' @param beta (`number`)\cr the second shape parameter of the Beta prior distribution
#' `rho0 = p(xmin)` the probability of a DLT at the minimum dose `xmin`.
#' @param shape (`number`)\cr the shape parameter of the Gamma prior distribution
#' `gamma` the dose with target toxicity probability `theta`.
#' @param rate (`number`)\cr the rate parameter of the Gamma prior distribution
#' `gamma` the dose with target toxicity probability `theta`.
#'
#' @export
#' @example examples/Model-class-LogisticKadaneBetaGamma.R
#'
LogisticKadaneBetaGamma <- function(theta, xmin, xmax, alpha, beta, shape, rate) {
model_lk <- LogisticKadane(theta = theta, xmin = xmin, xmax = xmax)
.LogisticKadaneBetaGamma(
model_lk,
alpha = alpha,
beta = beta,
shape = shape,
rate = rate,
priormodel = function() {
rho0 ~ dbeta(alpha, beta)
gamma ~ dgamma(shape, rate)
lowestdose <- xmin
highestdose <- xmax
DLTtarget <- theta
},
modelspecs = function() {
list(
theta = theta,
xmin = xmin,
xmax = xmax,
alpha = alpha,
beta = beta,
shape = shape,
rate = rate
)
}
)
}
## default constructor ----
#' @rdname LogisticKadaneBetaGamma-class
#' @note Typically, end users will not use the `.Default()` function.
#' @export
.DefaultLogisticKadaneBetaGamma <- function() {
LogisticKadaneBetaGamma(
theta = 0.3,
xmin = 0,
xmax = 7,
alpha = 1,
beta = 19,
shape = 0.5625,
rate = 0.125
)
}
# LogisticNormalMixture ----
## class ----
#' `LogisticNormalMixture`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`LogisticNormalMixture`] is the class for standard logistic regression model
#' with a mixture of two bivariate normal priors on the intercept and slope parameters.
#'
#' @details The covariate is the natural logarithm of the dose \eqn{x} divided by
#' the reference dose \eqn{x*}, i.e.:
#' \deqn{logit[p(x)] = alpha0 + alpha1 * log(x/x*),}
#' where \eqn{p(x)} is the probability of observing a DLT for a given dose \eqn{x}.
#' The prior
#' \deqn{(alpha0, alpha1) ~ w * Normal(mean1, cov1) + (1 - w) * Normal(mean2, cov2).}
#' The weight w for the first component is assigned a beta prior `B(a, b)`.
#'
#' @note The weight of the two normal priors is a model parameter, hence it is a
#' flexible mixture. This type of prior is often used with a mixture of a minimal
#' informative and an informative component, in order to make the CRM more robust
#' to data deviations from the informative component.
#'
#' @slot comp1 (`ModelParamsNormal`)\cr bivariate normal prior specification of
#' the first component.
#' @slot comp2 (`ModelParamsNormal`)\cr bivariate normal prior specification of
#' the second component.
#' @slot weightpar (`numeric`)\cr the beta parameters for the weight of the
#' first component. It must a be a named vector of length 2 with names `a` and
#' `b` and with strictly positive values.
#' @slot ref_dose (`positive_number`)\cr the reference dose.
#'
#' @seealso [`ModelParamsNormal`], [`ModelLogNormal`],
#' [`LogisticNormalFixedMixture`], [`LogisticLogNormalMixture`].
#'
#' @aliases LogisticNormalMixture
#' @export
#'
.LogisticNormalMixture <- setClass(
Class = "LogisticNormalMixture",
contains = "GeneralModel",
slots = c(
comp1 = "ModelParamsNormal",
comp2 = "ModelParamsNormal",
weightpar = "numeric",
ref_dose = "numeric"
),
prototype = prototype(
comp1 = ModelParamsNormal(mean = c(0, 1), cov = diag(2)),
comp2 = ModelParamsNormal(mean = c(-1, 1), cov = diag(2)),
weightpar = c(a = 1, b = 1),
ref_dose = 1
),
validity = v_model_logistic_normal_mix
)
## constructor ----
#' @rdname LogisticNormalMixture-class
#'
#' @param comp1 (`ModelParamsNormal`)\cr bivariate normal prior specification of
#' the first component. See [`ModelParamsNormal`] for more details.
#' @param comp2 (`ModelParamsNormal`)\cr bivariate normal prior specification of
#' the second component. See [`ModelParamsNormal`] for more details.
#' @param weightpar (`numeric`)\cr the beta parameters for the weight of the
#' first component. It must a be a named vector of length 2 with names `a` and
#' `b` and with strictly positive values.
#' @param ref_dose (`number`)\cr the reference dose \eqn{x*}
#' (strictly positive number).
#'
#' @export
#' @example examples/Model-class-LogisticNormalMixture.R
#'
LogisticNormalMixture <- function(comp1,
comp2,
weightpar,
ref_dose) {
assert_number(ref_dose)
.LogisticNormalMixture(
comp1 = comp1,
comp2 = comp2,
weightpar = weightpar,
ref_dose = ref_dose,
datamodel = function() {
# The logistic likelihood - the same as for non-mixture case.
for (i in 1:nObs) {
logit(p[i]) <- alpha0 + alpha1 * log(x[i] / ref_dose)
y[i] ~ dbern(p[i])
}
},
priormodel = function() {
w ~ dbeta(weightpar[1], weightpar[2])
wc <- 1 - w
comp0 ~ dbern(wc)
comp <- comp0 + 1
# Conditional on the component index "comp", which is 1 or 2.
# comp = 1 with probability "w" and comp = 2 with probability "1 - w".
theta ~ dmnorm(mean[1:2, comp], prec[1:2, 1:2, comp])
alpha0 <- theta[1]
alpha1 <- theta[2]
},
modelspecs = function(from_prior) {
ms <- list(
mean = cbind(comp1@mean, comp2@mean),
prec = array(data = c(comp1@prec, comp2@prec), dim = c(2, 2, 2)),
weightpar = weightpar
)
if (!from_prior) {
ms$ref_dose <- ref_dose
}
ms
},
init = function() {
list(theta = c(0, 1))
},
datanames = c("nObs", "y", "x"),
sample = c("alpha0", "alpha1", "w")
)
}
## default constructor ----
#' @rdname LogisticNormalMixture-class
#' @note Typically, end-users will not use the `.DefaultLogisticNormalMixture()` function.
#' @export
.DefaultLogisticNormalMixture <- function() { # nolint
LogisticNormalMixture(
comp1 = ModelParamsNormal(
mean = c(-0.85, 1),
cov = matrix(c(1, -0.5, -0.5, 1), nrow = 2)
),
comp2 = ModelParamsNormal(
mean = c(1, 1.5),
cov = matrix(c(1.2, -0.45, -0.45, 0.6), nrow = 2)
),
weightpar = c(a = 1, b = 1),
ref_dose = 50
)
}
# LogisticNormalFixedMixture ----
## class ----
#' `LogisticNormalFixedMixture`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`LogisticNormalFixedMixture`] is the class for standard logistic regression
#' model with fixed mixture of multiple bivariate (log) normal priors on the
#' intercept and slope parameters. The weights of the normal priors are fixed,
#' hence no additional model parameters are introduced. This type of prior is
#' often used to better approximate a given posterior distribution, or when the
#' information is given in terms of a mixture.
#'
#' @details The covariate is the natural logarithm of the dose \eqn{x} divided
#' by the reference dose \eqn{x*}, i.e.:
#' \deqn{logit[p(x)] = alpha0 + alpha1 * log(x/x*),}
#' where \eqn{p(x)} is the probability of observing a DLT for a given dose \eqn{x}.
#' The prior
#' \deqn{(alpha0, alpha1) ~ w1 * Normal(mean1, cov1) + ... + wK * Normal(meanK, covK),}
#' if a normal prior is used and
#' \deqn{(alpha0, log(alpha1)) ~ w1 * Normal(mean1, cov1) + ... + wK * Normal(meanK, covK),}
#' if a log normal prior is used.
#' The weights \eqn{w1, ..., wK} of the components are fixed and sum to 1.
#'
#' The slots of this class comprise a list with components parameters. Every
#' single component contains the mean vector and the covariance matrix of
#' bivariate normal distributions. Remaining slots are the weights of the
#' components as well as the reference dose. Moreover, a special indicator
#' slot specifies whether a log normal prior is used.
#'
#' @slot components (`list`)\cr the specifications of the mixture components,
#' a list with [`ModelParamsNormal`] objects for each bivariate (log) normal
#' prior.
#' @slot weights (`numeric`)\cr the weights of the components; these must be
#' positive and must sum to 1.
#' @slot ref_dose (`positive_number`)\cr the reference dose.
#' @slot log_normal (`flag`)\cr should a log normal prior be used, such
#' that the mean vectors and covariance matrices are valid for the intercept
#' and log slope?
#'
#' @seealso [`ModelParamsNormal`], [`ModelLogNormal`],
#' [`LogisticNormalMixture`], [`LogisticLogNormalMixture`].
#'
#' @aliases LogisticNormalFixedMixture
#' @export
#'
.LogisticNormalFixedMixture <- setClass(
Class = "LogisticNormalFixedMixture",
contains = "GeneralModel",
slots = c(
components = "list",
weights = "numeric",
ref_dose = "numeric",
log_normal = "logical"
),
prototype = prototype(
components = list(
comp1 = ModelParamsNormal(mean = c(0, 1), cov = diag(2)),
comp2 = ModelParamsNormal(mean = c(-1, 1), cov = diag(2))
),
weights = c(0.5, 0.5),
ref_dose = 1,
log_normal = FALSE
),
validity = v_model_logistic_normal_fixed_mix
)
## constructor ----
#' @rdname LogisticNormalFixedMixture-class
#'
#' @param components (`list`)\cr the specifications of the mixture components,
#' a list with [`ModelParamsNormal`] objects for each bivariate (log) normal
#' prior.
#' @param weights (`numeric`)\cr the weights of the components; these must be
#' positive and will be normalized to sum to 1.
#' @param ref_dose (`number`)\cr the reference dose \eqn{x*}
#' (strictly positive number).
#' @param log_normal (`flag`)\cr should a log normal prior be specified, such
#' that the mean vectors and covariance matrices are valid for the intercept
#' and log slope?
#'
#' @export
#' @example examples/Model-class-LogisticNormalFixedMixture.R
#'
LogisticNormalFixedMixture <- function(components,
weights,
ref_dose,
log_normal = FALSE) {
assert_numeric(weights)
assert_number(ref_dose)
assert_flag(log_normal)
# Normalize weights to sum to 1.
weights <- weights / sum(weights)
.LogisticNormalFixedMixture(
components = components,
weights = weights,
ref_dose = positive_number(ref_dose),
log_normal = log_normal,
datamodel = function() {
for (i in 1:nObs) {
logit(p[i]) <- alpha0 + alpha1 * log(x[i] / ref_dose)
y[i] ~ dbern(p[i])
}
},
priormodel = if (log_normal) {
function() {
comp ~ dcat(weights)
theta ~ dmnorm(mean[1:2, comp], prec[1:2, 1:2, comp])
alpha0 <- theta[1]
alpha1 <- exp(theta[2])
}
} else {
function() {
comp ~ dcat(weights)
theta ~ dmnorm(mean[1:2, comp], prec[1:2, 1:2, comp])
alpha0 <- theta[1]
alpha1 <- theta[2]
}
},
modelspecs = function(from_prior) {
ms <- list(
weights = weights,
mean = do.call(cbind, lapply(components, h_slots, "mean", simplify = TRUE)),
prec = array(
do.call(c, lapply(components, h_slots, "prec", simplify = TRUE)),
dim = c(2, 2, length(components))
)
)
if (!from_prior) {
ms$ref_dose <- ref_dose
}
ms
},
init = function() {
list(theta = c(0, 1))
},
datanames = c("nObs", "y", "x"),
sample = c("alpha0", "alpha1")
)
}
## default constructor ----
#' @rdname LogisticNormalFixedMixture-class
#' @note Typically, end-users will not use the `.DefaultLogisticNormalFixedMixture()`
#' function.
#' @export
.DefaultLogisticNormalFixedMixture <- function() { # nolint
LogisticNormalFixedMixture(
components = list(
comp1 = ModelParamsNormal(
mean = c(-0.85, 1),
cov = matrix(c(1, -0.5, -0.5, 1), nrow = 2)
),
comp2 = ModelParamsNormal(
mean = c(1, 1.5),
cov = matrix(c(1.2, -0.45, -0.45, 0.6), nrow = 2)
)
),
weights = c(0.3, 0.7),
ref_dose = 50
)
}
# LogisticLogNormalMixture ----
## class ----
#' `LogisticLogNormalMixture`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`LogisticLogNormalMixture`] is the class for standard logistic model with
#' online mixture of two bivariate log normal priors.
#'
#' @details This model can be used when data is arising online from the informative
#' component of the prior, at the same time with the data of the trial of
#' main interest. Formally, this is achieved by assuming that the probability
#' of a DLT at dose \eqn{x} is given by
#' \deqn{p(x) = \pi * p1(x) + (1 - \pi) * p2(x)}
#' where \eqn{\pi} is the probability for the model \eqn{p(x)} being the same
#' as the model \eqn{p1(x)}, which is the informative component of the prior.
#' From this model data arises in parallel: at doses `xshare`, DLT information
#' `yshare` is observed, in total `nObsshare` data points (see [`DataMixture`]).
#' On the other hand, \eqn{1 - \pi}, is the probability of a separate model
#' \eqn{p2(x)}. Both components have the same log normal prior distribution,
#' which can be specified by the user, and which is inherited from the
#' [`LogisticLogNormal`] class.
#'
#' @slot share_weight (`proportion`)\cr the prior weight for the share component
#' \eqn{p_{1}(x)}.
#'
#' @seealso [`ModelLogNormal`], [`LogisticNormalMixture`],
#' [`LogisticNormalFixedMixture`].
#'
#' @aliases LogisticLogNormalMixture
#' @export
#'
.LogisticLogNormalMixture <- setClass(
Class = "LogisticLogNormalMixture",
contains = "LogisticLogNormal",
slots = c(
share_weight = "numeric"
),
prototype = prototype(
share_weight = 0.1
),
validity = v_model_logistic_log_normal_mix
)
## constructor ----
#' @rdname LogisticLogNormalMixture-class
#'
#' @inheritParams ModelLogNormal
#' @param share_weight (`proportion`)\cr the prior weight for the share component.
#'
#' @export
#' @example examples/Model-class-LogisticLogNormalMixture.R
#'
LogisticLogNormalMixture <- function(mean,
cov,
ref_dose,
share_weight) {
assert_number(ref_dose)
params <- ModelParamsNormal(mean, cov)
.LogisticLogNormalMixture(
params = params,
ref_dose = positive_number(ref_dose),
share_weight = share_weight,
datamodel = function() {
for (i in 1:nObs) {
# comp gives the component: non-informative (1) or share (2) the two components.
stand_log_dose[i] <- log(x[i] / ref_dose)
logit(p[i]) <- alpha0[comp] + alpha1[comp] * stand_log_dose[i]
y[i] ~ dbern(p[i])
}
for (j in 1:nObsshare) {
stand_log_dose_share[j] <- log(xshare[j] / ref_dose)
logit(pshare[j]) <- alpha0[2] + alpha1[2] * stand_log_dose_share[j]
yshare[j] ~ dbern(pshare[j])
}
},
priormodel = function() {
for (k in 1:2) {
theta[k, 1:2] ~ dmnorm(mean, prec)
alpha0[k] <- theta[k, 1]
alpha1[k] <- exp(theta[k, 2])
}
# The component indicator.
comp ~ dcat(cat_probs)
},
modelspecs = function(from_prior) {
ms <- list(
cat_probs = c(1 - share_weight, share_weight),
mean = params@mean,
prec = params@prec
)
if (!from_prior) {
ms$ref_dose <- ref_dose
}
ms
},
init = function() {
list(theta = matrix(c(0, 0, 1, 1), nrow = 2))
},
datanames = c("nObs", "y", "x", "nObsshare", "yshare", "xshare"),
sample = c("alpha0", "alpha1", "comp")
)
}
## default constructor ----
#' @rdname LogisticLogNormalMixture-class
#' @note Typically, end users will not use the `.DefaultLogNormalMixture()` function.
#' @export
.DefaultLogisticLogNormalMixture <- function() { # nolint
LogisticLogNormalMixture(
share_weight = 0.1,
mean = c(-0.85, 1),
cov = matrix(c(1, -0.5, -0.5, 1), nrow = 2),
ref_dose = 50
)
}
# DualEndpoint ----
## class ----
#' `DualEndpoint`
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' [`DualEndpoint`] is the general class for the dual endpoint model.
#'
#' @details The idea of the dual-endpoint models is to model not only the
#' dose-toxicity relationship, but also to model, at the same time, the
#' relationship of a PD biomarker with the dose. The sub-classes of this class
#' define how the dose-biomarker relationship is parametrized. This class here
#' shall contain all the common features to reduce duplicate code.
#' (This class however, must not be virtual as we need to create objects
#' of it during the construction of subclass objects.)
#'
#' The dose-toxicity relationship is modeled with probit regression model
#' \deqn{probit[p(x)] = betaZ1 + betaZ2 * x/x*,}
#' or
#' \deqn{probit[p(x)] = betaZ1 + betaZ2 * log(x/x*),}
#' in case when the option `use_log_dose` is `TRUE`.
#' Here, \eqn{p(x)} is the probability of observing a DLT for a given
#' dose \eqn{x} and \eqn{x*} is the reference dose.
#' The prior \deqn{(betaZ1, log(betaZ2)) ~ Normal(mean, cov).}
#'
#' For the biomarker response \eqn{w} at a dose \eqn{x}, we assume
#' \deqn{w(x) ~ Normal(f(x), sigma2W),}
#' where \eqn{f(x)} is a function of the dose \eqn{x}, which is further
#' specified in sub-classes. The biomarker variance \eqn{sigma2W} can be fixed
#' or assigned an Inverse-Gamma prior distribution; see the details below under
#' slot `sigma2W`.
#'
#' Finally, the two endpoints \eqn{y} (the binary DLT variable) and \eqn{w}
#' (the biomarker) can be correlated, by assuming a correlation of level
#' \eqn{rho} between the underlying continuous latent toxicity variable \eqn{z}
#' and the biomarker \eqn{w}. Again, this correlation can be fixed or assigned
#' a prior distribution from the scaled Beta family; see the details below
#' under slot `rho`.
#'
#' Please see the example vignette by typing `crmPackExample()` for a full example.
#'
#' @slot betaZ_params (`ModelParamsNormal`)\cr for the probit toxicity model, it
#' contains the prior mean, covariance matrix and precision matrix which is
#' internally calculated as an inverse of the covariance matrix.
#' @slot ref_dose (`positive_number`)\cr for the probit toxicity model, the
#' reference dose.
#' @slot use_log_dose (`flag`)\cr for the probit toxicity model, whether a log
#' transformation of the (standardized) dose should be used?
#' @slot sigma2W (`numeric`)\cr the biomarker variance. Either a fixed value or
#' Inverse-Gamma distribution parameters, i.e. vector with two elements named
#' `a` and `b`.
#' @slot rho (`numeric`)\cr either a fixed value for the correlation
#' (between `-1` and `1`), or a named vector with two elements named `a` and `b`
#' for the Beta prior on the transformation `kappa = (rho + 1) / 2`, which is
#' in `(0, 1)`. For example, `a = 1, b = 1` leads to a uniform prior on `rho`.
#' @slot use_fixed (`logical`)\cr indicates whether a fixed value for `sigma2W`
#' or `rho` (for each parameter separately) is used or not. This slot is
#' needed for internal purposes and must not be touched by the user.
#'
#' @seealso [`DualEndpointRW`], [`DualEndpointBeta`], [`DualEndpointEmax`].
#'
#' @aliases DualEndpoint
#' @export
#'
.DualEndpoint <- setClass(
Class = "DualEndpoint",
slots = c(
betaZ_params = "ModelParamsNormal",
ref_dose = "positive_number",
use_log_dose = "logical",
sigma2W = "numeric",
rho = "numeric",
use_fixed = "logical"
),
prototype = prototype(
betaZ_params = ModelParamsNormal(mean = c(0, 1), cov = diag(2)),
ref_dose = positive_number(1),
use_log_dose = FALSE,
sigma2W = 1,
rho = 0,
use_fixed = c(sigma2W = TRUE, rho = TRUE)
),
contains = "GeneralModel",
validity = v_model_dual_endpoint
)
## constructor ----
#' @rdname DualEndpoint-class
#'
#' @param mean (`numeric`)\cr for the probit toxicity model, the prior mean vector.
#' @param cov (`matrix`)\cr for the probit toxicity model, the prior covariance
#' matrix. The precision matrix is internally calculated as an inverse of `cov`.
#' @param ref_dose (`number`)\cr for the probit toxicity model, the reference
#' dose \eqn{x*} (strictly positive number).
#' @param use_log_dose (`flag`)\cr for the probit toxicity model, whether a log
#' transformation of the (standardized) dose should be used?
#' @param sigma2W (`numeric`)\cr the biomarker variance. Either a fixed value or
#' Inverse-Gamma distribution parameters, i.e. vector with two elements named
#' `a` and `b`.
#' @param rho (`numeric`)\cr either a fixed value for the correlation
#' (between `-1` and `1`), or a named vector with two elements named `a` and `b`
#' for the Beta prior on the transformation `kappa = (rho + 1) / 2`, which is
#' in `(0, 1)`. For example, `a = 1, b = 1` leads to a uniform prior on `rho`.
#'
#' @export
#'
DualEndpoint <- function(mean,
cov,
ref_dose = 1,
use_log_dose = FALSE,
sigma2W,
rho) {
assert_number(ref_dose)
assert_numeric(sigma2W, min.len = 1, max.len = 2)
assert_numeric(rho, min.len = 1, max.len = 2)
use_fixed <- c(
sigma2W = test_number(sigma2W),
rho = test_number(rho)
)
beta_z_params <- ModelParamsNormal(mean, cov)
datamodel <- function() {
for (i in 1:nObs) {
# The toxicity model.
stand_dose_temp[i] <- x[i] / ref_dose
stand_dose[i] <- ifelse(use_log_dose, log(stand_dose_temp[i]), stand_dose_temp[i])
meanZ[i] <- betaZ[1] + betaZ[2] * stand_dose[i]
z[i] ~ dnorm(meanZ[i], 1)
y[i] ~ dinterval(z[i], 0)
# The conditional biomarker model; betaW defined in subclasses!
condMeanW[i] <- betaW[xLevel[i]] + rho / sqrt(precW) * (z[i] - meanZ[i])
w[i] ~ dnorm(condMeanW[i], condPrecW)
}
}
priormodel <- function() {
# Priors for betaW defined in subclasses!
theta ~ dmnorm(betaZ_mean, betaZ_prec)
betaZ[1] <- theta[1]
betaZ[2] <- exp(theta[2])
# Conditional precision for biomarker.
# Code for `precW` and `rho` will be added by
# `h_model_dual_endpoint_sigma2W()`, `h_model_dual_endpoint_rho()` helpers, below.
condPrecW <- precW / (1 - pow(rho, 2))
}
modelspecs_prior <- list(
betaZ_mean = beta_z_params@mean,
betaZ_prec = beta_z_params@prec
)
comp <- list(
priormodel = priormodel,
modelspecs = modelspecs_prior,
init = NULL,
sample = "betaZ"
)
# Update model components with regard to biomarker regression variance.
comp <- h_model_dual_endpoint_sigma2W(
use_fixed["sigma2W"],
sigma2W = sigma2W,
comp = comp
)
# Update model components with regard to DLT and biomarker correlation.
comp <- h_model_dual_endpoint_rho(
use_fixed["rho"],
rho = rho,
comp = comp
)
.DualEndpoint(
betaZ_params = beta_z_params,
ref_dose = positive_number(ref_dose),
use_log_dose = use_log_dose,
sigma2W = sigma2W,
rho = rho,
use_fixed = use_fixed,
datamodel = datamodel,
priormodel = comp$priormodel,
modelspecs = function(from_prior) {
if (!from_prior) {
comp$modelspecs$ref_dose <- ref_dose
comp$modelspecs$use_log_dose <- use_log_dose
}
comp$modelspecs
},
init = function(y) {
c(comp$init, list(z = ifelse(y == 0, -1, 1), theta = c(0, 1)))
},
datanames = c("nObs", "w", "x", "xLevel", "y"),
sample = comp$sample
)
}
## default constructor ----
#' @rdname DualEndpoint-class
#' @note Typically, end users will not use the `.DefaultDualEndpoint()` function.
#' @export
.DefaultDualEndpoint <- function() {
stop(paste0("Class DualEndpoint cannot be instantiated directly. Please use one of its subclasses instead."))
}
# DualEndpointRW ----
## class ----
#' `DualEndpointRW`
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' [`DualEndpointRW`] is the class for the dual endpoint model with random walk
#' prior for biomarker.
#'
#'
#' @details This class extends the [`DualEndpoint`] class so that the dose-biomarker
#' relationship \eqn{f(x)} is modelled by a non-parametric random walk of first
#' or second order. That means, for the first order random walk we assume
#' \deqn{betaW_i - betaW_i-1 ~ Normal(0, (x_i - x_i-1) * sigma2betaW),}
#' where \eqn{betaW_i = f(x_i)} is the biomarker mean at the \eqn{i}-th dose
#' gridpoint \eqn{x_i}.
#' For the second order random walk, the second-order differences instead of
#' the first-order differences of the biomarker means follow the normal distribution
#' with \eqn{0} mean and \eqn{2 * (x_i - x_i-2) * sigma2betaW} variance.
#'
#' The variance parameter \eqn{sigma2betaW} is important because it steers the
#' smoothness of the function \eqn{f(x)}, i.e.: if it is large, then \eqn{f(x)}
#' will be very wiggly; if it is small, then \eqn{f(x)} will be smooth.
#' This parameter can either be a fixed value or assigned an inverse gamma prior
#' distribution.
#'
#' @note Non-equidistant dose grids can be used now, because the difference
#' \eqn{x_i - x_i-1} is included in the modelling assumption above.
#' Please note that due to impropriety of the random walk prior distributions,
#' it is not possible to produce MCMC samples with empty data objects (i.e.,
#' sample from the prior). This is not a bug, but a theoretical feature of this
#' model.
#'
#' @slot sigma2betaW (`numeric`)\cr the prior variance factor of the random walk
#' prior for the biomarker model. Either a fixed value or Inverse-Gamma distribution
#' parameters, i.e. vector with two elements named `a` and `b`.
#' @slot rw1 (`flag`)\cr for specifying the random walk prior on the biomarker
#' level. When `TRUE`, random walk of first order is used. Otherwise, the
#' random walk of second order is used.
#'
#' @seealso [`DualEndpoint`], [`DualEndpointBeta`], [`DualEndpointEmax`].
#'
#' @aliases DualEndpointRW
#' @export
#'
.DualEndpointRW <- setClass(
Class = "DualEndpointRW",
slots = c(
sigma2betaW = "numeric",
rw1 = "logical"
),
prototype = prototype(
sigma2betaW = 1,
rw1 = TRUE,
use_fixed = c(
sigma2W = TRUE,
rho = TRUE,
sigma2betaW = TRUE
)
),
contains = "DualEndpoint",
validity = v_model_dual_endpoint_rw
)
## constructor ----
#' @rdname DualEndpointRW-class
#'
#' @param sigma2betaW (`numeric`)\cr the prior variance factor of the random walk
#' prior for the biomarker model. Either a fixed value or Inverse-Gamma distribution
#' parameters, i.e. vector with two elements named `a` and `b`.
#' @param rw1 (`flag`)\cr for specifying the random walk prior on the biomarker
#' level. When `TRUE`, random walk of first order is used. Otherwise, the
#' random walk of second order is used.
#' @param ... parameters passed to [DualEndpoint()].
#'
#' @export
#' @example examples/Model-class-DualEndpointRW.R
#'
DualEndpointRW <- function(sigma2betaW,
rw1 = TRUE,
...) {
assert_numeric(sigma2betaW, min.len = 1, max.len = 2)
assert_flag(rw1)
start <- DualEndpoint(...)
start@use_fixed["sigma2betaW"] <- length(sigma2betaW) == 1L
priormodel <- if (rw1) {
function() {
# The 1st order differences.
# Essentially dflat(), which is not available in JAGS.
betaW[1] ~ dnorm(0, 0.000001)
for (i in 2:nGrid) {
delta[i - 1] ~ dnorm(0, precBetaW / (doseGrid[i] - doseGrid[i - 1]))
betaW[i] <- betaW[i - 1] + delta[i - 1]
}
}
} else {
function() {
# The 2nd order differences.
delta[1] ~ dnorm(0, 0.000001)
betaW[1] ~ dnorm(0, 0.000001)
betaW[2] <- betaW[1] + delta[1]
for (i in 3:nGrid) {
# delta2: differences of the differences of betaW follow normal dist.
delta2[i - 2] ~ dnorm(0, 2 * precBetaW / (doseGrid[i] - doseGrid[i - 2]))
delta[i - 1] <- delta[i - 2] + delta2[i - 2]
betaW[i] <- betaW[i - 1] + delta[i - 1]
}
}
}
start@priormodel <- h_jags_join_models(start@priormodel, priormodel)
start@datanames_prior <- c("nGrid", "doseGrid")
start@sample <- c(start@sample, "betaW", "delta")
# Update model components with regard to biomarker regression variance.
start <- h_model_dual_endpoint_sigma2betaW(
start@use_fixed["sigma2betaW"],
sigma2betaW = sigma2betaW,
de = start
)
.DualEndpointRW(
start,
sigma2betaW = sigma2betaW,
rw1 = rw1
)
}
## default constructor ----
#' @rdname DualEndpointRW-class
#' @note Typically, end users will not use the `.DefaultDualEndpointRW()` function.
#' @export
.DefaultDualEndpointRW <- function() {
DualEndpointRW(
mean = c(0, 1),
cov = matrix(c(1, 0, 0, 1), nrow = 2),
sigma2W = c(a = 0.1, b = 0.1),
rho = c(a = 1, b = 1),
sigma2betaW = 0.01,
rw1 = TRUE
)
}
# DualEndpointBeta ----
## class ----
#' `DualEndpointBeta`
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' [`DualEndpointBeta`] is the class for the dual endpoint model with beta
#' function for dose-biomarker relationship.
#'
#' @details This class extends the [`DualEndpoint`] class so that the dose-biomarker
#' relationship \eqn{f(x)} is modelled by a parametric, rescaled beta density
#' function:
#' \deqn{f(x) = E0 + (Emax - E0) * Beta(delta1, delta2) * (x/x*)^{delta1} * (1 - x/x*)^{delta2},}
#' where \eqn{x*} is the maximum dose (end of the dose range to be considered),
#' \eqn{delta1} and \eqn{delta2} are the two beta function parameters, and
#' \eqn{E0}, \eqn{Emax} are the minimum and maximum levels, respectively.
#' For ease of interpretation, we use the parametrization based on \eqn{delta1}
#' and the mode, where
#' \deqn{mode = delta1 / (delta1 + delta2),}
#' so that multiplying this by \eqn{x*} gives the mode on the dose grid.
#'
#' All parameters can currently be assigned uniform distributions or be fixed
#' in advance. Note that \code{E0} and \code{Emax} can have negative values or
#' uniform distributions reaching into negative range, while \code{delta1} and
#' \code{mode} must be positive or have uniform distributions in the positive
#' range.
#'
#' @slot E0 (`numeric`)\cr either a fixed number or the two uniform distribution
#' parameters.
#' @slot Emax (`numeric`)\cr either a fixed number or the two uniform
#' distribution parameters.
#' @slot delta1 (`numeric`)\cr either a fixed positive number or the two
#' parameters of the uniform distribution, that can take only positive values.
#' @slot mode (`numeric`)\cr either a fixed positive number or the two
#' parameters of the uniform distribution, that can take only positive values.
#' @slot ref_dose_beta (`positive_number`)\cr the reference dose \eqn{x*} (note
#' that this is different from the `ref_dose` in the inherited [`DualEndpoint`]
#' model).
#'
#' @seealso [`DualEndpoint`], [`DualEndpointRW`], [`DualEndpointEmax`].
#'
#' @aliases DualEndpointBeta
#' @export
#'
.DualEndpointBeta <- setClass(
Class = "DualEndpointBeta",
slots = c(
E0 = "numeric",
Emax = "numeric",
delta1 = "numeric",
mode = "numeric",
ref_dose_beta = "positive_number"
),
prototype = prototype(
E0 = c(0, 100),
Emax = c(0, 500),
delta1 = c(0, 5),
mode = c(1, 15),
ref_dose_beta = positive_number(1),
use_fixed = c(
sigma2W = TRUE,
rho = TRUE,
E0 = FALSE,
Emax = FALSE,
delta1 = FALSE,
mode = FALSE
)
),
contains = "DualEndpoint",
validity = v_model_dual_endpoint_beta
)
## constructor ----
#' @rdname DualEndpointBeta-class
#'
#' @param E0 (`numeric`)\cr either a fixed number or the two uniform distribution
#' parameters.
#' @param Emax (`numeric`)\cr either a fixed number or the two uniform distribution
#' parameters.
#' @param delta1 (`numeric`)\cr either a fixed positive number or the two parameters
#' of the uniform distribution, that can take only positive values.
#' @param mode (`numeric`)\cr either a fixed positive number or the two parameters
#' of the uniform distribution, that can take only positive values.
#' @param ref_dose_beta (`number`)\cr the reference dose \eqn{x*} (strictly
#' positive number). Note that this is different from the `ref_dose` in the
#' inherited [`DualEndpoint`] model).
#' @param ... parameters passed to [DualEndpoint()].
#'
#' @export
#' @example examples/Model-class-DualEndpointBeta.R
#'
DualEndpointBeta <- function(E0,
Emax,
delta1,
mode,
ref_dose_beta = 1,
...) {
assert_numeric(E0, min.len = 1, max.len = 2)
assert_numeric(Emax, min.len = 1, max.len = 2)
assert_numeric(delta1, min.len = 1, max.len = 2)
assert_numeric(mode, min.len = 1, max.len = 2)
assert_number(ref_dose_beta)
start <- DualEndpoint(...)
ms <- start@modelspecs
start@modelspecs <- function(from_prior) {
c(list(ref_dose_beta = ref_dose_beta), ms(from_prior))
}
start@datanames_prior <- c("nGrid", "doseGrid")
start@sample <- c(start@sample, "betaW")
start <- h_model_dual_endpoint_beta(
param = E0,
param_name = "E0",
priormodel = function() {
E0 ~ dunif(E0_low, E0_high)
},
de = start
)
start <- h_model_dual_endpoint_beta(
param = Emax,
param_name = "Emax",
priormodel = function() {
Emax ~ dunif(Emax_low, Emax_high)
},
de = start
)
start <- h_model_dual_endpoint_beta(
param = delta1,
param_name = "delta1",
priormodel = function() {
delta1 ~ dunif(delta1_low, delta1_high)
},
de = start
)
start <- h_model_dual_endpoint_beta(
param = mode,
param_name = "mode",
priormodel = function() {
mode ~ dunif(mode_low, mode_high)
},
de = start
)
start@priormodel <- h_jags_join_models(
start@priormodel,
function() {
# delta2 <- delta1 * (1 - (mode/ref_dose_beta)) / (mode/ref_dose_beta) # nolint
delta2 <- delta1 * (ref_dose_beta / mode - 1)
# betafun <- (delta1 + delta2)^(delta1 + delta2) * delta1^(- delta1) * delta2^(- delta2) # nolint
betafun <- (1 + delta2 / delta1)^delta1 * (delta1 / delta2 + 1)^delta2
for (i in 1:nGrid) {
stand_dose_beta[i] <- doseGrid[i] / ref_dose_beta
betaW[i] <- E0 + (Emax - E0) * betafun * stand_dose_beta[i]^delta1 * (1 - stand_dose_beta[i])^delta2
}
}
)
.DualEndpointBeta(
start,
E0 = E0,
Emax = Emax,
delta1 = delta1,
mode = mode,
ref_dose_beta = positive_number(ref_dose_beta)
)
}
## default constructor ----
#' @rdname DualEndpointBeta-class
#' @note Typically, end users will not use the `.DefaultDualEndpointBeta()` function.
#' @export
.DefaultDualEndpointBeta <- function() {
DualEndpointBeta(
mean = c(0, 1),
cov = matrix(c(1, 0, 0, 1), nrow = 2),
ref_dose = 10,
use_log_dose = TRUE,
sigma2W = c(a = 0.1, b = 0.1),
rho = c(a = 1, b = 1),
E0 = c(0, 100),
Emax = c(0, 500),
delta1 = c(0, 5),
mode = c(1, 15),
ref_dose_beta = 1000
)
}
# DualEndpointEmax ----
## class ----
#' `DualEndpointEmax`
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' [`DualEndpointEmax`] is the class for the dual endpoint model with `Emax`
#' function for dose-biomarker relationship.
#'
#' @details This class extends the [`DualEndpoint`] class so that the dose-biomarker
#' relationship \eqn{f(x)} is modelled by a parametric `Emax` function:
#' \deqn{f(x) = E0 + [(Emax - E0) * (x/x*)]/[ED50 + (x/x*)],}
#' where \eqn{x*} is a reference dose, \eqn{E0} and \eqn{Emax} are the minimum
#' and maximum levels for the biomarker, and \eqn{ED50} is the dose achieving
#' half of the maximum effect \eqn{0.5 * Emax}.
#' All parameters can currently be assigned uniform distributions or be fixed.
#'
#' @slot E0 (`numeric`)\cr either a fixed number or the two uniform distribution
#' parameters.
#' @slot Emax (`numeric`)\cr either a fixed number or the two uniform
#' distribution parameters.
#' @slot ED50 (`numeric`)\cr either a fixed number or the two uniform
#' distribution parameters.
#' @slot ref_dose_emax (`positive_number`)\cr the reference dose \eqn{x*} (note
#' that this is different from the `ref_dose` in the inherited [`DualEndpoint`]
#' model).
#'
#' @seealso [`DualEndpoint`], [`DualEndpointRW`], [`DualEndpointBeta`].
#'
#' @aliases DualEndpointEmax
#' @export
#'
.DualEndpointEmax <- setClass(
Class = "DualEndpointEmax",
slots = c(
E0 = "numeric",
Emax = "numeric",
ED50 = "numeric",
ref_dose_emax = "numeric"
),
prototype = prototype(
E0 = c(0, 100),
Emax = c(0, 500),
ED50 = c(0, 500),
ref_dose_emax = positive_number(1),
use_fixed = c(
sigma2W = TRUE,
rho = TRUE,
E0 = FALSE,
Emax = FALSE,
ED50 = FALSE
)
),
contains = "DualEndpoint",
validity = v_model_dual_endpoint_emax
)
## constructor ----
#' @rdname DualEndpointEmax-class
#'
#' @param E0 (`numeric`)\cr either a fixed number or the two uniform distribution
#' parameters.
#' @param Emax (`numeric`)\cr either a fixed number or the two uniform distribution
#' parameters.
#' @param ED50 (`numeric`)\cr either a fixed number or the two uniform distribution
#' parameters.
#' @param ref_dose_emax (`number`)\cr the reference dose \eqn{x*} (strictly
#' positive number). Note that this is different from the `ref_dose` in the
#' inherited [`DualEndpoint`] model).
#' @param ... parameters passed to [DualEndpoint()].
#'
#' @export
#' @example examples/Model-class-DualEndpointEmax.R
#'
DualEndpointEmax <- function(E0,
Emax,
ED50,
ref_dose_emax = 1,
...) {
assert_numeric(E0, min.len = 1, max.len = 2)
assert_numeric(Emax, min.len = 1, max.len = 2)
assert_numeric(ED50, min.len = 1, max.len = 2)
assert_number(ref_dose_emax)
start <- DualEndpoint(...)
start@sample <- c(start@sample, "betaW")
start@datanames_prior <- c("nGrid", "doseGrid")
ms <- start@modelspecs
start@modelspecs <- function(from_prior) {
c(list(ref_dose_emax = ref_dose_emax), ms(from_prior))
}
start <- h_model_dual_endpoint_beta(
param = E0,
param_name = "E0",
priormodel = function() {
E0 ~ dunif(E0_low, E0_high)
},
de = start
)
start <- h_model_dual_endpoint_beta(
param = Emax,
param_name = "Emax",
priormodel = function() {
Emax ~ dunif(Emax_low, Emax_high)
},
de = start
)
start <- h_model_dual_endpoint_beta(
param = ED50,
param_name = "ED50",
priormodel = function() {
ED50 ~ dunif(ED50_low, ED50_high)
},
de = start
)
start@priormodel <- h_jags_join_models(
start@priormodel,
function() {
for (i in 1:nGrid) {
stand_dose_emax[i] <- doseGrid[i] / ref_dose_emax
betaW[i] <- E0 + (Emax - E0) * stand_dose_emax[i] / (ED50 + stand_dose_emax[i])
}
}
)
.DualEndpointEmax(
start,
E0 = E0,
Emax = Emax,
ED50 = ED50,
ref_dose_emax = positive_number(ref_dose_emax)
)
}
## default constructor ----
#' @rdname DualEndpointEmax-class
#' @note Typically, end users will not use the `.DefaultDualEndpointEmax()` function.
#' @export
.DefaultDualEndpointEmax <- function() {
DualEndpointEmax(
mean = c(0, 1),
cov = matrix(c(1, 0, 0, 1), nrow = 2),
sigma2W = c(a = 0.1, b = 0.1),
rho = c(a = 1, b = 1),
E0 = c(0, 100),
Emax = c(0, 500),
ED50 = c(10, 200),
ref_dose_emax = 1000
)
}
# ModelPseudo ----
## class ----
#' `ModelPseudo`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`ModelPseudo`] is the parent class for models that express their prior in
#' the form of pseudo data (as if there is some data before the trial starts).
#'
#' @seealso [`GeneralModel`].
#'
#' @aliases ModelPseudo
#' @export
#'
.ModelPseudo <- setClass(
Class = "ModelPseudo",
contains = "CrmPackClass"
)
## default constructor ----
#' @rdname ModelPseudo-class
#' @note Typically, end users will not use the `.DefaultModelPseudo()` function.
#' @export
.DefaultModelPseudo <- function() {
stop(paste0("Class ModelPseudo should not be instantiated directly. Please use one of its subclasses instead."))
}
# ModelTox ----
## class ----
#' `ModelTox`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`ModelTox`] is the parent class for DLE (dose-limiting events) models using
#' pseudo data prior. It is dedicated for DLE models or toxicity models that
#' have their prior specified in the form of pseudo data (as if there is some
#' data before the trial starts).
#'
#' The `data` must obey the convention of the [`Data`] class. This refers to any
#' observed DLE responses (`y` in [`Data`]), the dose levels (`x` in [`Data`])
#' at which these responses are observed, all dose levels considered in the
#' study (`doseGrid` in [`Data`]), and finally other specifications in [`Data`]
#' class that can be used to generate prior or posterior modal estimates or
#' samples estimates for model parameter(s).
#' If no responses are observed, at least `doseGrid` has to be specified
#' in `data` for which prior modal estimates or samples can be obtained for
#' model parameters based on the specified pseudo data.
#'
#' @slot data (`Data`)\cr observed data that is used to obtain model parameters
#' estimates or samples (see details above).
#'
#' @seealso [`ModelEff`].
#'
#' @aliases ModelTox
#' @export
#'
.ModelTox <- setClass(
Class = "ModelTox",
slots = c(
data = "Data"
),
contains = "ModelPseudo"
)
## default constructor ----
#' @rdname ModelTox-class
#' @note Typically, end users will not use the `.DefaultModelTox()` function.
#' @export
.DefaultModelTox <- function() {
stop(paste0("Class ModelTox should not be instantiated directly. Please use one of its subclasses instead."))
}
# ModelEff ----
## class ----
#' `ModelEff`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`ModelEff`] is the parent class for efficacy models using pseudo data prior.
#' It is dedicated all efficacy models that have their prior specified in the
#' form of pseudo data (as if there is some data before the trial starts).
#'
#' The `data` must obey the convention of the [`DataDual`] class. This refers to
#' any observed efficacy/biomarker responses (`w` in [`DataDual`]), the dose
#' levels at which these responses are observed (`x` in [`DataDual`]), all dose
#' levels considered in the study (`doseGrid` in [`DataDual`]), and finally
#' other specifications in [`DataDual`] class that can be used to generate prior
#' or posterior modal estimates or samples estimates for model parameter(s).
#' If no responses are observed, at least `doseGrid` has to be specified
#' in `data` for which prior modal estimates or samples can be obtained for
#' model parameters based on the specified pseudo data.
#'
#' @slot data (`DataDual`)\cr observed data that is used to obtain model
#' parameters estimates or samples (see details above).
#'
#' @seealso [`ModelTox`].
#'
#' @aliases ModelEff
#' @export
#'
.ModelEff <- setClass(
Class = "ModelEff",
slots = c(
data = "DataDual"
),
contains = "ModelPseudo"
)
## default constructor ----
#' @rdname ModelEff-class
#' @note Typically, end users will not use the `.DefaultModelEff()` function.
#' @export
.DefaultModelEff <- function() {
stop(paste0("Class ModelEff should not be instantiated directly. Please use one of its subclasses instead."))
}
# LogisticIndepBeta ----
## class ----
#' `LogisticIndepBeta`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`LogisticIndepBeta`] is the class for the two-parameters logistic regression
#' dose-limiting events (DLE) model with prior expressed in form of pseudo data.
#' This model describes the relationship between the binary DLE responses
#' and the dose levels. More specifically, it represents the relationship of the
#' probabilities of the occurrence of a DLE for corresponding dose levels in log
#' scale. This model is specified as
#' \deqn{p(x) = exp(phi1 + phi2 * log(x)) / (1 + exp(phi1 + phi2 * log(x)))}
#' where \eqn{p(x)} is the probability of the occurrence of a DLE at dose \eqn{x}.
#' The two parameters of this model are the intercept \eqn{phi1} and the slope
#' \eqn{phi2}. The `LogisticIndepBeta` inherits all slots from [`ModelTox`] class.
#'
#' In the context of pseudo data, the following three arguments are used,
#' `binDLE`, `DLEdose` and `DLEweights`. The `DLEdose` represents fixed dose
#' levels at which the pseudo DLE responses `binDLE` are observed. `DLEweights`
#' represents total number of subjects treated per each dose level in `DLEdose`.
#' The `binDLE` represents the number of subjects observed with DLE per each
#' dose level in `DLEdose`. Hence, all these three vectors must be of the same
#' length and the order of the elements in any of the vectors `binDLE`,
#' `DLEdose` and `DLEweights` must be kept, so that an element of a given vector
#' corresponds to the elements of the remaining two vectors (see the example for
#' more insight).
#' Finally, since at least two DLE pseudo responses are needed to
#' obtain prior modal estimates (same as the maximum likelihood estimates) for
#' the model parameters, the `binDLE`, `DLEdose` and `DLEweights` must all be
#' vectors of at least length 2.
#'
#' @details The pseudo data can be interpreted as if we obtain some observations
#' before the trial starts. It can be used to express our prior, i.e. the
#' initial beliefs for the model parameters. The pseudo data is expressed in
#' the following way. First, fix at least two dose levels, then ask for experts'
#' opinion on how many subjects are to be treated at each of these dose levels
#' and on the number of subjects observed with a DLE. At each dose level, the
#' number of subjects observed with a DLE, divided by the total number of
#' subjects treated, is the probability of the occurrence of a DLE at that
#' particular dose level. The probabilities of the occurrence of a DLE based
#' on this pseudo data are independent and they follow Beta distributions.
#' Therefore, the joint prior probability density function of all these
#' probabilities can be obtained. Hence, by a change of variable, the joint
#' prior probability density function of the two parameters in this model can
#' also be obtained. In addition, a conjugate joint prior density function of
#' the two parameters in the model is used. For details about the form of all
#' these joint prior and posterior probability density functions, please refer
#' to Whitehead and Willamson (1998).
#'
#' @slot binDLE (`numeric`)\cr a vector of total numbers of DLE responses.
#' It must be at least of length 2 and the order of its elements must
#' correspond to values specified in `DLEdose` and `DLEweights`.
#' @slot DLEdose (`numeric`)\cr a vector of the dose levels corresponding to
#' It must be at least of length 2 and the order of its elements must
#' correspond to values specified in `binDLE` and `DLEweights`.
#' @slot DLEweights (`integer`)\cr total number of subjects treated at each of
#' the pseudo dose level `DLEdose`.
#' It must be at least of length 2 and the order of its elements must
#' correspond to values specified in `binDLE` and `DLEdose`.
#' @slot phi1 (`number`)\cr the intercept of the model. This slot is used in
#' output to display the resulting prior or posterior modal estimate of the
#' intercept obtained based on the pseudo data and (if any) observed data/responses.
#' @slot phi2 (`number`)\cr the slope of the model. This slot is used in output
#' to display the resulting prior or posterior modal estimate of the slope
#' obtained based on the pseudo data and (if any) the observed data/responses.
#' @slot Pcov (`matrix`)\cr refers to the 2x2 covariance matrix of the intercept
#' (\eqn{phi1}) and the slope parameters (\eqn{phi2}) of the model.
#' This is used in output to display the resulting prior and posterior
#' covariance matrix of \eqn{phi1} and \eqn{phi2} obtained, based on the
#' pseudo data and (if any) the observed data and responses. This slot is
#' needed for internal purposes.
#'
#' @aliases LogisticIndepBeta
#' @export
#'
.LogisticIndepBeta <- setClass(
Class = "LogisticIndepBeta",
slots = c(
binDLE = "numeric",
DLEdose = "numeric",
DLEweights = "integer",
phi1 = "numeric",
phi2 = "numeric",
Pcov = "matrix"
),
prototype = prototype(
binDLE = c(0, 0),
DLEdose = c(1, 1),
DLEweights = c(1L, 1L)
),
contains = "ModelTox",
validity = v_model_logistic_indep_beta
)
## constructor ----
#' @rdname LogisticIndepBeta-class
#'
#' @param binDLE (`numeric`)\cr the number of subjects observed with a DLE, the
#' pseudo DLE responses, depending on dose levels `DLEdose`.
#' Elements of `binDLE` must correspond to the elements of `DLEdose` and
#' `DLEweights`.
#' @param DLEdose (`numeric`)\cr dose levels for the pseudo DLE responses.
#' Elements of `DLEdose` must correspond to the elements of `binDLE` and
#' `DLEweights`.
#' @param DLEweights (`numeric`)\cr the total number of subjects treated at each
#' of the dose levels `DLEdose`, pseudo weights.
#' Elements of `DLEweights` must correspond to the elements of `binDLE` and
#' `DLEdose`.
#' @param data (`Data`)\cr the input data to update estimates of the model
#' parameters.
#'
#' @export
#' @example examples/Model-class-LogisticIndepBeta.R
#'
LogisticIndepBeta <- function(binDLE,
DLEdose,
DLEweights,
data) {
assert_numeric(binDLE)
assert_numeric(DLEdose)
assert_integerish(DLEweights, lower = 0, any.missing = FALSE)
assert_class(data, "Data")
# Combine pseudo and observed data. It can also happen that data@nObs == 0.
y <- c(binDLE, data@y)
x <- c(DLEdose, data@x)
w <- c(DLEweights, rep(1, data@nObs))
fit_dle <- suppressWarnings(
glm(y / w ~ log(x), family = binomial(link = "logit"), weights = w)
)
phi1 <- coef(fit_dle)[["(Intercept)"]]
phi2 <- coef(fit_dle)[["log(x)"]]
Pcov <- vcov(fit_dle)
.LogisticIndepBeta(
binDLE = binDLE,
DLEdose = DLEdose,
DLEweights = as.integer(DLEweights),
phi1 = phi1,
phi2 = phi2,
Pcov = Pcov,
data = data
)
}
## default constructor ----
#' @rdname LogisticIndepBeta-class
#' @note Typically, end users will not use the `.DefaultLogisticIndepBeta()` function.
#' @export
.DefaultLogisticIndepBeta <- function() {
my_model <- LogisticIndepBeta(
binDLE = c(1.05, 1.8),
DLEweights = c(3L, 3L),
DLEdose = c(25, 300),
data = Data(doseGrid = seq(25, 300, 25))
)
}
# Effloglog ----
## class ----
#' `Effloglog`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`Effloglog`] is the class for the linear log-log efficacy model using pseudo
#' data prior. It describes the relationship between continuous efficacy
#' responses and corresponding dose levels in log-log scale. This efficacy
#' log-log model is given as
#' \deqn{y_i = theta1 + theta2 * log(log(x_i)) + epsilon_i,}
#' where \eqn{y_i} is the efficacy response for subject \eqn{i}, \eqn{x_i} is
#' the dose level treated for subject \eqn{i} and \eqn{epsilon_i} is the random
#' error term of efficacy model at subject \eqn{i}. The error term
#' \eqn{epsilon_i} is a random variable that follows normal distribution with
#' mean \eqn{0} and variance \eqn{nu^{-1}}, which is assumed to be the
#' same for all subjects.
#' There are three parameters in this model, the intercept \eqn{theta1}, the
#' slope \eqn{theta2} and the precision \eqn{nu} of the efficacy responses, also
#' known as the inverse of the variance of the pseudo efficacy responses. It can
#' be a fixed constant or having a gamma distribution. Therefore, a single scalar
#' value or a vector with two positive numbers values must be specified for `nu`
#' slot. If there are some observed efficacy responses available, in the output,
#' `nu` will display the updated value of the precision or the updated values
#' for the parameters of the gamma distribution.
#' The `Effloglog` inherits all slots from [`ModelEff`] class.
#'
#' @details The prior of this model is specified in form of pseudo data. First,
#' at least two dose levels are fixed. Then, using e.g. experts' opinion, the
#' efficacy values that correspond to these dose levels can be obtained,
#' The `eff` and `eff_dose` arguments represent the prior in form of the pseudo
#' data. The `eff` represents the pseudo efficacy values. The `eff_dose`
#' represents the dose levels at which these pseudo efficacy values are
#' observed. Hence, the positions of the elements specified in `eff` and
#' `eff_dose` must correspond to each other between these vectors.
#' Since at least 2 pseudo efficacy values are needed to obtain modal
#' estimates of the intercept and slope parameters, both `eff` and `eff_dose`
#' must be vectors of length at least 2.
#'
#' The joint prior distribution of the intercept \eqn{theta1} and the slope
#' \eqn{theta2} of this model follows bivariate normal distribution with mean
#' \eqn{mu} and covariance matrix \eqn{(nu * Q)^{-1}}.
#' The mean \eqn{mu} is a \eqn{2 x 1} column vector that contains the prior
#' modal estimates of the intercept and the slope.
#' Scalar \eqn{nu} is the precision of the pseudo efficacy responses and
#' \eqn{Q} is the prior or posterior (given that observed, no DLT data is
#' available) precision matrix.
#' It is specified as \eqn{Q = X0^T * X0 + X^T * X}, where \eqn{X0} is a
#' design matrix that is based on pseudo dose levels only, and \eqn{X} is a
#' design matrix that is based on dose levels corresponding to the no DLT
#' efficacy responses observed only (if any).
#' Hence, the \eqn{X0} (or \eqn{X}) will be of size \eqn{r x 2}, if
#' there are \eqn{r >= 2} pseudo efficacy responses specified (or
#' if there are \eqn{r} no DLT efficacy responses observed in the `data`).
#'
#' @slot eff (`numeric`)\cr the pseudo efficacy responses. Each element here
#' must represent responses treated based on one subject.
#' It must be a vector of length at least 2 and the order of its elements must
#' correspond to values specified in `eff_dose`.
#' @slot eff_dose (`numeric`)\cr the pseudo efficacy dose levels at which the
#' pseudo efficacy responses are observed.
#' It must be a vector of length at least 2 and the order of its elements must
#' correspond to values specified in `eff`.
#' @slot nu (`numeric`)\cr parameter of the prior precision of pseudo efficacy
#' responses. This is either a fixed value or a named vector with two positive
#' numbers, the shape (`a`), and the rate (`b`) parameters for the gamma
#' distribution.
#' @slot use_fixed (`flag`)\cr indicates whether `nu` specified is a fixed value
#' or a vector with two parameters for gamma distribution. This slot is for
#' internal purposes only and must not be used by the user.
#' @slot theta1 (`number`)\cr the intercept in this efficacy log-log model. This
#' slot is used in output to display the resulting prior or posterior modal
#' estimates obtained based on the pseudo and observed (if any) data.
#' @slot theta2 (`number`)\cr the slope in this efficacy log-log model. This
#' slot is used in output to display the resulting prior or posterior modal
#' estimates obtained based on the pseudo and observed (if any) data.
#' @slot Pcov (`matrix`)\cr refers to the \eqn{2 x 2} covariance matrix of the
#' estimators of the intercept \eqn{theta1} and the slope \eqn{theta2}
#' parameters in this model.
#' This is used in output to display the resulting prior and posterior
#' covariance matrix of \eqn{theta1} and \eqn{theta2} obtained, based on the
#' pseudo and observed (if any) data. This slot is needed for internal purposes.
#' @slot X (`matrix`)\cr is the design matrix that is based on either the pseudo
#' dose levels or observed dose levels (without DLT). This is used
#' in the output to display the design matrix for the pseudo or the observed
#' efficacy responses.
#' @slot Y (`numeric`)\cr is a vector that either contains the pseudo efficacy
#' responses or observed efficacy responses (without DLT).
#' @slot mu (`numeric`)\cr a vector of the prior or the posterior modal estimates
#' of the intercept (\eqn{theta1}) and the slope (\eqn{theta2}).
#' This slot is used in output to display as the mean of the prior or posterior
#' bivariate normal distribution for \eqn{theta1} and \eqn{theta2}.
#' @slot Q (`matrix`)\cr is the prior or posterior (given that observed, no DLT
#' data is available) precision matrix. It is specified as
#' \eqn{Q = X0^T * X0 + X^T * X}, where \eqn{X0} is a design matrix that is
#' based on pseudo dose levels only, and \eqn{X} is a design matrix that is
#' based on dose levels corresponding to the observed, no DLT efficacy values
#' only (if any).
#' @slot const (`number`)\cr a non-negative number (default to 0), leading to the
#' model form described above. In general, the model has the form
#' \eqn{y_i = theta1 + theta2 * log(log(x_i + const)) + epsilon_i}, such that
#' dose levels greater than \eqn{1 - const} can be considered as described in
#' Yeung et al. (2015).
#'
#' @aliases Effloglog
#' @export
#'
.Effloglog <- setClass(
Class = "Effloglog",
slots = c(
eff = "numeric",
eff_dose = "numeric",
nu = "numeric",
use_fixed = "logical",
theta1 = "numeric",
theta2 = "numeric",
Pcov = "matrix",
X = "matrix",
Y = "numeric",
mu = "numeric",
Q = "matrix",
const = "numeric"
),
prototype = prototype(
eff = c(0, 0),
eff_dose = c(1, 1),
nu = 1 / 0.025,
use_fixed = TRUE,
const = 0
),
contains = "ModelEff",
validity = v_model_eff_log_log
)
## constructor ----
#' @rdname Effloglog-class
#'
#' @param eff (`numeric`)\cr the pseudo efficacy responses.
#' Elements of `eff` must correspond to the elements of `eff_dose`.
#' @param eff_dose (`numeric`)\cr dose levels that correspond to pseudo efficacy
#' responses in `eff`.
#' @param nu (`numeric`)\cr the precision (inverse of the variance) of the
#' efficacy responses. This is either a fixed value or a named vector with two
#' positive numbers, the shape (`a`), and the rate (`b`) parameters for the
#' gamma distribution.
#' @param data (`DataDual`)\cr observed data to update estimates of the model
#' parameters.
#' @param const (`number`)\cr the constant value added to the dose level when
#' the dose level value is less than or equal to 1 and a special form of the
#' linear log-log has to applied (Yeung et al. (2015).).
#'
#' @export
#' @example examples/Model-class-Effloglog.R
#'
Effloglog <- function(eff,
eff_dose,
nu,
data,
const = 0) {
assert_numeric(eff)
assert_numeric(eff_dose, len = length(eff))
assert_numeric(nu, min.len = 1, max.len = 2)
assert_class(data, "Data")
assert_number(const, finite = TRUE)
use_fixed <- length(nu) == 1L
eff_dose <- eff_dose + const
# Get observed efficacy data without DLT (if any).
eff_obsrv_w_x <- getEff(data, no_dlt = TRUE)
eff_obsrv <- eff_obsrv_w_x$w_no_dlt
eff_obsrv_dose <- eff_obsrv_w_x$x_no_dlt + const
# Fit pseudo and observed (if any) efficacy.
w <- c(eff, eff_obsrv)
x <- c(eff_dose, eff_obsrv_dose)
fit_eff <- suppressWarnings(lm(w ~ log(log(x))))
X <- model.matrix(fit_eff)
Y <- w
mu <- coef(fit_eff) # This is [theta1, theta2]^T est.
Q <- crossprod(X)
Pcov <- vcov(fit_eff)
nobs_no_dlt <- length(eff_obsrv)
if (nobs_no_dlt > 0L) { # Observed data available.
# Set X, Y to observed data only.
X <- model.matrix(fit_eff)[-seq_along(eff), ]
Y <- eff_obsrv
fit_eff0 <- lm(eff ~ log(log(eff_dose))) # Pseudo only.
X0 <- model.matrix(fit_eff0)
mu0 <- coef(fit_eff0)
Q0 <- crossprod(X0)
# Note that mu = (Q0 + X^T * X)^{-1} * (Q0 * mu0 + X^T * X * (X^T * X)^{-1} X^T * Y),
# given that (X^T * X) is invertible and X, Y, mu0, Q0, are specified in this else block.
if (!use_fixed) {
nu["a"] <- nu["a"] + (nobs_no_dlt) / 2
nu["b"] <- nu["b"] + (crossprod(Y) + t(mu0) %*% Q0 %*% mu0 - t(mu) %*% Q %*% mu) / 2
}
}
.Effloglog(
eff = eff,
eff_dose = eff_dose,
nu = nu,
use_fixed = use_fixed,
theta1 = mu[["(Intercept)"]],
theta2 = mu[["log(log(x))"]],
Pcov = Pcov,
X = X,
Y = Y,
mu = as.vector(mu),
Q = Q,
const = const,
data = data
)
}
## default constructor ----
#' @rdname Effloglog-class
#' @note Typically, end users will not use the `.DefaultEffloglog()` function.
#' @export
.DefaultEffloglog <- function() {
emptydata <- DataDual(doseGrid = seq(25, 300, 25), placebo = FALSE)
my_data <- DataDual(
x = c(25, 50, 50, 75, 100, 100, 225, 300),
y = c(0, 0, 0, 0, 1, 1, 1, 1),
w = c(0.31, 0.42, 0.59, 0.45, 0.6, 0.7, 0.6, 0.52),
doseGrid = emptydata@doseGrid,
ID = 1L:8L,
cohort = as.integer(c(1, 2, 2, 3, 4, 4, 5, 6))
)
Effloglog(
eff = c(1.223, 2.513),
eff_dose = c(25, 300),
nu = c(a = 1, b = 0.025),
data = my_data
)
}
# EffFlexi ----
## class ----
#' `EffFlexi`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`EffFlexi`] is the class for the efficacy model in flexible form of prior
#' expressed in form of pseudo data. In this class, a flexible form is used to
#' describe the relationship between the efficacy responses and the dose levels
#' and it is specified as
#' \deqn{(W | betaW, sigma2W) ~ Normal(X * betaW, sigma2W * I),}
#' where \eqn{W} is a vector of the efficacy responses, \eqn{betaW} is a column
#' vector of the mean efficacy responses for all dose levels, and \eqn{X} is
#' the design matrix with entries \eqn{I_i,j} that are equal to 1 if subject
#' \eqn{i} is allocated to dose \eqn{j}, and \eqn{0} otherwise. The \eqn{sigma2W}
#' is the variance of the efficacy responses which can be either a fixed number
#' or a number from an inverse gamma distribution.
#' This flexible form aims to capture different shapes of the dose-efficacy
#' curve. In addition, the first (RW1) or second order (RW2) random walk model
#' can be used for smoothing data. That is the random walk model is used to model
#' the first or the second order differences of the mean efficacy responses to
#' its neighboring dose levels of their mean efficacy responses.
#'
#' The RW1 model is given as
#' \deqn{betaW_j - betaW_j-1) ~ Normal(0, sigma2betaW),}
#' and for RW2 as
#' \deqn{betaW_j-2 - 2 * betaW_j-1 + beta_j ~ Normal(0, sigma2betaW),}
#' where \eqn{betaW_j} is the vector of mean efficacy responses at dose j, and
#' the \eqn{sigma2betaW} is the prior variance which can be either a fixed
#' number or a number from an inverse gamma distribution.
#'
#' The `eff` and `eff_dose` are the pseudo efficacy responses and dose levels at
#' which these pseudo efficacy responses are observed. Both, `eff` and `eff_dose`
#' must be vectors of length at least 2. The positions of the elements specified
#' in `eff` and `eff_dose` must correspond to each other between these vectors.
#'
#' @details This model will output the updated value or the updated values of the
#' parameters of the inverse gamma distributions for \eqn{sigma2W} and
#' \eqn{sigma2betaW}. The `EffFlexi` inherits all slots from [`ModelEff`] class.
#'
#' @slot eff (`numeric`)\cr the pseudo efficacy responses. Each element here
#' must represent responses treated based on one subject.
#' It must be a vector of length at least 2 and the order of its elements must
#' correspond to values specified in `eff_dose`.
#' @slot eff_dose (`numeric`)\cr the pseudo efficacy dose levels at which the
#' pseudo efficacy responses are observed.
#' It must be a vector of length at least 2 and the order of its elements must
#' correspond to values specified in `eff`.
#' @slot sigma2W (`numeric`)\cr the prior variance of the flexible efficacy form.
#' This is either a fixed value or a named vector with two positive numbers,
#' the shape (`a`), and the rate (`b`) parameters for the gamma distribution.
#' @slot sigma2betaW (`numeric`)\cr the prior variance of the random walk model
#' for the mean efficacy responses. This is either a fixed value or a named
#' vector with two positive numbers, the shape (`a`), and the rate (`b`)
#' parameters for the gamma distribution.
#' @slot use_fixed (`logical`)\cr indicates whether a fixed value for
#' `sigma2W` and `sigma2betaW` (for each parameter separately) is used or not.
#' This slot is needed for internal purposes and must not be touched by the user.
#' @slot rw1 (`flag`)\cr used for smoothing data for this efficacy model. If it
#' is `TRUE`, the first-order random walk model is used for the mean efficacy
#' responses. Otherwise, the random walk of second order is used.
#' @slot X (`matrix`)\cr the design matrix for the efficacy responses. It is
#' based on both the pseudo and the observed efficacy responses.
#' @slot RW (`matrix`)\cr the difference matrix for the random walk model. This
#' slot is needed for internal purposes and must not be used by the user.
#' @slot RW_rank (`integer`)\cr is the rank of the difference matrix. This
#' slot is needed for internal purposes and must not be used by the user.
#'
#' @aliases EffFlexi
#' @export
#'
.EffFlexi <- setClass(
Class = "EffFlexi",
slots = c(
eff = "numeric",
eff_dose = "numeric",
sigma2W = "numeric",
sigma2betaW = "numeric",
use_fixed = "logical",
rw1 = "logical",
X = "matrix",
RW = "matrix",
RW_rank = "integer"
),
prototype = prototype(
eff = c(0, 0),
eff_dose = c(1, 1),
sigma2W = 0.025,
sigma2betaW = 1,
rw1 = TRUE,
use_fixed = c(sigma2W = TRUE, sigma2betaW = TRUE)
),
contains = "ModelEff",
validity = v_model_eff_flexi
)
## constructor ----
#' @rdname EffFlexi-class
#'
#' @param eff (`numeric`)\cr the pseudo efficacy responses.
#' Elements of `eff` must correspond to the elements of `eff_dose`.
#' @param eff_dose (`numeric`)\cr dose levels that correspond to pseudo efficacy
#' responses in `eff`.
#' @param sigma2W (`numeric`)\cr the prior variance of the efficacy responses.
#' This is either a fixed value or a named vector with two positive numbers,
#' the shape (`a`), and the rate (`b`) parameters for the inverse gamma
#' distribution.
#' @param sigma2betaW (`numeric`)\cr the prior variance of the random walk model
#' used for smoothing. This is either a fixed value or a named vector with two
#' positive numbers, the shape (`a`), and the rate (`b`) parameters for the
#' inverse gamma distribution.
#' @param rw1 (`flag`)\cr used for smoothing data for this efficacy model. If it
#' is `TRUE`, the first-order random walk model is used for the mean efficacy
#' responses. Otherwise, the random walk of second order is used.
#' @param data (`DataDual`)\cr observed data to update estimates of the model
#' parameters.
#'
#' @export
#' @example examples/Model-class-EffFlexi.R
#'
EffFlexi <- function(eff,
eff_dose,
sigma2W,
sigma2betaW,
rw1 = TRUE,
data) {
assert_numeric(eff)
assert_numeric(eff_dose)
assert_numeric(sigma2W, min.len = 1, max.len = 2)
assert_numeric(sigma2betaW, min.len = 1, max.len = 2)
assert_flag(rw1)
assert_class(data, "DataDual")
use_fixed <- c(
sigma2W = test_number(sigma2W),
sigma2betaW = test_number(sigma2betaW)
)
x <- c(eff_dose, getEff(data, no_dlt = TRUE)$x_no_dlt)
x_level <- match_within_tolerance(x, data@doseGrid)
X <- model.matrix(~ -1L + factor(x_level, levels = seq_len(data@nGrid)))
X <- matrix(as.integer(X), ncol = ncol(X)) # To remove some obsolete attributes.
# Set up the random walk penalty matrix and its rank.
# D1: difference matrix of order 1.
D1 <- cbind(0, diag(data@nGrid - 1)) - cbind(diag(data@nGrid - 1), 0)
if (rw1) { # the rank-deficient prior precision for the RW1 prior.
RW <- crossprod(D1)
RW_rank <- data@nGrid - 1L # rank = dimension - 1. # nolintr
} else { # Second-order difference.
D2 <- D1[-1, -1] %*% D1
RW <- crossprod(D2)
RW_rank <- data@nGrid - 2L # nolintr
}
.EffFlexi(
eff = eff,
eff_dose = eff_dose,
sigma2W = sigma2W,
sigma2betaW = sigma2betaW,
use_fixed = use_fixed,
rw1 = rw1,
X = X,
RW = RW,
RW_rank = RW_rank,
data = data
)
}
## default constructor ----
#' @rdname EffFlexi-class
#' @note Typically, end users will not use the `.DefaultEffFlexi()` function.
#' @export
.DefaultEffFlexi <- function() {
empty_data <- DataDual(doseGrid = seq(25, 300, 25))
EffFlexi(
eff = c(1.223, 2.513),
eff_dose = c(25, 300),
sigma2W = c(a = 0.1, b = 0.1),
sigma2betaW = c(a = 20, b = 50),
rw1 = FALSE,
data = empty_data
)
data <- DataDual(
x = c(25, 50, 50, 75, 100, 100, 225, 300),
y = c(0, 0, 0, 0, 1, 1, 1, 1),
w = c(0.31, 0.42, 0.59, 0.45, 0.6, 0.7, 0.6, 0.52),
doseGrid = empty_data@doseGrid,
ID = 1L:8L,
cohort = as.integer(c(1, 2, 2, 3, 4, 4, 5, 6))
)
EffFlexi(
eff = c(1.223, 2.513),
eff_dose = c(25, 300),
sigma2W = c(a = 0.1, b = 0.1),
sigma2betaW = c(a = 20, b = 50),
rw1 = FALSE,
data = data
)
}
# DALogisticLogNormal ----
## class ----
#' `DALogisticLogNormal`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`DALogisticLogNormal`] is the class for the logistic model with bivariate
#' (log) normal prior and data augmentation. This class inherits from the
#' [`LogisticLogNormal`] class.
#'
#' @note We still need to include here formula for the lambda prior.
#'
#' @slot npiece (`number`)\cr the number of pieces in the `PEM`.
#' @slot l (`numeric`)\cr a vector used in the lambda prior.
#' @slot c_par (`numeric`)\cr a parameter used in the lambda prior; according to
#' Liu's paper, `c_par = 2` is recommended.
#' @slot cond_pem (`flag`)\cr is a conditional piecewise-exponential model used?
#' (default). Otherwise an unconditional model is used.
#'
#' @seealso [`ModelLogNormal`], [`LogisticNormal`], [`LogisticLogNormal`].
#'
#' @aliases DALogisticLogNormal
#' @export
#'
.DALogisticLogNormal <- setClass(
Class = "DALogisticLogNormal",
slots = c(
npiece = "integer",
l = "numeric",
c_par = "numeric",
cond_pem = "logical"
),
prototype = prototype(
npiece = 3L,
l = 0.5,
c_par = 2,
cond_pem = TRUE
),
contains = "LogisticLogNormal",
validity = v_model_da_logistic_log_normal
)
## constructor ----
#' @rdname DALogisticLogNormal-class
#'
#' @param npiece (`number`)\cr the number of pieces in the `PEM`.
#' @param l (`numeric`)\cr a vector used in the lambda prior.
#' @param c_par (`numeric`)\cr a parameter used in the lambda prior; according to
#' Liu's paper, `c_par = 2` is recommended.
#' @param cond_pem (`flag`)\cr is a conditional piecewise-exponential model used?
#' (default). Otherwise an unconditional model is used.
#' @inheritDotParams LogisticLogNormal
#'
#' @export
#' @example examples/Model-class-DALogisticLogNormal.R
#'
DALogisticLogNormal <- function(npiece = 3,
l,
c_par = 2,
cond_pem = TRUE,
...) {
assert_flag(cond_pem)
start <- LogisticLogNormal(...)
datamodel <- function() {
for (i in 1:nObs) {
# Part I: describe the logistic model of DLTs vs dose.
logit(p[i]) <- alpha0 + alpha1 * log(x[i] / ref_dose)
# Part II: describe the piecewise exponential.
# Notice that:
# when y=1 -> DLT=1 and u=<T;
# when y=0 & T<t (u=T) -> DLT=0;
# when y=0 & T>t (u<T) -> DLT=NA/missing;
# when indx=0 -> censored, i.e u<T and event=0;
# when indx=1 -> not censored, i.e. u>=T or event=1;
indx[i] <- 1 - step(Tmax - u[i] - eps) * (1 - y[i])
for (j in 1:npiece) {
# When not censored, i.e DLT!=NA & t[i]=u[i];
# if t[i]<h[j], d[i,j]=0;
# if h[j]<t[i]=<h[j+1], d[i,j]=1
# if h[j+1]<t[i], d[i,j]=0
# When censored t[i]>u[i] -> d[i,j]=0
d[i, j] <- y[i] * step(u[i] - h[j] - eps) * step(h[j + 1] - u[i])
# DLT free survival(time) for patient i in interval I(j);
# if t[i]<h[j], s[i,j]=0;
# if h[j]<t[i]<=h[j+1], s[i,j]=t[i]-h[j]
# if h[j+1]<=t[i], s[i,j]=h[j+1]-h[j]
s[i, j] <- min(u[i] - h[j], h[j + 1] - h[j]) * step(u[i] - h[j])
# piecewise exponential hazard rate lambda[j];
mu_u[i, j] <- lambda[j] * s[i, j]
mu[i, j] <- d[i, j] * log(lambda[j]) - y[i] * mu_u[i, j]
}
# The likelihood function.
L_obs[i] <- exp(sum(mu[i, ])) * pow(p[i] / A, y[i]) * pow(1 - p[i], 1 - y[i]) # Not censored. # nolintr
L_cnsr[i] <- 1 - p[i] * (1 - exp(-sum(mu_u[i, ]))) / A # Censored. # nolintr
L[i] <- pow(L_obs[i], indx[i]) * pow(L_cnsr[i], 1 - indx[i])
# Apply zero trick in JAGS.
phi[i] <- -log(L[i]) + cadj
zeros[i] ~ dpois(phi[i])
}
}
priormodel <- h_jags_join_models(
start@priormodel,
function() {
g_beta <- 1 / c_par
for (j in 1:npiece) {
g_alpha[j] <- l[j] / c_par
lambda[j] ~ dgamma(g_alpha[j], g_beta)
mu_T[j] <- lambda[j] * (h[j + 1] - h[j]) # nolintr
}
# If cond = 1, then conditional PEM is used and A is defined as
# the probability to have DLT, i.e. t<T, otherwise
# cond = 0 and A is just 1 (so no impact in likelihood).
A <- cond * (1 - exp(-sum(mu_T))) + (1 - cond)
}
)
modelspecs <- function(nObs, Tmax, from_prior) {
ms <- list(
prec = start@params@prec,
mean = start@params@mean,
npiece = npiece,
l = l,
c_par = c_par,
h = seq(from = 0L, to = Tmax, length = npiece + 1),
cond = as.integer(cond_pem)
)
if (!from_prior) {
ms <- c(list(ref_dose = start@ref_dose, zeros = rep(0, nObs), eps = 1e-10, cadj = 1e10), ms)
}
ms
}
assert_integerish(npiece, lower = 1)
.DALogisticLogNormal(
start,
npiece = as.integer(npiece),
l = l,
c_par = c_par,
cond_pem = cond_pem,
datamodel = datamodel,
priormodel = priormodel,
modelspecs = modelspecs,
datanames = c("nObs", "y", "x", "u", "Tmax"),
sample = c("alpha0", "alpha1", "lambda")
)
}
## default constructor ----
#' @rdname DALogisticLogNormal-class
#' @note Typically, end users will not use the `.DefaultDALogisticLogNormal()` function.
#' @export
.DefaultDALogisticLogNormal <- function() {
npiece <- 10
Tmax <- 60
lambda_prior <- function(k) {
npiece / (Tmax * (npiece - k + 0.5))
}
DALogisticLogNormal(
mean = c(-0.85, 1),
cov = matrix(c(1, -0.5, -0.5, 1), nrow = 2),
ref_dose = 56,
npiece = npiece,
l = as.numeric(t(apply(as.matrix(c(1:npiece), 1, npiece), 2, lambda_prior))),
c_par = 2
)
}
# TITELogisticLogNormal ----
## class ----
#' `TITELogisticLogNormal`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`TITELogisticLogNormal`] is the class for TITE-CRM based on a logistic
#' regression model using a bivariate normal prior on the intercept and log
#' slope parameters.
#'
#' This class inherits from the [`LogisticLogNormal`].
#'
#' @slot weight_method (`string`)\cr the weight function method: either linear
#' or adaptive. This was used in Liu, Yin and Yuan's paper.
#'
#' @seealso [`DALogisticLogNormal`].
#'
#' @aliases TITELogisticLogNormal
#' @export
#'
.TITELogisticLogNormal <- setClass(
Class = "TITELogisticLogNormal",
slots = c(weight_method = "character"),
prototype = prototype(weight_method = "linear"),
contains = "LogisticLogNormal",
validity = v_model_tite_logistic_log_normal
)
## constructor ----
#' @rdname TITELogisticLogNormal-class
#'
#' @param weight_method (`string`)\cr the weight function method: either linear
#' or adaptive. This was used in Liu, Yin and Yuan's paper.
#' @inheritDotParams LogisticLogNormal
#'
#' @export
#' @example examples/Model-class-TITELogisticLogNormal.R
#'
TITELogisticLogNormal <- function(weight_method = "linear",
...) {
assert_character(weight_method, min.len = 1L, max.len = 2L, any.missing = FALSE)
start <- LogisticLogNormal(...)
datamodel <- function() {
for (i in 1:nObs) {
logit(p[i]) <- alpha0 + alpha1 * log(x[i] / ref_dose)
# The piecewise exponential likelihood. Notice that:
# when y=1 -> DLT=1 and u=<T;
# when y=0 & T<t (u=T) -> DLT=0;
# when y=0 & T>t (u<T) -> DLT=NA/missing;
# when indx=0 -> censored, i.e u<T and event=0;
# when indx=1 -> not censored, i.e. u>=T or event=1;
L[i] <- pow(p[i], y[i]) * pow((1 - w[i] * p[i]), (1 - y[i]))
# Apply zero trick in JAGS.
phi[i] <- -log(L[i]) + cadj
zeros[i] ~ dpois(phi[i])
}
}
modelspecs <- function(nObs, u, Tmax, y, from_prior) {
ms <- list(prec = start@params@prec, mean = start@params@mean)
# Calculate weights `w` based on the input data.
if (!from_prior && nObs > 0L) {
if (weight_method == "linear") {
w <- u / Tmax
} else if (weight_method == "adaptive") {
nDLT <- sum(y)
if (nDLT > 0) {
u_dlt <- sort(u[y == 1])
w <- sapply(u, function(u_i) {
m <- sum(u_i >= u_dlt)
w_i <- if (m == 0) {
u_i / u_dlt[1]
} else if (m < nDLT) {
m + (u_i - u_dlt[m]) / (u_dlt[m + 1] - u_dlt[m])
} else { # m == nDLT. nolintr
m + (u_i - u_dlt[m]) / (Tmax + 0.00000001 - u_dlt[m])
}
w_i / (nDLT + 1)
})
} else {
w <- u / Tmax
}
}
w[y == 1] <- 1
w[u == Tmax] <- 1
ms <- c(list(ref_dose = start@ref_dose, zeros = rep(0, nObs), cadj = 1e10, w = w), ms)
}
ms
}
.TITELogisticLogNormal(
start,
weight_method = weight_method,
datamodel = datamodel,
modelspecs = modelspecs,
datanames = c("nObs", "y", "x")
)
}
## default constructor ----
#' @rdname TITELogisticLogNormal-class
#' @note Typically, end users will not use the `.DefaultTITELogisticLogNormal()` function.
#' @export
.DefaultTITELogisticLogNormal <- function() {
TITELogisticLogNormal(
mean = c(0, 1),
cov = diag(2),
ref_dose = 1,
weight_method = "linear"
)
}
# OneParLogNormalPrior ----
## class ----
#' `OneParLogNormalPrior`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`OneParLogNormalPrior`] is the class for a standard CRM with a normal prior on
#' the log power parameter for the skeleton prior probabilities.
#'
#' @slot skel_fun (`function`)\cr function to calculate the prior DLT probabilities.
#' @slot skel_fun_inv (`function`)\cr inverse function of `skel_fun`.
#' @slot skel_probs (`numeric`)\cr skeleton prior probabilities. This is a vector
#' of unique and sorted probability values between 0 and 1.
#' @slot sigma2 (`number`)\cr prior variance of log power parameter alpha.
#'
#' @seealso [`ModelLogNormal`].
#'
#' @aliases OneParLogNormalPrior
#' @export
#'
.OneParLogNormalPrior <- setClass(
Class = "OneParLogNormalPrior",
slots = c(
skel_fun = "function",
skel_fun_inv = "function",
skel_probs = "numeric",
sigma2 = "numeric"
),
contains = "GeneralModel",
validity = v_model_one_par_exp_normal_prior
)
## constructor ----
#' @rdname OneParLogNormalPrior-class
#'
#' @param skel_probs (`numeric`)\cr skeleton prior probabilities. This is a vector
#' of unique and sorted probability values between 0 and 1.
#' @param dose_grid (`numeric`)\cr dose grid. It must be must be a sorted vector
#' of the same length as `skel_probs`.
#' @param sigma2 (`number`)\cr prior variance of log power parameter alpha.
#'
#' @export
#' @example examples/Model-class-OneParLogNormalPrior.R
#'
OneParLogNormalPrior <- function(skel_probs,
dose_grid,
sigma2) {
assert_probabilities(skel_probs, unique = TRUE, sorted = TRUE) # So that skel_fun_inv exists.
assert_numeric(dose_grid, len = length(skel_probs), any.missing = FALSE, unique = TRUE, sorted = TRUE)
skel_fun <- approxfun(x = dose_grid, y = skel_probs, rule = 2)
skel_fun_inv <- approxfun(x = skel_probs, y = dose_grid, rule = 2)
.OneParLogNormalPrior(
skel_fun = skel_fun,
skel_fun_inv = skel_fun_inv,
skel_probs = skel_probs,
sigma2 = sigma2,
datamodel = function() {
for (i in 1:nObs) {
p[i] <- skel_probs[xLevel[i]]^exp(alpha)
y[i] ~ dbern(p[i])
}
},
priormodel = function() {
alpha ~ dnorm(0, 1 / sigma2)
},
modelspecs = function(from_prior) {
ms <- list(sigma2 = sigma2)
if (!from_prior) {
ms$skel_probs <- skel_probs
}
ms
},
init = function() {
list(alpha = 1)
},
datanames = c("nObs", "y", "xLevel"),
sample = "alpha"
)
}
## default constructor ----
#' @rdname OneParLogNormalPrior-class
#' @return an instance of the `OneParLogNormalPrior` class
#' @export
.DefaultOneParLogNormalPrior <- function() {
OneParLogNormalPrior(
skel_probs = seq(from = 0.1, to = 0.9, length = 5),
dose_grid = 1:5,
sigma2 = 2
)
}
# OneParExpPrior ----
## class ----
#' `OneParExpPrior`
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' [`OneParExpPrior`] is the class for a standard CRM with an exponential prior
#' on the power parameter for the skeleton prior probabilities. It is an
#' implementation of a version of the one-parameter CRM (O’Quigley et al. 1990).
#'
#' @note Typically, end users will not use the `.DefaultOneparExpPrior()` function.
#'
#' @slot skel_fun (`function`)\cr function to calculate the prior DLT probabilities.
#' @slot skel_fun_inv (`function`)\cr inverse function of `skel_fun`.
#' @slot skel_probs (`numeric`)\cr skeleton prior probabilities. This is a vector
#' of unique and sorted probability values between 0 and 1.
#' @slot lambda (`number`)\cr rate parameter of prior exponential distribution
#' for theta.
#'
#' @aliases OneParExpPrior
#' @export
#'
.OneParExpPrior <- setClass(
Class = "OneParExpPrior",
slots = c(
skel_fun = "function",
skel_fun_inv = "function",
skel_probs = "numeric",
lambda = "numeric"
),
contains = "GeneralModel",
validity = v_model_one_par_exp_prior
)
## constructor ----
#' @rdname OneParExpPrior-class
#'
#' @param skel_probs see slot definition.
#' @param dose_grid (`numeric`)\cr dose grid. It must be must be a sorted vector
#' of the same length as `skel_probs`.
#' @param lambda see slot definition.
#'
#' @export
#' @example examples/Model-class-OneParExpPrior.R
#'
OneParExpPrior <- function(skel_probs,
dose_grid,
lambda) {
assert_probabilities(skel_probs, unique = TRUE, sorted = TRUE) # So that skel_fun_inv exists.
assert_numeric(dose_grid, len = length(skel_probs), any.missing = FALSE, unique = TRUE, sorted = TRUE)
skel_fun <- approxfun(x = dose_grid, y = skel_probs, rule = 2)
skel_fun_inv <- approxfun(x = skel_probs, y = dose_grid, rule = 2)
.OneParExpPrior(
skel_fun = skel_fun,
skel_fun_inv = skel_fun_inv,
skel_probs = skel_probs,
lambda = lambda,
datamodel = function() {
for (i in 1:nObs) {
p[i] <- skel_probs[xLevel[i]]^theta
y[i] ~ dbern(p[i])
}
},
priormodel = function() {
theta ~ dexp(lambda)
},
modelspecs = function(from_prior) {
ms <- list(lambda = lambda)
if (!from_prior) {
ms$skel_probs <- skel_probs
}
ms
},
init = function() {
list(theta = 1)
},
datanames = c("nObs", "y", "xLevel"),
sample = "theta"
)
}
## default constructor ----
#' @rdname OneParExpPrior-class
#' @note Typically, end users will not use the `.DefaultOneParLogNormalPrior()` function.
#' @export
.DefaultOneParExpPrior <- function() {
OneParExpPrior(
skel_probs = c(0.1, 0.3, 0.5, 0.7, 0.9),
dose_grid = 1:5,
lambda = 2
)
}
# FractionalCRM ----
## class ----
#' `FractionalCRM`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' [`FractionalCRM`] is the class for a fractional CRM model based on a one
#' parameter CRM (with normal prior on the log-power parameter) as well as
#' Kaplan-Meier based estimation of the conditional probability to experience a
#' DLT for non-complete observations.
#'
#' This fractional CRM model follows the paper and code by Guosheng Yin et al.
#'
#' @seealso [`TITELogisticLogNormal`].
#'
#' @aliases FractionalCRM
#' @export
#'
.FractionalCRM <- setClass(
Class = "FractionalCRM",
contains = "OneParLogNormalPrior"
)
## constructor ----
#' @rdname FractionalCRM-class
#'
#' @inheritDotParams OneParLogNormalPrior
#'
#' @export
#' @example examples/Model-class-FractionalCRM.R
#'
FractionalCRM <- function(...) {
start <- OneParLogNormalPrior(...)
# This is adapted from the TITELogisticLogNormal class.
datamodel <- function() {
for (i in 1:nObs) {
p[i] <- skel_probs[xLevel[i]]^exp(alpha)
# The piecewise exponential likelihood. Notice that:
# when y=1 -> DLT=1 and u=<T;
# when y=0 & T<t (u=T) -> DLT=0;
# when y=0 & T>t (u<T) -> DLT=NA/missing.
# Therefore, `yhat` is used instead of `y` for the likelihood f. (see `modelspecs`).
L[i] <- pow(p[i], yhat[i]) * pow((1 - p[i]), (1 - yhat[i]))
# Apply zero trick in JAGS.
phi[i] <- -log(L[i]) + cadj
zeros[i] ~ dpois(phi[i])
}
}
modelspecs <- function(nObs, u, Tmax, y, from_prior) {
ms <- list(sigma2 = start@sigma2)
if (!from_prior) {
# Calculate fractional contribution `yhat`
# based on the input data using the Kaplan-Meier method.
yhat <- if (nObs > 0) {
km <- survival::survfit(survival::Surv(u, y) ~ 1)
s_tau <- tail(km$surv[km$time <= Tmax], 1) # Survival probability = S(Tmax).
ifelse(
u < Tmax & y == 0L, # Within the assessment window and so far no DLT.
yes = 1 - s_tau / sapply(u, function(u_i) tail(km$surv[km$time <= u_i], 1)),
no = y
)
} else {
1L
}
ms <- c(
list(skel_probs = start@skel_probs, zeros = rep(0, nObs), cadj = 1e10, yhat = yhat),
ms
)
}
ms
}
.FractionalCRM(
start,
datamodel = datamodel,
modelspecs = modelspecs,
datanames = c("nObs", "xLevel")
)
}
## default constructor ----
#' @rdname FractionalCRM-class
#' @note Typically, end users will not use the `.DefaultTITELogisticLogNormal()` function.
#' @export
.DefaultFractionalCRM <- function() {
FractionalCRM(
skel_probs = c(0.1, 0.2, 0.3, 0.4),
dose_grid = c(10, 30, 50, 100),
sigma2 = 2
)
}
## class ----
#' `LogisticLogNormalOrdinal`
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' [`LogisticLogNormalOrdinal`] is the class for a logistic lognormal CRM model
#' using an ordinal toxicity scale.
#'
#' @aliases LogisticLogNormalOrdinal
#' @export
.LogisticLogNormalOrdinal <- setClass(
Class = "LogisticLogNormalOrdinal",
contains = "ModelLogNormal",
validity = v_logisticlognormalordinal
)
## constructor ----
#' @rdname LogisticLogNormalOrdinal-class
#' @inheritParams ModelLogNormal
#' @export
#' @example examples/Model-class-LogisticLogNormalOrdinal.R
LogisticLogNormalOrdinal <- function(mean, cov, ref_dose) {
params <- ModelParamsNormal(mean, cov)
.LogisticLogNormalOrdinal(
params = params,
ref_dose = positive_number(ref_dose),
priormodel = function() {
alpha[1] ~ dnorm(mean[1], prec[1, 1])
for (i in 2:(k - 1)) {
alpha[i] ~ dnorm(mean[i], prec[i, i]) %_% T(, alpha[i - 1])
}
gamma ~ dnorm(mean[k], prec[k, k])
beta <- exp(gamma)
},
datamodel = function() {
for (i in 1:nObs) {
xhat[i] <- log(x[i] / ref_dose)
for (j in 1:(k - 1)) {
z[i, j] <- alpha[j] + beta * xhat[i]
p[i, j] <- exp(z[i, j]) / (1 + exp(z[i, j]))
tox[i, j] ~ dbern(p[i, j])
}
}
},
modelspecs = function(y, from_prior) {
ms <- list(
mean = params@mean,
prec = params@prec,
k = length(mean),
tox = array(dim = c(length(y), length(mean) - 1))
)
if (!from_prior) {
for (i in seq_along(y)) {
for (j in 1:(ms$k - 1)) {
ms$tox[i, j] <- y[i] >= j
}
}
ms$ref_dose <- ref_dose
}
ms
},
init = function() {
list(alpha = sapply(1:(length(mean) - 1), function(x) -(x + 1)), gamma = 1)
},
datanames = c("nObs", "y", "x"),
# Need to provide JAGS column names here
sample = c(paste0("alpha[", 1:(length(mean) - 1), "]"), "beta")
)
}
## default constructor ----
#' @rdname LogisticLogNormalOrdinal-class
#' @note Typically, end users will not use the `.DefaultLogisticLogNormalOrdinal()` function.
#' @export
.DefaultLogisticLogNormalOrdinal <- function() {
LogisticLogNormalOrdinal(
mean = c(-3, -4, 1),
cov = diag(c(3, 4, 1)),
ref_dose = 50
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.