Nothing
########################################################################
# Functions for composite outcomes (under Gumbel model) ##
########################################################################
# The delta function using numerical integration
# the output is a vector of (delta_1,delta_2)
delta.fun = function(baseline){
lambda_D=baseline$lambda_D
lambda_H=baseline$lambda_H
kappa=baseline$kappa
tau_b=baseline$tau_b
tau=baseline$tau
lambda_L=baseline$lambda_L
# the input "x" is a vector containing two entries (s,t)
SD = function(x){
exp(-lambda_D*x[2])
}
g_tilde = function(x){
if(x[1] > tau - tau_b){
2/tau_b^2*exp(-2*lambda_L*x[1])*(tau-x[1])*(lambda_L*(tau-x[1])+1)
}else{
2*lambda_L*exp(-2*lambda_L*x[1])
}
}
delta1 = adaptIntegrate(function(x){
tem = ((lambda_D*x[1])^kappa+
(lambda_H*x[2])^kappa)^(1/kappa)
H0 = exp(-tem)
lambda1= tem^(1-2*kappa)*(1-kappa)*(lambda_D*lambda_H*x[1])^kappa*x[2]^(kappa-1)
return(ifelse(x[2]<=x[1],1,0)*(lambda_D*SD(x)^2 + H0^2*lambda1)*g_tilde(x))
},
c(0,0), c(tau,tau))
delta2 = adaptIntegrate(function(x){
tem = ((lambda_D*x[1])^kappa+
(lambda_H*x[2])^kappa)^(1/kappa)
H0 = exp(-tem)
lambda2 = (1-kappa)*tem^(1-2*kappa) *
(lambda_H)^(2*kappa)*x[2]^(2*kappa-1)+
kappa * tem^(1-kappa) *
x[2]^(kappa-1)*lambda_H^kappa
return(ifelse(x[2]<=x[1],1,0) * (H0^2*lambda2)*g_tilde(x))
},
c(0,0), c(tau,tau))
return(c(delta1 = delta1$integral,delta2 = delta2$integral))
}
###############################################################
# Win function for composite outcomes
# see argument "winfun" in WR.anal()
#############################################################
Wp = function(y1,y2){
ifelse(y2[1]<min(c(y1[1],y1[3],y2[3])),1,0) +
ifelse((min(c(y1[1],y2[1])) > min(c(y1[3], y2[3])) ) &
(y2[2] < min(c(y1[2],y1[3],y2[3]))),1,0)
}
## This is the win-loss difference function
win.comp = function(y1,y2){Wp(y1,y2) - Wp(y2,y1)}
# monte-carlo integration of zeta_0^2
zeta2.fun=function(baseline,N=2000,seed){
lambda_D=baseline$lambda_D
lambda_H=baseline$lambda_H
kappa=baseline$kappa
tau_b=baseline$tau_b
tau=baseline$tau
lambda_L=baseline$lambda_L
#simulate death, hosp, and censoring
# Death
set.seed(seed)
outcome=rgumbel(N,alpha=kappa,dim=2)
D=-log(outcome[,1])/lambda_D
H=-log(outcome[,2])/lambda_H
set.seed(seed)
Ca=tau_b*runif(N)+tau-tau_b
set.seed(seed)
L=rexp(N)/lambda_L
C=pmin(Ca,L)
dat=cbind(D,H,C)
tmp=combn(1:nrow(dat),2,function(x) win.comp(dat[x[1],],dat[x[2],]))
WL.matrix=matrix(0,N,N)
WL.matrix[lower.tri(WL.matrix)]=tmp
WL.matrix=t(WL.matrix)-WL.matrix
w0=sum(abs(WL.matrix))/(2*N*(N-1))
ranks=rowSums(WL.matrix)/(N-1)
zeta2=sum(ranks^2)/(N-3)
return(list(w0=w0,zeta2=zeta2))
}
#######################################################
# Function to compute the baseline parameters
# in the sample size formula of Pocock's win ratio
######################################################
#' Compute the baseline parameters needed for sample size calculation for
#' standard win ratio test
#'
#' @description Compute the baseline parameters \eqn{\zeta_0^2} and \eqn{\boldsymbol\delta_0}
#' needed for sample size calculation for standard win ratio test (see \code{\link{WRSS}}).
#' The calculation is based
#' on a Gumbel--Hougaard copula model for survival time \eqn{D^{(a)}} and nonfatal event
#' time \eqn{T^{(a)}} for group \eqn{a} (1: treatment; 0: control):
#'\deqn{{P}(D^{(a)}>s, T^{(a)}>t) =\exp\left(-\left[\left\{\exp(a\xi_1)\lambda_Ds\right\}^\kappa+
#' \left\{\exp(a\xi_2)\lambda_Ht\right\}^\kappa\right]^{1/\kappa}\right),}
#' where \eqn{\xi_1} and \eqn{\xi_2} are the component-wise log-hazard ratios to be used
#' as effect size in \code{\link{WRSS}}.
#' We also assume that patients are recruited uniformly over the period \eqn{[0, \tau_b]}
#' and followed until time \eqn{\tau} (\eqn{\tau\geq\tau_b}), with an exponential
#' loss-to-follow-up hazard \eqn{\lambda_L}.
#' @param lambda_D Baseline hazard \eqn{\lambda_D} for death.
#' @param lambda_H Baseline hazard \eqn{\lambda_H} for nonfatal event.
#' @param kappa Gumbel--Hougaard copula correlation parameter \eqn{\kappa}.
#' @param tau_b Length of the initial (uniform) accrual period \eqn{\tau_b}.
#' @param tau Total length of follow-up \eqn{\tau}.
#' @param lambda_L Exponential hazard rate \eqn{\lambda_L} for random loss to follow-up.
#' @param N Simulated sample size for monte-carlo integration.
#' @param seed Seed for monte-carlo simulation.
#' @return A list containing real number \code{zeta2} for \eqn{\zeta_0^2}
#' and bivariate vector \code{delta} for \eqn{\boldsymbol\delta_0}.
#' @seealso \code{\link{gumbel.est}}, \code{\link{WRSS}}
#' @import gumbel
#' @import cubature
#' @importFrom utils combn tail
#' @importFrom stats rexp runif
#' @keywords WRSS
#' @examples
#' # see the example for WRSS
#' @export
#' @references Mao, L., Kim, K. and Miao, X. (2021). Sample size formula for general win ratio analysis.
#' Biometrics, https://doi.org/10.1111/biom.13501.
base=function(lambda_D,lambda_H,kappa,tau_b,tau,lambda_L,N=1000,seed=12345){
baseline=list(lambda_D=lambda_D,lambda_H=lambda_H,kappa=kappa,tau_b=tau_b,
tau=tau,lambda_L=lambda_L)
obj=zeta2.fun(baseline,N,seed)
zeta2=obj$zeta2
w0=obj$w0
delta=delta.fun(baseline)
return(list(zeta2=zeta2,w0=w0,delta=delta))
}
#' Compute the sample size for standard win ratio test
#'
#' @description Compute the sample size for standard win ratio test.
#'
#' @param xi A bivariate vector of hypothesized component-wise (treatment-to-control) log-hazard
#' ratios under the Gumbel--Hougaard copula model described in \code{\link{base}}.
#' @param bparam A list containing baseline parameters \code{zeta2} for \eqn{\zeta_0^2}
#' and \code{delta} for \eqn{\boldsymbol\delta_0}; Can directly use the output of \link{base}.
#' @param q Proportion of patients assigned to treatment.
#' @param alpha Type I error rate.
#' @param side 2-sided or 1-sided test.
#' @param power Target power.
#' @return A list containing \code{n}, the computed sample size.
#' @seealso \code{\link{gumbel.est}}, \code{\link{base}}
#' @importFrom stats qnorm
#' @keywords WRSS
#' @export
#' @references Mao, L., Kim, K. and Miao, X. (2021). Sample size formula for general win ratio analysis.
#' Biometrics, https://doi.org/10.1111/biom.13501.
#' @export
#' @examples
#' # The following is not run in package checking to save time.
#' \dontrun{
#'## load the package and pilot dataset
#'library(WR)
#'head(hfaction_cpx9)
#'dat<-hfaction_cpx9
#'## subset to control group
#'pilot<-dat[dat$trt_ab==0,]
#'
#'## get the data ready for gumbel.est()
#'id<-pilot$patid
#'## convert time from month to year
#'time<-pilot$time/12
#'status<-pilot$status
#'## compute the baseline parameters for the Gumbel--Hougaard
#'## copula for death and hospitalization
#'gum<-gumbel.est(id, time, status)
#'
#'## get the baseline parameters
#'lambda_D<-gum$lambda_D
#'lambda_H<-gum$lambda_H
#'kappa<-gum$kappa
#'## set up design parameters and use base()
#'## to calculate bparam for WRSS()
#'# max follow-up 4 years
#'tau<-4
#'# 3 years of initial accrual
#'tau_b<-3
#'# loss to follow-up rate
#'lambda_L=0.05
#'# compute the baseline parameters
#'bparam<-base(lambda_D,lambda_H,kappa,tau_b,tau,lambda_L)
#'bparam
#'
#'## sample size with power=0.8 under hazard ratios
#'## 0.9 and 0.8 for death and hospitalization, respectively.
#'WRSS(xi=log(c(0.9,0.8)),bparam=bparam,q=0.5,alpha=0.05,
#' power=0.8)$n
#'## sample size under the same set-up but with power 0.9
#'WRSS(xi=log(c(0.9,0.8)),bparam=bparam,q=0.5,alpha=0.05,
#' power=0.9)$n
#' }
WRSS=function(xi,bparam,q=0.5,alpha=0.05,side=2,power=0.8){
xi=as.matrix(xi)
#sample size formula
za=qnorm(1-alpha/side)
zb=qnorm(power)
zeta2=bparam$zeta2
delta=bparam$delta
n=zeta2*(za+zb)^2/(q*(1-q)*(t(xi)%*%delta)^2)
return(list(zeta2=zeta2,delta=delta,n=n,xi=xi))
}
#===========================================================================================
# Estimating the baseline parameters lambda_D, lambda_H
# and kappa from pilot study data
#===========================================================================================
#' Estimate baseline parameters in the Gumbel--Hougaard model for sample size
#' calculation using pilot data
#'
#' @description Estimate baseline parameters in the Gumbel--Hougaard model
#' described in \code{\link{base}} for sample size calculation using pilot study data.
#'
#' @param id A vector of unique patient identifiers.
#' @param time A numeric vector of event times.
#' @param status A vector of event type variable; 2 = nonfatal event, 1 = death,
#' and 0 = censoring.
#' @return A list containing \code{lambda_D} for \eqn{\lambda_D},
#' \code{lambda_H} for \eqn{\lambda_H}, and \code{kappa} for \eqn{\kappa}
#' in the Gumbel--Hougaard model.
#' @seealso \code{\link{base}}, \code{\link{WRSS}}
#' @examples
#' # see the example for WRSS
#' @export
#' @keywords WRSS
#' @references Mao, L., Kim, K. and Miao, X. (2021). Sample size formula for general win ratio analysis.
#' Biometrics, https://doi.org/10.1111/biom.13501.
gumbel.est=function(id, time, status){
data=data.frame(id,time,status)
#Sort the data by time within each id
o=order(data$id,data$time)
data=data[o,]
######################################
#Composite (time-to-first) event rate#
######################################
#subset to first event data
#Sort the data by time within each id
# get the first row for each id
data.CE=data[!duplicated(data$id),]
# event rate
lambda_CE=sum(data.CE$status>0)/sum(data.CE$time)
#########################################
# Cause-specific hazard rate for hosp
##########################################
lambda_sH=sum(data.CE$status==2)/sum(data.CE$time)
######################
# Death rate ##
######################
data.D=data[data$status<2,]
# event rate
lambda_D=sum(data.D$status)/sum(data.D$time)
kappa=log(1-lambda_sH/lambda_CE)/log(lambda_D/lambda_CE)
lambda_H=(lambda_CE^kappa-lambda_D^kappa)^(1/kappa)
return(list(lambda_D=lambda_D,lambda_H=lambda_H,kappa=kappa))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.