Closed
Description
Hello,
I am trying to understand thetaH1 in getSimulationSurvival. I am interested in an adaptive design involving a thetaH0 bound of 0.9. Note that I am interested in testing H0: HR >= 0.9 vs. H1: HR < 0.9, that means the 0.9 bound is not a non-inferiority bound in the classical sense. In the first scenario I ran a simulation without specifying thetaH1, which should be using the test-statistics of stage 1 and 2 and the observed hazard ratio of stage 2 to reassess the required events for stage 3. I picked out an iteration which had a hazard ratio of 0.4503625 at stage 2. Then I re-ran the simulation, using the same seed, with thetaH1 = 0.4503625 and looked at the same iteration again. Given I used the same seed, it was expected that the results including stage 2 test-statistic and observed hazard ratio are the same. However, this time, the reassessement, which in my understanding only requires the test-statistics of the previous stages (which here are the same between the two scenarios) and an assumed hazard ratio, which I also made sure to be the same given my specification of thetaH1, reported an amount of events for stage 3 which differs from the first scenario.
My code is below. Unless my understanding of the subject is incorrect, I do suspect some kind of inconsistency in regards the thetah0 bound between the two scenarios.
Kind regards
Tim
library(dplyr)
library(rpact)
futilityBounds <- c(0, 0.688)
N <- 14000
informationRates <- c(0.4, 0.6, 1)
dIN <- getDesignInverseNormal(kMax = 3,
alpha = 0.025,
beta = 0.1,
sided = 1,
informationRates = informationRates,
futilityBounds = futilityBounds,
typeOfDesign = "asUser",
userAlphaSpending = c(0, 0.01, 0.025)
)
myCalcEventsFunction_with_division_by_thetaH0 <- function(...,
stage, thetaH0, conditionalPower, estimatedTheta,
plannedEvents, eventsOverStages,
minNumberOfEventsPerStage, maxNumberOfEventsPerStage, allocationRatioPlanned,
conditionalCriticalValue) {
theta <- max(1 + 1e-12, estimatedTheta)
cat(paste0("conditional critical value: ", conditionalCriticalValue, "\n"))
cat(paste0("estimated theta: ", estimatedTheta, "\n"))
requiredStageEvents <-
max(0, conditionalCriticalValue + qnorm(conditionalPower))^2 * (((1+allocationRatioPlanned)^2)/allocationRatioPlanned) / log(theta/thetaH0)^2
requiredStageEvents <- min(
max(minNumberOfEventsPerStage[stage], requiredStageEvents),
maxNumberOfEventsPerStage[stage]
) + eventsOverStages[stage - 1]
return(requiredStageEvents)
}
#wrapper function for convenience, allows us to create the same simulation but it can only differ in thetaH1 and the calceventsfunction
sim_wrapper <- function(thetaH1 = NA_real_, calcEventsFunction = NULL) {
dIN_sim <- getSimulationSurvival(design = dIN,
thetaH0 = 0.9,
directionUpper = F,
pi2 = 0.003,
hazardRatio = c(0.4),
eventTime = 12,
dropoutRate1 = 0.08,
dropoutRate2 = 0.08,
dropoutTime = 12,
accrualTime = c(0,0.01),
maxNumberOfSubjects = N,
plannedEvents = c(20, 30, 50),
maxNumberOfIterations = 75,
longTimeSimulationAllowed = T,
seed = 1910,
showStatistics = T,
maxNumberOfRawDatasetsPerStage = 2,
minNumberOfEventsPerStage = c(20, 10, 20),
maxNumberOfEventsPerStage = c(20, 10, 45),
conditionalPower = 0.9,
calcEventsFunction = calcEventsFunction,
thetaH1 = thetaH1
)}
#standard, use observed HR as thetaH1, interesting iteration 75 chosen
sim_wrapper(thetaH1 = NA_real_, calcEventsFunction = NULL) %>% getData() %>% filter(iterationNumber == 75) %>% print()
#at stage 2, some outputs:
#estimated HR: 0.4503625
#test statistic: 1.896057
#-> for stage 3: increase single events from 20 to 41
#using the same seed, the data looks the same up to IA2. Now using the HR that we have just seen from iteration 75 in the
#"use interim HR as theta case" (0.4503625) as a fixed thetaH1.
#-> if I understand the theory correctly, it should lead to the same SSR, as the test statistics of stage 1 and 2 are identical,
#and the theta for the SSR is the same
sim_wrapper(thetaH1 = 0.4503625, calcEventsFunction = NULL) %>% getData() %>% filter(iterationNumber == 75) %>% print()
#first of all consistency checks, estimated HR and test statistic of stage 2 are equal to the case above...
#estimated HR: 0.4503625
#test statistic: 1.896057
#BUT, the ssr now says we need 30 singular events for stage 3 instead of the 41 from above.
#now with custom function, translated from cpp code to the best of my ability, so it should be the same as rpact native if done correctly
#since we are using only 75 iterations, and we are interested in iteration 75, the last printed line in your console corresponds to the
#conditional critical value and estimated theta of that iteration 75 (didn't know how to do that more neatly)
sim_wrapper(thetaH1 = NA_real_, calcEventsFunction = myCalcEventsFunction_with_division_by_thetaH0) %>% getData() %>% filter(iterationNumber == 75) %>% print()
#conditional critical value: 0.924377935201586 (from console print...)
#estimated theta: 1.79855127135391 (from console print...)
#results look the same as the native rpact ssr function with 41 single events at stage 3
#but focus on estimatedTheta which gets passed internally
#estimated theta: 1.79855127135391
#while hazardRatioEstimateLR 0.4503625
#obviously there is a directionUpper mirroring
#but 1/1.79855127135391 != 0.4503625
#maybe estimatedTheta already the theta which would lead to the same test-statistic given thetaH0 = 0.9
#but 1/(0.4503625/0.9) != 1.79855127135391
#however 1/(0.4503625/0.9) == 1.79855127135391/0.9
#now with a fixed theta, again, we should not have differences here for iteration 75 as the
#chosen fixed theta is exactly the theta observed in stage 2 of iteration 75
sim_wrapper(thetaH1 = 0.4503625, calcEventsFunction = myCalcEventsFunction_with_division_by_thetaH0) %>% getData() %>% filter(iterationNumber == 75) %>% print()
#conditional critical value: 0.924377935201586 -> same as before, as expected
#estimated theta: 1.99839018568375 -> not the same
#this time this equation holds
#1/(0.4503625/0.9) == 1.99839018568375
#now I try this by hand, and on the negative side of the normal distribution first so i dont get confused by the directionUpper mirroring...
#max(0, conditionalCriticalValue + qnorm(conditionalPower))^2 * (((1+allocationRatioPlanned)^2)/allocationRatioPlanned) / log(theta/thetaH0)^2
min(0, -0.924377935201586 - qnorm(0.9))^2 * (((1+1)^2)/1) / log(0.4503625/0.9)^2
#40.6071 -> I guess this would be consistent with the 41 from scenario 1 with thetaH1=NA_real_ and rpact native ssr function
#however, I plugged in the true HR of 0.4503625
#if I want to do this the directionUpper way (applying thetaH0 first, then mirroring), trivially the same as log(1) is 0.
max(0, 0.924377935201586 + qnorm(0.9))^2 * (((1+1)^2)/1) / log(1/(0.4503625/0.9))^2
#40.6071
#what happens in the fixed thetah1 case where it said we need 30 events for stage 3?
#estimatedTheta internally is 1.99839018568375, i.e. 1/(0.4503625/0.9)
#but then gets plugged in
##max(0, conditionalCriticalValue + qnorm(conditionalPower))^2 * (((1+allocationRatioPlanned)^2)/allocationRatioPlanned) / log(theta/thetaH0)^2
#which divides it again by 0.9
#i.e.
max(0, 0.924377935201586 + qnorm(0.9))^2 * (((1+1)^2)/1) / log(((1/(0.4503625/0.9))/0.9))^2
#30.58873
#((1/(0.4503625/0.9))/0.9) = 2.220434 = 1/0.4503625