# How to pick the winner in the "Play the Winner" treatment assignment scheme (Urn model)

"Play the winner" is an intutitive procedure whereby successful treatments in earlier trials have a higher probability of being assigned in later trials. Quoting https://online.stat.psu.edu/stat509/lesson/8/8.6:

Suppose there is a two-armed clinical trial and the urn contains one "A" ball and one "B" ball for the first patient. Suppose that the patient randomly is assigned treatment A. Now you need to know if the treatment was successful with the patient that received this treatment. If the patient does well on treatment A, then the original "A" ball and another "A" ball are placed in the urn. If the patient fails on treatment A, then the original "A" ball and a "B" ball are placed in the urn. Thus, the second patient has probability of 1/3 or 2/3 of receiving treatment A depending on whether treatment A was a success or failure for the first patient. This process continues. If one treatment is more successful than the other, the odds are stacked in favor of that treatment.

While there are more complex adaptive procedures, this one can be explained by an Urn model. It also has nice statistical properties. However, my understanding is that the probability does not trend towards 100% for the superior treatment, rather a proportion that depends on the relative superiority of the best treatment. Thus, at some point you have to call the winner.

Should you set a total N for the cumulative number of trials ahead of time? If so, what would be the best way to choose N? Would you have to adjust in any way for the unequal probabilities?

Alternately, could you run this as a sequential procedure and call the winner using a given stopping rule? What would that look like?

I got my hands on the 1978 JASA paper The Randomized Play-the-Winner Rule in Medical Trials and saw that the authors analyzed an "inverse" stopping rule:

"We continue taking observations until either $$S_A(n) + F_B(n)$$ or $$S_B(n) + F_A(n)$$ is equal to some given value $$r$$."

Above, $$S_i(n)$$ and $$F_i(n)$$ are the number of sucesses and failures observed for treatment $$i$$ after trial $$n$$, for both $$i=A, B$$.

They talk about how you can view the sequential experiment on a 2d lattice, with state and transition probabilities, and they prove that the probability of "correct selection" (chosing the better treatment) tends towards 1 as r tends towards infinity. Still, I think they had to run a simulation with different values of $$r$$ and the sucess probabilities. (This was the late 70s. I'd be really curious to learn whether this was analytical or a Monte Carlo experiment using a computer.)

In the table they present below, $$p_A$$ and $$p_B$$ are the success probabilities of treatment A and B, respectively. $$Pr(CS)$$ is the probability of correct stopping (note A is always superior) and $$ASN$$ is the average sample number at which $$r$$ was reached. Ignore the "PW" columns, and note that when I said "Play the Winner" in my original post, I meant "Randomized Play the Winner" (RPW).

Below is code that implements a Monte Carlo experiment to replicate these numbers. Again, I'm curious if there's a more elegant way to get these. Before posting the code, I'll conclude with the procedure:

1. Choose r based on probability of correct selection and average sample number, given hypothetical $$p_A$$ and $$p_B$$.
2. Run the Randomized Play the Winner strategy until one treatment has $$r$$ balls in the urn.
simulate_trial <- function(p_a, p_b, stop_value) {
# Initial counts
r_a <- 0
r_b <- 0

# Keep track of the history
history_a <- c()
history_b <- c()

while (r_a < stop_value && r_b < stop_value) {
# Calculate the current probability to draw treatment A
prob_a <- r_a / (r_a + r_b)

# If r_a and r_b are both 0, we need to randomly select a treatment since prob_a will be NaN
if (is.nan(prob_a)) {
prob_a <- 0.5
}

# Randomly draw a treatment
treatment_drawn <- ifelse(runif(1) < prob_a, "A", "B")

if (treatment_drawn == "A") {
# Check for success or failure
if (runif(1) < p_a) {
# Success for treatment A
r_a <- r_a + 1
} else {
# Failure for treatment A
r_b <- r_b + 1
}
} else {
# Check for success or failure
if (runif(1) < p_b) {
# Success for treatment B
r_b <- r_b + 1
} else {
# Failure for treatment B
r_a <- r_a + 1
}
}

# Append to history
history_a <- c(history_a, r_a)
history_b <- c(history_b, r_b)
}

return(list(history_a = history_a, history_b = history_b))
}

# Simulation parameters
grid <- merge(
data.frame(r=c(5, 10, 15, 30, 50)),
data.frame(p_a=c(.3, .7, .9, .5, .8, .7), p_b = c(.1, .5, .7, .1, .4, .1)),
all=TRUE
)[, c(2, 3, 1)]

M <- 100000  # Monte Carlo reps
for (j in 1:nrow(grid)) {
p_a <- grid[j, "p_a"]
p_b <- grid[j, "p_b"]
r <- grid[j, "r"]
cat("running with p_a =", p_a, ", p_b =", p_b, ", r =", r, "for", M, "reps\n")
correct_stopping <- numeric(M)
sample_size <- numeric(M)
for (rep in 1:M) {
result <- simulate_trial(p_a, p_b, r)
correct_stopping[rep] <- max(result$$history_a) == r sample_size[rep] <- length(result$$history_a)
}
grid[j, "pr_correct_stopping"] <- mean(correct_stopping)
grid[j, "average_sample_number"] <- mean(sample_size)
}

print(grid, digits=3)


• Regarding the statement “the probability of "correct selection" (chosing the better treatment) tends towards 1 as r tends towards infinity”, I have asked this question. It considers the paths on the 2d lattice as a Pólya urn type of process with an absorbing boundary at $n_a=n_b$. If some fraction of the paths never end up absorbed, then the probability of that statement can not tend towards 1. (correct selection requires all paths that initial move into the wrong direction to eventually turn into the right direction and cross $n_a=n_b$) Commented Nov 27, 2023 at 22:27