Upload New File

power_biasadj_unadj_stratWL.R

+#####################################################################################################################################################################################################################################################################################################################################################
+
+#' Function to calculate the mean power of the bias-adjusted and bias-unadjusted stratified Wei-Lachin test regarding a sample of randomization lists of several randomization procedures
+
+#####################################################################################################################################################################################################################################################################################################################################################
+
+#Libraries
+
+library(randomizeR)
+library(sadists)
+library(latex2exp)
+
+
+
+
+#' @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 mu_E expected value of the experimental group
+#' @param mu_C expected value of the control group
+#' @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 power of the bias-adjusted and bias-undadjusted stratified Wei-Lachin test
+
+
+Power_strat_WL<-function(N, K , balancing, m, Sigma, mu_E, mu_C, randomization, eta, r_number, alpha=0.05, BSD_par=3, EBC_par=0.67, MP_par=3, PBR_par=4){
+  
+  set.seed(1)
+  
+  power_unadjusted_test<-c()
+  power_adjusted_test<-c()
+  
+  #' 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<-samplebias[[1]]
+  biasseq<-samplebias[[2]]
+  
+  for(j in 1:r_number){
+    
+    #' calculation of the doubly non-centrality t-distribution under H_0/H_1 of the stratified Wei-Lachin test statistic for biased conditions
+    
+    numerator_delta<-c()
+    numerator_lambda<-c()
+    
+    numerator_delta_crit<-c()
+    
+    
+    d<-matrix(1,nrow=m,ncol=1)
+    weights<-c()
+    
+    for (i in 1:K){
+      
+      n_E <- sum(sample[[i]][j,]=="A")
+      n_C <- sum(sample[[i]][j,]=="B")
+      
+      
+      weights<-c(weights,(n_E*n_C)/(n_E+n_C))
+      
+      
+      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)
+      
+      
+      
+      tau_E<-1/n_E*rowSums(tau[,sample[[i]][j,]=="A"])
+      tau_C<-1/n_C*rowSums(tau[,sample[[i]][j,]=="B"])
+      
+      
+      numerator_delta<-c(numerator_delta,(n_E*n_C)/(n_E+n_C)*t(d)%*%(mu_E-mu_C+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 )
+      
+      numerator_delta_crit<-c(numerator_delta_crit,(n_E*n_C)/(n_E+n_C)*t(d)%*%(tau_E-tau_C))
+    }
+    
+    
+    delta<-sum(numerator_delta)/sqrt(sum(weights)*t(d)%*%Sigma%*%d)
+    lambda<-sum(numerator_lambda)/(t(d)%*%Sigma%*%d)
+    
+    delta_crit<-sum(numerator_delta_crit)/sqrt(sum(weights)*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 power of the bias-unadjusted startified Wei-Lachin Test
+    
+    power_unadj<-1-doublyT(qt(1-alpha,N-2*K),N-2*K,delta, lambda,lb=lb_ac,ub=ub_ac)
+    
+    #' calculation of the power of the bias-adjusted startified Wei-Lachin Test
+    power_adj<-1-doublyT(qdnt(1-alpha, N-2*K, delta_crit, lambda),N-2*K,delta, lambda,lb=lb_ac,ub=ub_ac) 
+    
+    
+    power_adjusted_test<-c(power_adjusted_test,power_adj)
+    power_unadjusted_test<-c( power_unadjusted_test,power_unadj)
+    
+    
+  }
+  
+  #' calculating the mean power of the bias-adjusted and bias-unadjusted stratified Wei-Lachin test
+  mean_power_adj<-mean(power_adjusted_test)
+  mean_power_unadj<-mean(power_unadjusted_test)
+  
+  
+  #' computing the standard error of the Monte Carlo simulations
+  error_adj<-sd(power_adjusted_test)/r_number
+  error_unadj<-sd(power_unadjusted_test)/r_number
+  
+  out<-data.frame(samplesize=N,number_of_strata=K,balancing=I(list(balancing)),number_of_endpoints=m,biasfactor=I(list(eta)), rand=randomization, mean_power_adj=mean_power_adj, mean_power_unadj=mean_power_unadj)
+  
+  return(out)
+  
+}
+
+
+
+#Check the 80% power in the unbiased case
+
+s<-Power_strat_WL(80,  2, c(1,1),2,diag(2),c(0.4,0.4),rep(0,2),'PBR', matrix(rep(0,4),nrow=2,ncol=2),1)
+print(s)
+
+
+
+
+#######################################################################################################################################################################################################################
+
+#' 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"
+
+#######################################################################################################################################################################################################################
+
+#'Mean power across 10 000 randomization lists generated by different RPs in clinical trials with N = 32 patients, m = 2 endpoints, K = 2 strata and common allocation bias effect across endpoints and strata chosen as proportion ρ ∈ {0, 0.2, 0.4, 0.6, 0.8, 1} of the effect sizes ∆_{32,2,2} = 0.64.
+
+#set the file path to save the simulation output
+
+setwd("/filepath[...]")
+file.create("power_stratWeiLachin.csv")
+
+
+p=c(0,0.2,0.4,0.6,0.8,1,4,8)
+
+randomization=c( 'PBR', 'RAR','MP', 'BSD', 'EBC')
+
+
+for ( i in randomization){
+  for(j in p){
+    
+    write.table(Power_strat_WL(32,  2, c(1,1),2,diag(2),rep(0.64,2),rep(0,2),i, matrix(rep(j*0.64,4),nrow=2,ncol=2),10000),"power_stratWeiLachin.csv",,append=TRUE ,row.names=FALSE,col.names = FALSE)
+    
+    }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+