O problema

Dada a ideia geral, vamos definir os detalhes. Eu vou simplificar a vida colocando muita estrutura:

1. Existem dois tratamentos possíveis, tratamento zero e tratamento um
2. Cada tratamento tem distribuição binomial - ou o tratamento é um sucesso ou é um fracasso - com probabilidade de sucesso $p_0$ e $p_1$, respectivamente
3. As probabilidades são desconhecidas por quem vai escolher o tratamento
4. O objetivo de quem escolhe o tratamento é meramente maximizar o valor esperado
5. Nós vamos ter $N$ observações
6. Em cada período, a gente pode escolher entre alocar o sujeito pro tratamento 0 ou pro tratamento 1

Soluções?

Segue uma descrição da solução mais simples que eu consigo pensar pro problema: durante $n < N$ períodos, aloque alternadamente entre tratamento 0 ou tratamento 1. Calcule a média de sucessos de cada tratamento. Para períodos depois de $n$, aloque todo mundo para o tratamento com maior média.

Eu vou escrever uma função que implementa isso. A função retorna pra qual dos tratamentos o sujeito deve ser alocado:

set.seed(06012022)

simple_rule1 <- function(periodo,cutoff,media0,media1){
  
  if (periodo < cutoff){
    if(periodo %% 2 == 0){
      return(c(0,TRUE))
    } else{
      return(c(1,TRUE))
    }
    
  } else {
    if(media1 > media0){
      return(c(1,FALSE))
    } else{
      return(c(0,FALSE))
    }
  }
}

Veja que periodo %% 2 testa se o período é par (%% é a operação módulo, ou resto da divisão). Isso faz com que nos períodos par, a gente aloque para o tratamento zero e nos períodos ímpares a gente aloque pro tratamento um. A função retorna qual é o tratamento escolhido e também um TRUE ou FALSE. O motivo: eu quero que, uma vez que a função tome a decisão de qual é o melhor tratamento, ela nunca mais mude o valor da média. Isso não é necessariamente realista: você poderia sempre atualizar a sua média depois de tomar a decisão e no caso improvável que a média do tratamento que você escolheu fique pior que a média do outro tratamento, você muda o tratamento. Eu vou criar uma função é que o processo gerador de dados e nos dá tanto o sucesso ou fracasso nesta rodada quanto faz todas as atualizações de média necessárias:

dgp <- function(prob0,prob1,n0,n1,media0,media1,decision,update){
  
  if(decision == 0){
    
    suc <- rbinom(1,1,prob0)
    
    if(update){
    
      media0 <- (n0*media0 + suc)/(n0+1)
      n0 <- n0 + 1
    }
    
    return(c(suc,n0,n1,media0,media1))
    
  } else if(decision == 1){
    
    suc <- rbinom(1,1,prob1)
    
    if(update){
      
      media1 <- (n1*media1 + suc)/(n1+1)
      n1 <- n1 + 1
    }
    
    return(c(suc,n0,n1,media0,media1))
    
  }
  
}

Agora, vamos escrever a função que retorna uma simulação dado uma regra. Ela vai me retornar o payoff obtido pelo agente e qual é a decisão dele no último período:

simulacao <- function(N,prob0,prob1,rule,...){
  
  media0 <- 0
  media1 <- 0
  n0 <- 0
  n1 <- 0
  payoff <- 0
  
  for(i in 1:N){
    
    decision <- rule(periodo = i, media0 = media0, media1 = media1,...)
    
    if(decision[1] == 0){
      
      result <- dgp(prob0,prob1,n0,n1,media0,media1,decision[1],decision[2])
      n0 <- result[2]
      media0 <- result[4]
      payoff <- result[1] + payoff
      
    } else if(decision[1] == 1){
      
      result <- dgp(prob0,prob1,n0,n1,media0,media1,decision[1],decision[2])
      n1 <- result[3]
      media1 <- result[5]
      payoff <- result[1] + payoff
      
    } else{
      stop("Decisão deve ser 0 ou 1")
    }
    
  }
  return(c(payoff,decision[1]))
}

Vamos simular umas mil observações do nosso experimento e com a regra que nós temos. Eu vou deixar a probabilidade de sucesso do experimento 0 ser 0.5 e a probabilidade de sucesso do experimento 1 ser 0.7 e deixar ele alternar por metade dos 30 períodos:

library(purrr)

input <- replicate(1000,30,simplify = FALSE)

output <- map(input,simulacao, prob0 = 0.5,prob1 = 0.7,simple_rule1, cutoff = 15)

output <- do.call(rbind,output)
colnames(output) <- c("Payoff","Proporção de decisões 1")

colMeans(output)

##                  Payoff Proporção de decisões 1 
##                  18.589                   0.689

O próximo passo, naturalmente, é se perguntar quanto tempo você deveria alternar durante os dois tratamentos disponíveis. A próxima função vai permitir que a gente simule a regra simples com diferentes valores pro tempo que a gente passa testando as alternativas - ou explorando elas:

simulacao_cutoff <- function(N,prob0,prob1,cutoffs){
  
  cutoffs_list <- as.list(cutoffs)
  foo <- function(x){simulacao(N,prob0,prob1,simple_rule1,cutoff = x)}
  output <- map(cutoffs_list,foo)
  
  output <- do.call(rbind,output)
  
  return(output)
  
}

O código que faz a simulação, eu vou testar apenas 12 períodos de teste, 14 períodos de teste e 20 períodos de teste. É importante que o número de períodos de teste seja par para grantir que nós vamos ter o mesmo número de observações para cada tratamento. Eu vou fazer duas mil replicações:

input2 <- replicate(2000,c(12,14,20),simplify = FALSE)

output <- map(input2,simulacao_cutoff, prob0 = 0.5,prob1 = 0.7, N = 30)

O código a seguir só organiza a lista em um array com linhas, colunas e uma dimensão a mais e tira a média entre as 2000 simulações:

library(abind)

output <- do.call(abind,list(output,along = 3))

result <- apply(output,c(1,2),mean)
result <- cbind(c("12","14","20"),result)

colnames(result) <- c("n","Payoff","Proporção tratamento 1")

knitr::kable(result)

Esse resultado é bem interessante: note que o maior payoff, na média, é quando nós só gastamos 12 períodos explorando. Veja que a nossa taxa de escolher o tratamento com maior probabilidade é um pouco maior com 20 observações. Isso é totalmente natural, uma vez que a gente tem quase o dobro de observações. Mas isso também requer que a gente passe muito mais tempo jogando a alternativa ruim, o que reduz o payoff que a gente recebe.

Veja que a gente não explorou muita a estrutura de ser capaz de de decidir parar arbitrariamente os testes. De fato, a gente fez basicamente o que um Randomized Controlled Trial faria.

Uma regra simples que me ocorre que explora a possibilidade de parar a qualquer momento do tempo é a seguinte: alterne os tratamentos por $t_1$ períodos; de $t_1$ a $t_2$ continue alternando exceto se a diferença de probabilidades ultrapassar um limiar. Depois de $t_2$ ou de ultrapassar o limiar, jogue a estratégia com maior média.

Vamos escrever a regra como uma função no R. Dessa vez a função retorna tanto a decisão (como a primeira posição no vetor), se nós estamos atualizando a média ou não, e se nós paramos antes da hora::

simple_rule2 <- function(periodo,t_1,t_2,threshold,media0,media1){
  
  if(periodo < t_1){
    
    if(periodo %% 2 == 0){
      return(c(0,TRUE,FALSE))
    } else{
      return(c(1,TRUE,FALSE))
    }
    
  } else if (periodo < t_2){
    
    if(abs(media0 - media1) < threshold){
      
      if(periodo %% 2 == 0){
        return(c(0,TRUE,FALSE))
      } else{
        return(c(1,TRUE,FALSE))
      }
     
    } else{
      
      if(media0 > media1){
        return(c(0,FALSE,TRUE))
      } else{
        return(c(1,FALSE,TRUE))
      }
      
    }
    
    
  } else{
    
    if(media0 > media1){
      return(c(0,FALSE,FALSE))
    } else{
      return(c(1,FALSE,FALSE))
    }
  }
  
}

Eu vou querer saber quantas vezes eu paro antes de $n_2$, então eu vou escrever uma função de simulação nova (desculpa, don’t repeat yourself):

simulacao2 <- function(N,prob0,prob1,t_1,t_2,threshold){
  
  media0 <- 0
  media1 <- 0
  n0 <- 0
  n1 <- 0
  payoff <- 0
  early_stop <- FALSE
  
  for(i in 1:N){
    
    decision <- simple_rule2(i,t_1,t_2,threshold,media0,media1)
    
    if(decision[1] == 0){
      
      result <- dgp(prob0,prob1,n0,n1,media0,media1,decision[1],decision[2])
      n0 <- result[2]
      media0 <- result[4]
      payoff <- result[1] + payoff
      
      early_stop <- ifelse(early_stop,TRUE,ifelse(decision[3],TRUE,FALSE))
        
    } else if(decision[1] == 1){
      
      result <- dgp(prob0,prob1,n0,n1,media0,media1,decision[1],decision[2])
      n1 <- result[3]
      media1 <- result[5]
      payoff <- result[1] + payoff
      
      early_stop <- ifelse(early_stop,TRUE,ifelse(decision[3],TRUE,FALSE))
      
    } else{
      stop("Decisão deve ser 0 ou 1")
    }
    
  }
  return(c(payoff,decision[1],early_stop))
}

Vamos simular o caso em que nós observamos pelo menos 10 observações, e se a diferença não for grande o bastante, nós observamos até 15 períodos. Eu vou colocar o threshold em 0.5, o que é relativamente alto:

input <- replicate(1000,30,simplify = FALSE)

output <- map(input,simulacao2, prob0 = 0.5,prob1 = 0.7,t_1 = 10, t_2 = 15, threshold = 0.5)

output <- do.call(rbind,output)

colnames(output) <- c("Payoff","Proporção decisão 1", "Proporção Early Stop")

knitr::kable(colMeans(output))

Nós temos agora uns três parâmetros livres: quantas observações nós fazemos no mínimo antes de tomar a decisão, quantas observações nós fazemos no máximo antes de tomar a decisão, e qual a diferença mínima que nós exigimos para tomar uma decisão antes de ver o máximo de observações. Eu vou brincar de mudar só esse último parâmetro. Primeiro, uma função que permite a gente simular, para vários limiares, a escolha que a regra faria:

simulacao_threshold <- function(N,prob0,prob1,t_1,t_2,thresholds){
  
  thresholds_list <- as.list(thresholds)
  results <- map(thresholds_list,simulacao2,N = N,prob0 = prob0,prob1 = prob1,t_1 = t_1,t_2 = t_2)
  results <- do.call(rbind,results)
  rownames(results) <- thresholds
  colnames(results) <- c("Payoff","Proporção Tratamento 1","Proporção Early Stop")
  
  return(results)
  
}

A simulação em si. EU to mantendo tudo como nos problemas anteriores:

input <- replicate(2000,c(0.2,0.5,0.7),simplify = FALSE)
output <- map(input,simulacao_threshold,N = 30,prob0 = 0.5,prob1 = 0.7, t_1 = 10,t_2 = 15)

output <- do.call(abind,list(output,along = 3))
result <- apply(output,c(1,2),mean)

knitr::kable(result)

O resultado sugere que, depois de dez observações, mesmo se a diferença entre as probabilidades for só 0.2, o nosso payoff já é maior que com regras que exigem mais evidência de que um tratamento é melhor que o outro. Com 5 simulações em cada tratamento, isso significa que se um tratamento tiver um sucesso a mais que o outro, nós já deveríamos parar e escolher o tratamento com mais sucessos. A diferença de payoff entre a diferença de probabilidade ser 0.2 ou 0.5 ser tão pequena sugere que talvez algum valor intermediário seja melhor que 0.2 e 0.5.

Existem um milhão de maneiras de brincar com esse problema. No próximo post, eu vou tentar formular isso como um problema de otimização.