Nothing
R/utility_functions.R
In gsearly: Creates Group Sequential Trial Designs when Early Outcomes are Available
Defines functions
.vbetaUser
.vbetaUnif
.vbetaExp
.vbeta0
.vbeta
.userinformData
.tau0
.tau
tfuStandard
.rvarUnif
.rvarExp
.recruitData
.informData
.getRdata
corrUnif
corrExp
.bisection
Documented in
corrExp
corrUnif
tfuStandard
#-------------------------------------------------------------------------------
# Utility functions for gsearly
#-------------------------------------------------------------------------------
# 1st July 2026
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 2. .bisection
#-------------------------------------------------------------------------------
.bisection <- function(ffn, pow, largs, a, b, niter = 1000, tol = 1e-07) {
## Bisection algorithm to find target power Stop if the signs of
## the function differ
if (!(ffn(a, pow, largs = largs) < 0) && (ffn(b, pow, largs = largs) >
0)) {
stop("signs of f(a) and f(b) differ")
} else if ((ffn(a, pow, largs = largs) > 0) && (ffn(b, pow, largs = largs) <
0)) {
stop("signs of f(a) and f(b) differ")
}
for (i in 1:niter) {
## Midpoint c
c <- (a + b)/2
## Stop and return the root if criterion met
if ((ffn(c, pow, largs = largs) == 0) || ((b - a)/2) < tol) {
return(c)
}
## Start another iteration Check signs of the functions
ifelse(sign(ffn(c, pow, largs = largs)) == sign(ffn(a, pow, largs = largs)),
a <- c, b <- c)
}
## Max iterations
stop("Max iterations with no root: try increasing sopt$bisect$niter")
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 3. corrExp
#-------------------------------------------------------------------------------
corrExp <- function(rho = 0, tfu) {
## Construct exponential correlation matrix
if (rho >= 0 & rho < 1) {
stfu <- tfu$stfu
s <- length(stfu)
corrmat <- matrix(0, nrow = s, ncol = s)
for (i in 1:s) {
for (j in 1:s) {
diff <- abs(stfu[i] - stfu[j])
corrmat[i, j] <- rho^(diff)
}
}
} else {
corrmat <- diag(rep(1, s))
}
return(as.matrix(corrmat))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 4. corrUnif
#-------------------------------------------------------------------------------
corrUnif <- function(rho = 0, tfu) {
## Construct uniform correlation matrix
if (rho >= 0 & rho < 1) {
stfu <- as.numeric(tfu[["stfu"]])
s <- length(stfu)
corrmat <- matrix(c(rho), nrow = s, ncol = s)
diag(corrmat) <- rep(1, s)
} else {
corrmat <- diag(rep(1, s))
}
return(as.matrix(corrmat))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 10. .getRdata
#-------------------------------------------------------------------------------
.getRdata <- function(tt = NULL, tfu, n, trecruit, rmodel, tlag = NULL,
m) {
## Longitudinal data
longN <- function(n, trecruit, tt, tr, m, rmodel = "fix") {
mod <- .selectRmodel(rmodel = rmodel)
longitn <- mod$fk(n = n, trecruit = trecruit, m = m) * sapply(tt,
mod$fg, tr = tr, trecruit = trecruit, m = m)
return(longitn)
}
## Get recruitment data
alltfu <- c(0, as.numeric(tfu$tfu))
ss <- length(alltfu)
if (is.null(tt) == TRUE) {
tfuData <- function(tfu, trecruit, tlag) {
allfu <- c(0, as.numeric(tfu[["tfu"]]))
s0 <- length(allfu)
return(list(s = s0 - 1, data = seq(-tlag[1], allfu[s0] + trecruit +
tlag[2], 1)))
}
tt <- tfuData(tfu = tfu, trecruit = trecruit, tlag = tlag)[["data"]]
rdat <- lapply(alltfu, longN, n = n, trecruit = trecruit, tt = tt,
rmodel = rmodel, m = m)
rdat <- list(t = tt, n = rdat)
names(rdat[["n"]]) <- paste("n", alltfu, sep = "")
} else {
rdat <- sapply(alltfu, longN, n = n, trecruit = trecruit, tt = tt,
rmodel = rmodel, m = m)
if (length(tt) == 1) {
rdat <- matrix(rdat, nrow = 1, ncol = ss, byrow = TRUE)
}
colnames(rdat) <- as.character(alltfu)
rownames(rdat) <- as.character(tt)
}
return(rdat)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 16. .informData
#-------------------------------------------------------------------------------
.informData <- function(rdata, cmodel, sd, rho, intonly = FALSE) {
## Construct full information data
s <- rdata$s
vphi <- rdata$vphi
tint <- c(rdata$tinterims, rdata$trecruit + rdata$tfu$tfu[s])
iint <- sapply(tint, .tau, cmodel = cmodel, tfu = rdata$tfu, trecruit = rdata$trecruit,
sd = sd, rho = rho, vphi = vphi, n = rdata$n, rmodel = rdata$rmodel,
m = rdata$m)
varbint <- sapply(tint, .vbeta, cmodel = cmodel, tfu = rdata$tfu, n = rdata$n,
trecruit = rdata$trecruit, sd = sd, rho = rho, vphi = vphi, rmodel = rdata$rmodel,
m = rdata$m)
allinterims <- matrix(c(iint, varbint), nrow = 2, ncol = length(rdata$tinterims) +
1, byrow = TRUE)
colnames(allinterims) <- as.character(tint)
rownames(allinterims) <- c("tau", "vbeta")
if (intonly == FALSE) {
longt <- rdata$data$t
idata <- sapply(longt, .tau, cmodel = cmodel, tfu = rdata$tfu, trecruit = rdata$trecruit,
sd = sd, rho = rho, vphi = vphi, n = rdata$n, rmodel = rdata$rmodel,
m = rdata$m)
varb <- sapply(longt, .vbeta, cmodel = cmodel, tfu = rdata$tfu,
n = rdata$n, trecruit = rdata$trecruit, sd = sd, rho = rho,
vphi = vphi, rmodel = rdata$rmodel, m = rdata$m)
idata <- list(t = longt, inform = idata, var = varb)
}
if (intonly == FALSE) {
out <- list(cmodel = cmodel, sd = sd, data = idata, tlooks = tint,
inform = allinterims)
} else {
out <- list(cmodel = cmodel, sd = sd, tlooks = tint, inform = allinterims)
}
return(out)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 25. .recruitData
#-------------------------------------------------------------------------------
.recruitData <- function(rmodel, trecruit, s, tfu, tinterims, n, vphi,
contrat, tlag = NULL, m, intonly = FALSE) {
## Construct study full recruitment data
getinterims <- .getRdata(tt = tinterims, tfu = tfu, n = n, trecruit = trecruit,
rmodel = rmodel, tlag = tlag, m = m)
if (intonly == FALSE) {
rdata <- .getRdata(tfu = tfu, n = n, trecruit = trecruit, rmodel = rmodel,
tlag = tlag, m = m)
}
if (intonly == FALSE) {
out <- list(rmodel = rmodel, m = m, trecruit = trecruit, s = s, tfu = tfu,
n = n, vphi = vphi, contrat = contrat, tinterims = tinterims,
intnumbers = getinterims, tlag = tlag, data = rdata)
} else {
out <- list(rmodel = rmodel, m = m, trecruit = trecruit, s = s, tfu = tfu,
n = n, vphi = vphi, contrat = contrat, tinterims = tinterims,
intnumbers = getinterims)
}
return(out)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 27. .rvarExp
#-------------------------------------------------------------------------------
.rvarExp <- function(tfu, tt, alpha, trecruit, rmodel = "fix", m) {
## Relative efficiency of treatment effect estimate for exponential
## model
mod <- .selectRmodel(rmodel = rmodel)
stfu <- tfu$stfu
tfu <- tfu$tfu
s <- length(tfu)
ts <- tfu[s]
if (tt <= ts) {
rvar <- 0
} else if (tt >= as.numeric(trecruit + ts)) {
rvar <- 1
} else {
diffmat <- matrix(stfu, nrow = s, ncol = s, byrow = FALSE) - matrix(stfu,
nrow = s, ncol = s, byrow = TRUE)
term_s <- (1 - alpha^(2 * abs(diffmat[s, s - 1])))
if (s > 2) {
mid_theta <- 0
for (ms in 1:(s - 2)) {
aterm <- alpha^(2 * abs(diffmat[s, s - ms])) * ((1 - alpha^(2 *
abs(diffmat[s - ms, s - ms - 1]))))
fg_rat <- mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit,
m = m)/mod$fg(tt = tt, tr = tfu[s - ms], trecruit = trecruit,
m = m)
mid_theta[ms] <- fg_rat * aterm
}
}
aterm <- (alpha^(2 * abs(diffmat[s, 1])))
fg_rat <- mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit, m = m)/mod$fg(tt = tt,
tr = tfu[1], trecruit = trecruit, m = m)
term_1 <- fg_rat * aterm
if (s > 2) {
rvar <- term_s + sum(mid_theta) + term_1
} else {
rvar <- term_s + term_1
}
}
return(as.numeric(rvar))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 28. .rvarUnif
#-------------------------------------------------------------------------------
.rvarUnif <- function(tfu, tt, alpha, trecruit, rmodel = "fix", m) {
## Relative efficiency of treatment effect estimate for uniform
## model
mod <- .selectRmodel(rmodel = rmodel)
stfu <- tfu$stfu
tfu <- tfu$tfu
s <- length(tfu)
ts <- tfu[s]
if (tt <= ts) {
rvar <- 0
} else if (tt >= as.numeric(trecruit + ts)) {
rvar <- 1
} else {
mid_theta <- 0
for (ms in 1:(s - 1)) {
aterm <- ((1 - alpha) * (1 + ms * alpha))/(1 + (ms - 1) * alpha)
fg_rat1 <- mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit,
m = m)/mod$fg(tt = tt, tr = tfu[ms + 1], trecruit = trecruit,
m = m)
fg_rat2 <- mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit,
m = m)/mod$fg(tt = tt, tr = tfu[ms], trecruit = trecruit,
m = m)
mid_theta[ms] <- aterm * (fg_rat1 - fg_rat2)
}
term_1 <- mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit, m = m)/mod$fg(tt = tt,
tr = tfu[1], trecruit = trecruit, m = m)
rvar <- term_1 + sum(mid_theta)
}
return(as.numeric(rvar))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 33. tfuStandard
#-------------------------------------------------------------------------------
tfuStandard <- function(tfu, tref = c(1, 2)) {
## Standardise follow-up times using reference categories tref
tfu <- as.integer(tfu)
s <- length(tfu)
if (all(tfu >= 1) & all(diff(tfu) > 0)) {
if (is.element(tref[1], 1:s) == TRUE & is.element(tref[2], 1:s) ==
TRUE) {
tspan <- abs(tfu[tref[2]] - tfu[tref[1]])
stfu <- (tfu/tspan)
stfu <- stfu - stfu[1] + 1
} else {
tspan <- abs(tfu[2] - tfu[1])
stfu <- (tfu/tspan)
stfu <- stfu - stfu[1] + 1
}
} else {
tfu <- stfu <- 1:s
}
return(list(tfu = as.integer(tfu), stfu = as.numeric(stfu)))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 35. .tau
#-------------------------------------------------------------------------------
.tau <- function(tfu, tt, trecruit, sd, rho = 0.5, cmodel = "uniform", vphi = 0.5,
n, rmodel = "fix", m) {
## Information fraction for model cmodel
s <- length(tfu$tfu)
etau0 <- .tau0(tfu = tfu, tt = tt, trecruit = trecruit, rmodel = rmodel,
m = m)
if (tt <= tfu$tfu[s]) {
finfo <- 0
} else {
if (cmodel == "exponential") {
finfo <- etau0/.rvarExp(tfu = tfu, tt = tt, alpha = rho, trecruit = trecruit,
rmodel = rmodel, m = m)
} else if (cmodel == "uniform") {
finfo <- etau0/.rvarUnif(tfu = tfu, tt = tt, alpha = rho, trecruit = trecruit,
rmodel = rmodel, m = m)
}
}
return(finfo)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 36. .tau0
#-------------------------------------------------------------------------------
.tau0 <- function(tfu, tt, trecruit, rmodel = "fix", m) {
## Information fraction when correlation is zero
mod <- .selectRmodel(rmodel = rmodel)
tfu <- tfu$tfu
s <- length(tfu)
if (tt > (trecruit + tfu[s])) {
tt <- trecruit + tfu[s]
}
tau <- mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit, m = m)/mod$fg(t = tfu[s] +
trecruit, tr = tfu[s], trecruit = trecruit, m = m)
return(as.numeric(tau))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 37. .userinformData
#-------------------------------------------------------------------------------
.userinformData <- function(x, largs) {
## Calculate variance and information for user model
if (is.matrix(largs$ninterims) == TRUE) {
all_n <- rbind(largs$ninterims, matrix(rep(x, largs$s + 1), nrow = 1))
} else if (is.function(largs$ninterims) == TRUE) {
all_n <- rbind(largs$ninterim(x), matrix(rep(x, largs$s + 1), nrow = 1))
}
get_varb <- sapply(1:largs$nlooks, .vbetaUser, s = largs$s, vphi = largs$vphi,
alln = all_n, vcovmat = largs$vcovmat)
iint <- get_varb[largs$nlooks]/get_varb
sumint <- rbind(iint, get_varb)
alltfu <- c(0, as.numeric(largs$tfu$tfu))
rownames(sumint) <- c("tau", "vbeta")
colnames(sumint) <- paste("t", as.character(largs$tint), sep = "")
return(sumint)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 44. .vbeta
#-------------------------------------------------------------------------------
.vbeta <- function(tfu, tt, n, trecruit, sd, rho = 0.5, cmodel = "uniform",
vphi = 0.5, rmodel = "fix", m) {
## Variance of treatment effect estimate for model cmodel
s <- length(tfu$tfu)
if (tt <= tfu$tfu[s]) {
var_beta <- NA
} else {
if (cmodel == "exponential") {
var_beta <- .vbetaExp(tfu = tfu, tt = tt, alpha = rho, n = n,
sd = sd, trecruit = trecruit, vphi = vphi, rmodel = rmodel,
m = m)
} else if (cmodel == "uniform") {
var_beta <- .vbetaUnif(tfu = tfu, tt = tt, alpha = rho, n = n,
sd = sd, trecruit = trecruit, vphi = vphi, rmodel = rmodel,
m = m)
}
}
return(var_beta)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 45. .vbeta0
#-------------------------------------------------------------------------------
.vbeta0 <- function(tfu, tt, n, trecruit, sd, vphi, rmodel = "fix", m) {
## Variance of treatment effect estimate when correlation is zero
mod <- .selectRmodel(rmodel = rmodel)
tfu <- tfu$tfu
s <- length(tfu)
if (tt <= tfu[s]) {
var_beta <- NA
} else {
kfg <- mod$fk(n = n, trecruit = trecruit, m = m) * mod$fg(tt = tt,
tr = tfu[s], trecruit = trecruit, m = m)
var_beta <- (sd^2)/(kfg * vphi * (1 - vphi))
}
return(as.numeric(var_beta))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 46. .vbetaExp
#-------------------------------------------------------------------------------
.vbetaExp <- function(tfu, tt, sd, alpha, n, trecruit, vphi, rmodel = "fix",
m) {
## Variance of treatment effect estimate for exponential model
mod <- .selectRmodel(rmodel = rmodel)
stfu <- tfu$stfu
tfu <- tfu$tfu
s <- length(tfu)
if (alpha == 0) {
var_beta <- (sd^2)/(mod$fk(n = n, trecruit = trecruit, m = m) *
vphi * (1 - vphi) * mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit,
m = m))
} else {
diffmat <- matrix(stfu, nrow = s, ncol = s, byrow = FALSE) - matrix(stfu,
nrow = s, ncol = s, byrow = TRUE)
# term s
term_s <- (1 - alpha^(2 * abs(diffmat[s, s - 1])))/mod$fg(tt = tt,
tr = tfu[s], trecruit = trecruit, m = m)
# mid terms
if (s > 2) {
mid_theta <- 0
for (ms in 1:(s - 2)) {
mid_term <- alpha^(2 * abs(diffmat[s, s - ms]))
mid_theta[ms] <- ((1 - alpha^(2 * abs(diffmat[s - ms, s -
ms - 1])))) * mid_term/mod$fg(tt = tt, tr = tfu[s - ms],
trecruit = trecruit, m = m)
}
}
# term 1
term_1 <- (alpha^(2 * abs(diffmat[s, 1])))/mod$fg(tt = tt, tr = tfu[1],
trecruit = trecruit, m = m)
# var
kconst <- (sd^2)/(mod$fk(n = n, trecruit = trecruit, m = m) * vphi *
(1 - vphi))
if (s > 2) {
var_beta <- kconst * (term_s + sum(mid_theta) + term_1)
} else {
var_beta <- kconst * (term_s + term_1)
}
}
return(as.numeric(var_beta))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 47. .vbetaUnif
#-------------------------------------------------------------------------------
.vbetaUnif <- function(tfu, tt, sd, alpha, n, trecruit, vphi, rmodel = "fix",
m) {
## Variance of treatment effect estimate for uniform model
mod <- .selectRmodel(rmodel = rmodel)
stfu <- tfu$stfu
tfu <- tfu$tfu
s <- length(tfu)
if (alpha == 0) {
var_beta <- (sd^2)/(mod$fk(n = n, trecruit = trecruit, m = m) *
vphi * (1 - vphi) * mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit,
m = m))
} else {
sum_term <- 0
for (ms in 1:(s - 1)) {
term_m1 <- 1/mod$fg(tt = tt, tr = tfu[ms + 1], trecruit = trecruit,
m = m)
term_m <- 1/mod$fg(tt = tt, tr = tfu[ms], trecruit = trecruit,
m = m)
sum_term[ms] <- ((1 - alpha) * (1 + ms * alpha)/(1 + (ms -
1) * alpha)) * (term_m1 - term_m)
}
term_1 <- 1/mod$fg(tt = tt, tr = tfu[1], trecruit = trecruit, m = m)
var_beta <- ((sd^2)/(mod$fk(n = n, trecruit = trecruit, m = m) *
vphi * (1 - vphi))) * (term_1 + sum(sum_term))
}
return(as.numeric(var_beta))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 48. .vbetaUser
#-------------------------------------------------------------------------------
.vbetaUser <- function(x, s, vphi, alln, vcovmat) {
## Get variance for each interim analyses and study end
all_nx <- alln[x, 2:(s + 1)]
n_0 <- vphi * all_nx
n_1 <- (1 - vphi) * all_nx
## Calculate variance for known sample sizes and covariance matrix
n_0 <- c(n_0, 0)
n_1 <- c(n_1, 0)
V0 <- V1 <- matrix(0, nrow = s, ncol = s)
for (r in 1:s) {
for (c in 1:s) {
mx <- max(r, c)
for (mx in max(r, c):s) {
iest_vcov <- solve(vcovmat[1:mx, 1:mx])
V1[r, c] <- V1[r, c] + (n_1[mx] - n_1[mx + 1]) * iest_vcov[r,
c]
V0[r, c] <- V0[r, c] + (n_0[mx] - n_0[mx + 1]) * iest_vcov[r,
c]
}
}
}
mvcov <- solve(V0) + solve(V1)
var_beta <- mvcov[s, s]
return(var_beta)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# end
#-------------------------------------------------------------------------------
Try the gsearly package in your browser
Any scripts or data that you put into this service are public.
gsearly documentation built on July 10, 2026, 5:09 p.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 .vbetaUser .vbetaUnif .vbetaExp .vbeta0 .vbeta .userinformData .tau0 .tau tfuStandard .rvarUnif .rvarExp .recruitData .informData .getRdata corrUnif corrExp .bisection
Documented in corrExp corrUnif tfuStandard
#-------------------------------------------------------------------------------
# Utility functions for gsearly
#-------------------------------------------------------------------------------
# 1st July 2026
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 2. .bisection
#-------------------------------------------------------------------------------
.bisection <- function(ffn, pow, largs, a, b, niter = 1000, tol = 1e-07) {
## Bisection algorithm to find target power Stop if the signs of
## the function differ
if (!(ffn(a, pow, largs = largs) < 0) && (ffn(b, pow, largs = largs) >
0)) {
stop("signs of f(a) and f(b) differ")
} else if ((ffn(a, pow, largs = largs) > 0) && (ffn(b, pow, largs = largs) <
0)) {
stop("signs of f(a) and f(b) differ")
}
for (i in 1:niter) {
## Midpoint c
c <- (a + b)/2
## Stop and return the root if criterion met
if ((ffn(c, pow, largs = largs) == 0) || ((b - a)/2) < tol) {
return(c)
}
## Start another iteration Check signs of the functions
ifelse(sign(ffn(c, pow, largs = largs)) == sign(ffn(a, pow, largs = largs)),
a <- c, b <- c)
}
## Max iterations
stop("Max iterations with no root: try increasing sopt$bisect$niter")
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 3. corrExp
#-------------------------------------------------------------------------------
corrExp <- function(rho = 0, tfu) {
## Construct exponential correlation matrix
if (rho >= 0 & rho < 1) {
stfu <- tfu$stfu
s <- length(stfu)
corrmat <- matrix(0, nrow = s, ncol = s)
for (i in 1:s) {
for (j in 1:s) {
diff <- abs(stfu[i] - stfu[j])
corrmat[i, j] <- rho^(diff)
}
}
} else {
corrmat <- diag(rep(1, s))
}
return(as.matrix(corrmat))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 4. corrUnif
#-------------------------------------------------------------------------------
corrUnif <- function(rho = 0, tfu) {
## Construct uniform correlation matrix
if (rho >= 0 & rho < 1) {
stfu <- as.numeric(tfu[["stfu"]])
s <- length(stfu)
corrmat <- matrix(c(rho), nrow = s, ncol = s)
diag(corrmat) <- rep(1, s)
} else {
corrmat <- diag(rep(1, s))
}
return(as.matrix(corrmat))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 10. .getRdata
#-------------------------------------------------------------------------------
.getRdata <- function(tt = NULL, tfu, n, trecruit, rmodel, tlag = NULL,
m) {
## Longitudinal data
longN <- function(n, trecruit, tt, tr, m, rmodel = "fix") {
mod <- .selectRmodel(rmodel = rmodel)
longitn <- mod$fk(n = n, trecruit = trecruit, m = m) * sapply(tt,
mod$fg, tr = tr, trecruit = trecruit, m = m)
return(longitn)
}
## Get recruitment data
alltfu <- c(0, as.numeric(tfu$tfu))
ss <- length(alltfu)
if (is.null(tt) == TRUE) {
tfuData <- function(tfu, trecruit, tlag) {
allfu <- c(0, as.numeric(tfu[["tfu"]]))
s0 <- length(allfu)
return(list(s = s0 - 1, data = seq(-tlag[1], allfu[s0] + trecruit +
tlag[2], 1)))
}
tt <- tfuData(tfu = tfu, trecruit = trecruit, tlag = tlag)[["data"]]
rdat <- lapply(alltfu, longN, n = n, trecruit = trecruit, tt = tt,
rmodel = rmodel, m = m)
rdat <- list(t = tt, n = rdat)
names(rdat[["n"]]) <- paste("n", alltfu, sep = "")
} else {
rdat <- sapply(alltfu, longN, n = n, trecruit = trecruit, tt = tt,
rmodel = rmodel, m = m)
if (length(tt) == 1) {
rdat <- matrix(rdat, nrow = 1, ncol = ss, byrow = TRUE)
}
colnames(rdat) <- as.character(alltfu)
rownames(rdat) <- as.character(tt)
}
return(rdat)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 16. .informData
#-------------------------------------------------------------------------------
.informData <- function(rdata, cmodel, sd, rho, intonly = FALSE) {
## Construct full information data
s <- rdata$s
vphi <- rdata$vphi
tint <- c(rdata$tinterims, rdata$trecruit + rdata$tfu$tfu[s])
iint <- sapply(tint, .tau, cmodel = cmodel, tfu = rdata$tfu, trecruit = rdata$trecruit,
sd = sd, rho = rho, vphi = vphi, n = rdata$n, rmodel = rdata$rmodel,
m = rdata$m)
varbint <- sapply(tint, .vbeta, cmodel = cmodel, tfu = rdata$tfu, n = rdata$n,
trecruit = rdata$trecruit, sd = sd, rho = rho, vphi = vphi, rmodel = rdata$rmodel,
m = rdata$m)
allinterims <- matrix(c(iint, varbint), nrow = 2, ncol = length(rdata$tinterims) +
1, byrow = TRUE)
colnames(allinterims) <- as.character(tint)
rownames(allinterims) <- c("tau", "vbeta")
if (intonly == FALSE) {
longt <- rdata$data$t
idata <- sapply(longt, .tau, cmodel = cmodel, tfu = rdata$tfu, trecruit = rdata$trecruit,
sd = sd, rho = rho, vphi = vphi, n = rdata$n, rmodel = rdata$rmodel,
m = rdata$m)
varb <- sapply(longt, .vbeta, cmodel = cmodel, tfu = rdata$tfu,
n = rdata$n, trecruit = rdata$trecruit, sd = sd, rho = rho,
vphi = vphi, rmodel = rdata$rmodel, m = rdata$m)
idata <- list(t = longt, inform = idata, var = varb)
}
if (intonly == FALSE) {
out <- list(cmodel = cmodel, sd = sd, data = idata, tlooks = tint,
inform = allinterims)
} else {
out <- list(cmodel = cmodel, sd = sd, tlooks = tint, inform = allinterims)
}
return(out)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 25. .recruitData
#-------------------------------------------------------------------------------
.recruitData <- function(rmodel, trecruit, s, tfu, tinterims, n, vphi,
contrat, tlag = NULL, m, intonly = FALSE) {
## Construct study full recruitment data
getinterims <- .getRdata(tt = tinterims, tfu = tfu, n = n, trecruit = trecruit,
rmodel = rmodel, tlag = tlag, m = m)
if (intonly == FALSE) {
rdata <- .getRdata(tfu = tfu, n = n, trecruit = trecruit, rmodel = rmodel,
tlag = tlag, m = m)
}
if (intonly == FALSE) {
out <- list(rmodel = rmodel, m = m, trecruit = trecruit, s = s, tfu = tfu,
n = n, vphi = vphi, contrat = contrat, tinterims = tinterims,
intnumbers = getinterims, tlag = tlag, data = rdata)
} else {
out <- list(rmodel = rmodel, m = m, trecruit = trecruit, s = s, tfu = tfu,
n = n, vphi = vphi, contrat = contrat, tinterims = tinterims,
intnumbers = getinterims)
}
return(out)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 27. .rvarExp
#-------------------------------------------------------------------------------
.rvarExp <- function(tfu, tt, alpha, trecruit, rmodel = "fix", m) {
## Relative efficiency of treatment effect estimate for exponential
## model
mod <- .selectRmodel(rmodel = rmodel)
stfu <- tfu$stfu
tfu <- tfu$tfu
s <- length(tfu)
ts <- tfu[s]
if (tt <= ts) {
rvar <- 0
} else if (tt >= as.numeric(trecruit + ts)) {
rvar <- 1
} else {
diffmat <- matrix(stfu, nrow = s, ncol = s, byrow = FALSE) - matrix(stfu,
nrow = s, ncol = s, byrow = TRUE)
term_s <- (1 - alpha^(2 * abs(diffmat[s, s - 1])))
if (s > 2) {
mid_theta <- 0
for (ms in 1:(s - 2)) {
aterm <- alpha^(2 * abs(diffmat[s, s - ms])) * ((1 - alpha^(2 *
abs(diffmat[s - ms, s - ms - 1]))))
fg_rat <- mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit,
m = m)/mod$fg(tt = tt, tr = tfu[s - ms], trecruit = trecruit,
m = m)
mid_theta[ms] <- fg_rat * aterm
}
}
aterm <- (alpha^(2 * abs(diffmat[s, 1])))
fg_rat <- mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit, m = m)/mod$fg(tt = tt,
tr = tfu[1], trecruit = trecruit, m = m)
term_1 <- fg_rat * aterm
if (s > 2) {
rvar <- term_s + sum(mid_theta) + term_1
} else {
rvar <- term_s + term_1
}
}
return(as.numeric(rvar))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 28. .rvarUnif
#-------------------------------------------------------------------------------
.rvarUnif <- function(tfu, tt, alpha, trecruit, rmodel = "fix", m) {
## Relative efficiency of treatment effect estimate for uniform
## model
mod <- .selectRmodel(rmodel = rmodel)
stfu <- tfu$stfu
tfu <- tfu$tfu
s <- length(tfu)
ts <- tfu[s]
if (tt <= ts) {
rvar <- 0
} else if (tt >= as.numeric(trecruit + ts)) {
rvar <- 1
} else {
mid_theta <- 0
for (ms in 1:(s - 1)) {
aterm <- ((1 - alpha) * (1 + ms * alpha))/(1 + (ms - 1) * alpha)
fg_rat1 <- mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit,
m = m)/mod$fg(tt = tt, tr = tfu[ms + 1], trecruit = trecruit,
m = m)
fg_rat2 <- mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit,
m = m)/mod$fg(tt = tt, tr = tfu[ms], trecruit = trecruit,
m = m)
mid_theta[ms] <- aterm * (fg_rat1 - fg_rat2)
}
term_1 <- mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit, m = m)/mod$fg(tt = tt,
tr = tfu[1], trecruit = trecruit, m = m)
rvar <- term_1 + sum(mid_theta)
}
return(as.numeric(rvar))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 33. tfuStandard
#-------------------------------------------------------------------------------
tfuStandard <- function(tfu, tref = c(1, 2)) {
## Standardise follow-up times using reference categories tref
tfu <- as.integer(tfu)
s <- length(tfu)
if (all(tfu >= 1) & all(diff(tfu) > 0)) {
if (is.element(tref[1], 1:s) == TRUE & is.element(tref[2], 1:s) ==
TRUE) {
tspan <- abs(tfu[tref[2]] - tfu[tref[1]])
stfu <- (tfu/tspan)
stfu <- stfu - stfu[1] + 1
} else {
tspan <- abs(tfu[2] - tfu[1])
stfu <- (tfu/tspan)
stfu <- stfu - stfu[1] + 1
}
} else {
tfu <- stfu <- 1:s
}
return(list(tfu = as.integer(tfu), stfu = as.numeric(stfu)))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 35. .tau
#-------------------------------------------------------------------------------
.tau <- function(tfu, tt, trecruit, sd, rho = 0.5, cmodel = "uniform", vphi = 0.5,
n, rmodel = "fix", m) {
## Information fraction for model cmodel
s <- length(tfu$tfu)
etau0 <- .tau0(tfu = tfu, tt = tt, trecruit = trecruit, rmodel = rmodel,
m = m)
if (tt <= tfu$tfu[s]) {
finfo <- 0
} else {
if (cmodel == "exponential") {
finfo <- etau0/.rvarExp(tfu = tfu, tt = tt, alpha = rho, trecruit = trecruit,
rmodel = rmodel, m = m)
} else if (cmodel == "uniform") {
finfo <- etau0/.rvarUnif(tfu = tfu, tt = tt, alpha = rho, trecruit = trecruit,
rmodel = rmodel, m = m)
}
}
return(finfo)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 36. .tau0
#-------------------------------------------------------------------------------
.tau0 <- function(tfu, tt, trecruit, rmodel = "fix", m) {
## Information fraction when correlation is zero
mod <- .selectRmodel(rmodel = rmodel)
tfu <- tfu$tfu
s <- length(tfu)
if (tt > (trecruit + tfu[s])) {
tt <- trecruit + tfu[s]
}
tau <- mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit, m = m)/mod$fg(t = tfu[s] +
trecruit, tr = tfu[s], trecruit = trecruit, m = m)
return(as.numeric(tau))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 37. .userinformData
#-------------------------------------------------------------------------------
.userinformData <- function(x, largs) {
## Calculate variance and information for user model
if (is.matrix(largs$ninterims) == TRUE) {
all_n <- rbind(largs$ninterims, matrix(rep(x, largs$s + 1), nrow = 1))
} else if (is.function(largs$ninterims) == TRUE) {
all_n <- rbind(largs$ninterim(x), matrix(rep(x, largs$s + 1), nrow = 1))
}
get_varb <- sapply(1:largs$nlooks, .vbetaUser, s = largs$s, vphi = largs$vphi,
alln = all_n, vcovmat = largs$vcovmat)
iint <- get_varb[largs$nlooks]/get_varb
sumint <- rbind(iint, get_varb)
alltfu <- c(0, as.numeric(largs$tfu$tfu))
rownames(sumint) <- c("tau", "vbeta")
colnames(sumint) <- paste("t", as.character(largs$tint), sep = "")
return(sumint)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 44. .vbeta
#-------------------------------------------------------------------------------
.vbeta <- function(tfu, tt, n, trecruit, sd, rho = 0.5, cmodel = "uniform",
vphi = 0.5, rmodel = "fix", m) {
## Variance of treatment effect estimate for model cmodel
s <- length(tfu$tfu)
if (tt <= tfu$tfu[s]) {
var_beta <- NA
} else {
if (cmodel == "exponential") {
var_beta <- .vbetaExp(tfu = tfu, tt = tt, alpha = rho, n = n,
sd = sd, trecruit = trecruit, vphi = vphi, rmodel = rmodel,
m = m)
} else if (cmodel == "uniform") {
var_beta <- .vbetaUnif(tfu = tfu, tt = tt, alpha = rho, n = n,
sd = sd, trecruit = trecruit, vphi = vphi, rmodel = rmodel,
m = m)
}
}
return(var_beta)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 45. .vbeta0
#-------------------------------------------------------------------------------
.vbeta0 <- function(tfu, tt, n, trecruit, sd, vphi, rmodel = "fix", m) {
## Variance of treatment effect estimate when correlation is zero
mod <- .selectRmodel(rmodel = rmodel)
tfu <- tfu$tfu
s <- length(tfu)
if (tt <= tfu[s]) {
var_beta <- NA
} else {
kfg <- mod$fk(n = n, trecruit = trecruit, m = m) * mod$fg(tt = tt,
tr = tfu[s], trecruit = trecruit, m = m)
var_beta <- (sd^2)/(kfg * vphi * (1 - vphi))
}
return(as.numeric(var_beta))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 46. .vbetaExp
#-------------------------------------------------------------------------------
.vbetaExp <- function(tfu, tt, sd, alpha, n, trecruit, vphi, rmodel = "fix",
m) {
## Variance of treatment effect estimate for exponential model
mod <- .selectRmodel(rmodel = rmodel)
stfu <- tfu$stfu
tfu <- tfu$tfu
s <- length(tfu)
if (alpha == 0) {
var_beta <- (sd^2)/(mod$fk(n = n, trecruit = trecruit, m = m) *
vphi * (1 - vphi) * mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit,
m = m))
} else {
diffmat <- matrix(stfu, nrow = s, ncol = s, byrow = FALSE) - matrix(stfu,
nrow = s, ncol = s, byrow = TRUE)
# term s
term_s <- (1 - alpha^(2 * abs(diffmat[s, s - 1])))/mod$fg(tt = tt,
tr = tfu[s], trecruit = trecruit, m = m)
# mid terms
if (s > 2) {
mid_theta <- 0
for (ms in 1:(s - 2)) {
mid_term <- alpha^(2 * abs(diffmat[s, s - ms]))
mid_theta[ms] <- ((1 - alpha^(2 * abs(diffmat[s - ms, s -
ms - 1])))) * mid_term/mod$fg(tt = tt, tr = tfu[s - ms],
trecruit = trecruit, m = m)
}
}
# term 1
term_1 <- (alpha^(2 * abs(diffmat[s, 1])))/mod$fg(tt = tt, tr = tfu[1],
trecruit = trecruit, m = m)
# var
kconst <- (sd^2)/(mod$fk(n = n, trecruit = trecruit, m = m) * vphi *
(1 - vphi))
if (s > 2) {
var_beta <- kconst * (term_s + sum(mid_theta) + term_1)
} else {
var_beta <- kconst * (term_s + term_1)
}
}
return(as.numeric(var_beta))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 47. .vbetaUnif
#-------------------------------------------------------------------------------
.vbetaUnif <- function(tfu, tt, sd, alpha, n, trecruit, vphi, rmodel = "fix",
m) {
## Variance of treatment effect estimate for uniform model
mod <- .selectRmodel(rmodel = rmodel)
stfu <- tfu$stfu
tfu <- tfu$tfu
s <- length(tfu)
if (alpha == 0) {
var_beta <- (sd^2)/(mod$fk(n = n, trecruit = trecruit, m = m) *
vphi * (1 - vphi) * mod$fg(tt = tt, tr = tfu[s], trecruit = trecruit,
m = m))
} else {
sum_term <- 0
for (ms in 1:(s - 1)) {
term_m1 <- 1/mod$fg(tt = tt, tr = tfu[ms + 1], trecruit = trecruit,
m = m)
term_m <- 1/mod$fg(tt = tt, tr = tfu[ms], trecruit = trecruit,
m = m)
sum_term[ms] <- ((1 - alpha) * (1 + ms * alpha)/(1 + (ms -
1) * alpha)) * (term_m1 - term_m)
}
term_1 <- 1/mod$fg(tt = tt, tr = tfu[1], trecruit = trecruit, m = m)
var_beta <- ((sd^2)/(mod$fk(n = n, trecruit = trecruit, m = m) *
vphi * (1 - vphi))) * (term_1 + sum(sum_term))
}
return(as.numeric(var_beta))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 48. .vbetaUser
#-------------------------------------------------------------------------------
.vbetaUser <- function(x, s, vphi, alln, vcovmat) {
## Get variance for each interim analyses and study end
all_nx <- alln[x, 2:(s + 1)]
n_0 <- vphi * all_nx
n_1 <- (1 - vphi) * all_nx
## Calculate variance for known sample sizes and covariance matrix
n_0 <- c(n_0, 0)
n_1 <- c(n_1, 0)
V0 <- V1 <- matrix(0, nrow = s, ncol = s)
for (r in 1:s) {
for (c in 1:s) {
mx <- max(r, c)
for (mx in max(r, c):s) {
iest_vcov <- solve(vcovmat[1:mx, 1:mx])
V1[r, c] <- V1[r, c] + (n_1[mx] - n_1[mx + 1]) * iest_vcov[r,
c]
V0[r, c] <- V0[r, c] + (n_0[mx] - n_0[mx + 1]) * iest_vcov[r,
c]
}
}
}
mvcov <- solve(V0) + solve(V1)
var_beta <- mvcov[s, s]
return(var_beta)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# end
#-------------------------------------------------------------------------------
Try the gsearly package in your browser
Any scripts or data that you put into this service are public.
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.