*! Date : 08 May 2018 *! Version : 1.0 /* 08/05/18 v1.0 Initial version complete */ program define des_curtailed, rclass version 11.0 syntax, [J(integer 2) pi0(real 0.1) pi1(real 0.3) Alpha(real 0.05) /// Beta(real 0.2) thetaf(numlist) thetae(numlist) nmin(integer 1) /// nmax(integer 30) FUTility(integer 1) EFFicacy(integer 0) /// optimality(string) EQual_n(integer 0) ENSign(integer 0) /// SUMmary(integer 0) PLot *] local xopt `"`options'"' preserve if ("`optimality'" == "") { local optimality "null_ess" } ///// Perform checks on input variables //////////////////////////////////////// if (`j' <= 0) { di "{error} J must be an integer greater than or equal to one." exit(198) } if ((`pi0' <= 0) | (`pi0' >= 1)) { di "{error} pi0 must be a real strictly between 0 and 1." exit(198) } if ((`pi1' <= 0) | (`pi1' >= 1)) { di "{error} pi1 must be a real strictly between 0 and 1." exit(198) } if (`pi0' >= `pi1') { di "{error} pi0 must be a real strictly less than the real pi1." exit(198) } if ((`alpha' <= 0) | (`alpha' >= 1)) { di "{error} alpha must be a real strictly between 0 and 1." exit(198) } if ((`beta' <= 0) | (`beta' >= 1)) { di "{error} beta must be a real strictly between 0 and 1." exit(198) } if ("`thetaf'" ~= "") { local lenthetaf:list sizeof thetaf if (`lenthetaf' ~= `j') { di "{error} thetaf must be a numlist of length J." exit(198) } if ("`thetae'" == "") { forvalues i = 1/`lenthetaf' { local thetafi:word `i' of `thetaf' if (`thetafi' > 1) { di "{error} Elements in thetaf must be less than or equal to 1." exit(198) } } } } if ("`thetae'" ~= "") { local lenthetae:list sizeof thetae if (`lenthetae' ~= `j') { di "{error} thetae must be a numlist of length J." exit(198) } if ("`thetaf'" == "") { forvalues i = 1/`lenthetae' { local thetaei:word `i' of `thetae' if (`thetaei' < 0) { di "{error} Elements in thetae must be greater than or equal to 0." exit(198) } } } } if ("`thetaf'" ~= "" & "`thetae'" ~= "") { forvalues i = 1/`lenthetaf' { local thetafi:word `i' of `thetaf' local thetaei:word `i' of `thetae' if (`thetafi' > `thetaei') { di "{error} Elements in thetaf must be less than or equal to their corresponding element in thetae." exit(198) } } } if (`nmin' <= 0) { di "{error} nmin must be a strictly positive integer." exit(198) } if (`nmax' <= 0) { di "{error} nmax must be a strictly positive integer." exit(198) } if (`nmin' > `nmax') { di "{error} nmin must be an integer less than or equal to the integer nmax." exit(198) } if (`futility' ~= 1 & `efficacy' ~= 1) { di "{error} At least one of futility and efficacy must be equal to 1." exit(198) } if ("`optimality'" ~= "minimax" & "`optimality'" ~= "null_ess" & "`optimality'" ~= "alt_ess" & "`optimality'" ~= "null_med" & "`optimality'" ~= "alt_med") { di "{error} optimality must be one of: minimax, null_ess, alt_ess, null_med, and alt_med." exit(198) } if (`j' > 2) { di "{error} J > 2 is not currently supported." exit(198) } // Set up matrices to pass to mata if ("`thetaf'" ~= "") { local matathetaf "" foreach i of local thetaf{ if "`matathetaf'" == "" local matathetaf "`i'" else local matathetaf "`matathetaf',`i'" } mat thetaf = (`matathetaf') } else { mat thetaf = . } if ("`thetae'" ~= "") { local matathetae "" foreach i of local thetae{ if "`matathetae'" == "" local matathetae "`i'" else local matathetae "`matathetae',`i'" } mat thetae = (`matathetae') } else { mat thetae = . } ///// Compute and Output /////////////////////////////////////////////////////// mata: DesCurtailed(`j', `pi0', `pi1', `alpha', `beta', `nmin', `nmax', `futility', `efficacy', "`optimality'", `equal_n', `ensign', `summary', "`plot'", `"`xopt'"') return mat des = des matrix colnames opchar = "P(pi0)" "P(pi1)" "ESS(pi0)" "ESS(pi1)" "Med(pi0)" "Med(pi1)" "VSS(pi0)" "VSS(pi1)" "max(N)" return mat opchar = opchar return mat J = J return mat n = n return mat a = a return mat r = r return mat J_curt = J_curt return mat n_curt = n_curt return mat a_curt = a_curt return mat r_curt = r_curt if (`j' == 1) { matrix colnames feasible = n1 a1 r1 "P(pi0)" "P(pi1)" } else if (`j' == 2) { matrix colnames feasible = n1 n2 a1 a2 r1 r2 "P(pi0)" "P(pi1)" "PET1(pi0)" "PET1(pi1)" "ESS(pi0)" "ESS(pi1)" "Med(pi0)" "Med(pi1)" "VSS(pi0)" "VSS(pi1)" "max(N)" } else { matrix colnames feasible = n1 n2 n3 a1 a2 a3 r1 r2 r3 "P(pi0)" "P(pi1)" "PET1(pi0)" "PET1(pi1)" "PET2(pi0)" "PET2(pi1)" "ESS(pi0)" "ESS(pi1)" "Med(pi0)" "Med(pi1)" "VSS(pi0)" "VSS(pi1)" "max(N)" } return mat feasible = feasible return mat pi0 = pi0 return mat pi1 = pi1 return mat alpha = alpha return mat beta = beta return mat thetaf = thetaF return mat thetae = thetaE return mat nmin = Nmin return mat nmax = Nmax return mat futility = futility return mat efficacy = efficacy return local optimality "`optimality'" return mat equal_n = equal_n return mat ensign = ensign return mat summary = summary restore end ///// Mata //////////////////////////////////////////////////////////////////// mata: void DesCurtailed(J, pi0, pi1, alpha, beta, Nmin, Nmax, futility, efficacy, optimality, equal_n, ensign, summary, plot, xopt) { thetaF = st_matrix("thetaf") thetaE = st_matrix("thetae") if (thetaF == .) { thetaF = J(1, J, 0) } if (thetaE == .) { thetaE = J(1, J, 1) } ///// Print Summary ////////////////////////////////////////////////////////// ///// Main Computations ////////////////////////////////////////////////////// if (J == 1) { dbinomial_pi0 = J(Nmax + 1, Nmax, 0) dbinomial_pi1 = J(Nmax + 1, Nmax, 0) for (i = 1; i <= Nmax; i++) { dbinomial_pi0[1::(i + 1), i] = binomialp(i, 0::i, pi0) dbinomial_pi1[1::(i + 1), i] = binomialp(i, 0::i, pi1) } if (exact == 1) { possible = permutations(0::Nmax, 2) possible = possible[mm_which(possible[, 1] :> possible[, 2] :& 0 :< possible[, 1]), ] } else { possible = Nmin::Nmax possible = (possible, round(possible*pi0 + invnormal(1 - alpha)*sqrt(possible*pi0*(1 - pi0)))) possible = possible[mm_which(possible[, 1] :> possible[, 2] :& 0 :<= possible[, 2]), ] } possible = (possible, possible[, 2] + J(rows(possible), 1, 1), J(rows(possible), 2, 0)) for (i = 1; i <= rows(possible); i++) { possible[i, (4, 5)] = (sum(dbinomial_pi0[(possible[i, 3] + 1)::possible[i, 1], possible[i, 1]]), sum(dbinomial_pi1[(possible[i, 3] + 1)::possible[i, 1], possible[i, 1]])) if (allof((summary == 1, mod(i, 1000) == 0), 1)) { printf("{txt}...{res}%g{txt} designs evaluated...\n", i) } } if (summary == 1) { printf("...determining feasible designs, and then the optimal design...") } feasible = possible[mm_which(possible[, 4] :<= alpha :& possible[, 5] :>= 1 - beta), ] if (rows(feasible) > 0) { feasible = sort(feasible, (1, -5)) des = feasible[1, ] n = feasible[1, 1] a = feasible[1, 2] r = feasible[1, 3] if (plot ~= "") { // Need plot of state-space here (use xopt) } } else { feasible = . n = . a = . r = . des = . } } else { feasible = IntDesGS(J, pi0, pi1, alpha, beta, Nmin, Nmax, futility, efficacy, equal_n, ensign, summary) if (feasible[1, 1] > 0) { if (summary == 1) { printf("...feasible designs in range of considered maximal allowed sample size identified....") printf("...now identifying optimal design(s)...") } feasible = (feasible[, 1::(2*J + J - 1)], feasible[, 2*J] + J(rows(feasible), 1, 1), feasible[, (2*J + J)::cols(feasible)], rowsum(feasible[, 1::J])) if (futility ~= 1) { feasible[, (J + 1)::(2*J - 1)] = J(rows(feasible), J - 1, .) } if (efficacy ~= 1) { feasible[, (2*J + 1)::(3*J - 1)] = J(rows(feasible), J - 1, .) } if (optimality == "null-ess") { feasible = sort(feasible, (5*J + 1, 5*J + 7)) } else if (optimality == "alt-ess") { feasible = sort(feasible, (5*J + 2, 5*J + 7)) } else if (optimality == "null-med") { feasible = sort(feasible, (5*J + 3, 5*J + 7)) } else if (optimality == "alt-med") { feasible = sort(feasible, (5*J + 4, 5*J + 7)) } else if (optimality == "minimax") { feasible = sort(feasible, (5*J + 7, 5*J + 1)) } des = feasible[1, ] n = feasible[1, 1::J] a = feasible[1, (J + 1)::(2*J)] r = feasible[1, (2*J + 1)::(3*J)] } else { feasible = des = . if (summary) { printf("...no feasible designs found in range of considered maximal allowed sample size. Consider decreasing nmin and increasing nmax...") } } } if (des ~= .) { cp_mat = J(sum(n) + 1, sum(n), 0) for (m = 1; m <= sum(n); m++) { for (s = 0; s <= m; s++) { cp_mat[s + 1, m] = cr_gs(pi1, s, m, J, a, r, n) } } a_curt = J(1, sum(n), .) r_curt = J(1, sum(n), .) n_curt = J(1, sum(n), 1) for (i = 1; i <= sum(n); i++) { j = min(mm_which(i :<= runningsum(n))) if (any(cp_mat[1::(i + 1), i] :<= thetaF[j])) { a_curt[i] = max(mm_which(cp_mat[1::(i + 1), i] :<= thetaF[j])) - 1 } if (any(cp_mat[, i] :>= thetaE[j])) { r_curt[i] = min(mm_which(cp_mat[, i] :>= thetaE[j])) - 1 } } opchar_H0 = int_opchar_curtailed(pi0, J, sum(n), a_curt, r_curt, n, n_curt, runningsum(n), runningsum(n_curt), (1::J)') opchar_H1 = int_opchar_curtailed(pi1, J, sum(n), a_curt, r_curt, n, n_curt, runningsum(n), runningsum(n_curt), (1::J)') a_curt[mm_which(a_curt :== -1)] = J(1, length(mm_which(a_curt :== -1)), .) r_curt[mm_which(r_curt :== sum(n) + 1)] = J(1, length(mm_which(r_curt :== sum(n) + 1)), .) maxN = max(mm_which(a_curt :+ 1 :< r_curt)) + 1 J_curt = sum(n) des = (J, n, a, r, n_curt, a_curt, r_curt, opchar_H0[2], opchar_H1[2], opchar_H0[(6 + 4*J_curt + 2*J)::(5 + 4*J_curt + 3*J)], opchar_H1[(6 + 4*J_curt + 2*J)::(5 + 4*J_curt + 3*J)], opchar_H0[3], opchar_H1[3], opchar_H0[5], opchar_H1[5], opchar_H0[4], opchar_H1[4], maxN) opchar = (opchar_H0[2], opchar_H1[2], opchar_H0[3], opchar_H1[3], opchar_H0[5], opchar_H1[5], opchar_H0[4], opchar_H1[4], maxN) if (plot ~= "") { a_curt[mm_which(a_curt :== .)] = J(1, length(mm_which(a_curt :== .)), -1) r_curt[mm_which(r_curt :== .)] = J(1, length(mm_which(r_curt :== .)), sum(n) + 1) states = permutations((0::n_curt[1]), 2) states = states[mm_which((states[, 1] :<= states[, 2]) :& (states[, 2] :> 0)), ] outcome = J(rows(states), 1, 0) for (i = 1; i <= rows(states); i++) { if (states[i, 1] <= a_curt[1] & states[i, 2] == n_curt[1]) { outcome[i] = 1 } else if (states[i, 1] >= r_curt[1] & states[i, 2] == n_curt[1]) { outcome[i] = 2 } } states = (states, outcome) cont = (max((0, a_curt[1] + 1)), min((r_curt[1] - 1, n_curt[1]))) for (j = 2; j <= J_curt; j++) { vals_j = permutations((0::n_curt[j]), 2) vals_j = vals_j[mm_which((vals_j[, 1] :<= vals_j[, 2]) :& (vals_j[, 2] :> 0)), ] states_j = J(1, 2, .) for (sj = cont[1]; sj <= cont[2]; sj++) { states_j = (states_j \ (vals_j[, 1] :+ sj, vals_j[, 2] :+ sum(n_curt[1::(j - 1)]))) } states_j = states_j[2::rows(states_j), ] outcome = J(rows(states_j), 1, 0) for (i = 1; i <= rows(states_j); i++) { if (states_j[i, 1] <= a_curt[j] & states_j[i, 2] == sum(n_curt[1::j])) { outcome[i] = 1 } else if (states_j[i, 2] >= r_curt[j] & states_j[i, 2] == sum(n_curt[1::j])) { outcome[i] = 2 } } states_j = (states_j, outcome) cont = (min(states_j[mm_which(states_j[, 1] :> a_curt[j]), 1]), max(states_j[mm_which(states_j[, 1] :< r_curt[j]), 1])) states = (states \ states_j) } states_cont = states[mm_which(states[, 3] :== 0), 1::2] states_acc = states[mm_which(states[, 3] :== 1), 1::2] states_rej = states[mm_which(states[, 3] :== 2), 1::2] st_matrix("states_cont", states_cont) st_matrix("states_acc", states_acc) st_matrix("states_rej", states_rej) stata("qui svmat states_cont") stata("qui svmat states_acc") stata("qui svmat states_rej") stata(`"twoway (scatter states_cont1 states_cont2, msymbol(oh) color(gs13)) (scatter states_acc1 states_acc2, msymbol(X) color(red)) (scatter states_rej1 states_rej2, msymbol(+) color(green)), xtitle({it:s}) ytitle({it:m}) legend(lab(1 "Continue") lab(2 "Do not reject {it:H}{subscript:0}") lab(3 "Reject {it:H}{subscript:0}"))"'+ xopt) } } else { des = . opchar = J(1, 9, .) cols(feasible) feasible = J(1, (J == 1 ? 5 : (J == 2 ? 17 : 22)), .) J_curt = a_curt = n_curt = r_curt = . } ///// Return ///////////////////////////////////////////////////////////////// if (summary == 1) { printf("...outputting") } st_matrix("des", des) st_matrix("opchar", opchar) st_matrix("J", J) st_matrix("n", n) st_matrix("a", a) st_matrix("r", r) st_matrix("J_curt", J_curt) st_matrix("n_curt", n_curt) st_matrix("a_curt", a_curt) st_matrix("r_curt", r_curt) st_matrix("feasible", feasible) st_matrix("pi0", pi0) st_matrix("pi1", pi1) st_matrix("alpha", alpha) st_matrix("beta", beta) st_matrix("thetaF", thetaF) st_matrix("thetaE", thetaE) st_matrix("Nmin", Nmin) st_matrix("Nmax", Nmax) st_matrix("futility", futility) st_matrix("efficacy", efficacy) st_local("optimality", optimality) st_matrix("equal_n", equal_n) st_matrix("ensign", ensign) st_matrix("summary", summary) } // Function to determine feasible group sequential designs real matrix IntDesGS(real scalar J, real scalar pi0, real scalar pi1, real scalar alpha, real scalar beta, real scalar Nmin, real scalar Nmax, real scalar futility, real scalar efficacy, real scalar equal_n, real scalar ensign, real scalar summary) { dbinomial_pi0 = J(Nmax, Nmax - 1, 0) dbinomial_pi1 = J(Nmax, Nmax - 1, 0) for (i = 1; i <= Nmax - 1; i++) { dbinomial_pi0[1::(i + 1), i] = binomialp(i, 0::i, pi0) dbinomial_pi1[1::(i + 1), i] = binomialp(i, 0::i, pi1) } A_pi0 = A_pi1 = R_pi0 = R_pi1 = J(1, J, 0) C_pi0 = C_pi1 = J(Nmax + 1, J, 0) feasible_designs = J(500000, (J == 2 ? 15 : 20), 0) counter = 1 counter2 = 0 if (J == 2) { for (n1 = (futility + efficacy == 2 ? 2 : 1); n1 <= (equal_n == 1 ? floor(Nmax/2) : (Nmax - 1)); n1++) { C_pi0[, 1] = C_pi1[, 1] = J(Nmax + 1, 1, 0) C_pi0[1::(n1 + 1), 1] = dbinomial_pi0[1::(n1 + 1), n1] C_pi1[1::(n1 + 1), 1] = dbinomial_pi1[1::(n1 + 1), n1] for (a1 = (futility == 1 ? 0 : -1); a1 <= (futility == 1 ? (ensign == 1 ? 0 : (efficacy == 1 ? n1 - 2 : n1 - 1)) : -1); a1++) { if (a1 >= 0) { A_pi0[1] = sum(C_pi0[1::(a1 + 1), 1]) A_pi1[1] = sum(C_pi1[1::(a1 + 1), 1]) } else { A_pi0[1] = A_pi1[1] = 0 } if (A_pi1[1] > beta) { break } for (r1 = (efficacy == 1 ? n1 : n1 + 1); r1 >= (efficacy == 1 ? a1 + 2 : n1 + 1); r1--) { if (r1 <= n1) { R_pi0[1] = sum(C_pi0[(r1 + 1)::(n1 + 1), 1]) R_pi1[1] = sum(C_pi1[(r1 + 1)::(n1 + 1), 1]) } else { R_pi0[1] = R_pi1[1] = 0 } if (R_pi0[1] > alpha) { break } for (n2 = (equal_n == 1 ? n1 : max((1, Nmin - n1))); n2 <= (equal_n == 1 ? n1 : Nmax - n1); n2++) { C_pi0[, 2] = J(Nmax + 1, 1, 0) C_pi1[, 2] = J(Nmax + 1, 1, 0) for (s = max((a1 + 1, 0)); s <= min((r1 - 1, n1)); s++) { C_pi0[(s + 1)::(s + 1 + n2), 2]= C_pi0[(s + 1)::(s + 1 + n2), 2] + C_pi0[s + 1, 1]*dbinomial_pi0[1::(n2 + 1), n2] C_pi1[(s + 1)::(s + 1 + n2), 2]= C_pi1[(s + 1)::(s + 1 + n2), 2] + C_pi1[s + 1, 1]*dbinomial_pi1[1::(n2 + 1), n2] } for (a2 = max((0, a1 + 1)); a2 <= min((r1 + n2 - 2, n1 + n2 - 1)); a2++) { R_pi1[2] = sum(C_pi1[(a2 + 2)::(n1 + n2 + 1), 2]) counter2++ if (allof((summary == 1, mod(counter2, 10000) == 0), 1)) { printf("{txt}...over {res}%g{txt} designs evaluated...\n", counter2) } if (sum(R_pi1) < 1 - beta) { break } R_pi0[2] = sum(C_pi0[(a2 + 2)::(n1 + n2 + 1), 2]) if (sum(R_pi0) <= alpha) { A_pi0[2] = sum(C_pi0[1::(a2 + 1), 2]) A_pi1[2] = sum(C_pi1[1::(a2 + 1), 2]) PET1_pi0 = A_pi0[1] + R_pi0[1] PET1_pi1 = A_pi1[1] + R_pi1[1] ESS_pi0 = n1 + (1 - PET1_pi0)*n2 ESS_pi1 = n1 + (1 - PET1_pi1)*n2 Med_pi0 = (PET1_pi0 < 0.5 ? n1 + n2 : (PET1_pi0 > 0.5 ? n1 : n1 + 0.5*n2)) Med_pi1 = (PET1_pi1 < 0.5 ? n1 + n2 : (PET1_pi1 > 0.5 ? n1 : n1 + 0.5*n2)) Var_pi0 = PET1_pi0*n1^2 + (1 - PET1_pi0)*(n1 + n2)^2 - ESS_pi0^2 Var_pi1 = PET1_pi1*n1^2 + (1 - PET1_pi1)*(n1 + n2)^2 - ESS_pi1^2 feasible_designs[counter, ] = (n1, n2, a1, a2, r1, sum(R_pi0), sum(R_pi1), PET1_pi0, PET1_pi1, ESS_pi0, ESS_pi1, Med_pi0, Med_pi1, Var_pi0, Var_pi1) counter++ } } } } } } } else { for (n1 = (futility + efficacy == 2 ? 2 : 1); n1 <= (equal_n == 1 ? floor(Nmax/3) : (Nmax - 3)); n1++) { C_pi0[, 1] = C_pi1[, 1] = J(Nmax + 1, 1, 0) C_pi0[1::(n1 + 1), 1] = dbinomial_pi0[1::(n1 + 1), n1] C_pi1[1::(n1 + 1), 1] = dbinomial_pi1[1::(n1 + 1), n1] for (a1 = (futility == 1 ? 0 : -1); a1 <= (futility == 1 ? (ensign == 1 ? 0 : (efficacy == 1 ? n1 - 2 : n1 - 1)) : -1); a1++) { if (a1 >= 0) { A_pi0[1] = sum(C_pi0[1::(a1 + 1), 1]) A_pi1[1] = sum(C_pi1[1::(a1 + 1), 1]) } else { A_pi0[1] = A_pi1[1] = 0 } if (A_pi1[1] > beta) { break } for (r1 = (efficacy == 1 ? n1 : n1 + 1); r1 >= (efficacy == 1 ? a1 + 2 : n1 + 1); r1--) { if (r1 <= n1) { R_pi0[1] = sum(C_pi0[(r1 + 1)::(n1 + 1), 1]) R_pi1[1] = sum(C_pi1[(r1 + 1)::(n1 + 1), 1]) } else { R_pi0[1] = R_pi1[1] = 0 } if (R_pi0[1] > alpha) { break } for (n2 = (equal_n == 1 ? n1 : 2); n2 <= (equal_n == 1 ? n1 : Nmax - n1 - 1); n2++) { C_pi0[, 2] = J(Nmax + 1, 1, 0) C_pi1[, 2] = J(Nmax + 1, 1, 0) for (s = max((a1 + 1, 0)); s <= min((r1 - 1, n1)); s++) { C_pi0[(s + 1)::(s + 1 + n2), 2]= C_pi0[(s + 1)::(s + 1 + n2), 2] + C_pi0[s + 1, 1]*dbinomial_pi0[1::(n2 + 1), n2] C_pi1[(s + 1)::(s + 1 + n2), 2]= C_pi1[(s + 1)::(s + 1 + n2), 2] + C_pi1[s + 1, 1]*dbinomial_pi1[1::(n2 + 1), n2] } for (a2 = (futility == 1 ? max((0, a1 + 1)) : -1); a2 <= (futility == 1 ? min((r1 + n2 - 2, n1 + n2 - 1)) : -1); a2++) { if (a2 >= 0) { A_pi0[2] = sum(C_pi0[1::(a2 + 1), 2]) A_pi1[2] = sum(C_pi1[1::(a2 + 1), 2]) } else { A_pi0[2] = A_pi1[2] = 0 } if (A_pi1[1] + A_pi1[2] > beta){ break } for (r2 = (efficacy == 1 ? r1 + n2 - 1 : n1 + n2 + 1); r2 >= (efficacy == 1 ? max((a2 + 2, r1)) : n1 + n2 + 1); r2--) { if (r2 <= n1 + n2) { R_pi0[2] = sum(C_pi0[(r2 + 1)::(n1 + n2 + 1), 1]) R_pi1[2] = sum(C_pi1[(r2 + 1)::(n1 + n2 + 1), 1]) } else { R_pi0[2] = R_pi1[2] = 0 } if (R_pi0[1] + R_pi0[2] > alpha){ break } for (n3 = (equal_n == 1 ? n1 : max((1, Nmin - n1 - n2))); n3 <= (equal_n == 1 ? n1 : Nmax - n1 - n2); n3++) { for (s = max((a2 + 1, 0)); s <= min((r2 - 1, n1 + n2)); s++) { C_pi0[(s + 1)::(s + 1 + n3), 3]= C_pi0[(s + 1)::(s + 1 + n3), 3] + C_pi0[s + 1, 2]*dbinomial_pi0[1::(n3 + 1), n3] C_pi1[(s + 1)::(s + 1 + n3), 3]= C_pi1[(s + 1)::(s + 1 + n3), 3] + C_pi1[s + 1, 2]*dbinomial_pi1[1::(n3 + 1), n3] } for (a3 = max((0, a2 + 1)); a3 <= min((r2 + n3 - 2, n1 + n2 + n3 - 1)); a3++) { R_pi1[3] = sum(C_pi1[(a3 + 2)::(n1 + n2 + n3 + 1), 1]) counter2++ if (allof((summary == 1, mod(counter2, 10000) == 0), 1)) { printf("{txt}...over {res}%g{txt} designs evaluated...\n", counter2) } if (sum(R_pi1) < 1 - beta) { break } R_pi0[3] = sum(C_pi0[(a3 + 2)::(n1 + n2 + n3 + 1), 1]) if (sum(R_pi0) <= alpha) { A_pi0[3] = sum(C_pi0[1::(a3 + 1), 3]) A_pi1[3] = sum(C_pi1[1::(a3 + 1), 3]) PET1_pi0 = A_pi0[1] + R_pi0[1] PET1_pi1 = A_pi1[1] + R_pi1[1] PET2_pi0 = A_pi0[2] + R_pi0[2] PET2_pi1 = A_pi1[2] + R_pi1[2] ESS_pi0 = PET1_pi0*n1 + PET2_pi0*(n1 + n2) + (1 - PET1_pi0 - PET2_pi0)*(n1 + n2 + n3) ESS_pi1 = PET1_pi1*n1 + PET2_pi1*(n1 + n2) + (1 - PET1_pi1 - PET2_pi1)*(n1 + n2 + n3) if (PET1_pi0 == 0.5) { Med_pi0 = n1 + 0.5*n2 } else if (PET1_pi0 + PET2_pi0 == 0.5){ Med_pi0 = n1 + n2 + 0.5*n3 } else { if (PET1_pi0 > 0.5){ Med_pi0 = n1 } else if ((PET1_pi0 < 0.5) & (PET1_pi0 + PET2_pi0 > 0.5)) { Med_pi0 = n1 + n2 } else { Med_pi0 = n1 + n2 + n3 } } if (PET1_pi1 == 0.5) { Med_pi1 = n1 + 0.5*n2 } else if (PET1_pi1 + PET2_pi1 == 0.5) { Med_pi1 = n1 + n2 + 0.5*n3 } else { if (PET1_pi1 > 0.5) { Med_pi1 = n1 } else if ((PET1_pi1 < 0.5) & (PET1_pi1 + PET2_pi1 > 0.5)) { Med_pi1 = n1 + n2 } else { Med_pi1 = n1 + n2 + n3 } } Var_pi0 = PET1_pi0*n1^2 + PET2_pi0*(n1 + n2)^2 + (1 - PET1_pi0 - PET2_pi0)*(n1 + n2 + n3)^2 - ESS_pi0^2 Var_pi1 = PET1_pi1*n1^2 + PET2_pi1*(n1 + n2)^2 + (1 - PET1_pi1 - PET2_pi1)*(n1 + n2 + n3)^2 - ESS_pi1^2 feasible_designs[counter, ] = (n1, n2, n3, a1, a2, a3, r1, r2, sum(R_pi0), sum(R_pi1), PET1_pi0, PET1_pi1, PET2_pi0, PET2_pi1, ESS_pi0, ESS_pi1, Med_pi0, Med_pi1, Var_pi0, Var_pi1) counter++ } } } } } } } } } } return(feasible_designs[1::(counter > 1 ? counter - 1 : 1), ]) } // Function to return conditional rejection probability in a group sequential // design real scalar cr_gs(real scalar pi, real scalar s, real scalar m, real scalar J, real rowvector a, real rowvector r, real rowvector n) { if (J > 1) { a[mm_which(a :== .)] = J(1, length(mm_which(a :== .)), -1) r[mm_which(r :== .)] = J(1, length(mm_which(r :== .)), sum(n) + 1) } if (J == 1) { cr = binomialtail(n - m, r - s, pi) } else { if (m <= n[1]) { if (s >= r[1]) { cr = 1 } else { poss_s = permutations(uniqrows(s :+ 0::(n[1] - m) \ 0::n[2]), 2) poss_s = select(poss_s, poss_s[, 1] :>= s :& poss_s[, 1] :<= s + n[1] - m :& poss_s[, 2] :>= 0 :& poss_s[, 2] :<= n[2]) poss_s = select(poss_s, poss_s[, 1] :>= r[1] :| (poss_s[, 1] :>= max((0, a[1] + 1)) :& poss_s[, 1] :<= min((n[1], r[1] - 1)) :& (poss_s[, 1] + poss_s[, 2] :>= r[2]))) if (rows(poss_s) > 0) { if (m < n[1]) { prob_poss_s = (binomialp(n[1] - m, poss_s[, 1] :- s, pi), binomialp(n[2], poss_s[, 2], pi)) } else { prob_poss_s = (J(rows(poss_s), 1, 1), binomialp(n[2], poss_s[, 2], pi)) } cr = sum(prob_poss_s[, 1]:*prob_poss_s[, 2]) } else { cr = 0 } } } else { cr = binomialtail(sum(n) - m, r[2] - s, pi) } } return(cr) } // Function for determining pmf of group sequential design real matrix pmf_curtailed(real scalar pi, real scalar J, real scalar J_curt, real rowvector a_curt, real rowvector r_curt, real rowvector n, real rowvector n_curt, real rowvector k,| real colvector dbinom_pi) { a_curt[mm_which(a_curt :== .)] = J(1, length(mm_which(a_curt :== .)), -1) r_curt[mm_which(r_curt :== .)] = J(1, length(mm_which(r_curt :== .)), sum(n) + 1) if (args() < 9) { dbinom_pi = J(max(n_curt) + 1, J_curt, 0) for (j = 1; j <= J_curt; j++) { dbinom_pi[1::(n_curt[j] + 1), j] = binomialp(n_curt[j], 0::n_curt[j], pi) } } pmf_mat = J(sum(n) + 1, J_curt, 0) pmf_mat[1::(n_curt[1] + 1), 1] = dbinom_pi[1::(n_curt[1] + 1), 1] cont = (max((0, a_curt[1] + 1)), min((r_curt[1] - 1, n_curt[1]))) for (j = 1; j <= J_curt - 1; j++) { for (i = cont[1]; i <= cont[2]; i++) { pmf_mat[(i + 1)::(i + 1 + n_curt[j + 1]), j + 1] = pmf_mat[(i + 1)::(i + 1 + n_curt[j + 1]), j + 1] :+ pmf_mat[i + 1, j]*dbinom_pi[1::(n_curt[j + 1] + 1), j + 1] } pmf_mat[(cont[1] + 1)::(cont[2] + 1), j] = J(cont[2] - cont[1] + 1, 1, 0) upd = cont[1]::(cont[2] + n_curt[j + 1]) cont = (min(select(upd, upd :> a_curt[j + 1])), max(select(upd, upd :< r_curt[j + 1]))) } N = (0, runningsum(n)) full_k = 0 for (i = 1; i <= length(k); i++) { full_k = (full_k \ (N[k[i]] + 1)::N[k[i] + 1]) } full_k = full_k[2::length(full_k)] if (length(full_k) < J_curt) { for (stage = 1; stage <= J_curt; stage++) { if (!any(stage :== full_k)) { pmf_mat[, stage] = J(sum(n) + 1, 1, 0) } } pmf_mat = pmf_mat/sum(pmf_mat) } terminal = terminal_states_gs(J_curt, a_curt, r_curt, n_curt, full_k') if (sum(pmf_mat :> 0) > 1) { true_k = J(rows(terminal), 1, 0) for (i = 1; i <= rows(terminal); i++) { true_k[i] = min(mm_which(N :>= terminal[i, 2])) - 1 } f = . for (j = 1; j <= J_curt; j++) { f = (f \ select(pmf_mat[, j], pmf_mat[, j] :> 0)) } pmf = (J(rows(terminal), 1, pi), terminal[, 1], terminal[, 2], true_k, f[2::length(f)]) } else { true_k = J(rows(terminal), 1, 0) for (i = 1; i <= rows(terminal); i++) { true_k[i] = min(mm_which(N :>= terminal[i, 3])) - 1 } pmf = (J(rows(terminal), 1, pi), terminal[, 1], terminal[, 2], terminal[, 3], J(rows(terminal), 1, 0)) non_zero = (mm_which(rowsum(pmf_mat) :> 0) - 1, mm_which(colsum(pmf_mat :> 0))) pmf[mm_which(pmf[, 2] :== non_zero[1] :& pmf[, 4] :== non_zero[2]), 5] = 1 } return(pmf) } // Function for determining operating characteristics of group sequential design real rowvector int_opchar_curtailed(real scalar pi, real scalar J, real scalar J_curt, real rowvector a_curt, real rowvector r_curt, real rowvector n, real rowvector n_curt, real rowvector N, real rowvector N_curt, real rowvector k,| real matrix pmf_pi) { a_curt[mm_which(a_curt :== .)] = J(1, length(mm_which(a_curt :== .)), -1) r_curt[mm_which(r_curt :== .)] = J(1, length(mm_which(r_curt :== .)), sum(n) + 1) if (args() < 11) { pmf_pi = pmf_curtailed(pi, J, J_curt, a_curt, r_curt, n, n_curt, k) } A = J(1, J_curt, 0) R = J(1, J_curt, 0) for (j = 1; j <= J_curt; j++) { A[j] = sum(select(pmf_pi, pmf_pi[, 2] :<= a_curt[j] :& pmf_pi[, 3] :== j)[, 5]) R[j] = sum(select(pmf_pi, pmf_pi[, 2] :>= r_curt[j] :& pmf_pi[, 3] :== j)[, 5]) } S = A + R cum_S = runningsum(S) N = (0, N) Atilde = J(1, J, 0) Rtilde = J(1, J, 0) for (j = 1; j <= J; j++) { Atilde[j] = sum(A[(1 + N[j])::N[j + 1]]) Rtilde[j] = sum(R[(1 + N[j])::N[j + 1]]) } N = N[2::length(N)] Stilde = Atilde + Rtilde cum_Stilde = runningsum(Stilde) if (any(cum_S :== 0.5)) { Med = 0.5*(N_curt[mm_which(cum_S :== 0.5)] + N_curt[mm_which(cum_S :== 0.5) + 1]) } else { Med = N_curt[mm_which(cum_S :> 0.5)[1]] } return((pi, sum(R), sum(N_curt:*S), sum(N_curt:^2:*S) - sum(N_curt:*S)^2, Med, A, R, S, cum_S, Atilde, Rtilde, Stilde, cum_Stilde)) } // Function to determine all permutations of length q from a vector vec with // replacement real matrix permutations(real colvector vec, real scalar q) { if (q == 1) { return(vec) } lvec = length(vec) if (lvec == 1) { return(J(1, q, vec)) } l = lvec^q allperms = J(l, q, 0) d = vec[2::lvec] :- vec[1::(lvec - 1)] ld = length(d) vl = (-sum(d) \ d) tmp = J(1, l/lvec, vl) allperms[1::l, q] = vec(tmp) elts = 1 count = 1 while (elts[count] + lvec^(q - 1) <= l) { elts = (elts \ elts[count] + lvec^(q - 1)) count = count + 1 } allperms[elts, 1] = vl if (q > 2) { for (i = 2; i <= q - 1; i++) { elts = 1 count = 1 while (elts[count] + lvec^(i - 1) <= l) { elts = (elts \ elts[count] + lvec^(i - 1)) count = count + 1 } tmp = J(1, length(elts)/(ld + 1), vl) allperms[elts, q - i + 1] = vec(tmp) } } allperms[1, 1::q] = J(1, q, vec[1]) origallperms = allperms for (i = 2; i <= l; i++) { for (j = 1; j <= q; j++) { allperms[i, j] = sum(origallperms[1::i, j]) } } return(allperms) } // Function to determine terminal states in a group sequential design real matrix terminal_states_gs(real scalar J, real rowvector a, real rowvector r, real rowvector n, real rowvector k) { a[mm_which(a :== .)] = J(1, length(mm_which(a :== .)), -1) r[mm_which(r :== .)] = J(1, length(mm_which(r :== .)), sum(n) + 1) terminal = J(1, 3, 0) if (a[1] >= 0) { s1 = 0::a[1] terminal = (terminal \ (s1, J(a[1] + 1, 1, n[1]), J(a[1] + 1, 1, 1))) } if (r[1] <= n[1]) { s1 = r[1]::n[1] terminal = (terminal \ (s1, J(n[1] - r[1] + 1, 1, n[1]), J(n[1] - r[1] + 1, 1, 1))) } cont = (max((0, a[1] + 1)), min((r[1] - 1, n[1]))) if (J >= 3) { for (j = 1; j <= J - 2; j++) { num_rows = cont[2] + n[j + 1] - cont[1] + 1 vals = cont[1]::(cont[2] + n[j + 1]) upd = (vals, J(num_rows, 1, sum(n[1::(j + 1)])), J(num_rows, 1, j + 1)) terminal = (terminal \ select(upd, upd[, 1] :<= a[j + 1] :| upd[, 1] :>= r[j + 1])) cont = (min(select(upd, upd[, 1] :> a[j + 1])[, 1]), max(select(upd, upd[, 1] :< r[j + 1])[, 1])) } } num_rows = cont[2] + n[J] - cont[1] + 1 vals = cont[1]::(cont[2] + n[J]) terminal = (terminal \ (vals, J(num_rows, 1, sum(n)), J(num_rows, 1, J))) keep_rows = mm_which(terminal[, 3] :== k[1]) if (length(k) > 1) { for (els = 2; els <= length(k); els++) { keep_rows = (keep_rows \ mm_which(terminal[, 3] :== k[els])) } } terminal = terminal[keep_rows, ] return(terminal) } end