Gândindu -mă la covariate într -o analiză a unui RCT

URMĂREȘTE-NE
16,065FaniÎmi place
1,142CititoriConectați-vă

Am discutat recent despre planul analitic pentru un studiu controlat randomizat (RCT) cu un colaborator clinic când a întrebat dacă este potrivit să se ajusteze pentru covariatele de bază pre-specificate. Această întrebare este atât de interesantă, deoarece abordează problemele fundamentale ale inferenței – atât cauzale, cât și statistice. Care este estimarea țintă într -un RCT – adică ce efect măsurăm de fapt? Ce sperăm să învățăm din eșantionul specific recrutat pentru proces (adică cum pot fi analizate rezultatele într -un mod care să îmbunătățească generalizarea)? Ce presupuneri de bază despre replicabilitate, eșantionare și incertitudine informează argumentele pentru și împotriva ajustării covariate? Acestea sunt mari întrebări, la care nu se vor răspunde neapărat aici, dar trebuie să se țină cont atunci când ne gândim la meritele ajustării covariate

Unii cercetători rezistă la ajustarea covariatei în analiza primară, îngrijorată de faptul că ar putea complica interpretabilitatea sau poate limita transparența. Alții ar putea să -i placă claritatea și simplitatea simplă a comparației randomizate. Dar poate cea mai mare problemă pe care oamenii o au cu ajustarea covariatului este o preocupare de lungă durată că modelarea flexibilă s -ar putea transforma într -o expediție de pescuit – căutarea covariatelor care obțin cea mai favorabilă estimare a efectului.

După acea conversație cu colegul meu, am revizuit o lucrare din 1994 a lui Stephen Senn, care susține că, mai degrabă decât să verific dezechilibrele covariate ale șanselor înainte de a face ajustări, „Statisticistul practic va face bine să stabilească în prealabil o listă limitată de covariate considerate utile și să le potrivească Indiferent. O astfel de strategie va de obicei duce la un câștig de puterenu are efect negativ asupra dimensiunii necondiționate și controlează dimensiunea condițională în raport cu covariatele identificate. ” O lucrare ulterioară de Pocock și colab. consolidează această perspectivă. Deși observă că „experiența arată că pentru majoritatea studiilor clinice, analizele care se ajustează pentru covariatele de bază sunt în strânsă acorduri cu comparațiile de tratament mai simple neajustate”, susțin că ajustarea pentru covariate poate fi justificată dacă ajută: (1) să obțină cel mai potrivit cel mai potrivit Valoarea P pentru diferențele de tratament, (2) oferă estimări nepărtinitoare și (3) îmbunătățesc precizia.

Motivat de Pocock și colab., Am creat câteva simulări pentru a explora caracteristicile operaționale ale ajustării covariate pe care le împărtășesc aici. Am fost distras mai recent cu scrierea pe hârtie și editarea manuscrisului, așa că sunt fericit să mă întorc la un pic R Codificare.

Simulări

Pentru a începe lucrurile, iată R pachete utilizate în simulări.

library(simstudy)
library(data.table)
library(stargazer)
library(parallel)

Definiții de date

Folosesc aici două seturi de definiții de date, împărțind crearea covariatelor de bază ( (x_1 ) şi (x_2 )) și alocarea grupului din rezultat (y ). Presupunem că (y ) Are o distribuție gaussiană (normală).

def_c <- 
  defData(varname = "x1", formula = 0, variance = "..s_x1^2", dist = "normal") |>
  defData(varname = "x2", formula = 0, variance = "..s_x2^2", dist = "normal") |>
  defData(varname = "A", formula = "1;1", dist = "trtAssign")

def_y <- defDataAdd(
  varname = "y", 
  formula = "5 + ..delta * A + ..b1 * x1 + ..b2 * x2", 
  variance = "..s_y^2", 
  dist = "normal"
)

Parametre inițiale

Iată parametrii folosiți în generarea de date. (x_1 ) este foarte corelat cu rezultatul (y )întrucât (x_2 ) nu este. În primul set de simulări, presupunem o dimensiune a efectului adevărat ( delta = 5 ).

s_x1 <- 8
s_x2 <- 9
s_y <- 12
b1 <- 1.50
b2 <- 0.0
delta <- 5

Generare de seturi de date unice

Generatorul de numere recursiv al lui L’Ecuyer (CMRG) combinat multiplu recursiv (CMRG) este utilizat aici, deoarece replicarea datelor (și analizele) se face folosind un proces paralel pentru a accelera puțin lucrurile.

RNGkind("L'Ecuyer-CMRG")
set.seed(55)

dc <- genData(250, def_c) 
dd <- addColumns(def_y, dc)

head(dd)
## Key: 
##       id         x1          x2     A          y
##                        
## 1:     1  3.4075117 -1.66327988     0  8.1461927
## 2:     2 -0.3040474 -5.60073657     0 -1.7577859
## 3:     3 -4.4460516  1.20189340     1 -1.6324336
## 4:     4  0.6834332  0.09974478     1  9.5938736
## 5:     5 -4.6324773  7.85745373     0 -0.4782144
## 6:     6 -6.5650815  0.49812462     0 -3.4082489

Pentru acest set de date unic, putem vedea că mijloacele pentru (x_1 ) în cadrul fiecărui grup sunt ușor diferite, în timp ce mijloacele pentru (x_2 ) sunt mai asemănătoare.

dd(, .(mu_x1 = mean(x1), mu_x2 = mean(x2)), keyby = A)
## Key: 
##        A      mu_x1       mu_x2
##                 
## 1:     0  0.3960784 -0.36510184
## 2:     1 -0.4821337  0.08564994

Aceste diferențe sunt confirmate prin calcularea dezechilibrului standardizat (Z_x ) și diferența standardizată (d ). (Diferența dintre (Z_x ) şi (d ) este asta (Z_x ) are o ajustare pentru dimensiunile eșantionului de grup.)

calc_diff <- function(dx, rx, v) {
  
  mean_diff <- dx(get(rx)==1, mean(get(v))) - dx(get(rx)==0, mean(get(v)))
  s_pooled <- sqrt(
    (dx(get(rx)==1, (.N - 1) * var(get(v))) + 
       dx(get(rx)==0, (.N - 1) * var(get(v))) ) / dx(, .N - 2))
  
  Z_x <- mean_diff / ( s_pooled * sqrt(1/dx(get(rx)==0, .N) + 1/dx(get(rx)==1, .N)) )
  d <- mean_diff / s_pooled
  
  return(list(Z_x = Z_x, d = d))
  
  }

calc_diff(dd, "A", "x1")
## $Z_x
## (1) -0.9424664
## 
## $d
## (1) -0.1192136
calc_diff(dd, "A", "x2")
## $Z_x
## (1) 0.3862489
## 
## $d
## (1) 0.04885705

După cum este proiectat, (x_1 ) este puternic corelat cu rezultatul (y )întrucât (x_2 ) nu este.

dd(, .(rho_x1.y = cor(x1, y), rho_x2.y = cor(x2, y)))
##     rho_x1.y     rho_x2.y
##                
## 1: 0.6807215 -0.003174757

Estimarea modelului

Încasăm patru modele la aceste date: (1) fără ajustare pentru covariate, (2) ajustarea pentru (x_1 ) singur, (3) ajustarea pentru (x_2 ) Singur, și (4) ajustarea pentru ambele covariate.

model_1 <- lm(data = dd, formula = y ~ A)
model_2 <- lm(data = dd, formula = y ~ A + x1)
model_3 <- lm(data = dd, formula = y ~ A + x2)
model_4 <- lm(data = dd, formula = y ~ A + x1 + x2)

Două preluări cheie din acest set de date unice sunt că (1) de când (x_1 ) este un confuzor (deși slab), care nu reușește să se ajusteze pentru covariate duce la o subestimare a efectului de tratament (datorită corelației negative (mici) a (x_1 ) şi (O)) și (2) de când (x_1 ) este atât de foarte corelat cu (y )modelele pentru care se reglează (x_1 ) au erori standard mai mici pentru estimarea efectului de tratament (în jur de 2,0 pentru modelele 1 și 3 și mai aproape de 1,5 pentru modelele 2 și 4).

## 
## ============================================
##            (1)      (2)      (3)      (4)   
## --------------------------------------------
## A         3.040   4.367***  3.044   4.367***
##          (2.039)  (1.480)  (2.044)  (1.484) 
##                                             
## x1                1.512***          1.512***
##                   (0.101)           (0.101) 
##                                             
## x2                          -0.010   0.001  
##                            (0.111)  (0.081) 
##                                             
## Constant 6.237*** 5.638*** 6.234*** 5.639***
##          (1.442)  (1.046)  (1.446)  (1.048) 
##                                             
## ============================================
## ============================================
## 

Caracteristici de operare (bazate pe seturi de date replicate)

Pentru a înțelege meritele relative ale diferitelor abordări de modelare, trebuie să reproducem mai multe seturi de date sub același set de presupuneri utilizate pentru a genera setul de date unic. Vom genera 2000 de seturi de date și vom estima toate cele patru modele pentru fiecare set de date. Pentru fiecare replicare, folosim funcția est_ancova Pentru a calcula o valoare P unilaterală. Vom urmări estimarea punctului, estimarea erorilor standard și valoarea p pentru fiecare iterație.

est_ancova <- function(dx, vars) {

  formula <- as.formula(paste("y ~", paste(vars, collapse = " + ")))
  model <- lm(data = dx, formula = formula)
  
  coef_summary <- summary(model)$coefficients("A", )
  t_stat <- coef_summary("t value")
  
  p_value <- pt(t_stat, df = model$df.residual, lower.tail = FALSE)
  ests <- data.table(t(coef_summary(1:2)), p_value)
  setnames(ests, c("est", "se", "pval"))
  
  return(ests)

}

replicate <- function() {
  
  dc <- genData(250, def_c) 
  dd <- addColumns(def_y, dc)
  
  est_1 <- est_ancova(dd, vars = "A")
  est_2 <- est_ancova(dd, vars = c("A", "x1"))
  est_3 <- est_ancova(dd, vars = c("A", "x2"))
  est_4 <- est_ancova(dd, vars = c("A", "x1", "x2"))
  
  return(list(est_1 = est_1, est_2 = est_2, est_3 = est_3, est_4 = est_4))
  
}
res <- mclapply(1:2000, function(x) replicate())

Toate cele patru modele produc estimări relativ nepărtinitoare, deși modelele care se ajustează (x_1 ) (potențialul confuzor) are ca rezultat o prejudecată redusă în raport cu cele care nu. Cu toate acestea, avantajul clar al modelelor 2 și 4 (cele care se adaptează (x_1 )) este variația redusă a estimatorului efectului de tratament:

get.field <- function(x, field) {
  data.table(t(sapply(x, function(x) x((field))) ))
}

ests <- rbindlist(lapply(res, function(x) get.field(x, "est")))
sapply(ests, function(x) c(bias = mean(x) - delta, var = var(x)))
##            est_1       est_2       est_3         est_4
## bias -0.05513569 -0.00157402 -0.05263476 -0.0002470811
## var   4.51844704  2.33883113  4.55788040  2.3507163302

Reducerea variației se traduce direct la o putere crescută pentru modelele care se ajustează (x_1 )de la aproximativ 63% la 90%. Acest lucru pare un motiv destul de bun pentru a vă ajusta pentru o covariate de bază pe care (a) o colectați și (b) este foarte corelată cu rezultatul.

pvals <- rbindlist(lapply(res, function(x) get.field(x, "pval")))
sapply(pvals, function(x) c(mean(x < 0.025))) 
##  est_1  est_2  est_3  est_4 
## 0.6255 0.9055 0.6275 0.9035

Explorarea ratelor de eroare de tip 1

Partea flip a puterii statistice este eroarea de tip 1 – probabilitatea de a concluziona că există un efect de tratament atunci când, de fapt, nu există niciun efect de tratament. Putem evalua acest lucru prin stabilirea ( delta = 0 ) și rularea unui alt număr mare de replici. Dacă facem acest lucru de 2.000 de ori generând un set de date complet nou de fiecare dată, vedem că ratele de eroare observate sunt apropiate de 0,025 pentru toate modelele, deși modelele care se adaptează (x_1 ) sunt mai aproape de valoarea teoretică.

delta <- 0

res <- mclapply(1:2000, function(x) replicate())

pvals <- rbindlist(lapply(res, function(x) get.field(x, "pval")))
sapply(pvals, function(x) c(mean(x < 0.025))) 
##  est_1  est_2  est_3  est_4 
## 0.0180 0.0265 0.0185 0.0285

Atât Senn, cât și Pocock și colab. Sugerați-vă că un avantaj esențial al ajustării pentru covariatele de bază este faptul că ajută la obținerea ratelor de eroare dorite, în special pentru testele unilaterale. Presupunând că Toate RCT -urile posibile sunt realizate cu același nivel de dezechilibru covariatmodelele care includ ajustări de covariate de bază vor produce rate de eroare exacte. În schimb, modelele care nu se ajustează pentru covariate importante (extrem de corelate) vor produce rate de eroare dezumfliate. Acest lucru apare în principal deoarece erorile standard ale estimărilor efectului sunt supraestimate sistematic, reducând probabilitatea de a respinge vreodată ipoteza nulă.

Pentru a imita cerința ca setul de date să fie eșantionat condiționat de un nivel fix de dezechilibru covariat, generăm covariatele de bază și alocarea tratamentului o singură dată, în timp ce rezultatul este generat din nou pentru fiecare set de date. În cadrul acestei abordări, covariatele și alocarea tratamentului sunt fixate și numai rezultatul pentru o anumită unitate variază de -a lungul iterațiilor. O abordare alternativă ar fi generarea unui număr mare de seturi de date folosind procesul complet de randomizare – crearea de noi valori covariate, alocări de tratament și rezultate pentru fiecare iterație. Seturile de date ar fi analizate numai dacă ar corespunde nivelului de dezechilibru covariat pre-specificat. Deși această abordare produce aceleași rezultate ca metoda ales (am confirmat cu simulările care nu sunt prezentate aici), procesul de eșantionare apare excesiv de artificial, complicând în continuare interpretarea valorii p.

replicate_2 <- function() {
  
  dd <- addColumns(def_y, dc)
  
  est_1 <- est_ancova(dd, vars = "A")
  est_2 <- est_ancova(dd, vars = c("A", "x1"))
  est_3 <- est_ancova(dd, vars = c("A", "x2"))
  est_4 <- est_ancova(dd, vars = c("A", "x1", "x2"))
  
  return(list(est_1 = est_1, est_2 = est_2, est_3 = est_3, est_4 = est_4))
  
}

dc <- genData(250, def_c) 

res <- mclapply(1:2000, function(x) replicate_2())

pvals <- rbindlist(lapply(res, function(x) get.field(x, "pval")))
sapply(pvals, function(x) c(mean(x < 0.025))) 
##  est_1  est_2  est_3  est_4 
## 0.0075 0.0285 0.0070 0.0280

Metode de inferență cauzală pentru echilibrare

Pentru mine, cel mai puternic argument împotriva ajustării pentru covariatele de bază în analiză este riscul ca anchetatorii să poată părea excesiv de dornici să demonstreze succesul intervenției. Pre-specificarea planului de analiză merge mult spre atenuarea acestor preocupări. În plus, abordările alternative din metodele de inferență cauzală pot reduce și mai mult dependența de presupunerile modelului de rezultat. În special, metodele de echilibrare, cum ar fi ponderarea probabilității inverse (IPW) și greutățile suprapuse (OW) pot aborda dezechilibrele covariate, păstrând în același timp estimarea inițială. Aceste tehnici din nou re-înregistrează eșantionul pentru a crea pseudo-populații echilibrate fără a modifica direct modelul de rezultat, oferind o alternativă viabilă la ajustările bazate pe regresie. Acestea au avantajul de a separa modelul de proiectare de modelul de rezultat (deoarece modelele de expunere și rezultat sunt două etape distincte). Echilibrarea se poate face înainte de a analiza datele rezultatelor – deci niciun risc de pescuit pentru rezultate. Planific să împărtășesc simulări folosind aceste abordări cândva în viitor.

Referințe:

Senn, Ștefan. „Testarea echilibrului de bază în studiile clinice.” Statistici în medicină 13, nr. 17 (1994): 1715-1726.

Pocock, Stuart J., Susan E. Assmann, Laura E. Enos și Linda E. Kasten. „Analiza subgrupurilor, ajustarea covariate și comparațiile de bază în raportarea studiilor clinice: practică și probleme curente.” Statistici în medicină 21, nr. 19 (2002): 2917-2930.

Dominic Botezariu
Dominic Botezariuhttps://www.noobz.ro/
Creator de site și redactor-șef.

Cele mai noi știri

Pe același subiect

LĂSAȚI UN MESAJ

Vă rugăm să introduceți comentariul dvs.!
Introduceți aici numele dvs.