##################################################################
# This code generates the numerical results in chapter 2 ##
##################################################################
# load the survival package
library(survival)
# install and load the WR package for
# 1. dataset hfaction_cpx9;
# 2. function WRrec() for win ratio test (of recurrent events and death)
# 3. functions base() and WRSS() for sample size calculation
# install.packages("WR")
library(WR)
library(tidyverse) # for data wrangling
##### Read in HF-ACTION DATA########
# same as rmt::hfaction used in chap 1
# (except for status coding)
data(hfaction_cpx9)
hfaction <- hfaction_cpx9
head(hfaction)
# count unique patients in each arm
hfaction |>
group_by(trt_ab) |>
distinct(patid) |>
count(trt_ab)
#### demo ############
obj <- WRrec(ID = hfaction$patid, time = hfaction$time,
status = hfaction$status, trt = hfaction$trt_ab,
strata = hfaction$age60, naive = TRUE)
# summary results
obj
# LWR
beta <- obj$log.WR
se <- obj$se
# test
pval <- 2*(1-pnorm(abs(beta/se)))
pval
# NWR
beta.naive<-obj$log.WR.naive
se.naive<-obj$se.naive
# test
pval.naive<-2*(1-pnorm(abs(beta.naive/se.naive)))
pval.naive
# FWR
beta.FI<-obj$log.WR.FI
se.FI<-obj$se.FI
# test
pval.FI<-2*(1-pnorm(abs(beta.FI/se.FI)))
pval.FI
#################################
# Win ratio analyses: tabulate #
#################################
data <- hfaction
### create a dataset with only the first hospitalization data.H1
# hospitalization data
tmpH <- data[data$status==2,]
# get the first record of each id
o <- order(tmpH$patid,tmpH$time)
tmpH <- tmpH[o,]
tmpFH <- tmpH[!duplicated(tmpH$patid),]
# combine it with mortality data
data.H1 <- rbind(tmpFH,data[data$status!=2,])
o <- order(data.H1$patid,data.H1$time)
data.H1 <- data.H1[o,]
# Function to create table for estimates of
# PWR, NWR, FWR, and LWR, 95% CI nad p-values
# ind: index for data
# ind1: index for data.H1
# r: number of decimal point
gwr.fun=function(ind,ind1,r=2){
# fit NWR, FWR, and LWR
obj <- WRrec(ID=data$patid[ind],time=data$time[ind],status=data$status[ind],
trt=data$trt_ab[ind],strata=data$age60[ind],naive=T)
# fit sWR
obj1 <- WRrec(ID=data.H1$patid[ind1],time=data.H1$time[ind1],status=data.H1$status[ind1],
trt=data.H1$trt_ab[ind1],strata=data.H1$age60[ind1],naive=F)
# critical value of p=0.05
za <- qnorm(0.975)
## LWR
beta <- obj$log.WR
se <- obj$se
theta <- obj$theta
r4 <- c(paste0(round(100*theta[1],1),"%"),
paste0(round(100*theta[2],1),"%"),
paste0(round(exp(beta),r)," (",round(exp(beta-za*se),r),", ",round(exp(beta+za*se),r),")"),
round(1-pchisq((beta/se)^2,1),3)
)
## PWR
beta1=obj1$log.WR
se1=obj1$se
theta1=obj1$theta
r1=c(paste0(round(100*theta1[1],1),"%"),
paste0(round(100*theta1[2],1),"%"),
paste0(round(exp(beta1),r)," (",round(exp(beta1-za*se1),r),", ",round(exp(beta1+za*se1),r),")"),
round(1-pchisq((beta1/se1)^2,1),3)
)
## NWR
beta.naive=obj$log.WR.naive
se.naive=obj$se.naive
theta.naive=obj$theta.naive
r2=c(paste0(round(100*theta.naive[1],1),"%"),
paste0(round(100*theta.naive[2],1),"%"),
paste0(round(exp(beta.naive),r)," (",round(exp(beta.naive-za*se.naive),r),", ",
round(exp(beta.naive+za*se.naive),r),")"),
round(1-pchisq((beta.naive/se.naive)^2,1),3)
)
## FWR
beta.FI=obj$log.WR.FI
se.FI=obj$se.FI
theta.FI=obj$theta.FI
r3=c(paste0(round(100*theta.FI[1],1),"%"),
paste0(round(100*theta.FI[2],1),"%"),
paste0(round(exp(beta.FI),r)," (",round(exp(beta.FI-za*se.FI),r),", ",
round(exp(beta.FI+za*se.FI),r),")"),
round(1-pchisq((beta.FI/se.FI)^2,1),3)
)
result=rbind(r1,r2,r3,r4)
rownames(result)=c("PWR","NWR","FWR","LWR")
return(result)
}
# Create table
## Age <= 60 years
ind=(data$age60==0)
ind1=(data.H1$age60==0)
result.lt60=gwr.fun(ind,ind1,r=2)
## Age > 60 years
ind=(data$age60==1)
ind1=(data.H1$age60==1)
result.ge60=gwr.fun(ind,ind1,r=2)
## overall
ind=rep(T,nrow(data))
ind1=rep(T,nrow(data.H1))
result.all=gwr.fun(ind,ind1,r=2)
# combine results
results=rbind(result.lt60,result.ge60,result.all)
colnames(results)=c("Win", "Loss", "Win ratio (95% CI)", "p-value")
noquote(results)
############################################################################
# Sample size calculation #
############################################################################
# get training arm data
pilot <- hfaction |>
filter(trt_ab == 1)
# number of subjects
pilot |> distinct(patid) |>
count()
############## estimate parameters ##############
# Get the variables from pilot dataset
# to estimate baseline parameters
# lambda_D, lambda_H, kappa
outcome_base <- gumbel.est(pilot$patid, pilot$time / 12, pilot$status)
lambda_D <- outcome_base$lambda_D
lambda_H <- outcome_base$lambda_H
kappa <- outcome_base$kappa
lambda_D
lambda_H
kappa
## Kendall's rank correlation
1 - 1/kappa
#> [1] 0.360812
### demo ###################
# set design parameters
tau_b <- 3
tau <- 4
lambda_L <- 0.01
# use base() function to compute zeta2 and delta
## may take up to 30s
bparam <- base(lambda_D, lambda_H, kappa, tau_b, tau, lambda_L)
# compute sample size under HRs 0.8 and 0.9
# for death and nonfatal event, respectively
obj <- WRSS(xi = log(c(0.9,0.8)), bparam = bparam, q = 0.5, alpha = 0.05,
power = 0.8)
obj$n
## effect size specification
thetaD <- seq(0.6, 0.95,by = 0.05) ## hazard ratio for death
thetaH <- seq(0.6, 0.95,by = 0.05) ## hazard ratio for hospitalization
## create a matrix "SS08" for sample size powered at 80%
## under each combination of thetaD and thetaH
mD <- length(thetaD)
mH <- length(thetaH)
SS08 <- matrix(NA, mD, mH)
rownames(SS08) <- thetaD
colnames(SS08) <- thetaH
## fill in the computed sample size values
for (i in 1:mD){
for (j in 1:mH){
## sample size under hazard ratios thetaD[i] for death and thetaH[j] for hospitalization
SS08[i,j]<-WRSS(xi=log(c(thetaD[i],thetaH[j])),bparam=bparam,q=0.5,alpha=0.05,
power=0.8)$n
}
}
## print the calculated sample sizes
print(SS08)
# 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95
# 0.6 198.1636 261.3528 351.2579 484.2468 691.0512 1034.785 1661.812 2976.570
# 0.65 209.5088 278.6457 378.4143 528.6465 767.7521 1177.940 1961.252 3727.580
# 0.7 220.9045 296.2334 406.4659 575.4320 850.7474 1338.718 2316.976 4708.455
# 0.75 232.3703 314.1480 435.4823 624.8038 940.7326 1519.930 2742.885 6016.294
# 0.8 243.9237 332.4188 465.5327 676.9740 1038.4859 1724.935 3257.279 7803.776
# 0.85 255.5795 351.0732 496.6864 732.1690 1144.8828 1957.769 3884.618 10321.618
# 0.9 267.3513 370.1374 529.0141 790.6317 1260.9129 2223.322 4658.150 14003.925
# 0.95 279.2513 389.6366 562.5886 852.6243 1387.6996 2527.561 5623.948 19653.870
## repeating the same calculation for power = 90%
SS09 <- matrix(NA, mD, mH)
rownames(SS09) <- thetaD
colnames(SS09) < -thetaH
## fill in the computed sample size values
for (i in 1:mD){
for (j in 1:mH){
## sample size under hazard ratios thetaD[i] for death and thetaH[j] for hospitalization
SS09[i,j]<-WRSS(xi=log(c(thetaD[i],thetaH[j])),bparam=bparam,q=0.5,alpha=0.05,
power=0.9)$n
}
}
## print the calculated sample sizes
print(SS09)
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
# 0.6 265.2849 349.8773 470.2346 648.2691 925.1215 1385.283 2224.695 3984.783
# 0.65 280.4729 373.0275 506.5893 707.7076 1027.8022 1576.928 2625.560 4990.172
# 0.7 295.7284 396.5725 544.1425 770.3402 1138.9094 1792.164 3101.773 6303.285
# 0.75 311.0780 420.5551 582.9873 836.4350 1259.3740 2034.755 3671.945 8054.111
# 0.8 326.5446 445.0145 623.2161 906.2762 1390.2379 2309.198 4360.573 10447.042
# 0.85 342.1484 469.9874 664.9221 980.1666 1532.6732 2620.897 5200.401 13817.718
# 0.9 357.9075 495.5089 708.1999 1058.4316 1688.0046 2976.398 6235.942 18747.282
# 0.95 373.8383 521.6128 753.1465 1141.4221 1857.7360 3383.687 7528.871 26310.956
oldpar <- par(mfrow = par("mfrow"))
par(mfrow=c(1,2))
persp(thetaD, thetaH, SS08/1000, theta = 50, phi = 15, expand = 0.8, col = "gray",
ltheta = 180, lphi=180, shade = 0.75,
ticktype = "detailed",
xlab = "\n HR on Death", ylab = "\n HR on Hospitalization",
zlab=paste0("\n Sample Size (10e3)"),
main="Power = 80%",
zlim=c(0, 26)
)
persp(thetaD, thetaH, SS09/1000, theta = 50, phi = 15, expand = 0.8, col = "gray",
ltheta = 180, lphi=180, shade = 0.75,
ticktype = "detailed",
xlab = "\nHR on Death", ylab = "\nHR on Hospitalization",
zlab=paste0("\n Sample Size (10e3)"),
main="Power = 90%",
zlim=c(0, 26)
)