### Install RARtrials from GitHub Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/README.md Install the RARtrials package from GitHub using devtools. The commented line is a prerequisite for installing devtools if it is not already installed. ```r #install.packages('devtools') devtools::install_github("yayayaoyaoyao/RARtrials") ``` -------------------------------- ### pmax_beta example with 5 arms and upper side Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/pmax_beta.Rd Example call for a trial with 5 arms, specifying Beta parameters for each arm and side='upper' to compute the probability that arm 1 is the best (maximal). ```R pmax_beta(armn=5,a1=8,b1=10,a2=5,b2=19,a3=8,b3=21, a4=6, b4=35, a5=15, b5=4, side='upper') ``` -------------------------------- ### pmax_beta example with 4 arms and lower side Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/pmax_beta.Rd Example call for a trial with 4 arms, specifying Beta parameters for each arm and side='lower' to compute the probability that arm 1 is the worst (minimal). ```R pmax_beta(armn=4,a1=56,b1=98,a2=25,b2=70,a3=87,b3=107, a4=106, b4=202, side='lower') ``` -------------------------------- ### pmax_beta example with 3 arms and upper side Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/pmax_beta.Rd Example call for a trial with 3 arms, specifying Beta parameters for each arm and side='upper' to compute the probability that arm 1 is the best (maximal). ```R pmax_beta(armn=3,a1=60,b1=46,a2=55,b2=46,a3=35,b3=36,side='upper') ``` -------------------------------- ### Install RARtrials from CRAN Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/README.md Install the RARtrials package from CRAN. ```r install.packages('RARtrials') ``` -------------------------------- ### pmax_normal examples for 3, 4, and 5 arms Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/pmax_normal.Rd Three example calls to pmax_normal with different arm counts (5, 4, 3) and side values ('upper', 'lower', 'upper'). Each call specifies the required mean and sd arguments for the given number of arms. ```r pmax_normal(armn=5,mean1=0.8,sd1=0.2,mean2=0.5,sd2=0.1,mean3=0.8, sd3=0.5,mean4=0.6,sd4=0.2,mean5=0.6,sd5=0.2,side='upper') ``` ```r pmax_normal(armn=4,mean1=8,sd1=2,mean2=8.5,sd2=2,mean3=8.3, sd3=1.8,mean4=8.7,sd4=2,side='lower') ``` ```r pmax_normal(armn=3,mean1=80,sd1=20,mean2=50,sd2=10,mean3=80, sd3=15,side='upper') ``` -------------------------------- ### Integration-based posterior probability calculation with pgreater_beta and pmax_beta Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/vignettes/RARtrials.Rmd Uses direct integration functions pgreater_beta and pmax_beta to compute the same posterior probabilities as the simulation example, but faster. Requires the same alpha and beta vectors; delta=0.1 with side='upper' for larger effects, delta=-0.1 with side='lower' for smaller effects. Results are printed as comments. ```r #Example: Sample code using Integration #Expect the treatment group to have a larger treatment effect compared to the control group #Calculate results of Pr(p_k>p_{control}+\delta|data_{t-1}) pgreater_beta(a1=alpha[1],b1=beta[1],a2=alpha[2], b2=beta[2],delta=0.1,side='upper') pgreater_beta(a1=alpha[1],b1=beta[1],a2=alpha[3], b2=beta[3],delta=0.1,side='upper') #Return results: 0.7951487 0.3477606 #Calculate results of Pr(p_k=max\{p_1,...,p_K\}|data_{t-1}) pmax_beta(armn=3,a2=alpha[3],b2=beta[3],a3=alpha[2],b3=beta[2],a1=alpha[1], b1=beta[1],side='upper') pmax_beta(armn=3,a2=alpha[1],b2=beta[1],a3=alpha[3],b3=beta[3],a1=alpha[2], b1=beta[2],side='upper') pmax_beta(armn=3,a2=alpha[1],b2=beta[1],a3=alpha[2],b3=beta[2],a1=alpha[3], b1=beta[3],side='upper') #Return results: 0.01796526 0.8788907 0.1031441 #Expect the treatment group to have a smaller treatment effect compared to the control group #Calculate results of Pr(p_{control}>p_k+\delta|data_{t-1}) pgreater_beta(a1=alpha[1],b1=beta[1],a2=alpha[2],b2=beta[2],delta=-0.1,side='lower') pgreater_beta(a1=alpha[1],b1=beta[1],a2=alpha[3],b2=beta[3],delta=-0.1,side='lower') #Return results: 0.001093548 0.03348547 #Calculate results of Pr(p_k=min\{p_1,...,p_K\}|data_{t-1}) pmax_beta(armn=3,a2=alpha[3],b2=beta[3],a3=alpha[2],b3=beta[2],a1=alpha[1], b1=beta[1],side='lower') pmax_beta(armn=3,a2=alpha[1],b2=beta[1],a3=alpha[3],b3=beta[3],a1=alpha[2], b1=beta[2],side='lower') pmax_beta(armn=3,a2=alpha[1],b2=beta[1],a3=alpha[2],b3=beta[2],a1=alpha[3], b1=beta[3],side='lower') #Return results: 0.7560864 0.01230027 0.2316133 ``` -------------------------------- ### Simulate A-Optimal Allocation Trial with sim_A_optimal_known_var Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/sim_A_optimal_known_var.Rd Calls sim_A_optimal_known_var with patient accrual, outcome delay, enrollment rate, sample sizes, arm count, means, standard deviations, alpha, arm labels, and side. The example shows the required argument order and typical values for a two-arm trial with known variances. ```R sim_A_optimal_known_var( Pats = 10, nMax = 100, TimeToOutcome = expression(distribution(sample = 1, rate = 0.5)), enrollrate = 1, N2 = 50, armn = 2, mean = c(0, 0), sd = c(1, 1), alphaa = 0.025, armlabel = c("Control", "Treatment"), side = 1 ) ``` -------------------------------- ### sim_dabcd_max_power Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/sim_dabcd_max_power.Rd Documentation for the sim_dabcd_max_power function, including usage, arguments, return value, and examples. ```APIDOC ## sim_dabcd_max_power ### Description Simulates a trial for doubly adaptive biased coin design with maximal power strategy targeting generalized Neyman allocation (2-5 arms) or generalized RSIHR allocation (2-3 arms) for binary outcomes. The function returns an object of class "dabcd" containing the final decision, Z test statistics, the simulated data set, and participants accrued for each arm. ### Usage ```r sim_dabcd_max_power(Pats, nMax, TimeToOutcome, enrollrate, N1, N2, armn, armlabel, h, BB, type, side) ``` ### Arguments - **Pats** (numeric) - Number of patients accrued within a certain time frame. - **nMax** (numeric) - Assumed maximum accrued number of patients with the disease in the population. - **TimeToOutcome** (expression) - Distribution of delayed response times or a fixed delay time for responses. - **enrollrate** (numeric) - Patient enrollment rate. - **N1** (numeric) - First stage sample size. - **N2** (numeric) - Second stage sample size. - **armn** (numeric) - Number of arms (2-5 for Neyman, 2-3 for RSIHR). - **armlabel** (vector) - Labels for the arms. - **h** (numeric vector) - Allocation probabilities or related parameters. - **BB** (numeric) - Biased coin parameter. - **type** (character) - Allocation type: 'Neyman' or 'RSIHR'. - **side** (character) - Direction of one-sided test: 'upper' or 'lower'. ### Value Returns an object of class "dabcd", a list containing: - Final decision based on Z test statistics (1 = selected, 0 = not selected). - Z test statistics. - Simulated data set with columns: participant ID, enrollment time, observed time of results, allocated arm, and participant's result. - Participants accrued for each arm at termination time. ### Details - For more than two arms, one-sided nominal level is alphaa divided by arm*(arm-1)/2 (Bonferroni correction). - Methods are not smoothed and based on Hu & Zhang's formula. - Delay mechanism parameters adapted from IntroBayesianSimulation. ### Examples ```r sim_dabcd_max_power(Pats=10,nMax=50000,TimeToOutcome=expression(rnorm(length( vStartTime ),30, 3)),enrollrate=0.9,N1=30,N2=300,armn=3,armlabel=c(1,2,3),h=c(0.2,0.3,0.2),BB=0.1,type='Neyman',side='upper') sim_dabcd_max_power(Pats=10,nMax=50000,TimeToOutcome=expression(rnorm(length( vStartTime ),60, 3)),enrollrate=0.1,N1=50,N2=500,armn=3,armlabel=c(1,2,3),h=c(0.2,0.3,0.3),BB=0.15,type='RSIHR',side='upper') ``` ``` -------------------------------- ### sim_RSIHR_optimal_known_var Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/sim_RSIHR_optimal_known_var.Rd Documentation for the sim_RSIHR_optimal_known_var function, including usage, arguments, return value, details, examples, and references. ```APIDOC ## sim_RSIHR_optimal_known_var ### Description Simulates a trial for continuous endpoints with known variances and fixed allocation probabilities, aiming to minimize the criteria \(\sum_{i=1}^{K}n_i\Psi_i\) with constraints \(\frac{\sigma_1^2}{n_1}+\frac{\sigma_k^2}{n_k}\leq C\). ### Usage ```r sim_RSIHR_optimal_known_var(Pats, nMax, TimeToOutcome, enrollrate, N1, N2, armn, mean, sd, alphaa, armlabel, cc, side) ``` ### Arguments - **Pats** (numeric) - Number of patients accrued within a certain time frame. - **nMax** (numeric) - Assumed maximum accrued number of patients with the disease in the population. - **TimeToOutcome** (expression) - Distribution of delayed response times or a fixed delay time for responses. - **enrollrate** (numeric) - Enrollment rate. - **N1** (numeric) - First interim sample size. - **N2** (numeric) - Second interim sample size. - **armn** (numeric) - Number of arms. - **mean** (numeric vector) - Mean for each arm. - **sd** (numeric vector) - Standard deviation for each arm. - **alphaa** (numeric) - One-sided nominal level for each test (with Bonferroni correction for more than two arms). - **armlabel** (numeric vector) - Labels for each arm. - **cc** (numeric) - Fixed constant C in the constraint. - **side** (character) - Direction of a one-sided test, with values 'upper' or 'lower'. ### Value Returns an object of class "RSIHRoptimal", a list containing: - Final decision based on the Z test statistics (1 for selected, 0 for not selected). - Z test statistics. - The simulated data set and participants accrued for each arm at the time of termination of that group in one trial. The simulated data set includes 5 columns: participant ID number, enrollment time, observed time of results, allocated arm, and participants' result. ### Details This function aims to minimize the criteria \(\sum_{i=1}^{K}n_i\Psi_i\) with constraints \(\frac{\sigma_1^2}{n_1}+\frac{\sigma_k^2}{n_k}\leq C\), where \(k=2,...,K\) for some fixed C. It is equivalent to generalized RSIHR allocation for continuous endpoints with known variances. With more than two arms, the one-sided nominal level of each test is \(\alpha_a\) divided by \(arm \times (arm-1)/2\); a Bonferroni correction. Considering the delay mechanism, parameters \(Pats\), \(nMax\), and \(TimeToOutcome\) are adapted from \url{https://github.com/kwathen/IntroBayesianSimulation}. ### Examples ```r # Run the function with delayed responses follow a normal distribution with # a mean of 30 days and a standard deviation of 3 days under null hypothesis # in a three-armed trial sim_RSIHR_optimal_known_var(Pats=10,nMax=50000,TimeToOutcome=expression( rnorm(length( vStartTime ),30, 3)),enrollrate=0.9,N1=9,N2=132,armn=3, mean=c(9.1/100,9.1/100,9.1/100),sd=c(0.009,0.009,0.009),alphaa=0.025, armlabel=c(1,2,3),cc=mean(c(9.1/100,9.1/100,9.1/100)),side='lower') # Run the function with delayed responses follow a normal distribution with # a mean of 30 days and a standard deviation of 3 days under alternative hypothesis # in a three-armed trial sim_RSIHR_optimal_known_var(Pats=10,nMax=50000,TimeToOutcome=expression( rnorm(length( vStartTime ),30, 3)),enrollrate=0.9,N1=9,N2=132,armn=3, mean=c(9.1/100,8.47/100,8.47/100),sd=c(0.009,0.009,0.009),alphaa=0.025, armlabel=c(1,2,3),cc=mean(c(9.1/100,8.47/100,8.47/100)),side='lower') ``` ### References Biswas (2011) as cited in the RARtrials package. ``` -------------------------------- ### Select a_U by calling brar_select_au_binary 2000 times Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/vignettes/RARtrials.Rmd Runs the brar_select_au_binary function 2000 times with a fixed seed to generate the test statistic distribution under the intermediate hypothesis c(0.3,0.3). The example sets the alternative hypothesis to c(0.2,0.4) and uses the intermediate hypothesis to avoid excessively large or small power. ```R set.seed(12345) #Example: Select a_U by calling brar_au_binary 2000 times simnull3<-lapply(1:2000,function(x){ set.seed(x) brar_select_au_binary(Pats=10,nMax=50000,TimeToOutcome=0,enrollrate=0.9,N1=24,armn=2, h=c(0.3,0.3),N2=224,tp=1,armlabel=c(1, 2),blocksize=4,alpha1=1,beta1=1, alpha2=1,beta2=1,minstart=24,deltaa=-0.07,tpp=0,deltaa1=0.1,side='upper') }) ``` -------------------------------- ### Simulate cut-off and power for FLGI with known variances Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/vignettes/RARtrials.Rmd R code that first selects a critical cut-off value under the null hypothesis using flgi_cut_off_known_var(), then calculates power using sim_flgi_known_var(). The example uses continuous outcomes with known variances and the 'FLGI PM' rule. It recommends running more simulations for accuracy and notes the heavy computational burden. ```r #Example: The forward-looking Gittins Index rule applied in continuous outcomes with known #variances #Select cut-off values under null hypothesis for continuous outcomes with known variances. #The delayed responses follow a normal distribution with a mean of 30 days, sd of 3 days #and enrollment rate of 1. set.seed(123452) stopbound5<-lapply(1:1000,function(x){flgi_cut_off_known_var( Gittinstype='KV',df=0.995,Pats=10,nMax=50000,TimeToOutcome=0,enrollrate=1,armn=3, noRuns2=100,N2=120,block=8,rule='FLGI PM',prior_n=rep(1,3),prior_mean=rep(0,3), mean=c(-0.05,-0.05,-0.05),sd=c(0.346,0.346,0.346),side='upper')}) stopbound5a<-do.call(rbind,stopbound5) sum(stopbound5a[,1]>2.074 | stopbound5a[,2]> 2.074)/1000 #0.05 #The selected cut-off value is 2.074 with an overall upper one-sided type I error of 0.05. #It is recommended to run more simulations to obtain a more accurate cut-off value. #Calculate power based on 1000 simulations set.seed(123452) sim5<-lapply(1:1000,function(x){sim_flgi_known_var(Gittinstype='KV', df=0.995,Pats=10,nMax=50000,TimeToOutcome=0,enrollrate=1,armn=3, noRuns2=100,N2=120,block=8,rule='FLGI PM', prior_n=rep(1,3), prior_mean=rep(0,3), mean=c(-0.05,0.07,0.13),sd=c(0.346,0.346,0.346), stopbound =2.074,side='upper')}) h1decision<-t(sapply(sim5, "[[", 1)) sum(h1decision[,1]==1)/1000 #0.074 sum(h1decision[,2]==1)/1000 #0.244 sum(h1decision[,1]==1 | h1decision[,2]==1)/1000 #0.267 #Marginal power of rejecting H02 is 0.074 #Marginal power of rejecting H03 is 0.244 #Power of rejecting H02 or H03 is 0.267 ``` -------------------------------- ### Convert Normal-Inverse-Gamma parameters to Normal-Inverse-Chi-Squared Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/convert_gamma_to_chisq.Rd Converts a Normal-Inverse-Gamma parameter list (m, V, a, b) to a Normal-Inverse-Chi-Squared parameter list (mu, kappa, nu, sigsq). The example uses V=1/2, a=0.5, m=9.1/100, and b=0.00002. ```r convert_gamma_to_chisq(list(V=1/2,a=0.5,m=9.1/100,b=0.00002)) ``` -------------------------------- ### Convert Normal-Inverse-Chi-Squared parameters to Normal-Inverse-Gamma Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/convert_chisq_to_gamma.Rd Converts a Normal-Inverse-Chi-Squared parameter list (mu, kappa, nu, sigsq) to a Normal-Inverse-Gamma parameter list (m, V, a, b). The example uses mu=0.091, kappa=2, nu=1, and sigsq=4e-05. ```r convert_chisq_to_gamma(list(mu=0.091,kappa=2,nu=1,sigsq=4e-05)) ``` -------------------------------- ### Run flgi_cut_off_binary 20000 times to obtain cut-off value Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/flgi_cut_off_binary.Rd Simulates the forward-looking Gittins Index rule with delayed responses following a normal distribution (mean 60 days, sd 3 days). Run this command 20000 times to obtain a cut-off value around 1.9991 with an upper one-sided type I error of 0.025. The example is wrapped in \donttest{} and uses the 'FLGI PM' rule with unpooled z-statistics. ```r stopbound1<-lapply(1:20000,function(x){ flgi_cut_off_binary(Gittinstype='Binary',df=0.5,Pats=10, nMax=50000,TimeToOutcome=expression(rnorm( length( vStartTime ),60, 3)), enrollrate=0.9,I0= matrix(1,nrow=2,ncol=2),armn=2,N2=992,h=c(0.6,0.6),block=20, rule='FLGI PM',ztype='unpooled')}) stopbound1a<-do.call(rbind,stopbound1) sum(stopbound1a>(1.9991) )/20000 ``` -------------------------------- ### Simulate FLGI rule for continuous outcomes with known variances Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/vignettes/RARtrials.pdf Simulates the forward-looking Gittins Index rule for continuous outcomes with known variances. First selects the critical cut-off value under the null hypothesis using flgi_cut_off_known_var() with 1000 runs, then calculates power using sim_flgi_known_var() with the selected stopbound of 2.074. The example uses 3 arms, a normal delayed response distribution (mean 30 days, sd 3 days), and an enrollment rate of 1. The output shows an overall type I error of 0.05 and marginal powers of 0.074 and 0.244 for rejecting H02 and H03, respectively, with a combined power of 0.267. ```R #Example: The forward-looking Gittins Index rule applied in continuous outcomes with known #variances #Select cut-off values under null hypothesis for continuous outcomes with known variances. #The delayed responses follow a normal distribution with a mean of 30 days, sd of 3 days #and enrollment rate of 1. set.seed(123452) stopbound5<-lapply(1: 1000,function(x){flgi_cut_off_known_var( Gittinstype='KV',df=0.995,Pats=10,nMax=50000,TimeToOutcome=0,enrollrate=1,armn=3, noRuns2=100,N2=120,block=8,rule='FLGI PM',prior_n=rep(1,3),prior_mean=rep(0,3), mean=c(-0.05,-0.05,-0.05),sd=c(0.346,0.346,0.346),side='upper')}) stopbound5a<-do.call(rbind,stopbound5) sum(stopbound5a[,1]>2.074 | stopbound5a[,2]> 2.074)/1000 #0.05 #The selected cut-off value is 2.074 with an overall upper one-sided type I error of 0.05. #It is recommended to run more simulations to obtain a more accurate cut-off value. #Calculate power based on 1000 simulations set.seed(123452) sim5<-lapply(1: 1000,function(x){sim_flgi_known_var(Gittinstype='KV', df=0.995,Pats=10,nMax=50000,TimeToOutcome=0,enrollrate=1,armn=3, noRuns2=100,N2=120,block=8,rule='FLGI PM',prior_n=rep(1,3), prior_mean=rep(0,3),mean=c(-0.05,0.07,0.13),sd=c(0.346,0.346,0.346), stopbound =2.074,side='upper')}) h1decision<-t(sapply(sim5,"[[",1)) sum(h1decision[,1]==1)/1000 #0.074 sum(h1decision[,2]==1)/1000 #0.244 sum(h1decision[,1]==1 | h1decision[,2]==1)/1000 #0.267 #Marginal power of rejecting H02 is 0.074 #Marginal power of rejecting H03 is 0.244 #Power of rejecting H02 or H03 is 0.267 ``` -------------------------------- ### Simulate brar_select_au_unknown_var with delayed responses and evaluate stopping boundary Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/brar_select_au_unknown_var.Rd Runs 5 simulations of brar_select_au_unknown_var with delayed responses following a normal distribution (mean 30 days, sd 3 days) and identical prior distributions across arms. The example then computes futility statistics and evaluates the selected stopping boundary of 0.85, yielding an overall lower one-sided type I error of 0.2. It is recommended to run more simulations (e.g., 1000) for accurate deltaa and au. ```R set.seed(789) stopbound1<-lapply(1:5,function(x){ brar_select_au_unknown_var(Pats=10,nMax=50000,TimeToOutcome=expression(rnorm( length( vStartTime ),30, 3)), enrollrate=0.1, N1=48, armn=3, N2=480, tp=1, armlabel=c(1,2,3), blocksize=6, mean=c((9.1/100+8.74/100+8.74/100)/3, (9.1/100+8.74/100+8.74/100)/3,(9.1/100+8.74/100+8.74/100)/3) , sd=c(0.009,0.009,0.009),minstart=48, deltaa=c(-0.0003,-0.00035), tpp=0, deltaa1=c(0,0),V01=1/2,a01=0.3,m01=9/100,b01=0.00001,side='lower') }) simf<-list() simf1<-list() for (xx in 1:5){ if (any(stopbound1[[xx]][48:479,2]<0.01)){ simf[[xx]]<-NA } else{ simf[[xx]]<-stopbound1[[xx]][480,2] } if (any(stopbound1[[xx]][48:479,3]<0.01)){ simf1[[xx]]<-NA } else{ simf1[[xx]]<-stopbound1[[xx]][480,3] } } simf2<-do.call(rbind,simf) sum(is.na(simf2)) #1, achieve around 20\% futility simf3<-do.call(rbind,simf1) sum(is.na(simf3)) #1, achieve around 20\% futility stopbound1a<-cbind(simf2,simf3) stopbound1a[is.na(stopbound1a)] <- 0 sum(stopbound1a[,1]>0.85 | stopbound1a[,2]>0.85)/5 #0.2 ``` -------------------------------- ### Run sim_Aa_optimal_known_var with normal delayed responses Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/sim_Aa_optimal_known_var.Rd Executes the sim_Aa_optimal_known_var function for a two-armed trial with known variance. The TimeToOutcome argument uses rnorm to generate delayed responses with mean 30 and standard deviation 3 days. The example sets Pats=10, nMax=50000, enrollrate=0.9, N2=88, armn=2, means 9.1/100 and 8.47/100, sd=0.009, alphaa=0.025, arm labels 1 and 2, and side='lower'. ```R #Run the function with delayed responses follow a normal distribution with #a mean of 30 days and a standard deviation of 3 days under alternative hypothesis #in a two-armed trial sim_Aa_optimal_known_var(Pats=10,nMax=50000,TimeToOutcome=expression( rnorm(length( vStartTime ),30, 3)),enrollrate=0.9,N2=88,armn=2, mean=c(9.1/100,8.47/100),sd=c(0.009,0.009),alphaa=0.025,armlabel = c(1,2),side='lower') ``` -------------------------------- ### Simulate Bayesian response-adaptive trial with binary outcomes using sim_brar_binary Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/sim_brar_binary.Rd Simulates a two-arm Bayesian response-adaptive randomized trial with a control group for binary outcomes. This example uses delayed responses following a normal distribution with mean 30 and standard deviation 3 days, with h=c(0.2,0.4) and tp=0.5. The function returns an object of class 'brar' containing the final decision, test statistics, simulated data, and participants accrued per arm. ```R sim_brar_binary(Pats=10,nMax=50000,TimeToOutcome=expression(rnorm( length( vStartTime ),30, 3)), enrollrate=0.1,N1=24,armn=2,h=c(0.2,0.4),au=0.36,N2=224,tp=0.5,armlabel=c(1,2),blocksize=4, alpha1=1,beta1=1,alpha2=1,beta2=1,minstart=24,deltaa=-0.01,tpp=0,deltaa1=0.1,side='upper') ``` -------------------------------- ### sim_brar_unknown_var example with delayed normal responses Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/sim_brar_unknown_var.Rd Runs a three-arm trial with delayed responses following a normal distribution with mean 30 days and standard deviation 3 days under the null hypothesis. Uses mean=c(9.19/100,8.74/100,8.74/100), sd=c(0.009,0.009,0.009), tp=1, and minimal effect size 0. The side parameter is set to 'lower'. ```R sim_brar_unknown_var(Pats=10,nMax=50000,TimeToOutcome=expression(rnorm( length(vStartTime ),30,3)),enrollrate=0.1, N1=48,armn=3,au=c(0.85,0.85), N2=480,tp=1,armlabel=c(1, 2,3),blocksize=6,mean=c(9.19/100,8.74/100,8.74/100), sd=c(0.009,0.009,0.009), minstart=48,deltaa=c(-0.000325,-0.000325), tpp=0,deltaa1=c(0,0),V01=1/2,a01=0.3,m01=9/100,b01=0.00001,side='lower') ``` -------------------------------- ### Select a_U and estimate power with brar_select_au_unknown_var Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/vignettes/RARtrials.Rmd Selects the futility threshold a_U by calling brar_select_au_unknown_var 2000 times under the null, then uses the chosen a_U (0.9836) in sim_brar_unknown_var to estimate power. The first block runs the null simulations and extracts test statistics, ensuring about 1% futility stops and a type I error near 0.025. The second block runs 1000 power simulations and reports the proportion of 'Superiorityfinal' decisions (0.882). ```R set.seed(12345) simnull4<-lapply(1:2000,function(x){ brar_select_au_unknown_var(Pats=10,nMax=50000,TimeToOutcome=expression( rnorm(length( vStartTime ),30, 3)),enrollrate=0.1, N1=192,armn=3, N2=1920,tp=1,armlabel=c(1,2,3),blocksize=6, mean=c((9.1/100+8.92/100+8.92/100)/3,(9.1/100+8.92/100+8.92/100)/3, (9.1/100+8.92/100+8.92/100)/3),sd=c(0.009,0.009,0.009),minstart=192, deltaa=c(0.00077,0.00077),tpp=1,deltaa1=c(0,0),V01=1/2,a01=0.3,m01=9/100, b01=0.00001,side='lower') }) #Obtain the data set of test statistics simf<-list() for (xx in 1:2000){ if (any(simnull4[[xx]][192:1919,2]<0.01)){ simf[[xx]]<-NA } else{ simf[[xx]]<-simnull4[[xx]][1920,2] } } simf4a<-do.call(rbind,simf) simf<-list() for (xx in 1:2000){ if (any(simnull4[[xx]][192:1919,3]<0.01)){ simf[[xx]]<-NA } else{ simf[[xx]]<-simnull4[[xx]][1920,3] } } simf4b<-do.call(rbind,simf) #Ensure that around 1% of the trials stop for futility sum(is.na(simf4a)) #21 sum(is.na(simf4b)) #22 #Select a_U to make sure that the overall type I error is around 0.025 sum((simf4a[,1]>0.9836| simf4b[,1]>0.9836),na.rm=T )/2000#0.025 #The selected a_U is 0.9836. ``` ```R sim4<-lapply(1:1000,function(x) { sim_brar_unknown_var(Pats=10,nMax=50000,TimeToOutcome=expression(rnorm( length(vStartTime ),30, 3)),enrollrate=0.1,N1=192,armn=3, a_U=c(0.9836,0.9836),N2=1920,tp=1,armlabel=c(1,2,3),blocksize=6, mean=c(9.1/100,8.92/100,8.92/100),sd=c(0.009,0.009,0.009), minstart=192,deltaa=c(0.00077,0.00077),tpp=1,deltaa1=c(0,0), V01=1/2,a01=0.3,m01=9/100,b01=0.00001,side='lower') }) decision<-t(sapply(sim4, "[[", 1)) sum(decision[,2]=='Superiorityfinal')/1000 #0.882 #The simulated power from 1000 simulations is 0.882. ``` -------------------------------- ### Simulate NIG parameters and estimate posterior probabilities Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/vignettes/RARtrials.Rmd Simulates Normal-Inverse-Gamma parameters from chi-square hyper-parameters, then uses Monte Carlo draws to estimate posterior probabilities that each arm is the minimum (or maximum) response rate. The first block defines the conversion and sampling functions and generates 10,000 NIG draws per arm. The second and third blocks compute the probabilities for the lower and upper comparisons, respectively, with the commented results shown. ```R convert_chisq_to_gamma<-function(cpar){ list( m=cpar$mu, V=1/cpar$kappa, a=cpar$nu/2, b=cpar$nu*cpar$sigsq/2 ) } rnigamma<-function(nsim,par){ sigma2<-1/rgamma(nsim, shape=par$a,rate=par$b) mu<-rnorm(nsim, par$m, sqrt(sigma2*par$V)) cbind(mu,sigma2) } set.seed(12341) NIG_par1<-rnigamma(10000,convert_chisq_to_gamma(par1)) set.seed(12342) NIG_par2<-rnigamma(10000,convert_chisq_to_gamma(par2)) set.seed(12343) NIG_par3<-rnigamma(10000,convert_chisq_to_gamma(par3)) ``` ```R dat<-matrix(NA,10000,5) for (i in 1:10000){ dat1<-rnorm(10000,NIG_par1[i,1],NIG_par1[i,2]) dat2<-rnorm(10000,NIG_par2[i,1],NIG_par2[i,2]) dat3<-rnorm(10000,NIG_par3[i,1],NIG_par3[i,2]) dat[i,1]<-sum(dat1>dat2)/10000 dat[i,2]<-sum(dat1>dat3)/10000 minimal<-base::pmin(dat1,dat2,dat3) dat[i,3]<-sum(minimal==dat1)/10000 dat[i,4]<-sum(minimal==dat2)/10000 dat[i,5]<-sum(minimal==dat3)/10000 } mean(dat[,1]) mean(dat[,2]) #Return results: 0.1994863 0.8113217 mean(dat[,3]) mean(dat[,4]) mean(dat[,5]) #Return results: 0.1802267 0.0285785 0.7911948 ``` ```R dat<-matrix(NA,10000,5) for (i in 1:10000){ dat1<-rnorm(10000,NIG_par1[i,1],NIG_par1[i,2]) dat2<-rnorm(10000,NIG_par2[i,1],NIG_par2[i,2]) dat3<-rnorm(10000,NIG_par3[i,1],NIG_par3[i,2]) dat[i,1]<-sum(dat11.988,na.rm=T)/5000 #0.025 #Select a cut-off value to attain type I error 0.025 using the sample size 192 sim1b<-lapply(1: 5000, function(x){ sim_RPTW(Pats=10,nMax=50000,TimeToOutcome=0,enrollrate=1,na0=1,nb0=1,na1=1,nb1=1, h=c(0.5,0.7),N2=192,side='upper',Z=1.988)}) sum(sapply(sim1b,"[[",1)==1)/5000 #Using the selected cut-off value 1.988, we obtain power 0.7938, which is close to 0.8 #Example: RPTW(1,1) with the first 1 represents the initial number of balls for #each treatment group in the urn and the second 1 represents the number of balls #added to the urn when result of each participant becomes available. #Directly using asymptotic chi-square test which is equivalent to Z test statistics set.seed(12345) sim1<-lapply(1: 5000, function(x){ sim_RPTW(Pats=10,nMax=50000,TimeToOutcome=0,enrollrate=1,na0=1,nb0=1,na1=1,nb1=1, h=c(0.5,0.7),N2=192,side='upper')}) sum(sapply(sim1,"[[",1)==1)/5000 #Using Z test statistics from normal distribution, we obtain power of 0.8038 ``` -------------------------------- ### Gittins function usage Source: https://github.com/yayayaoyaoyao/rartrials/blob/main/man/Gittins.Rd Shows the function signature and argument documentation for Gittins. The function accepts a Gittinstype and a discount factor df, returning indices for binary or normal reward processes. ```r Gittins(Gittinstype, df) ```