Simularea unei simple răspunsuri randomizare adaptativă - trebuie să o văd pentru a -l crede

În simulările mele de randomizare adaptativă a răspunsului, am descoperit că se comportă comparabil la alocarea fixă de 50-50 în identificarea efectelor tratamentului. Abordarea adaptativă pare să funcționeze! Cu toate acestea, cu doar 10 încercări, am zgâriat doar suprafața. Există limitări importante – riscuri de prejudecăți temporale, ineficiență statistică și ajustări complexe de multiplicitate în cadrele bayesiene.

Obiective

Care este randomizarea adaptativă a răspunsului?

Randomizarea adaptivă a răspunsului (RAR) este o tehnică utilizată în studiile clinice în care alocarea pacienților la diferite modificări ale brațelor de tratament pe baza rezultatelor provizorii colectate în timpul studiului. Spre deosebire de randomizarea fixă tradițională, unde raporturile de alocare a tratamentului rămân constante pe parcursul studiului, RAR ajustează probabilitatea alocării grupurilor de tratament pe măsură ce datele se acumulează, de obicei pentru a atribui mai mulți pacienți la tratamente care par să fie mai bune. Această abordare este concepută pentru a atinge două obiective principale: maximizarea informațiilor obținute despre tratamente superioare, reducând în același timp numărul de pacienți expuși la tratamente inferioare în timpul studiului în sine.

Motivația principală pentru utilizarea RAR este etică – își propune să trateze mai mulți participanți la încercări în mod eficient, în timp ce încă adună date suficiente pentru concluzii științifice. Este luat în considerare în special în studiile care implică condiții grave cu rate mari de mortalitate sau în studiile cu mai multe etaje în care sunt comparate mai multe tratamente experimentale. RAR a primit o atenție teoretică considerabilă încă din anii 1930, dar rămâne controversat în practică, susținătorii citând potențialele avantaje etice ale acestuia, în timp ce criticii indică îngrijorări cu privire la validitatea statistică, prejudecățile potențiale din tendințele temporale și o complexitate crescută în proiectarea și analiza încercărilor.

Remap-Cap a fost primul studiu pe care l-am întâlnit care m-a prezentat la acest design și am găsit teoria destul de interesantă și interesantă! Pare prea bine pentru a fi adevărat, dar pentru mine nu pot obține intuiția în spatele modului în care randomizarea adaptativă poate să o identifice. Din fericire, există un alt mod în care ne putem convinge că acest lucru funcționează și că este să o simulăm pentru noi înșine, stilul bayesian!

Formula adaptivă pentru randomizarea adaptativă a răspunsului

Abordarea Thall and Wathen

Aceasta este o formulă utilizată frecvent care ajustează probabilitățile de randomizare bazate pe probabilitățile posterioare. Pentru un studiu cu două brațe, probabilitatea alocării unui pacient la tratamentul 1 este:

begin{gather} π₁,ᵢ = frac{(P(p₁ > p₀|data))^c}{(P(p₁ > p₀|data))^c + (1-P(p₁ > p₀|data))^c} end{gather}

În cazul în care p (p₁> p₀ | date) este probabilitatea posterioară ca tratamentul 1 să fie mai bun decât tratamentul 0, iar C este un parametru de reglare care controlează viteza de adaptare. Când C = 0, acest lucru se reduce la randomizarea egală; Când C = 1, devine eșantionare Thompson.

Pot exista și alte abordări ale RAR, dar acesta este cel care are sens și ușor de implementat. Scufundați -vă mai adânc

Plan de simulare

Sperăm că simularea acestui lucru de mai multe ori ar putea să -mi ofere o idee mai bună despre cum va funcționa în practică. Cum facem asta?

Alocare adaptivă față de 50-50

Stabiliți o lume în care avem întreaga populație și știm cum răspund cu și fără tratament.
Vom stabili un efect de tratament ca -1.09ceea ce înseamnă că tratamentul este mai bun la reducerea unei proporții de eveniment decât la control. Formula ar fi log odds = b0 + b1 * treatmentunde b0 este 0, și b1 este efectul nostru de tratament.
Vom scrie apoi un cod pentru a simula ambele Response Adaptive Randomization şi 50-50 allocationcu 50 de seturi. Fiecare set va avea apoi 4 analize secvențiale și eșantionare. Prima analiză va avea 50 de pacienți cu alocare aleatorie cu 50%. Pentru analizele provizorii ulterioare, grupul de randomizare adaptativ de răspuns va depinde de răspuns, în timp ce 50-50 grup de alocare va rămâne la 50%. Pentru randomizarea adatpivă, vom folosi abordarea Thall și Wathen cu c între 0 şi 0.5. Vom folosi formula Thall and Warten pentru actualizarea c cu n / (2*N)unde n este numărul de pacienți eșantionați până acum, în plus, cei vor fi prelevați acest lucru în studiul provizoriu și N este numărul maxim de pacienți în studiu. După cum puteți vedea, atunci când este eșantionat în total 200 de pacienți, C va fi 0.5 la ultima analiză provizorie.
Fiecare studiu pentru ambele grupuri este setat cu aceeași sămânță pentru a asigura reproductibilitatea.
După fiecare studiu, extragem coeficientul de interes (în cazul nostru Beta1) și, de asemenea, evaluăm care este probabilitatea ca Beta1 să fie mai mică de 0, ceea ce înseamnă că tratamentul este util la reducerea oricărui rezultat pe care dorim să îl reducem.

Întrebări pe care le am pentru mine

Care este diferența dintre Response Adaptive Randomization şi 50-50 allocation În ceea ce privește estimările efectului de tratament? RAR va avea aceeași precizie la identificarea efectului tratamentului? Va identifica RAR REAL să identifice efectul tratamentului mai repede de 50-50 de alocare?
Dacă nu există niciun efect, Rar va putea să -l tachineze? Și cum ar arăta asta?

Cod

library(tidyverse)
library(cmdstanr)
library(glue)

## comparing 50-50 allocation vs Adaptive Randomization
method_list <- c("50_50","adaptive")
seeds <- sample(1:100000, size = 10, replace = F)


## set up empty dataframe
result <- tibble(
  b0 = numeric(),
  b0_lower = numeric(),
  b0_upper = numeric(),
  b1 = numeric(),
  b1_lower = numeric(),
  b1_upper = numeric(),
  prob = numeric(),
  treatment_num = numeric(),
  num_sample = numeric(),
  treatment_num_cum = numeric(),
  num = numeric(),
  method = character(),
  seed = numeric()
)

for (seed in seeds) {
for (method in method_list){
treatment_effect <- -1.09

### Set Up Entire Population, knowing both effects of placebo and treatment of each patient
set.seed(seed)
N <- 100000
x0 <- replicate(N, 0)
y0 <- rbinom(N, 1, plogis(treatment_effect*x0))
x1 <- replicate(N,1)
y1 <- rbinom(N, 1, plogis(treatment_effect*x1))
df <- tibble(y0,y1) |>
  mutate(id = row_number())


### Max sample
max_n <- 200

### How many sampling trials
n <- 50
sample_number <- max_n / n
b0 = b0_lower = b0_upper = b1 = b1_lower = b1_upper = prob_vec = treatment_num = num_vec = vector(mode = "numeric", length = sample_number)
x_vec <- c()

### Changing c parameter
n_sum <- 0
c_param <- function(x) {
  value <- x / (2*max_n)
  return(value)
}

for (i in 0:sample_number) {
  ## Set initial parameters and don't log it
  if (i == 0) {
    b0mu <- 0
    b0sd <- 10
    b1mu <- 0
    b1sd <- 2.5
    diff <- 0.5
    x <- rbinom(n, 1, 0.5)
    num_vec(i+1) <- n
    n_sum <- n
  } else {
  ## Update prior
b0mu <- fit$summary("beta0")(('median'))
b0sd <- fit$summary("beta0")(('sd'))
b1mu <- fit$summary("beta1")(('median'))
b1sd <- fit$summary("beta1")(('sd'))
b0(i) <- b0mu
b0_lower(i) <- fit$summary("beta0")(("q5"))
b0_upper(i) <- fit$summary("beta0")(("q95"))
b1(i) <- b1mu
b1_lower(i) <- b1_lower_i
b1_upper(i) <- b1_upper_i
n_sum <- n_sum + n


## Assigment of sampling proportion; 50% allocation at all times vs Adaptive
if (method == "50_50") {
  diff <- 0.5
  } else {
    
    prob <- beta1 |>
      mutate(treatment_benefit = ifelse(beta1 < 0, 1, 0)) |>
      summarize(prop = mean(treatment_benefit)) |>
      pull()
    
    c <- c_param(n_sum)
    # diff <- min(max(prob, 0.1), 0.9)
    diff <- prob^c / (prob^c + (1-prob)^c)
  }

prob_vec(i) <- diff
treatment_num(i) = sum(x)

## Each Sampling is 50 patients
num_vec(i+1) <- n
x <- rbinom(n, 1, diff)
}
  

## Sampling from population
df_list <- df |>
  slice_sample(n = n) |>
  bind_cols(x=x) |>
  mutate(y = case_when(
    x == 1 ~ y1,
    x == 0 ~ y0
  ))

## Bayesian Model

## main model
stan_model <- glue("
data {{
  int N;  
  array(N) int x;  
  array(N) int y;  
}}
parameters {{
  real beta0;  
  real beta1;  
}}
model {{
  beta0 ~ normal({b0mu},{b0sd});
  beta1 ~ normal({b1mu},{b1sd});
  y ~ bernoulli_logit(beta0 + beta1*to_vector(x));
}}
")
  
## compile model
mod <- write_stan_file(stan_model)
model <- cmdstan_model(mod)
fit <- model$sample(
  data = list(N=nrow(df_list), x=df_list$x, y=df_list$y),
  chains = 4, 
  iter_sampling = 2000,
  iter_warmup = 1000,
  seed = 1,
  parallel_chains = 4
)

## Remove patients who are already sampled from the population
sample <- df_list |> pull(id)
df <- df |> filter(!id %in% sample)
x_vec <- c(x_vec,x)
print(i)

## Extract MCMC
mcmc <- as.data.frame(fit$draws(inc_warmup = F))

## Assess probability that treatment effect is < 0 
beta1 <- mcmc |>
  select(contains("beta1")) |>
  pivot_longer(cols = everything(), names_to = "column", values_to = "beta1") 

b1_lower_i <- beta1 |>
  summarize(lower = quantile(beta1, 0.025)) |>
  pull(lower)

b1_upper_i <- beta1 |>
  summarize(upper = quantile(beta1, 0.975)) |>
  pull(upper)


}


## log results
result <- result |>
  bind_rows(tibble(b0=b0,b0_lower=b0_lower,b0_upper=b0_upper,b1=b1,b1_lower=b1_lower,b1_upper=b1_upper,prob=prob_vec,treatment_num=treatment_num, num_sample=num_vec(1:(length(num_vec)-1))) |> 
  mutate(treatment_num_cum = cumsum(treatment_num),
         num = cumsum(num_sample),
         method = method,
         seed = seed) |>
  slice_head(n=30))


}
}

Interpretare

Vizualizare

## Visualize the result
result |>
  ggplot(aes(x=num,y=b1)) +
  geom_point() +
  geom_line() +
  geom_ribbon(aes(ymin=b1_lower,ymax=b1_upper), alpha=0.5) +
  geom_hline(yintercept = 0) +
  geom_hline(yintercept = treatment_effect, color = "blue") +
  geom_label(aes(x=num,y=b1,label=treatment_num_cum), size=3) +
  ggtitle("Reponse Adaptive Randomization vs 50-50 Fixed Randomization", subtitle = "Blue line = True value, Black line = No different between treatment & placebo, Numbers labeled = Number of treatment") +
  xlab("Patients") +
  theme_bw() +
  xlim(40,210) +
  facet_grid(method~seed)

Graficul de mai sus prezintă parametrul coeficientului beta1 din distribuția posterioară. Cele 10 coloane sunt semințele pentru fiecare set de încercare. Rândul reprezintă metoda de alocare, fie 50-50, fie randomizare adaptativă. Linia albastră este adevăratul efect de tratament, care este -1.09. Linia gri este atunci când nu există niciun efect (0), iar zona umbrită este intervalul credibil de 95%. Simplu euristic este că, dacă zona gri coboară sub linia neagră, atunci înseamnă că există 95% probabilitatea de a exista un efect de tratament. Acum să vedem dacă putem răspunde la întrebările pe care le -am avut.

Va avea nevoie de mai puțini pacienți pentru a arăta un efect de tratament pozitiv?

Nu pare ca, atât 50-50, cât și RAR par să aibă un număr similar de pacienți care să arate un efect de tratament pozitiv.

RAR va avea aceeași precizie la identificarea efectului tratamentului?

Da. Modelele de pe ambele metode par a fi destul de similare.

Dacă nu există niciun efect, Rar va putea să -l tachineze? Și cum ar arăta asta?

Pot fi? Atât 50-50, cât și randomizarea adpativă par să fie capabile să identifice că nu există niciun efect de tratament. Cu toate acestea, alocarea 50-50 pare a fi mai stabilă decât randomizarea adaptativă. Acest lucru este doar din 10 încercări, greu de spus. De asemenea, nu am stabilit criterii pentru oprirea inutilității.

Pe o notă laterală, observați cum cu sămânță 26561 şi 87907vom concluziona în mod fals că există un efect de tratament atunci când nu există.

Limitări

Conform rezistenței tentației randomizării adaptive a răspunsului, randomizarea adaptativă a răspunsului (RAR) provoacă multe probleme, inclusiv:

Prejudiciul de la tendințele temporale în studiile clinice atunci când pacienții înscriși în diferite momente se confruntă cu condiții sistematic diferite.
Ineficiența estimării efectului de tratament care necesită adesea dimensiuni mai mari ale eșantionului pentru a menține puterea statistică.
Volatilitatea distribuțiilor de mărime a eșantionului care, în mod paradoxal, pot atribui mai mulți pacienți tratamentelor inferioare datorită variației aleatorii.
Dificultatea de a analiza în mod valabil rezultatele atunci când dimensiunile eșantionului în sine conțin informații despre eficacitatea tratamentului.
Potențialul prejudecății de selecție atunci când cercetătorii devin neintenționat neinpozate la rezultatele provizorii pe măsură ce probabilitățile de alocare se schimbă.

Ce zici de ajustarea multiplicității și eroarea de tip 1?

Proiectele adaptive bayesiene se confruntă cu o dilemă practică în ceea ce privește ajustarea multiplicității. În timp ce teoretic, inferența bayesiană nu necesită corecții pentru mai multe analize ale datelor, cerințele de reglementare de obicei solicită controlul de eroare de tip I pentru studiile de confirmare, indiferent de abordarea statistică. Studiile de caz arată că analizele provizorii neajustate, cu eficacitatea precoce, oprirea ratelor de eroare de tip I, în timp ce oprirea doar inutilitatea scade eroarea de tip I, dar reduce puterea. Pentru cercetătorii care implementează proiectele adaptative bayesiene cu oprirea eficacității timpurii, ajustările la granițe devin necesare pe măsură ce analizele cresc, creând tensiune între filozofia bayesiană și cerințele de reglementare, în special în studiile care vizează schimbarea practicii clinice.

Citește mai mult: Trebuie să ne adaptăm pentru analize provizorii într -un design de încercare adaptiv Bayesian?

Oportunitate de îmbunătățire

Am făcut doar 10 încercări în total, ceea ce nu este suficient pentru a trage concluzii concrete.
Distribuție Cauchy pentru coeficientul de regresie logistică Citește mai mult
Setați frânghia (regiunea de echivalență practică) pentru efectul de tratament și vedeți dacă putem face o idee mai bună despre câți pacienți trebuie să probăm înainte de a putea concluziona că nu există niciun efect de tratament.

Îl trimisem mai devreme pe prietenul meu Alec Wong scenariul. Iată comentariile pentru care să mă îmbunătățesc pentru referințe viitoare.

Extrageți funcționalitatea în funcții dedicate, în mod ideal, mici, concise și extrem de descriptive a ceea ce face. Scrieți -le pentru a evita cât mai mult variabilele globale și, în schimb, să fie explicit despre ce variabile depinde funcția și să le transmiteți.
O mulțime de linii care sunt independente de bucla „metode”, dar totuși apar în interiorul buclei. Orice lucru care nu depinde de buclă nu ar trebui să aparțină în buclă. Modificați DF de date simulate în timpul buclei, da, dar acest lucru nu trebuie să fie cazul, puteți modifica la o nouă variabilă în loc să suprascrieți DF.
Bucla „Metoda” este de fapt destul de mică și nu face prea mult. Se pare că nu există prea multe motive pentru a scrie scenariul, astfel încât totul să fie rulat de două ori.
Vă reutilizați de multe ori în acest script, ceea ce face dificil să argumentați ce valoarea lui N este la un moment dat. Dacă urmați sfaturile din (1), puteți reutiliza N în funcții în interior fără ambiguitate, deoarece acestea sunt obținute la funcție.

Mulțumesc Alec! Voi lucra la acestea!

Lecții învățate

Reîmprospătat pe statisticile noastre Bayesiene și Codul Stan
Am învățat despre RAR și cum arată în comparație cu alocarea 50-50
Căutare a literaturii cu privire la întrebarea privind multiplicitatea și statisticile bayesiene
Distribuție Cauchy pentru coeficientul de regresie logistică Citește mai mult

Dacă vă place acest articol:

Simularea unei simple răspunsuri randomizare adaptativă – trebuie să o văd pentru a -l crede