Tratamento Sequencial I
Neste post, eu vou fingir que eu sou um microeconomista aplicado ou um médico e a minha pergunta é se um dado tratamento - no sentido literal se você for um médico - é mais eficiente que o outro.
Nós estamos acostumados a pensar testes de hipótese da seguinte forma: nós temos \(N\) observações, e algumas foram tratadas - talvez aleatoriamente, talvez selecionado em observáveis ou não observáveis. Nós então formulamos o teste de hipótese e buscamos cometer um erro do tipo I com uma probabilidade no máximo \(\alpha\).
Mas nós podemos ter um outro cenário: primeiro, uma unidade chega na nossa mão e nós escolhemos se alocamos essa unidade para o tratamento A ou B (o tratamento B pode muito bem ser um placebo). Nós observamos o resultado da unidade. Chega uma nova unidade e nós precisamos decidir se alocamos ela para tratamento A e B.
Além de médicos e microeconomistas, esse cenário é interessante em vários outros contextos. Por exemplo, você é dono de um site que faz vendas online e precisa decidir onde vai colocar o botão de comprar de forma a maximizar a chance do usuário comprar o produto.
Este post não vai ser nada estruturado, eu vou simplesmente fazer umas simulações com ideias razoáveis e deixar qualquer formalização da solução do problema pro post seguinte.
O problema
Dada a ideia geral, vamos definir os detalhes. Eu vou simplificar a vida colocando muita estrutura:
- Existem dois tratamentos possíveis, tratamento zero e tratamento um
- Cada tratamento tem distribuição binomial - ou o tratamento é um sucesso ou é um fracasso - com probabilidade de sucesso \(p_0\) e \(p_1\), respectivamente
- As probabilidades são desconhecidas por quem vai escolher o tratamento
- O objetivo de quem escolhe o tratamento é meramente maximizar o valor esperado
- Nós vamos ter \(N\) observações
- 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.689O 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)| n | Payoff | Proporção tratamento 1 |
|---|---|---|
| 12 | 19.1205 | 0.77 |
| 14 | 19.0655 | 0.7875 |
| 20 | 18.8195 | 0.8385 |
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))| x | |
|---|---|
| Payoff | 19.306 |
| Proporção decisão 1 | 0.836 |
| Proporção Early Stop | 0.310 |
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)| Payoff | Proporção Tratamento 1 | Proporção Early Stop | |
|---|---|---|---|
| 0.2 | 19.2420 | 0.8205 | 0.8075 |
| 0.5 | 19.3550 | 0.8525 | 0.3395 |
| 0.7 | 19.2185 | 0.8510 | 0.0710 |
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.