Clinical trials simple randomization scheme simulation [closed]

Question

I'm trying to do 1000 simulations for how to randomize 100 patients to two treatment groups (A and B) using simple randomization scheme (with probability 0.5) and figure out what is the chance of a greater than 55-45% split (i.e.45 or fewer subjects will be assigned one of the treatments and 55 or more on the other treatment). Assuming the response rates are 0.6 and 0.5 for treatment A and B, respectively.

I know the basic statistic idea of how to do this: Simple Randomization scheme is equivalent with tossing a fair coin for each subject that enters the trial ( Heads = Treatment A; Tails = Treatment B)

so we can use rbinom (1,100,0.5) to generate 100 random numbers --> then we have our two treatment groups ( 1 = group A; 0 = group B) and count the number of 1s as sample size in group A, and 0s as sample size in group B --> if the split is greater than 55-45% we count as 1; else count as 0. Then we repeat this 1000 times, count how many times the split is greater than 55-45 and get the probability as n/1000.

But I dont know how to translate these steps into R code?

Answer 1

I will first assume that you are required to perform all these steps, as in an R exercise.

Let's start with

rbinom(100, 1, 0.5) 

which gets you one assignment vector. The ones are patients assigned to treatment A, the zeroes those assigned to treatment B. Notice here that the first argument of the rbinom function is the number of draws from the binomial distribution. Here, we draw 100 independent Binomial(n=1, p=0.5), which, since n=1, are Bernulli random variables. You can take the sum of these 100 Bernulli random variables, and call them sum1: this corresponds to the number of people assigned to treatment A in this first assignment.

sum1 <- sum(rbinom(100,1,0.5))

As you said, you want to check whether you are outside the 45-55% split range. This means that you want to check whether sum1 is greater than 55, or lower than 45. (By symmetry, this also checks this range for those assigned to treatment B, because they must sum up to 100.)

If we are outside this range for this one assignment vector, we want to return 1. We use a logical test for that:

(sum1 > 55)  || (sum1 < 45).

If either side of || is TRUE, then this will be TRUE.

Or, alternatively, you could check whether the difference to 50 (a balanced assignment) is greater than 5, in absolute value:

abs(sum1 - 50) > 5

You're done for 1 assignment. Now you have to do it for N=1000, say, and compute the proportion of assignments that fall outside the pre-specified range. You could write a for loop:

# Option 1: `naive` simulation
N <- 1000
n_out <- 0
for(i in 1:N){
  sum1 <- sum(rbinom(100, 1, 0.5))
  outside <- abs(sum1 - 50) > 5
  n_out <- n_out + outside
}
proportion <- n_out/N

Now you certainly know that a sum of i.i.d Bernulli is a Binomial, and so sum1 is simply one draw from a Binomial(n=100, p=0.5)

sum2 <- rbinom(1, 100, 0.5)

We want N=1000 such samples, and count how many of them are outside the "balancing" range. This can be written sunccinctly in one single line:

# Option 2: better simulation
N <- 1000
n_out <- sum(abs(rbinom(N, 100, 0.5) - 50) > 5)
proportion <- n_out/N

Finally, denote by $X$ the random variable "number of patients assigned to treatment A out of 100 patients". Recall that $X$ is a Binomial(n=100, p=0.5) random variable. You can compute the exact probability of the event "the sample is outside the "balancing" range", which we call event $A$:

$P(A) = P(X < 45) + P (X > 55) = P(X \leq 44) + (1- P(X \leq 55))$

with the corresponding Rcode:

# Option 3: true proportion
proportion <- pbinom(44, 100, 0.5) + (1 - pbinom(55, 100, 0.5))

You'll see, as you increase N in your simulations, that you will come closer and closer to this number.