R/functions_bias.R
In Sterniii3/drugdevelopR: Utility-Based Optimal Phase II/III Drug Development Planning
Defines functions
utility_R2
EPsProg_R2
Ed3_R2
Epgo_R2
utility_R
EPsProg_R
Ed3_R
utility_L2
EPsProg_L2
Ed3_L2
Epgo_L2
utility_L
EPsProg_L
Ed3_L
Documented in
Ed3_L
Ed3_L2
Ed3_R
Ed3_R2
Epgo_L2
Epgo_R2
EPsProg_L
EPsProg_L2
EPsProg_R
EPsProg_R2
utility_L
utility_L2
utility_R
utility_R2
# 1.1. conservative sample size calculation: use lower bound of one-sided confidence intervall
##############################################################################################
# prior distribution
# as above
# expected probability to go to phase III
# as above
# expected number of events for phase III when going to phase III
#' Expected sample size for phase III for bias adjustment programs and time-to-event outcomes
#'
#' To discount for overoptimistic results in phase II when calculating the optimal sample size in phase III,
#' it is necessary to use the functions `Ed3_L()`, `Ed3_L2()`, `Ed3_R()` and `Ed3_R2()`.
#' Each function describes a specific case:
#' - `Ed3_L()`: calculates the optimal sample size for an additive adjustment factor (i.e. adjust the lower bound of the one-sided confidence interval),
#' however the go-decision is not affected by the bias adjustment
#' - `Ed3_L2()`: calculates the optimal sample size for an additive adjustment factor (i.e. adjust the lower bound of the one-sided confidence interval)
#' when the go-decision is also affected by the bias adjustment
#' - `Ed3_R()`: calculates the optimal sample size for a multiplicative adjustment factor (i.e. use estimate with a retention factor),
#' however the go-decision is not affected by the bias adjustment
#' - `Ed3_R2()`: calculates the optimal sample size for a multiplicative adjustment factor (i.e. use estimate with a retention factor)
#' when the go-decision is also affected by the bias adjustment
#' @param HRgo threshold value for the go/no-go decision rule
#' @param d2 total events for phase II; must be even number
#' @param Adj adjustment parameter
#' @param alpha significance level
#' @param beta `1-beta` power for calculation of sample size for phase III
#' @param w weight for mixture prior distribution
#' @param hr1 first assumed true treatment effect on HR scale for prior distribution
#' @param hr2 second assumed true treatment effect on HR scale for prior distribution
#' @param id1 amount of information for `hr1` in terms of number of events
#' @param id2 amount of information for `hr2` in terms of number of events
#' @param fixed choose if true treatment effects are fixed or random, if TRUE `hr1` is used as fixed effect
#' @return The output of the functions `Ed3_L`, `Ed3_L2`, `Ed3_R` and `Ed3_R2` is the expected number of participants in phase III.
#' @examples res <- Ed3_L(HRgo = 0.8, d2 = 50, Adj = 0.4,
#' alpha = 0.025, beta = 0.1, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- Ed3_L2(HRgo = 0.8, d2 = 50, Adj = 0.4,
#' alpha = 0.025, beta = 0.1, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- Ed3_R(HRgo = 0.8, d2 = 50, Adj = 0.9,
#' alpha = 0.025, beta = 0.1, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- Ed3_R2(HRgo = 0.8, d2 = 50, Adj = 0.9,
#' alpha = 0.025, beta = 0.1, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#'
#' @name Ed3_bias
#' @export
#' @keywords internal
Ed3_L<-function(HRgo, d2, Adj, alpha, beta, w, hr1, hr2, id1, id2, fixed){
if(fixed){
theta <- -log(hr1)
int = try(integrate(function(y){
sapply(y,function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y-qnorm(1-Adj)*sqrt(4/d2))^2))*
dnorm(y,
mean=theta,
sd=sqrt(4/d2))
})
}, -log(HRgo),Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}else{
int = try(integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y-qnorm(1-Adj)*sqrt(4/d2))^2))*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo),Inf)$value
})
}, -Inf, Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}
}
#' Expected probability of a successful program for bias adjustment programs with time-to-event outcomes
#'
#' To discount for overoptimistic results in phase II when calculating the optimal sample size in phase III,
#' it is necessary to use the following functions, which each describe a specific case:
#' - `EPsProg_L()`: calculates the expected probability of a successful for an additive adjustment factor (i.e. adjust the lower bound of the one-sided confidence interval),
#' however the go-decision is not affected by the bias adjustment
#' - `EPsProg_L2()`: calculates the expected probability of a successful for an additive adjustment factor (i.e. adjust the lower bound of the one-sided confidence interval)
#' when the go-decision is also affected by the bias adjustment
#' - `EPsProg_R()`: calculates the expected probability of a successful for a multiplicative adjustment factor (i.e. use estimate with a retention factor),
#' however the go-decision is not affected by the bias adjustment
#' - `EPsProg_R2()`: calculates the expected probability of a successful for a multiplicative adjustment factor (i.e. use estimate with a retention factor)
#' when the go-decision is also affected by the bias adjustment
#' @param HRgo threshold value for the go/no-go decision rule
#' @param d2 total events for phase II; must be even number
#' @param Adj adjustment parameter
#' @param alpha significance level
#' @param beta `1-beta` power for calculation of sample size for phase III
#' @param step1 lower boundary for effect size
#' @param step2 upper boundary for effect size
#' @param w weight for mixture prior distribution
#' @param hr1 first assumed true treatment effect on HR scale for prior distribution
#' @param hr2 second assumed true treatment effect on HR scale for prior distribution
#' @param id1 amount of information for `hr1` in terms of number of events
#' @param id2 amount of information for `hr2` in terms of number of events
#' @param fixed choose if true treatment effects are fixed or random, if TRUE `hr1` is used as fixed effect
#' @return The output of the functions `EPsProg_L()`, `EPsProg_L2()`, `EPsProg_R()` and `EPsProg_R2()` is the expected probability of a successful program.
#' @examples res <- EPsProg_L(HRgo = 0.8, d2 = 50, Adj = 0.4,
#' alpha = 0.025, beta = 0.1,
#' step1 = 1, step2 = 0.95,
#' w = 0.3, hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- EPsProg_L2(HRgo = 0.8, d2 = 50, Adj = 0.4,
#' alpha = 0.025, beta = 0.1,
#' step1 = 1, step2 = 0.95,
#' w = 0.3, hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- EPsProg_R(HRgo = 0.8, d2 = 50, Adj = 0.9,
#' alpha = 0.025, beta = 0.1,
#' step1 = 1, step2 = 0.95,
#' w = 0.3, hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- EPsProg_R2(HRgo = 0.8, d2 = 50, Adj = 0.9,
#' alpha = 0.025, beta = 0.1,
#' step1 = 1, step2 = 0.95,
#' w = 0.3, hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' @name EPsProg_bias
#' @export
#' @keywords internal
EPsProg_L<-function(HRgo, d2, Adj, alpha, beta, step1, step2, w, hr1, hr2, id1, id2, fixed){
c=(qnorm(1-alpha)+qnorm(1-beta))^2
if(fixed){
theta <- -log(hr1)
return(integrate(function(y){
sapply(y,function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=theta/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=theta/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1) )*
dnorm(y,
mean=theta,
sd=sqrt(4/d2))
})
}, -log(HRgo),Inf)$value)
}else{
return(
integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=x/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=x/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1) )*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo),Inf)$value
})
}, -Inf, Inf)$value
)
}
}
#' Utility function for bias adjustment programs with time-to-event outcomes.
#'
#' The utility function calculates the expected utility of our drug development program and is given as gains minus costs and depends on the parameters and the expected probability of a successful program.
#' The utility is in a further step maximized by the `optimal_bias()` function.
#' @param HRgo threshold value for the go/no-go decision rule
#' @param d2 total events for phase II; must be even number
#' @param Adj adjustment parameter
#' @param alpha significance level
#' @param beta `1-beta` power for calculation of sample size for phase III
#' @param xi2 event rate for phase II
#' @param xi3 event rate for phase III
#' @param w weight for mixture prior distribution
#' @param hr1 first assumed true treatment effect on HR scale for prior distribution
#' @param hr2 second assumed true treatment effect on HR scale for prior distribution
#' @param id1 amount of information for `hr1` in terms of number of events
#' @param id2 amount of information for `hr2` in terms of number of events
#' @param c2 variable per-patient cost for phase II
#' @param c3 variable per-patient cost for phase III
#' @param c02 fixed cost for phase II
#' @param c03 fixed cost for phase III
#' @param K constraint on the costs of the program, default: Inf, e.g. no constraint
#' @param N constraint on the total expected sample size of the program, default: Inf, e.g. no constraint
#' @param S constraint on the expected probability of a successful program, default: -Inf, e.g. no constraint
#' @param steps1 lower boundary for effect size category `"small"` in RR scale, default: 1
#' @param stepm1 lower boundary for effect size category `"medium"` in RR scale = upper boundary for effect size category "small" in RR scale, default: 0.95
#' @param stepl1 lower boundary for effect size category `"large"` in RR scale = upper boundary for effect size category "medium" in RR scale, default: 0.85
#' @param b1 expected gain for effect size category `"small"`
#' @param b2 expected gain for effect size category `"medium"`
#' @param b3 expected gain for effect size category `"large"`
#' @param fixed choose if true treatment effects are fixed or random, if TRUE `hr1` is used as fixed effect
#' @return The output of the functions `utility_L()`, `utility_L2()`, `utility_R()` and `utility_R2()` is the expected utility of the program.
#' @examples res <- utility_L(d2 = 50, HRgo = 0.8, Adj = 0.4, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, xi2 = 0.7, xi3 = 0.7,
#' alpha = 0.025, beta = 0.1,
#' c2 = 0.75, c3 = 1, c02 = 100, c03 = 150,
#' K = Inf, N = Inf, S = -Inf,
#' steps1 = 1, stepm1 = 0.95, stepl1 = 0.85,
#' b1 = 1000, b2 = 2000, b3 = 3000,
#' fixed = TRUE)
#' res <- utility_L2(d2 = 50, HRgo = 0.8, Adj = 0.4, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, xi2 = 0.7, xi3 = 0.7,
#' alpha = 0.025, beta = 0.1,
#' c2 = 0.75, c3 = 1, c02 = 100, c03 = 150,
#' K = Inf, N = Inf, S = -Inf,
#' steps1 = 1, stepm1 = 0.95, stepl1 = 0.85,
#' b1 = 1000, b2 = 2000, b3 = 3000,
#' fixed = TRUE)
#' res <- utility_R(d2 = 50, HRgo = 0.8, Adj = 0.9, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, xi2 = 0.7, xi3 = 0.7,
#' alpha = 0.025, beta = 0.1,
#' c2 = 0.75, c3 = 1, c02 = 100, c03 = 150,
#' K = Inf, N = Inf, S = -Inf,
#' steps1 = 1, stepm1 = 0.95, stepl1 = 0.85,
#' b1 = 1000, b2 = 2000, b3 = 3000,
#' fixed = TRUE)
#' res <- utility_R2(d2 = 50, HRgo = 0.8, Adj = 0.9, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, xi2 = 0.7, xi3 = 0.7,
#' alpha = 0.025, beta = 0.1,
#' c2 = 0.75, c3 = 1, c02 = 100, c03 = 150,
#' K = Inf, N = Inf, S = -Inf,
#' steps1 = 1, stepm1 = 0.95, stepl1 = 0.85,
#' b1 = 1000, b2 = 2000, b3 = 3000,
#' fixed = TRUE)
#' @name utility_bias
#' @export
#' @keywords internal
utility_L <- function(d2, HRgo, Adj, w, hr1, hr2, id1, id2,
alpha, beta, xi2, xi3,
c2, c3, c02, c03,
K, N, S,
steps1, stepm1, stepl1,
b1, b2, b3,
fixed){
steps2 <- stepm1
stepm2 <- stepl1
stepl2 <- 0
d3 <- Ed3_L(HRgo = HRgo, d2 = d2, Adj = Adj,
alpha = alpha, beta = beta, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
if(is.na(d3)){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# sample size is rounded up to next even natural number
n2 <- ceiling(d2 * (1/xi2))
if(round(n2/2) != n2 / 2) {n2 <- n2 + 1}
n3 <- ceiling(d3 * (1/xi3))
if(round(n3/2) != n3 / 2) {n3 <- n3 + 1}
# expected number of events is rounded to natural number for presentation
d3 <- ceiling(d3)
if(n2+n3>N){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
pg <- Epgo_tte(HRgo = HRgo, d2 = d2, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
K2 <- c02 + c2 * n2 # cost phase II
K3 <- c03 * pg + c3 * n3 # cost phase III
if(K2+K3>K){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# probability of a successful program; small, medium, large effect size
prob1 <- EPsProg_L(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = steps1, step2 = steps2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob2 <- EPsProg_L(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepm1, step2 = stepm2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob3 <- EPsProg_L(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepl1, step2 = stepl2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
SP <- prob1 + prob2 + prob3
if(SP<S){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
G <- b1 * prob1 + b2 * prob2 + b3 * prob3 #gain
EU <- - K2 - K3 + G
return(c(EU, d3, SP, pg, K2, K3, prob1, prob2, prob3, n2, n3))
#output: expected utility Eud, En3, EsP, Epgo, cost phase II and III
}
}
}
}
}
# 1.2. conservative decision rule and sample size calculation:
# use lower bound of one-sided confidence intervall
##############################################################################################
# prior distribution
# as above
#' Expected probability to go to phase III for bias adjustment programs with time-to-event outcomes
#'
#' In the case we do not only want do discount for overoptimistic results in phase II when calculating the sample size in phase III,
#' but also when deciding whether to go to phase III or not the functions `Epgo_L2` and `Epgo_R2` are necessary.
#' The function `Epgo_L2` uses an additive adjustment parameter (i.e. adjust the lower bound of the one-sided confidence interval),
#' the function `Epgo_R2` uses a multiplicative adjustment parameter (i.e. use estimate with a retention factor)
#' @param HRgo threshold value for the go/no-go decision rule
#' @param d2 total number of events for phase II; must be even number
#' @param Adj adjustment parameter
#' @param w weight for mixture prior distribution
#' @param hr1 first assumed true treatment effect on HR scale for prior distribution
#' @param hr2 second assumed true treatment effect on HR scale for prior distribution
#' @param id1 amount of information for `hr1` in terms of number of events
#' @param id2 amount of information for `hr2` in terms of number of events
#' @param fixed choose if true treatment effects are fixed or random, if TRUE `hr1` is used as fixed effect
#' @return The output of the functions `Epgo_L2` and `Epgo_R2` is the expected probability to go to phase III with conservative decision rule and sample size calculation.
#' @examples res <- Epgo_L2(HRgo = 0.8, d2 = 50, Adj = 0.4,
#' w = 0.3, hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- Epgo_R2(HRgo = 0.8, d2 = 50, Adj = 0.9,
#' w = 0.3, hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' @name Epgo_bias
#' @export
#' @keywords internal
Epgo_L2<-function(HRgo, d2, Adj, w, hr1, hr2, id1, id2, fixed){
if(fixed){
return(
pnorm((-log(hr1)+log(HRgo)-qnorm(1-Adj)*sqrt(4/d2))/sqrt(4/d2))
)
}else(
return(
integrate(function(x){
sapply(x,function(x){
pnorm((x+log(HRgo)-qnorm(1-Adj)*sqrt(4/d2))/sqrt(4/d2))*prior_tte(x,w,hr1,hr2,id1,id2)
})
}, -Inf, Inf)$value
)
)
}
# expected number of events for phase III when going to phase III
#' @rdname Ed3_bias
#' @export
Ed3_L2<-function(HRgo, d2, Adj, alpha, beta, w, hr1, hr2, id1, id2, fixed){
if(fixed){
int = try(
integrate(function(y){
sapply(y,function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y-qnorm(1-Adj)*sqrt(4/d2))^2))*
dnorm(y,
mean=-log(hr1),
sd=sqrt(4/d2))
})
}, -log(HRgo)+qnorm(1-Adj)*sqrt(4/d2),Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}else{
int = try(integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y-qnorm(1-Adj)*sqrt(4/d2))^2))*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo)+qnorm(1-Adj)*sqrt(4/d2),Inf)$value
})
}, -Inf, Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}
}
# expected probability of a successful program
#' @rdname EPsProg_bias
#' @keywords internal
#' @export
EPsProg_L2<-function(HRgo, d2, Adj, alpha, beta, step1, step2, w, hr1, hr2, id1, id2, fixed){
c=(qnorm(1-alpha)+qnorm(1-beta))^2
if(fixed){
return(
integrate(function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=-log(hr1)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=-log(hr1)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1) )*
dnorm(y,
mean=-log(hr1),
sd=sqrt(4/d2))
}, -log(HRgo)+qnorm(1-Adj)*sqrt(4/d2),Inf)$value
)
}else{
return(
integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=x/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=x/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1) )*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo)+qnorm(1-Adj)*sqrt(4/d2),Inf)$value
})
}, -Inf, Inf)$value
)
}
}
# Utility function
#' @rdname utility_bias
#' @keywords internal
#' @export
utility_L2 <- function(d2, HRgo, Adj, w, hr1, hr2, id1, id2,
alpha, beta, xi2, xi3,
c2, c3, c02, c03,
K, N, S,
steps1, stepm1, stepl1,
b1, b2, b3,
fixed){
steps2 <- stepm1
stepm2 <- stepl1
stepl2 <- 0
d3 <- Ed3_L2(HRgo = HRgo, d2 = d2, Adj = Adj,
alpha = alpha, beta = beta, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
if(is.na(d3)){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# sample size is rounded up to next even natural number
n2 <- ceiling(d2 * (1/xi2))
if(round(n2/2) != n2 / 2) {n2 <- n2 + 1}
n3 <- ceiling(d3 * (1/xi3))
if(round(n3/2) != n3 / 2) {n3 <- n3 + 1}
# expected number of events is rounded to natural number for presentation
d3 <- ceiling(d3)
if(n2+n3>N){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
pg <- Epgo_L2(HRgo = HRgo, d2 = d2, Adj = Adj, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
K2 <- c02 + c2 * n2 # cost phase II
K3 <- c03 * pg + c3 * n3 # cost phase III
if(K2+K3>K){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# probability of a successful program; small, medium, large effect size
prob1 <- EPsProg_L2(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = steps1, step2 = steps2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob2 <- EPsProg_L2(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepm1, step2 = stepm2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob3 <- EPsProg_L2(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepl1, step2 = stepl2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
SP <- prob1 + prob2 + prob3
if(SP<S){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
G <- b1 * prob1 + b2 * prob2 + b3 * prob3 #gain
EU <- - K2 - K3 + G
return(c(EU, d3, SP, pg, K2, K3, prob1, prob2, prob3, n2, n3))
#output: expected utility Eud, En3, EsP, Epgo, cost phase II and III
}
}
}
}
}
# 2.1. conservative sample size calculation: use estimate with retetion factor
##############################################################################################
# prior distribution
# as above
# expected probability to go to phase III
# as above
# expected number of events for phase III when going to phase III
#' @rdname Ed3_bias
#' @keywords internal
#' @export
Ed3_R<-function(HRgo, d2, Adj, alpha, beta, w, hr1, hr2, id1, id2, fixed){
if(fixed){
theta <- -log(hr1)
int = try(integrate(function(y){
sapply(y,function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y*Adj)^2))*
dnorm(y,
mean=theta,
sd=sqrt(4/d2))
})
}, -log(HRgo),Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}else{
int = try(integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y*Adj)^2))*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo),Inf)$value
})
}, -Inf, Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}
}
# expected probability of a successful program
#' @rdname EPsProg_bias
#' @keywords internal
#' @export
EPsProg_R<-function(HRgo, d2, Adj, alpha, beta, step1, step2, w, hr1, hr2, id1, id2, fixed){
c=(qnorm(1-alpha)+qnorm(1-beta))^2
if(fixed){
theta <- -log(hr1)
return( integrate(function(y){
sapply(y,function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y*Adj)^2/c)),
mean=theta/(sqrt((y*Adj)^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y*Adj)^2/c)),
mean=theta/(sqrt((y*Adj)^2/c)),
sd=1) )*
dnorm(y,
mean=theta,
sd=sqrt(4/d2))
})
}, -log(HRgo),Inf)$value)
}else{
return(
integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y*Adj)^2/c)),
mean=x/(sqrt((y*Adj)^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y*Adj)^2/c)),
mean=x/(sqrt((y*Adj)^2/c)),
sd=1) )*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo),Inf)$value
})
}, -Inf, Inf)$value
)
}
}
# Utility function
#' @rdname utility_bias
#' @keywords internal
#' @export
utility_R <- function(d2, HRgo, Adj, w, hr1, hr2, id1, id2,
alpha, beta, xi2, xi3,
c2, c3, c02, c03,
K, N, S,
steps1, stepm1, stepl1,
b1, b2, b3,
fixed){
steps2 <- stepm1
stepm2 <- stepl1
stepl2 <- 0
d3 <- Ed3_R(HRgo = HRgo, d2 = d2, Adj = Adj,
alpha = alpha, beta = beta, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
if(is.na(d3)){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# sample size is rounded up to next even natural number
n2 <- ceiling(d2 * (1/xi2))
if(round(n2/2) != n2 / 2) {n2 <- n2 + 1}
n3 <- ceiling(d3 * (1/xi3))
if(round(n3/2) != n3 / 2) {n3 <- n3 + 1}
# expected number of events is rounded to natural number for presentation
d3 <- ceiling(d3)
if(n2+n3>N){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
pg <- Epgo_tte(HRgo = HRgo, d2 = d2, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
K2 <- c02 + c2 * n2 # cost phase II
K3 <- c03 * pg + c3 * n3 # cost phase III
if(K2+K3>K){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# probability of a successful program; small, medium, large effect size
prob1 <- EPsProg_R(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = steps1, step2 = steps2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob2 <- EPsProg_R(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepm1, step2 = stepm2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob3 <- EPsProg_R(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepl1, step2 = stepl2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
SP <- prob1 + prob2 + prob3
if(SP<S){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
G <- b1 * prob1 + b2 * prob2 + b3 * prob3 #gain
EU <- - K2 - K3 + G
return(c(EU, d3, SP, pg, K2, K3, prob1, prob2, prob3, n2, n3))
#output: expected utility Eud, En3, EsP, Epgo, cost phase II and III
}
}
}
}
}
# 2.2. conservative decision rule and sample size calculation:
# use estimate with retetion factor
##############################################################################################
# prior distribution
# as above
# expected probability to go to phase III
#' @rdname Epgo_bias
#' @keywords internal
#' @export
Epgo_R2<-function(HRgo, d2, Adj, w, hr1, hr2, id1, id2, fixed){
if(fixed){
return(
pnorm((-log(hr1)+(log(HRgo)/Adj))/sqrt(4/d2))
)
}else{
return(
integrate(function(x){
sapply(x,function(x){
pnorm((x+(log(HRgo)/Adj))/sqrt(4/d2))*prior_tte(x,w,hr1,hr2,id1,id2)
})
}, -Inf, Inf)$value
)
}
}
# expected number of events for phase III when going to phase III
#' @rdname Ed3_bias
#' @keywords internal
#' @export
Ed3_R2<-function(HRgo, d2, Adj, alpha, beta, w, hr1, hr2, id1, id2, fixed){
if(fixed){
theta <- -log(hr1)
int = try(integrate(function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y*Adj)^2))*
dnorm(y,
mean=theta,
sd=sqrt(4/d2))
}, -log(HRgo)/Adj,Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}else{
int = try(integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y*Adj)^2))*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo)/Adj,Inf)$value
})
}, -Inf, Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}
}
# expected probability of a successful program
#' @rdname EPsProg_bias
#' @keywords internal
#' @export
EPsProg_R2<-function(HRgo, d2, Adj, alpha, beta, step1, step2, w, hr1, hr2, id1, id2, fixed){
c=(qnorm(1-alpha)+qnorm(1-beta))^2
if(fixed){
return(
integrate(function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y*Adj)^2/c)),
mean=-log(hr1)/(sqrt((y*Adj)^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y*Adj)^2/c)),
mean=-log(hr1)/(sqrt((y*Adj)^2/c)),
sd=1) )*
dnorm(y,
mean=-log(hr1),
sd=sqrt(4/d2))
}, -log(HRgo)/Adj,Inf)$value
)
}else{
return(
integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y*Adj)^2/c)),
mean=x/(sqrt((y*Adj)^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y*Adj)^2/c)),
mean=x/(sqrt((y*Adj)^2/c)),
sd=1) )*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo)/Adj,Inf)$value
})
}, -Inf, Inf)$value
)
}
}
# Utility function
#' @rdname utility_bias
#' @keywords internal
#' @export
utility_R2 <- function(d2, HRgo, Adj, w, hr1, hr2, id1, id2,
alpha, beta, xi2, xi3,
c2, c3, c02, c03,
K, N, S,
steps1, stepm1, stepl1,
b1, b2, b3,
fixed){
steps2 <- stepm1
stepm2 <- stepl1
stepl2 <- 0
d3 <- Ed3_R2(HRgo = HRgo, d2 = d2, Adj = Adj,
alpha = alpha, beta = beta, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
if(is.na(d3)){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# sample size is rounded up to next even natural number
n2 <- ceiling(d2 * (1/xi2))
if(round(n2/2) != n2 / 2) {n2 <- n2 + 1}
n3 <- ceiling(d3 * (1/xi3))
if(round(n3/2) != n3 / 2) {n3 <- n3 + 1}
# expected number of events is rounded to natural number for presentation
d3 <- ceiling(d3)
if(n2+n3>N){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
pg <- Epgo_R2(HRgo = HRgo, d2 = d2, Adj = Adj, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
K2 <- c02 + c2 * n2 # cost phase II
K3 <- c03 * pg + c3 * n3 # cost phase III
if(K2+K3>K){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# probability of a successful program; small, medium, large effect size
prob1 <- EPsProg_R2(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = steps1, step2 = steps2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob2 <- EPsProg_R2(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepm1, step2 = stepm2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob3 <- EPsProg_R2(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepl1, step2 = stepl2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
SP <- prob1 + prob2 + prob3
if(SP<S){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
G <- b1 * prob1 + b2 * prob2 + b3 * prob3 #gain
EU <- - K2 - K3 + G
return(c(EU, d3, SP, pg, K2, K3, prob1, prob2, prob3, n2, n3))
#output: expected utility Eud, En3, EsP, Epgo, cost phase II and III
}
}
}
}
}
Sterniii3/drugdevelopR documentation built on Jan. 26, 2024, 6:17 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions utility_R2 EPsProg_R2 Ed3_R2 Epgo_R2 utility_R EPsProg_R Ed3_R utility_L2 EPsProg_L2 Ed3_L2 Epgo_L2 utility_L EPsProg_L Ed3_L
Documented in Ed3_L Ed3_L2 Ed3_R Ed3_R2 Epgo_L2 Epgo_R2 EPsProg_L EPsProg_L2 EPsProg_R EPsProg_R2 utility_L utility_L2 utility_R utility_R2
# 1.1. conservative sample size calculation: use lower bound of one-sided confidence intervall
##############################################################################################
# prior distribution
# as above
# expected probability to go to phase III
# as above
# expected number of events for phase III when going to phase III
#' Expected sample size for phase III for bias adjustment programs and time-to-event outcomes
#'
#' To discount for overoptimistic results in phase II when calculating the optimal sample size in phase III,
#' it is necessary to use the functions `Ed3_L()`, `Ed3_L2()`, `Ed3_R()` and `Ed3_R2()`.
#' Each function describes a specific case:
#' - `Ed3_L()`: calculates the optimal sample size for an additive adjustment factor (i.e. adjust the lower bound of the one-sided confidence interval),
#' however the go-decision is not affected by the bias adjustment
#' - `Ed3_L2()`: calculates the optimal sample size for an additive adjustment factor (i.e. adjust the lower bound of the one-sided confidence interval)
#' when the go-decision is also affected by the bias adjustment
#' - `Ed3_R()`: calculates the optimal sample size for a multiplicative adjustment factor (i.e. use estimate with a retention factor),
#' however the go-decision is not affected by the bias adjustment
#' - `Ed3_R2()`: calculates the optimal sample size for a multiplicative adjustment factor (i.e. use estimate with a retention factor)
#' when the go-decision is also affected by the bias adjustment
#' @param HRgo threshold value for the go/no-go decision rule
#' @param d2 total events for phase II; must be even number
#' @param Adj adjustment parameter
#' @param alpha significance level
#' @param beta `1-beta` power for calculation of sample size for phase III
#' @param w weight for mixture prior distribution
#' @param hr1 first assumed true treatment effect on HR scale for prior distribution
#' @param hr2 second assumed true treatment effect on HR scale for prior distribution
#' @param id1 amount of information for `hr1` in terms of number of events
#' @param id2 amount of information for `hr2` in terms of number of events
#' @param fixed choose if true treatment effects are fixed or random, if TRUE `hr1` is used as fixed effect
#' @return The output of the functions `Ed3_L`, `Ed3_L2`, `Ed3_R` and `Ed3_R2` is the expected number of participants in phase III.
#' @examples res <- Ed3_L(HRgo = 0.8, d2 = 50, Adj = 0.4,
#' alpha = 0.025, beta = 0.1, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- Ed3_L2(HRgo = 0.8, d2 = 50, Adj = 0.4,
#' alpha = 0.025, beta = 0.1, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- Ed3_R(HRgo = 0.8, d2 = 50, Adj = 0.9,
#' alpha = 0.025, beta = 0.1, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- Ed3_R2(HRgo = 0.8, d2 = 50, Adj = 0.9,
#' alpha = 0.025, beta = 0.1, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#'
#' @name Ed3_bias
#' @export
#' @keywords internal
Ed3_L<-function(HRgo, d2, Adj, alpha, beta, w, hr1, hr2, id1, id2, fixed){
if(fixed){
theta <- -log(hr1)
int = try(integrate(function(y){
sapply(y,function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y-qnorm(1-Adj)*sqrt(4/d2))^2))*
dnorm(y,
mean=theta,
sd=sqrt(4/d2))
})
}, -log(HRgo),Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}else{
int = try(integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y-qnorm(1-Adj)*sqrt(4/d2))^2))*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo),Inf)$value
})
}, -Inf, Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}
}
#' Expected probability of a successful program for bias adjustment programs with time-to-event outcomes
#'
#' To discount for overoptimistic results in phase II when calculating the optimal sample size in phase III,
#' it is necessary to use the following functions, which each describe a specific case:
#' - `EPsProg_L()`: calculates the expected probability of a successful for an additive adjustment factor (i.e. adjust the lower bound of the one-sided confidence interval),
#' however the go-decision is not affected by the bias adjustment
#' - `EPsProg_L2()`: calculates the expected probability of a successful for an additive adjustment factor (i.e. adjust the lower bound of the one-sided confidence interval)
#' when the go-decision is also affected by the bias adjustment
#' - `EPsProg_R()`: calculates the expected probability of a successful for a multiplicative adjustment factor (i.e. use estimate with a retention factor),
#' however the go-decision is not affected by the bias adjustment
#' - `EPsProg_R2()`: calculates the expected probability of a successful for a multiplicative adjustment factor (i.e. use estimate with a retention factor)
#' when the go-decision is also affected by the bias adjustment
#' @param HRgo threshold value for the go/no-go decision rule
#' @param d2 total events for phase II; must be even number
#' @param Adj adjustment parameter
#' @param alpha significance level
#' @param beta `1-beta` power for calculation of sample size for phase III
#' @param step1 lower boundary for effect size
#' @param step2 upper boundary for effect size
#' @param w weight for mixture prior distribution
#' @param hr1 first assumed true treatment effect on HR scale for prior distribution
#' @param hr2 second assumed true treatment effect on HR scale for prior distribution
#' @param id1 amount of information for `hr1` in terms of number of events
#' @param id2 amount of information for `hr2` in terms of number of events
#' @param fixed choose if true treatment effects are fixed or random, if TRUE `hr1` is used as fixed effect
#' @return The output of the functions `EPsProg_L()`, `EPsProg_L2()`, `EPsProg_R()` and `EPsProg_R2()` is the expected probability of a successful program.
#' @examples res <- EPsProg_L(HRgo = 0.8, d2 = 50, Adj = 0.4,
#' alpha = 0.025, beta = 0.1,
#' step1 = 1, step2 = 0.95,
#' w = 0.3, hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- EPsProg_L2(HRgo = 0.8, d2 = 50, Adj = 0.4,
#' alpha = 0.025, beta = 0.1,
#' step1 = 1, step2 = 0.95,
#' w = 0.3, hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- EPsProg_R(HRgo = 0.8, d2 = 50, Adj = 0.9,
#' alpha = 0.025, beta = 0.1,
#' step1 = 1, step2 = 0.95,
#' w = 0.3, hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- EPsProg_R2(HRgo = 0.8, d2 = 50, Adj = 0.9,
#' alpha = 0.025, beta = 0.1,
#' step1 = 1, step2 = 0.95,
#' w = 0.3, hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' @name EPsProg_bias
#' @export
#' @keywords internal
EPsProg_L<-function(HRgo, d2, Adj, alpha, beta, step1, step2, w, hr1, hr2, id1, id2, fixed){
c=(qnorm(1-alpha)+qnorm(1-beta))^2
if(fixed){
theta <- -log(hr1)
return(integrate(function(y){
sapply(y,function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=theta/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=theta/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1) )*
dnorm(y,
mean=theta,
sd=sqrt(4/d2))
})
}, -log(HRgo),Inf)$value)
}else{
return(
integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=x/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=x/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1) )*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo),Inf)$value
})
}, -Inf, Inf)$value
)
}
}
#' Utility function for bias adjustment programs with time-to-event outcomes.
#'
#' The utility function calculates the expected utility of our drug development program and is given as gains minus costs and depends on the parameters and the expected probability of a successful program.
#' The utility is in a further step maximized by the `optimal_bias()` function.
#' @param HRgo threshold value for the go/no-go decision rule
#' @param d2 total events for phase II; must be even number
#' @param Adj adjustment parameter
#' @param alpha significance level
#' @param beta `1-beta` power for calculation of sample size for phase III
#' @param xi2 event rate for phase II
#' @param xi3 event rate for phase III
#' @param w weight for mixture prior distribution
#' @param hr1 first assumed true treatment effect on HR scale for prior distribution
#' @param hr2 second assumed true treatment effect on HR scale for prior distribution
#' @param id1 amount of information for `hr1` in terms of number of events
#' @param id2 amount of information for `hr2` in terms of number of events
#' @param c2 variable per-patient cost for phase II
#' @param c3 variable per-patient cost for phase III
#' @param c02 fixed cost for phase II
#' @param c03 fixed cost for phase III
#' @param K constraint on the costs of the program, default: Inf, e.g. no constraint
#' @param N constraint on the total expected sample size of the program, default: Inf, e.g. no constraint
#' @param S constraint on the expected probability of a successful program, default: -Inf, e.g. no constraint
#' @param steps1 lower boundary for effect size category `"small"` in RR scale, default: 1
#' @param stepm1 lower boundary for effect size category `"medium"` in RR scale = upper boundary for effect size category "small" in RR scale, default: 0.95
#' @param stepl1 lower boundary for effect size category `"large"` in RR scale = upper boundary for effect size category "medium" in RR scale, default: 0.85
#' @param b1 expected gain for effect size category `"small"`
#' @param b2 expected gain for effect size category `"medium"`
#' @param b3 expected gain for effect size category `"large"`
#' @param fixed choose if true treatment effects are fixed or random, if TRUE `hr1` is used as fixed effect
#' @return The output of the functions `utility_L()`, `utility_L2()`, `utility_R()` and `utility_R2()` is the expected utility of the program.
#' @examples res <- utility_L(d2 = 50, HRgo = 0.8, Adj = 0.4, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, xi2 = 0.7, xi3 = 0.7,
#' alpha = 0.025, beta = 0.1,
#' c2 = 0.75, c3 = 1, c02 = 100, c03 = 150,
#' K = Inf, N = Inf, S = -Inf,
#' steps1 = 1, stepm1 = 0.95, stepl1 = 0.85,
#' b1 = 1000, b2 = 2000, b3 = 3000,
#' fixed = TRUE)
#' res <- utility_L2(d2 = 50, HRgo = 0.8, Adj = 0.4, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, xi2 = 0.7, xi3 = 0.7,
#' alpha = 0.025, beta = 0.1,
#' c2 = 0.75, c3 = 1, c02 = 100, c03 = 150,
#' K = Inf, N = Inf, S = -Inf,
#' steps1 = 1, stepm1 = 0.95, stepl1 = 0.85,
#' b1 = 1000, b2 = 2000, b3 = 3000,
#' fixed = TRUE)
#' res <- utility_R(d2 = 50, HRgo = 0.8, Adj = 0.9, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, xi2 = 0.7, xi3 = 0.7,
#' alpha = 0.025, beta = 0.1,
#' c2 = 0.75, c3 = 1, c02 = 100, c03 = 150,
#' K = Inf, N = Inf, S = -Inf,
#' steps1 = 1, stepm1 = 0.95, stepl1 = 0.85,
#' b1 = 1000, b2 = 2000, b3 = 3000,
#' fixed = TRUE)
#' res <- utility_R2(d2 = 50, HRgo = 0.8, Adj = 0.9, w = 0.3,
#' hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, xi2 = 0.7, xi3 = 0.7,
#' alpha = 0.025, beta = 0.1,
#' c2 = 0.75, c3 = 1, c02 = 100, c03 = 150,
#' K = Inf, N = Inf, S = -Inf,
#' steps1 = 1, stepm1 = 0.95, stepl1 = 0.85,
#' b1 = 1000, b2 = 2000, b3 = 3000,
#' fixed = TRUE)
#' @name utility_bias
#' @export
#' @keywords internal
utility_L <- function(d2, HRgo, Adj, w, hr1, hr2, id1, id2,
alpha, beta, xi2, xi3,
c2, c3, c02, c03,
K, N, S,
steps1, stepm1, stepl1,
b1, b2, b3,
fixed){
steps2 <- stepm1
stepm2 <- stepl1
stepl2 <- 0
d3 <- Ed3_L(HRgo = HRgo, d2 = d2, Adj = Adj,
alpha = alpha, beta = beta, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
if(is.na(d3)){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# sample size is rounded up to next even natural number
n2 <- ceiling(d2 * (1/xi2))
if(round(n2/2) != n2 / 2) {n2 <- n2 + 1}
n3 <- ceiling(d3 * (1/xi3))
if(round(n3/2) != n3 / 2) {n3 <- n3 + 1}
# expected number of events is rounded to natural number for presentation
d3 <- ceiling(d3)
if(n2+n3>N){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
pg <- Epgo_tte(HRgo = HRgo, d2 = d2, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
K2 <- c02 + c2 * n2 # cost phase II
K3 <- c03 * pg + c3 * n3 # cost phase III
if(K2+K3>K){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# probability of a successful program; small, medium, large effect size
prob1 <- EPsProg_L(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = steps1, step2 = steps2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob2 <- EPsProg_L(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepm1, step2 = stepm2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob3 <- EPsProg_L(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepl1, step2 = stepl2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
SP <- prob1 + prob2 + prob3
if(SP<S){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
G <- b1 * prob1 + b2 * prob2 + b3 * prob3 #gain
EU <- - K2 - K3 + G
return(c(EU, d3, SP, pg, K2, K3, prob1, prob2, prob3, n2, n3))
#output: expected utility Eud, En3, EsP, Epgo, cost phase II and III
}
}
}
}
}
# 1.2. conservative decision rule and sample size calculation:
# use lower bound of one-sided confidence intervall
##############################################################################################
# prior distribution
# as above
#' Expected probability to go to phase III for bias adjustment programs with time-to-event outcomes
#'
#' In the case we do not only want do discount for overoptimistic results in phase II when calculating the sample size in phase III,
#' but also when deciding whether to go to phase III or not the functions `Epgo_L2` and `Epgo_R2` are necessary.
#' The function `Epgo_L2` uses an additive adjustment parameter (i.e. adjust the lower bound of the one-sided confidence interval),
#' the function `Epgo_R2` uses a multiplicative adjustment parameter (i.e. use estimate with a retention factor)
#' @param HRgo threshold value for the go/no-go decision rule
#' @param d2 total number of events for phase II; must be even number
#' @param Adj adjustment parameter
#' @param w weight for mixture prior distribution
#' @param hr1 first assumed true treatment effect on HR scale for prior distribution
#' @param hr2 second assumed true treatment effect on HR scale for prior distribution
#' @param id1 amount of information for `hr1` in terms of number of events
#' @param id2 amount of information for `hr2` in terms of number of events
#' @param fixed choose if true treatment effects are fixed or random, if TRUE `hr1` is used as fixed effect
#' @return The output of the functions `Epgo_L2` and `Epgo_R2` is the expected probability to go to phase III with conservative decision rule and sample size calculation.
#' @examples res <- Epgo_L2(HRgo = 0.8, d2 = 50, Adj = 0.4,
#' w = 0.3, hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' res <- Epgo_R2(HRgo = 0.8, d2 = 50, Adj = 0.9,
#' w = 0.3, hr1 = 0.69, hr2 = 0.81,
#' id1 = 280, id2 = 420, fixed = FALSE)
#' @name Epgo_bias
#' @export
#' @keywords internal
Epgo_L2<-function(HRgo, d2, Adj, w, hr1, hr2, id1, id2, fixed){
if(fixed){
return(
pnorm((-log(hr1)+log(HRgo)-qnorm(1-Adj)*sqrt(4/d2))/sqrt(4/d2))
)
}else(
return(
integrate(function(x){
sapply(x,function(x){
pnorm((x+log(HRgo)-qnorm(1-Adj)*sqrt(4/d2))/sqrt(4/d2))*prior_tte(x,w,hr1,hr2,id1,id2)
})
}, -Inf, Inf)$value
)
)
}
# expected number of events for phase III when going to phase III
#' @rdname Ed3_bias
#' @export
Ed3_L2<-function(HRgo, d2, Adj, alpha, beta, w, hr1, hr2, id1, id2, fixed){
if(fixed){
int = try(
integrate(function(y){
sapply(y,function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y-qnorm(1-Adj)*sqrt(4/d2))^2))*
dnorm(y,
mean=-log(hr1),
sd=sqrt(4/d2))
})
}, -log(HRgo)+qnorm(1-Adj)*sqrt(4/d2),Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}else{
int = try(integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y-qnorm(1-Adj)*sqrt(4/d2))^2))*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo)+qnorm(1-Adj)*sqrt(4/d2),Inf)$value
})
}, -Inf, Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}
}
# expected probability of a successful program
#' @rdname EPsProg_bias
#' @keywords internal
#' @export
EPsProg_L2<-function(HRgo, d2, Adj, alpha, beta, step1, step2, w, hr1, hr2, id1, id2, fixed){
c=(qnorm(1-alpha)+qnorm(1-beta))^2
if(fixed){
return(
integrate(function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=-log(hr1)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=-log(hr1)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1) )*
dnorm(y,
mean=-log(hr1),
sd=sqrt(4/d2))
}, -log(HRgo)+qnorm(1-Adj)*sqrt(4/d2),Inf)$value
)
}else{
return(
integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=x/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
mean=x/(sqrt((y-qnorm(1-Adj)*sqrt(4/d2))^2/c)),
sd=1) )*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo)+qnorm(1-Adj)*sqrt(4/d2),Inf)$value
})
}, -Inf, Inf)$value
)
}
}
# Utility function
#' @rdname utility_bias
#' @keywords internal
#' @export
utility_L2 <- function(d2, HRgo, Adj, w, hr1, hr2, id1, id2,
alpha, beta, xi2, xi3,
c2, c3, c02, c03,
K, N, S,
steps1, stepm1, stepl1,
b1, b2, b3,
fixed){
steps2 <- stepm1
stepm2 <- stepl1
stepl2 <- 0
d3 <- Ed3_L2(HRgo = HRgo, d2 = d2, Adj = Adj,
alpha = alpha, beta = beta, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
if(is.na(d3)){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# sample size is rounded up to next even natural number
n2 <- ceiling(d2 * (1/xi2))
if(round(n2/2) != n2 / 2) {n2 <- n2 + 1}
n3 <- ceiling(d3 * (1/xi3))
if(round(n3/2) != n3 / 2) {n3 <- n3 + 1}
# expected number of events is rounded to natural number for presentation
d3 <- ceiling(d3)
if(n2+n3>N){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
pg <- Epgo_L2(HRgo = HRgo, d2 = d2, Adj = Adj, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
K2 <- c02 + c2 * n2 # cost phase II
K3 <- c03 * pg + c3 * n3 # cost phase III
if(K2+K3>K){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# probability of a successful program; small, medium, large effect size
prob1 <- EPsProg_L2(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = steps1, step2 = steps2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob2 <- EPsProg_L2(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepm1, step2 = stepm2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob3 <- EPsProg_L2(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepl1, step2 = stepl2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
SP <- prob1 + prob2 + prob3
if(SP<S){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
G <- b1 * prob1 + b2 * prob2 + b3 * prob3 #gain
EU <- - K2 - K3 + G
return(c(EU, d3, SP, pg, K2, K3, prob1, prob2, prob3, n2, n3))
#output: expected utility Eud, En3, EsP, Epgo, cost phase II and III
}
}
}
}
}
# 2.1. conservative sample size calculation: use estimate with retetion factor
##############################################################################################
# prior distribution
# as above
# expected probability to go to phase III
# as above
# expected number of events for phase III when going to phase III
#' @rdname Ed3_bias
#' @keywords internal
#' @export
Ed3_R<-function(HRgo, d2, Adj, alpha, beta, w, hr1, hr2, id1, id2, fixed){
if(fixed){
theta <- -log(hr1)
int = try(integrate(function(y){
sapply(y,function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y*Adj)^2))*
dnorm(y,
mean=theta,
sd=sqrt(4/d2))
})
}, -log(HRgo),Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}else{
int = try(integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y*Adj)^2))*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo),Inf)$value
})
}, -Inf, Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}
}
# expected probability of a successful program
#' @rdname EPsProg_bias
#' @keywords internal
#' @export
EPsProg_R<-function(HRgo, d2, Adj, alpha, beta, step1, step2, w, hr1, hr2, id1, id2, fixed){
c=(qnorm(1-alpha)+qnorm(1-beta))^2
if(fixed){
theta <- -log(hr1)
return( integrate(function(y){
sapply(y,function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y*Adj)^2/c)),
mean=theta/(sqrt((y*Adj)^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y*Adj)^2/c)),
mean=theta/(sqrt((y*Adj)^2/c)),
sd=1) )*
dnorm(y,
mean=theta,
sd=sqrt(4/d2))
})
}, -log(HRgo),Inf)$value)
}else{
return(
integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y*Adj)^2/c)),
mean=x/(sqrt((y*Adj)^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y*Adj)^2/c)),
mean=x/(sqrt((y*Adj)^2/c)),
sd=1) )*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo),Inf)$value
})
}, -Inf, Inf)$value
)
}
}
# Utility function
#' @rdname utility_bias
#' @keywords internal
#' @export
utility_R <- function(d2, HRgo, Adj, w, hr1, hr2, id1, id2,
alpha, beta, xi2, xi3,
c2, c3, c02, c03,
K, N, S,
steps1, stepm1, stepl1,
b1, b2, b3,
fixed){
steps2 <- stepm1
stepm2 <- stepl1
stepl2 <- 0
d3 <- Ed3_R(HRgo = HRgo, d2 = d2, Adj = Adj,
alpha = alpha, beta = beta, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
if(is.na(d3)){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# sample size is rounded up to next even natural number
n2 <- ceiling(d2 * (1/xi2))
if(round(n2/2) != n2 / 2) {n2 <- n2 + 1}
n3 <- ceiling(d3 * (1/xi3))
if(round(n3/2) != n3 / 2) {n3 <- n3 + 1}
# expected number of events is rounded to natural number for presentation
d3 <- ceiling(d3)
if(n2+n3>N){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
pg <- Epgo_tte(HRgo = HRgo, d2 = d2, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
K2 <- c02 + c2 * n2 # cost phase II
K3 <- c03 * pg + c3 * n3 # cost phase III
if(K2+K3>K){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# probability of a successful program; small, medium, large effect size
prob1 <- EPsProg_R(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = steps1, step2 = steps2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob2 <- EPsProg_R(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepm1, step2 = stepm2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob3 <- EPsProg_R(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepl1, step2 = stepl2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
SP <- prob1 + prob2 + prob3
if(SP<S){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
G <- b1 * prob1 + b2 * prob2 + b3 * prob3 #gain
EU <- - K2 - K3 + G
return(c(EU, d3, SP, pg, K2, K3, prob1, prob2, prob3, n2, n3))
#output: expected utility Eud, En3, EsP, Epgo, cost phase II and III
}
}
}
}
}
# 2.2. conservative decision rule and sample size calculation:
# use estimate with retetion factor
##############################################################################################
# prior distribution
# as above
# expected probability to go to phase III
#' @rdname Epgo_bias
#' @keywords internal
#' @export
Epgo_R2<-function(HRgo, d2, Adj, w, hr1, hr2, id1, id2, fixed){
if(fixed){
return(
pnorm((-log(hr1)+(log(HRgo)/Adj))/sqrt(4/d2))
)
}else{
return(
integrate(function(x){
sapply(x,function(x){
pnorm((x+(log(HRgo)/Adj))/sqrt(4/d2))*prior_tte(x,w,hr1,hr2,id1,id2)
})
}, -Inf, Inf)$value
)
}
}
# expected number of events for phase III when going to phase III
#' @rdname Ed3_bias
#' @keywords internal
#' @export
Ed3_R2<-function(HRgo, d2, Adj, alpha, beta, w, hr1, hr2, id1, id2, fixed){
if(fixed){
theta <- -log(hr1)
int = try(integrate(function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y*Adj)^2))*
dnorm(y,
mean=theta,
sd=sqrt(4/d2))
}, -log(HRgo)/Adj,Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}else{
int = try(integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( (4*(qnorm(1-alpha)+qnorm(1-beta))^2)/((y*Adj)^2))*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo)/Adj,Inf)$value
})
}, -Inf, Inf),silent=TRUE)
if(inherits(int ,'try-error')){
warning(as.vector(int))
integrated <- NA_real_
} else {
integrated <- int$value
}
return(integrated)
}
}
# expected probability of a successful program
#' @rdname EPsProg_bias
#' @keywords internal
#' @export
EPsProg_R2<-function(HRgo, d2, Adj, alpha, beta, step1, step2, w, hr1, hr2, id1, id2, fixed){
c=(qnorm(1-alpha)+qnorm(1-beta))^2
if(fixed){
return(
integrate(function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y*Adj)^2/c)),
mean=-log(hr1)/(sqrt((y*Adj)^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y*Adj)^2/c)),
mean=-log(hr1)/(sqrt((y*Adj)^2/c)),
sd=1) )*
dnorm(y,
mean=-log(hr1),
sd=sqrt(4/d2))
}, -log(HRgo)/Adj,Inf)$value
)
}else{
return(
integrate(function(x){
sapply(x,function(x){
integrate(function(y){
( pnorm(qnorm(1-alpha)-log(step2)/(sqrt((y*Adj)^2/c)),
mean=x/(sqrt((y*Adj)^2/c)),
sd=1)-
pnorm(qnorm(1-alpha)-log(step1)/(sqrt((y*Adj)^2/c)),
mean=x/(sqrt((y*Adj)^2/c)),
sd=1) )*
dnorm(y,
mean=x,
sd=sqrt(4/d2))*
prior_tte(x,w,hr1,hr2,id1,id2)
}, -log(HRgo)/Adj,Inf)$value
})
}, -Inf, Inf)$value
)
}
}
# Utility function
#' @rdname utility_bias
#' @keywords internal
#' @export
utility_R2 <- function(d2, HRgo, Adj, w, hr1, hr2, id1, id2,
alpha, beta, xi2, xi3,
c2, c3, c02, c03,
K, N, S,
steps1, stepm1, stepl1,
b1, b2, b3,
fixed){
steps2 <- stepm1
stepm2 <- stepl1
stepl2 <- 0
d3 <- Ed3_R2(HRgo = HRgo, d2 = d2, Adj = Adj,
alpha = alpha, beta = beta, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
if(is.na(d3)){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# sample size is rounded up to next even natural number
n2 <- ceiling(d2 * (1/xi2))
if(round(n2/2) != n2 / 2) {n2 <- n2 + 1}
n3 <- ceiling(d3 * (1/xi3))
if(round(n3/2) != n3 / 2) {n3 <- n3 + 1}
# expected number of events is rounded to natural number for presentation
d3 <- ceiling(d3)
if(n2+n3>N){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
pg <- Epgo_R2(HRgo = HRgo, d2 = d2, Adj = Adj, w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
K2 <- c02 + c2 * n2 # cost phase II
K3 <- c03 * pg + c3 * n3 # cost phase III
if(K2+K3>K){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
# probability of a successful program; small, medium, large effect size
prob1 <- EPsProg_R2(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = steps1, step2 = steps2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob2 <- EPsProg_R2(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepm1, step2 = stepm2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
prob3 <- EPsProg_R2(HRgo = HRgo, d2 = d2, Adj = Adj, alpha = alpha, beta = beta,
step1 = stepl1, step2 = stepl2,
w = w, hr1 = hr1, hr2 = hr2, id1 = id1, id2 = id2,
fixed = fixed)
SP <- prob1 + prob2 + prob3
if(SP<S){
return(c(-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999))
}else{
G <- b1 * prob1 + b2 * prob2 + b3 * prob3 #gain
EU <- - K2 - K3 + G
return(c(EU, d3, SP, pg, K2, K3, prob1, prob2, prob3, n2, n3))
#output: expected utility Eud, En3, EsP, Epgo, cost phase II and III
}
}
}
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.