Data Generator

Description

It generates data for a dose-finding trial.

Usage

data_generator(doses, sample.size, distr, parm, 
	censoring.rate = NULL)

Arguments

Details

If distr = "weibull", the list parm should contain two components - lambda and sigma - that are the scale and shape parameters in the following parametrization of the Weibull distribution:

f(t;\lambda,\sigma) = \frac{1}{\sigma\lambda^{1/\sigma}} t^{1/\sigma-1} \exp\left\{ -\left( t/\lambda \right)^{1/\sigma} \right\}, \ t > 0,

with hazard rate given by

h(t) = \frac{1}{\lambda^{1/\sigma}\sigma}t^{1/\sigma - 1}

and regression structure

\log(\lambda_i) = d_i\beta_i.

where \log(\lambda_i) represents the model effect for dose i, doses[i].

Value

a data frame of dimension [length(doses)\timessample.size] \times 3 when distr = "weibull" containing time-to-event, censoring indicator and dose.

References

Diniz, Márcio A. and Gallardo, Diego I. and Magalhães, Tiago M. (2023). Improved inference for MCP-Mod approach for time-to-event endpoints with small sample sizes. arXiv <doi.org/10.48550/arXiv.2301.00325>

Examples

library(DoseFinding)
library(MCPModBC)

## doses scenarios 
doses <- c(0, 5, 25, 50, 100)
nd <- length(doses)

# median survival time for placebo dose
mst.control <- 4 

# shape parameter
sigma.true <- 0.5

# maximum hazard ratio between active dose and placebo dose 
hr.ratio <- 4  
# minimum hazard ratio between active dose and placebo dose
hr.Delta <- 2 

# hazard rate for placebo dose
placEff <- log(mst.control/(log(2)^sigma.true)) 

# maximum hazard rate for active dose
maxEff <- log((mst.control*(hr.ratio^sigma.true))/(log(2)^sigma.true))

# minimum hazard rate for active dose
minEff.Delta <- log((mst.control*(hr.Delta^sigma.true))/(log(2)^sigma.true))
Delta <- (minEff.Delta - placEff)

## MCP Parameters  
emax <- guesst(d = doses[4], p = 0.5, model="emax")
exp <- guesst(d = doses[4], p = 0.1, model="exponential", Maxd = doses[nd])
logit <- guesst(d = c(doses[3], doses[4]), p = c(0.1,0.8), "logistic",  
	Maxd= doses[nd])
betam <- guesst(d = doses[2], p = 0.3, "betaMod", scal=120, dMax=50, 
	Maxd= doses[nd])

models.candidate <- Mods(emax = emax, linear = NULL,
                         exponential = exp, logistic = logit,
                         betaMod = betam, doses = doses,
                         placEff = placEff, maxEff = (maxEff- placEff))
plot(models.candidate)

## True Model
model.true <- "emax"
response <- model_response(doses = doses,
                           distr = "weibull", 
                           model.true = model.true, 
                           models.candidate = models.candidate) 
lambda.true <- response$lambda
parm <- list(lambda = lambda.true, sigma = sigma.true)

## Scenario: Censoring 10%
censoring.rate <- 0.1

dt <- data_generator(doses = doses,
                     sample.size = 20,
                     distr = "weibull",
                     parm = parm,                    
                     censoring.rate = censoring.rate)
## Print data
#dt

Simulation to obtain operating characteristics for MCP-Mod design

Description

It simulates dose-finding trials using MCP-Mod design with Maximum Likelihood Estimator and Fisher Information (MLE), Maximum Likelihood Estimator and second-order Fisher Information (MLE2), Cox and Snell's Bias-Corrected Estimator and Fisher Information (BCE), Cox and Snell's Bias-Corrected Estimator and second-order Fisher Information (BCE2), and Firth Bias-Corrected estimators (Firth) as discussed in Diniz, Magalhães and Gallardo.

Usage

mcpmod_simulation(doses, parm, sample.size, model.true,
  models.candidate, selModel = "AIC",significance.level = 0.025,
  Delta, distr = "weibull", censoring.rate = NULL,
  sigma.estimator = NULL, n.cores, seed,n.sim)

Arguments

Value

An object of class mcpmod_simulation with the following components:

mle: a matrix of dimension n.sim \times 4 with results when using the MCP-Mod approach with MLE;

mle2: a matrix of dimension n.sim \times 4 with results when using the MCP-Mod approach with MLE2;

bce: a matrix of dimension n.sim \times 4 with results when using the MCP-Mod approach with BCE;

bce2: a matrix of dimension n.sim \times 4 with results when using the MCP-Mod approach with BCE2;

firth: a matrix of dimension n.sim \times 4 with results when using the MCP-Mod approach with Firth's estimator;

All matrices contain the following columns: (1) the first column indicates whether proof-of-concept (1 = "yes", 0 = "no"), in other words, the p-value of Wald test was statistically significant; (2) the second column indicates whether the true model was selected to estimate the dose-response curve (1 = "yes", 0 = "no") when proof-of-concept is demonstrated; (3) the third column contains the estimated target dose; (4) the fourth column contains the sample size considered in the trial.

conditions: a list containing the conditions of the simulation.

Author(s)

Diniz, M.A., Gallardo D.I., Magalhaes, T.M.

References

Bretz F, Pinheiro JC, Branson M. Combining multiple comparisons and modeling techniques in dose‐response studies. Biometrics. 2005 Sep;61(3):738-48.

Bornkamp B, Pinheiro J, Bretz F. MCPMod: An R package for the design and analysis of dose-finding studies. Journal of Statistical Software. 2009 Feb 20;29:1-23.

Diniz, Márcio A. and Gallardo, Diego I. and Magalhães, Tiago M. (2023). Improved inference for MCP-Mod approach for time-to-event endpoints with small sample sizes. arXiv <doi.org/10.48550/arXiv.2301.00325>

Pinheiro J, Bornkamp B, Glimm E, Bretz F. Model‐based dose finding under model uncertainty using general parametric models. Statistics in medicine. 2014 May 10;33(10):1646-61.

Examples

library(DoseFinding)
library(MCPModBC)

## doses scenarios 
doses <- c(0, 5, 25, 50, 100)
nd <- length(doses)
sample.size <- 25

# shape parameter
sigma.true <- 0.5

# median survival time for placebo dose
mst.control <- 4 

# maximum hazard ratio between active dose and placebo dose 
hr.ratio <- 4  
# minimum hazard ratio between active dose and placebo dose
hr.Delta <- 2 

# hazard rate for placebo dose
placEff <- log(mst.control/(log(2)^sigma.true)) 

# maximum hazard rate for active dose
maxEff <- log((mst.control*(hr.ratio^sigma.true))/(log(2)^sigma.true))

# minimum hazard rate for active dose
minEff.Delta <- log((mst.control*(hr.Delta^sigma.true))/(log(2)^sigma.true))
Delta <- (minEff.Delta - placEff)
	
## MCP Parameters 
significance.level <- 0.05
selModel <- "AIC"

emax <- guesst(d = doses[4], p = 0.5, model="emax")
exp <- guesst(d = doses[4], p = 0.1, model="exponential", Maxd = doses[nd])
logit <- guesst(d = c(doses[3], doses[4]), p = c(0.1,0.8), "logistic",  
	Maxd= doses[nd])
betam <- guesst(d = doses[2], p = 0.3, "betaMod", scal=120, dMax=50, 
	Maxd= doses[nd])

models.candidate <- Mods(emax = emax, linear = NULL,
                         exponential = exp, logistic = logit,
                         betaMod = betam, doses = doses,
                         placEff = placEff, maxEff = (maxEff- placEff))
plot(models.candidate)

## Simulation Parameters
n.sim <- 10
seed <- 1234 
n.cores <- 1

## True Model
model.true <- "emax"
response <- model_response(doses = doses,
                           distr = "weibull", 
                           model.true = model.true, 
                           models.candidate = models.candidate) 
lambda.true <- response$lambda
parm <- list(lambda = lambda.true, sigma = sigma.true)

## Scenario: Censoring 10%
censoring.rate <- 0.1

test <- mcpmod_simulation(doses = doses,
           parm = parm, sample.size = sample.size,
           model.true = model.true,
           models.candidate = models.candidate,
           selModel = selModel,
           significance.level = significance.level,
           Delta = Delta, distr = "weibull",
           censoring.rate = censoring.rate,
           sigma.estimator = "jackknife",
           n.cores = n.cores, seed = seed, n.sim = n.sim)

test
summary(test)

Model Responser

Description

It calculates the model response and parameters of interest for a given distribution.

Usage

model_response(doses, distr, model.true, models.candidate)

Arguments

Value

a data frame with dimension length(doses) \times 3 with the following columns: (1) model response (2) model parameter and (3) doses

Author(s)

Diniz, M.A., Gallardo D.I., Magalhaes, T.M.

References

Diniz, Márcio A. and Gallardo, Diego I. and Magalhães, Tiago M. (2023). Improved inference for MCP-Mod approach for time-to-event endpoints with small sample sizes. arXiv <doi.org/10.48550/arXiv.2301.00325>

Examples

library(DoseFinding)
library(MCPModBC)

## doses scenarios 
doses <- c(0, 5, 25, 50, 100)
nd <- length(doses)

# median survival time for placebo dose
mst.control <- 4 

# shape parameter
sigma.true <- 0.5

# maximum hazard ratio between active dose and placebo dose 
hr.ratio <- 4  
# minimum hazard ratio between active dose and placebo dose
hr.Delta <- 2 

# hazard rate for placebo dose
placEff <- log(mst.control/(log(2)^sigma.true)) 

# maximum hazard rate for active dose
maxEff <- log((mst.control*(hr.ratio^sigma.true))/(log(2)^sigma.true))

# minimum hazard rate for active dose
minEff.Delta <- log((mst.control*(hr.Delta^sigma.true))/(log(2)^sigma.true))
Delta <- (minEff.Delta - placEff)

## MCP Parameters  
emax <- guesst(d = doses[4], p = 0.5, model="emax")
exp <- guesst(d = doses[4], p = 0.1, model="exponential", Maxd = doses[nd])
logit <- guesst(d = c(doses[3], doses[4]), p = c(0.1,0.8), "logistic",  Maxd= doses[nd])
betam <- guesst(d = doses[2], p = 0.3, "betaMod", scal=120, dMax=50, Maxd= doses[nd])

models.candidate <- Mods(emax = emax, linear = NULL,
                         exponential = exp, logistic = logit,
                         betaMod = betam, doses = doses,
                         placEff = placEff, maxEff = (maxEff- placEff))
plot(models.candidate)

## True Model
model.true <- "emax"
response <- model_response(doses = doses,
                           distr = "weibull", 
                           model.true = model.true, 
                           models.candidate = models.candidate) 
response

lambda.true <- response$lambda
parm <- list(lambda = lambda.true, sigma = sigma.true)

Summary of simulation results

Description

It summarizes results of simulations of dose-finding trials following the MCP-Mod approach with bias-corrected and second-order covariance matrices.

Usage

## S3 method for class 'mcpmod_simulation'
summary(object, ...)

Arguments

Value

A data frame with a summary with the information provided by mcpmod_simulation.

References

Diniz, Márcio A. and Gallardo, Diego I. and Magalhães, Tiago M. (2023). Improved inference for MCP-Mod approach for time-to-event endpoints with small sample sizes. arXiv <doi.org/10.48550/arXiv.2301.00325>

Examples

library(DoseFinding)
library(MCPModBC)

## doses scenarios 
doses <- c(0, 5, 25, 50, 100)
nd <- length(doses)
sample.size <- 25

# shape parameter
sigma.true <- 0.5

# median survival time for placebo dose
mst.control <- 4 

# maximum hazard ratio between active dose and placebo dose 
hr.ratio <- 4  
# minimum hazard ratio between active dose and placebo dose
hr.Delta <- 2 

# hazard rate for placebo dose
placEff <- log(mst.control/(log(2)^sigma.true)) 

# maximum hazard rate for active dose
maxEff <- log((mst.control*(hr.ratio^sigma.true))/(log(2)^sigma.true))

# minimum hazard rate for active dose
minEff.Delta <- log((mst.control*(hr.Delta^sigma.true))/(log(2)^sigma.true))
Delta <- (minEff.Delta - placEff)
	
## MCP Parameters 
significance.level <- 0.05
selModel <- "AIC"

emax <- guesst(d = doses[4], p = 0.5, model="emax")
exp <- guesst(d = doses[4], p = 0.1, model="exponential", Maxd = doses[nd])
logit <- guesst(d = c(doses[3], doses[4]), p = c(0.1,0.8), "logistic",  
	Maxd= doses[nd])
betam <- guesst(d = doses[2], p = 0.3, "betaMod", scal=120, dMax=50, 
	Maxd= doses[nd])

models.candidate <- Mods(emax = emax, linear = NULL,
                         exponential = exp, logistic = logit,
                         betaMod = betam, doses = doses,
                         placEff = placEff, maxEff = (maxEff- placEff))
plot(models.candidate)

## Simulation Parameters
n.sim <- 10
seed <- 1234 
n.cores <- 1

## True Model
model.true <- "emax"
response <- model_response(doses = doses,
                           distr = "weibull", 
                           model.true = model.true, 
                           models.candidate = models.candidate) 
lambda.true <- response$lambda
parm <- list(lambda = lambda.true, sigma = sigma.true)

## Scenario: Censoring 10%
censoring.rate <- 0.1

test <- mcpmod_simulation(doses = doses,
           parm = parm,
           sample.size = sample.size,
           model.true = model.true,
           models.candidate = models.candidate,
           selModel = selModel,
           significance.level = significance.level,
           Delta = Delta,
           distr = "weibull",
           censoring.rate = censoring.rate,
           sigma.estimator = "jackknife",
           n.cores = n.cores, seed = seed, n.sim = n.sim)
summary(test)

Print a summary for a object of the weibreg class.

Description

Summarizes the results for a object of the weibreg class.

Usage

## S3 method for class 'weibreg'
summary(object, ...)

## S3 method for class 'weibreg'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

Arguments

Value

A complete summary for the coefficients extracted from a weibreg object.

Functions

• print(weibreg):

References

Diniz, Márcio A. and Gallardo, Diego I. and Magalhães, Tiago M. (2023). Improved inference for MCP-Mod approach for time-to-event endpoints with small sample sizes. arXiv <doi.org/10.48550/arXiv.2301.00325>

Examples

require(survival)
set.seed(2100)

##Generating covariates
n=20 
x<-runif(n, max=10)
lambda<-exp(1.2-0.5*x)
sigma<-1.5

##Drawing T from Weibull model and fixing censoring at 1.5
T<-rweibull(n, shape=1/sigma, scale=lambda)
L<-rep(1.5, n)

##Defining the observed times and indicators of failure
t<-pmin(T,L)
delta<-ifelse(T<=L, 1, 0)
data=data.frame(t=t, delta=delta, x=x)

##Fitting for Weibull regression model

##Traditional MLE with corrected variance
ex1=weibfit(Surv(t,delta)~x, data=data, L=L, estimator="MLE", 
	corrected.var=TRUE)
summary(ex1)

##BCE without corrected variance
ex2=weibfit(Surv(t,delta)~x, data=data, L=L, estimator="Firth", 
	corrected.var=FALSE)
summary(ex2)

##BCE with corrected variance
ex3=weibfit(Surv(t,delta)~x, data=data, L=L, estimator="BCE", 
	corrected.var=TRUE)
summary(ex3)

##Firth's correction without corrected variance
ex4=weibfit(Surv(t,delta)~x, data=data, L=L, estimator="BCE", 
	corrected.var=FALSE)
summary(ex4)

Computes different estimators for the censored Weibull regression model

Description

Computes the maximum likelihood estimators (MLE) for the censored Weibull regression model. The bias-corrected estimators based on the Cox and Snell's and Firth's methods also are available. In addition, for the covariance matrix the corrected estimators discussed in Magalhaes et al. 2021 also are available.

Usage

weibfit(formula, data, L = Inf, estimator = "MLE", 
	corrected.var = FALSE)

Arguments

Details

The Weibull distribution considered here has probability density function

f(t;\lambda, \sigma)=\frac{1}{\sigma \lambda^{1/\sigma}}t^{1/\sigma-1}\exp\left\{-\left(\frac{t}{\lambda}\right)^{1/\sigma}\right\}, \quad t, \sigma, \lambda>0.

The regression structure is incorporated as

\log(\lambda_i)={\bm x}_i^\top {\bm \beta}, \quad i=1,\ldots,n.

For the computation of the bias-corrected estimators, \sigma is assumed as fixed in the jackknife estimator based on the traditional MLE.

The Fisher information matrix for \bm \beta is given by {\bm K}=\sigma^{-2} {\bm X}^\top {\bm W} {\bm X}, where {\bm X}=({\bm x}_1,\ldots,{\bm x}_n)^\top, {\bm W}=\mbox{diag}(w_1,\ldots,w_n), and

w_i=E\left[\exp\left(\frac{y_i-\log \lambda_i}{\sigma}\right)\right]=q \times \left\{ 1 - \exp\left[ -L_i^{1/\sigma} \exp(-\mu_i/\sigma) \right] \right\} + (1-q)\times \left(r/n\right),

with q = P\left(W_{(r)}\leq \log L_i\right) and W_{(r)} denoting the rth order statistic from W_1, \ldots, W_n, with q=1 and q=0 for types I and II censoring, respectively. (See Magalhaes et al. 2019 for details).

The bias-corrected maximum likelihood estimator based on the Cox and Snell's method (say \widetilde{\bm \beta}) is based on a corrective approach given by \widetilde{\beta}=\widehat{\beta}-B(\widehat{\beta}), where

B({\bm \beta})= - \frac{1}{2 \sigma^3} {\bm P} {\bm Z}_d \left({\bm W} + 2 \sigma {\bm W}^{\prime}\right) {\bm 1},

with {\bm P} = {\bm K}^{-1} {\bm X}^{\top}, {\bm Z} = {\bm X} {\bm K}^{-1} {\bm X}^{\top}, {\bm Z}_d is a diagonal matrix with diagonal given by the diagonal of {\bm Z}, {\bm W}^{\prime} = diag(w_1^{\prime}, \ldots, w_n^{\prime}), w_i^{\prime} = - \sigma^{-1} L_i^{1/\sigma} \exp\{ -L_i^{1/\sigma} \exp(-\mu_i/\sigma) - \mu_i/\sigma \} and {\bm 1} is a n-dimensional vector of ones.

The bias-corrected maximum likelihood estimator based on the Firth's method (say \check{\bm \beta}) is based on a preventive approach, which is the solution for the equation {\bm U}_{{\bm \beta}}^{\star} = {\bm 0}, where

{\bm U}_{{\bm \beta}}^{\star} = {\bm U}_{{\bm \beta}} - {\bm K}_{{\bm \beta} {\bm \beta}} B({\bm \beta}).

The covariance correction is based on the general result of Magalhaes et al. 2021 given by

\mbox{{\bf Cov}}_{\bm 2}^{\bm \tau}({\bm \beta}^{\star}) = {\bm K}^{-1} + {\bm K}^{-1} \left\{ {\bm \Delta} + {\bm \Delta}^{\top} \right\} {\bm K}^{-1} + \mathcal{O}(n^{-3})

where {\Delta} = -0.5 {\Delta}^{(1)} + 0.25 {\Delta}^{(2)} + 0.5 \tau_2 {\Delta}^{(3)}, with

\Delta^{(1)} = \frac{1}{\sigma^4} {\bm X}^{\top} {\bm W}^{\star} {\bm Z}_{d} {\bm X},

\Delta^{(2)} = - \frac{1}{\sigma^6} {\bm X}^{\top} \left[ {\bm W} {\bm Z}^{(2)} {\bm W} - 2 \sigma {\bm W} {\bm Z}^{(2)} {\bm W}^{\prime} - 6 \sigma^2 {\bm W}^{\prime} {\bm Z}^{(2)} {\bm W}^{\prime} \right] {\bm X},

and

\Delta^{(3)} = \frac{1}{\sigma^5} {\bm X}^{\top} {\bm W}^{\prime} {\bm W}^{\star\star} {\bm X},

where {\bm W}^{\star} = diag(w_1^{\star}, \ldots, w_n^{\star}), w_i^{\star} = w_i (w_i -2) - 2 \sigma w_i^{\prime} + \sigma \tau_1 (w_i^{\prime} + 2 \sigma w_i^{\prime\prime}), {\bm Z}^{(2)} = {\bm Z} \odot {\bm Z}, with \odot representing a direct product of matrices (Hadamard product), {\bm W}^{\star\star} is a diagonal matrix, with {\bm Z} ( {\bm W} + 2 \sigma {\bm W}^{\prime}) {\bm Z}_{d} {\bm 1} as its diagonal, {\bm W}^{\prime\prime} = diag(w_1^{\prime\prime}, \ldots, w_n^{\prime\prime}), w_i^{\prime\prime} = - \sigma^{-1} w_i^{\prime} \left[ L_i^{1/\sigma} \exp(-\mu_i/\sigma) - 1 \right], {\bm \tau} = (\tau_1, \tau_2) = (1, 1) indicating the second-order covariance matrix of the MLE {\bm \beta}^{\star} = \widehat{{\bm \beta}} denoted by \mbox{{\bf Cov}}_{\bm 2}(\widehat{{\bm \beta}}) and {\bm \tau} = (0, -1) indicating the second-order covariance matrix of the BCE {\bm \beta}^{\star} = \widetilde{{\bm \beta}} denoted by \mbox{{\bf Cov}}_{\bm 2}(\widetilde{{\bm \beta}}).

Value

Author(s)

Gallardo D.I., Diniz, M.A., Magalhaes, T.M.

References

Cox, D.R., Snell E.J. A general definition of residuals Journal of the Royal Statistical Society. Series B (Methodological). 1968;30:248-275.

Diniz, Márcio A. and Gallardo, Diego I. and Magalhães, Tiago M. (2023). Improved inference for MCP-Mod approach for time-to-event endpoints with small sample sizes. arXiv <doi.org/10.48550/arXiv.2301.00325>

Firth, D. Bias reduction of maximum likelihood estimates Biometrika. 1993;80:27-38.

Magalhaes Tiago M., Botter Denise A., Sandoval Monica C. A general expression for second- order covariance matrices - an application to dispersion models Brazilian Journal of Probability and Statistics. 2021;35:37-49.

Examples

require(survival)
set.seed(2100)

##Generating covariates
n=20; 
x<-runif(n, max=10)
lambda<-exp(1.2-0.5*x); sigma<-1.5

##Drawing T from Weibull model and fixing censoring at 1.5
T<-rweibull(n, shape=1/sigma, scale=lambda); L<-rep(1.5, n)

##Defining the observed times and indicators of failure
y<-pmin(T,L); 
delta<-ifelse(T<=L, 1, 0)
data=data.frame(y=y, delta=delta, x=x)

##Fitting for Weibull regression model

##Traditional MLE with corrected variance
ex1=weibfit(Surv(y,delta)~x, data=data, L=L, estimator="MLE", 
	corrected.var=TRUE)
summary(ex1)

##BCE without corrected variance
ex2=weibfit(Surv(y,delta)~x, data=data, L=L, estimator="BCE", 
	corrected.var=FALSE)
summary(ex2)

##BCE with corrected variance
ex3=weibfit(Surv(y,delta)~x, data=data, L=L, estimator="BCE", 
	corrected.var=TRUE)
summary(ex3)

##Firth's correction without corrected variance
ex4=weibfit(Surv(y,delta)~x, data=data, L=L, estimator="BCE", 
	corrected.var=FALSE)
summary(ex4)