Upload New File

T1EstratWL.R

+#####################################################################################################################################################################################################################################################################################################################################################
+
+#' Function to calculate the mean T1E under misspecification and the probabilities of obtaining an T1E under misspecification that not exceeds the significance level [P(T1E<=0.05)] regarding different randomization procedures under allocation bias in clinical trials that are evaluated with the stratified Wei-Lachin test.
+
+#####################################################################################################################################################################################################################################################################################################################################################
+
+#Libraries
+library(randomizeR)
+
+
+
+
+#' @param N total sample size
+#' @param K number of strata
+#' @param balancing the balancing of the total sample across the strata, i.e. in the case of balanced strata balancing=rep(1,K)
+#' @param m number of endpoints
+#' @param randomization one of the randomization procedures "RAR", "PBR",  "BSD", "MP", "EBC" (default parameter for this randomization procedures: PBR(4), BSD(3), EBC(0.67), MP(3))
+#' @param r_number number of randomization lists that contained in the sample
+#' @return A list that consists of two lists containing matrices representing the samples of randomization lists per strata and the matrices containing the sign of the biasing factor corresponding to these randomization lists
+
+stratRandsamples<-function(N, K, balancing, randomization, r_number,BSD_par=3,EBC_par=0.67,MP_par=3,PBR_par=4){
+  
+  #set seed for replication of simulations
+  set.seed(1)
+  
+  #check input variables
+  if(N<=0 | K<=0 | r_number<=0){
+    stop('Invalid number for N,m or r_number')
+  }
+  if(length(balancing)!=K)
+  {
+    stop('Invalid length of balancing vector')
+  }
+  
+    biasseq<-list()
+    sample_Randlist<-list()
+  
+  
+  # determination of the sample size per strata
+  sizepercenter=N/sum(balancing)*balancing
+  
+  #generating strata-wise randomization lists
+  for (i in 1:K){
+    if(randomization=='BSD'){
+      Rand<- bsdPar(sizepercenter[i], BSD_par, groups = LETTERS[1:2])
+      myPar<- genSeq(Rand,r_number)
+      Rand_list<- getRandList(myPar)
+    }else if(randomization=='EBC'){
+      Rand<-ebcPar(sizepercenter[i], EBC_par, groups = LETTERS[1:2])
+      myPar<- genSeq(Rand,r_number)
+      Rand_list<- getRandList(myPar)
+    }else if(randomization=='MP'){
+      Rand<-mpPar(sizepercenter[i] ,mti=MP_par, ratio = rep(1, 2), groups = LETTERS[1:2])
+      myPar<- genSeq(Rand,r_number)
+      Rand_list<- getRandList(myPar)
+    }else if(randomization=='PBR'){
+      if(sizepercenter[i] %% PBR_par==0){
+        Rand<-pbrPar(rep(PBR_par, times= sizepercenter[i]/PBR_par), K = 2, groups = LETTERS[1:2])
+        myPar<- genSeq(Rand,r_number)
+        Rand_list<- getRandList(myPar)
+      }else{
+        stop("Size of one center cannot be divided by PBR_par")
+      }
+    }else if(randomization=='RAR'){
+      Rand<-rarPar(sizepercenter[i],K=2,groups = LETTERS[1:2])
+      myPar<- genSeq(Rand,r_number)
+      Rand_list<- getRandList(myPar)
+    }else{
+      stop("Randomization procedure is not defined")
+    }
+    
+    sample_Randlist<-append(sample_Randlist,list(Rand_list))
+    biasseq<-append(biasseq,list(getExpectation(myPar, selBias("CS", 1, "exact"))))
+    
+  }
+  samplebias<-list(sample_Randlist,biasseq)
+  
+  return(samplebias)
+  
+}
+
+
+
+
+
+
+#' @param N total sample size
+#' @param K number of strata
+#' @param balancing the balancing of the total sample across the strata, i.e. in the case of balanced strata balancing=rep(1,K)
+#' @param m number of endpoints
+#' @param Sigma covariance matrix of the patient responses regarding the m endpoints
+#' @param randomization one of the randomization procedures "RAR", "PBR",  "BSD", "MP", "EBC" (default parameter for this randomization procedures: PBR(4), BSD(3), EBC(0.67), MP(3))
+#' @param eta Kxm matrix containing the endpoint and strata specific biasing factors
+#' @param r_number number of randomization lists that be used to calculate mean T1E under misspecification and the probabilities of obtaining an a T1E under misspecification that exceeds the 5% significance level
+#' @param alpha significance level (default: 0.05)
+#' @return a dataframe with variables @param N, @param K, @param m, @param eta, @param randomization,  mean T1E under misspecification, P(T1E<=0.05), standard error of the Monte Carlo simulations
+
+
+stratWeiLachin<-function(N, K , balancing, m, Sigma, randomization, eta, r_number,alpha=0.05,BSD_par=3,EBC_par=0.67,MP_par=3,PBR_par=4){
+  
+  #' generation of the sample of stratified randomization lists and the corresponding matrices containing the signs of the biasing factors
+  
+  samplebias<-stratRandsamples(N,K,balancing, randomization,r_number, BSD_par=BSD_par,EBC_par=EBC_par,MP_par=MP_par,PBR_par=PBR_par)
+  
+  sample_Randlist<-samplebias[[1]]
+  biasseq<-samplebias[[2]]
+  
+  
+  typeIerror<-c()
+  counter<-0
+  
+  for(j in 1:r_number){
+    
+    #' calculation of the doubly non-centrality t-distribution of the stratified Wei-Lachin test statistic under biased conditions
+    
+    numerator_delta<-c()
+    numerator_lambda<-c()
+    
+    d<-matrix(1,nrow=m,ncol=1)
+    weights<-c()
+    
+    for (i in 1:K){
+      
+      n_E <- sum(sample_Randlist[[i]][j,]=="A")
+      n_C <- sum(sample_Randlist[[i]][j,]=="B")
+      
+      
+      weights<-c(weights,(n_E*n_C)/(n_E+n_C))
+      
+      print(matrix(rep(eta[i,], times = n_E+n_C), ncol  =  n_E+n_C, byrow = FALSE))
+      tau <-  matrix(rep(eta[i,], times = n_E+n_C), ncol  =  n_E+n_C, byrow = FALSE)* matrix(rep(biasseq[[i]][j,], times = m), nrow  = m, byrow = TRUE) 
+      
+      if(sum(sample_Randlist[[i]][j,]=="A")>1){
+        tau_E<-1/n_E*rowSums(tau[,sample_Randlist[[i]][j,]=="A"])
+      }else if(sum(sample_Randlist[[i]][j,]=="A")==1){
+        tau_E<-1/n_E*tau[,sample_Randlist[[i]][j,]=="A"] 
+      }
+      else{tau_E=rep(0,m)}
+      
+      
+      
+      if(sum(sample_Randlist[[i]][j,]=="B")>1){
+        tau_C<-1/n_C*rowSums(tau[,sample_Randlist[[i]][j,]=="B"])
+      }else if(sum(sample_Randlist[[i]][j,]=="B")==1){
+        tau_C<-1/n_C*tau[,sample_Randlist[[i]][j,]=="B"] 
+      }
+      else{tau_C=rep(0,m)}
+      
+      
+      
+      numerator_delta<-c(numerator_delta,(n_E*n_C)/(n_E+n_C)*t(d)%*%(tau_E-tau_C))
+      numerator_lambda<-c(numerator_lambda, sum((t(d)%*%tau)^2)-n_E*(t(d)%*%tau_E)^2-n_C*(t(d)%*%tau_C)^2)
+      
+    }
+    
+    delta<-sum(numerator_delta)/sqrt(sum(weights)*t(d)%*%Sigma%*%d)
+    lambda<-sum(numerator_lambda)/(t(d)%*%Sigma%*%d)
+    
+    ub_ac<-ceiling(lambda/2 + qpois(.995, lambda/2))
+    lb_ac<-max(floor(lambda/2 - qpois(.995, lambda/2)), 0)
+    
+    #' calculation of the misspecified T1E of the stratified Wei-Lachin Test
+    
+    T1E<-1-doublyT(qt(1-alpha,N-2*K),N-2*K,delta, lambda,lb=lb_ac,ub=ub_ac)
+    
+    #'computing the number of T1E's under misspecification that less or equal to @param alpha
+    
+    if(isTRUE(all.equal(T1E, alpha))){
+      counter=counter+1
+    }else if(T1E<alpha){
+      counter=counter+1
+    }
+    
+    typeIerror<-c(typeIerror,T1E)
+  }
+  
+  #' calculating the mean T1E under misspecification regarding the sample of stratified randomization lists
+  mean_typeIerror<- mean(typeIerror)
+  
+  #' calculating the go_nogo criteria, the probability of obtaining a T1E under misspecification less or equal to @param alpha (P(T1E<=0.05))  
+  go_nogo<-counter/(r_number)*100
+  
+  #' computing the standard error of the Monte Carlo simulations
+  sd<-sd(typeIerror)
+  stderror<-sd/sqrt(r_number)
+  
+  #' table with summary measures
+  out<-data.frame(samplesize=N,number_of_strata=K,balancing=I(list(balancing)),number_of_endpoints=m,biasfactor=I(list(eta)), rand=randomization, mean_typeIerror=mean_typeIerror , go_nogo=go_nogo,stderror=stderror)
+  
+  #' table with all calculated error rates
+  
+  out1<-data.frame(typeIerror=typeIerror,samplesize=rep(N,r_number),number_of_strata=rep(K,r_number), balancing=I(replicate(r_number, balancing, simplify = FALSE)), number_of_endpoints=rep(m,r_number),biasfactor=I(replicate(r_number, eta, simplify = FALSE)), rand=rep(randomization,r_number))
+  
+  
+  return(out) 
+}
+
+
+
+
+#######################################################################################################################################################################################################################
+
+#' Simulations of the simulation study outlined in the manuscript "A bias-adjusted analysis of stratified clinical trials with multi-component endpoints using the Wei-Lachin test"
+
+#######################################################################################################################################################################################################################
+
+#set the file path to save the simulation output
+
+setwd("/filepath[...]")
+file.create("stratWeiLachin.csv")
+
+
+#' Note ∆_{N,K,m} denotes the effect size that corresponds to an 80% power of the stratified Wei-Lachin test in trials with N patients, m endpoints and K strata (∆_{32,2,2} = 0.64, ∆_{32,2,4} = 0.46, ∆_{32,4,2} = 0.66 and ∆_{32,4,4} = 0.47)
+effectsizes=c(0.64,0.46,0.66,0.47)
+effectsizes=matrix(effectsizes,ncol=2,nrow=2, byrow=T)
+
+
+#' number of endpoint components
+m=c(2,4)
+
+#' number of strata
+K=c(2,4)
+
+
+#' the common biasing factor for each strata and endpoint component is chosen as proportion (5%,10%) of the effect size that corresponds to a nominal power of 80% of the stratified Wei-Lachin
+p=c(0.05,0.1)
+
+#' analyzed stratified randomization procedures
+rand=c('RAR','BSD','MP','EBC','PBR')
+
+#' run the simulations across the different simulation settings
+
+for(e in 1:length(m)){  
+  
+  for (a in rand){
+    
+    for(b in 1:length(K)){
+      
+      for (z in p) {
+        
+        write.table(stratWeiLachin(N=32,K=K[b],balancing=rep(1,K[b]),m=m[e],Sigma=diag(m[e]), randomization=a, eta= matrix(rep(z*effectsizes[b,e],m[e]*K[b]),nrow=K[b],ncol=m[e]), r_number=10000),"stratWeiLachin.csv")
+      }
+    }
+  }
+}
+
+
+