Jaki jest najlepszy sposób oszacowania średniego efektu leczenia w badaniu podłużnym?

W badaniu podłużnym wyniki $Y_{it}$ jednostek $i$ są wielokrotnie mierzone w punktach czasowych $t$ z łącznie $m$ ustalone okazje pomiarowe (ustalone = pomiary jednostek są wykonywane w tym samym czasie).

Jednostki są losowo przypisywane do leczenia, $G=1$ lub do grupy kontrolnej, $G=0$ . Chcę oszacować i przetestować średni efekt leczenia, tj

ZA T. mi = mi (Y | sol = 1) - mi (Y | sol = 0),

$ATE=E(Y | G=1) - E(Y | G=0),$ gdzie oczekiwania są uwzględniane w czasie i osobach. Rozważam zastosowanie do tego celu modelu wielopoziomowego (mieszanego) o ustalonej okazji:

Y_{ja t} = α + β {sol}_{ja} + u_{0 ja} + {mi}_{ja t}

$Y_{it} = \alpha + \beta G_i + u_{0i} + e_{it}$

z $\alpha$ przechwycenie, $\beta$ $ATE$ , $u$ losowe przechwytywanie między jednostkami, oraz $e$ resztkowe.

Teraz rozważam alternatywny model

Y_{ja t} = \tilde{β} {sol}_{ja} + \sum_{jot = 1}^{m} κ_{jot} {re}_{ja jot} + \sum_{jot = 1}^{m} γ_{jot} {re}_{ja jot} {sol}_{ja} + {\tilde{u}}_{0 ja} + {\tilde{mi}}_{ja t}

$Y_{it} = \tilde{\beta} G_i + \sum_{j=1}^m \kappa_j d_{ij} + \sum_{j=1}^m \gamma_j d_{ij} G_i + \tilde{u}_{0i} + \tilde{e}_{it}$

który zawiera ustalone efekty $\kappa_j$ na każdą okazję $t$ gdzie manekin $d_t=1$ gdyby $j=t$ i $0$ jeszcze. Ponadto model ten zawiera interakcję między leczeniem a czasem z parametrami $\gamma$ . Tak więc ten model bierze pod uwagę efekt $G$ mogą się różnić w czasie. Jest to samo w sobie pouczające, ale uważam, że powinno to również zwiększyć precyzję szacowania parametrów, ponieważ niejednorodność w $Y$ jest brany pod uwagę.

Jednak w tym modelu $\tilde{\beta}$ współczynnik wydaje się nie być równy $ATE$ już. Zamiast tego reprezentuje ATE przy pierwszej okazji ( $t=1$ ). Więc oszacowanie $\tilde{\beta}$ może być bardziej wydajny niż $\beta$ ale to nie reprezentuje $ATE$ już.

Moje pytania to :

Jaki jest najlepszy sposób oszacowania efektu leczenia w tym projekcie badania podłużnego?
Czy muszę korzystać z modelu 1, czy jest sposób na wykorzystanie (być może bardziej wydajnego) modelu 2?
Czy jest na to sposób $\tilde{\beta}$ mieć interpretację $ATE$ i $\gamma$ odchylenie specyficzne dla okazji (np. przy użyciu kodowania efektów)?

— tomka
źródło

W modelu 2 ATE nie jest równe

\tilde{β}

$\tilde{\beta}$ plus średnia z

γ_{j}

$\gamma_j$ ?

— jujae

Jeśli Twoim celem jest wyłącznie oszacowanie ATE, wystarczy model 1, ponieważ będzie on bezstronny. Dodanie okresu lub interakcji w modelu zmniejszy wariancję twojego oszacowania. I myślę, że możesz spróbować kodować

γ

$\gamma$ jako kodowanie odchyleń (odchylenie od średniej)?

— jujae

@jujae Głównym powodem dla modelu 2 jest redukcja wariancji, tak. Ale zastanawiam się, jak wyciągnąć ATE z modelu 2. Twój pierwszy komentarz wydaje się być wskaźnikiem. Możesz to pokazać lub rozwinąć? To byłaby odpowiedź na moje pytanie!

— tomka

Kiedy dopasujesz model 2,

\tilde{β}

$\tilde{\beta}$ interpretuje ATE w okresie 1. Współczynniki terminu interakcji, dla celów identyfikacji, będą kodowane ATE w okresie 1 jako poziom odniesienia. W związku z tym

γ_{j}

$\gamma_j$ jest faktycznie różnicą między leczeniem w danym okresie

j

$j$ i leczenie w okresie 1 od produkcji oprogramowania. Tak więc w każdym okresie

j

$j$ , ATE jest

\tilde{β} + γ_{j}

$\tilde{\beta} + \gamma_j$ a gdy uśredni ATE dla danego okresu, doprowadzi to do wielkiej średniej ATE, którą jest

β

$\beta$ w twoim modelu 1.

— jujae

Odpowiedzi:

Odpowiadając na twoje pytanie „Zastanawiam się, jak wyciągnąć ATE z modelu 2” w komentarzach:

Po pierwsze, w twoim modelu 2 nie wszystkie $\gamma_j$ jest identyfikowalny, co prowadzi do problemu niedoboru rang w matrycy projektowej. Konieczne jest opuszczenie jednego poziomu, na przykład zakładanie $\gamma_j =0$ dla $j=1$ . Oznacza to, że stosując kodowanie kontrastu i zakładając, że efekt leczenia w okresie 1 wynosi 0. W R koduje termin interakcji z efektem leczenia w okresie 1 jako poziom odniesienia, i to jest również powód, dla którego $\tilde{\beta}$ ma interpretację efektu leczenia w okresie 1. W SAS koduje efekt leczenia w okresie $m$ jako poziom odniesienia $\tilde{\beta}$ ma interpretację efektu leczenia w danym okresie $m$ , nie kropka już 1.

Zakładając, że kontrast jest tworzony w sposób R, wówczas współczynniki szacowane dla każdego składnika interakcji (nadal będę to oznaczać przez $\gamma_j$ , chociaż nie jest to dokładnie to, co zdefiniowałeś w swoim modelu) ma interpretację różnicy efektu leczenia między okresem $j$ i okres 1. Oznacz ATE w każdym okresie $\mathrm{ATE}_j$ , następnie $\gamma_j= \mathrm{ATE}_j - \mathrm{ATE}_1$ dla $j=2,\dots, m$ . Dlatego estymator dla $\mathrm{ATE}_j$ jest $\tilde{\beta} + \gamma_j$ . (ignorując różnicę notacji między prawdziwym parametrem a samym estymatorem, ponieważ lenistwo) I oczywiście twoje $\mathrm{ATE}=\beta=\frac{1}{m} \sum_{j=1}^m \mathrm{ATE}_j=\frac{\tilde{\beta}+(\tilde{\beta}+\gamma_2)+\cdots+(\tilde{\beta}+\gamma_m)}{m} = \tilde{\beta}+\frac{1}{m}(\gamma_2 + \cdots + \gamma_m)$ .

Wykonałem prostą symulację w R, aby to sprawdzić:

set.seed(1234)
time <- 4
n <-2000
trt.period <- c(2,3,4,5) #ATE=3.5
kj <- c(1,2,3,4)
intercept <- rep(rnorm(n, 1, 1), each=time)
eij <- rnorm(n*time, 0, 1.5)
trt <- rep(c(rep(0,n/2),rep(1,n/2)), each=time)
y <- intercept + trt*(rep(trt.period, n))+rep(kj,n)+eij
sim.data <- data.frame(id=rep(1:n, each=time), period=factor(rep(1:time, n)), y=y, trt=factor(trt))

library(lme4)
fit.model1 <- lmer(y~trt+(1|id), data=sim.data)
beta <- getME(fit.model1, "fixef")["trt1"]

fit.model2 <- lmer(y~trt*period + (1|id), data=sim.data)
beta_t <- getME(fit.model2, "fixef")["trt1"]
gamma_j <- getME(fit.model2, "fixef")[c("trt1:period2","trt1:period3","trt1:period4")]

results <-c(beta, beta_t+sum(gamma_j)/time)
names(results)<-c("ATE.m1", "ATE.m2")
print(results)

Wyniki potwierdzają to:

  ATE.m1   ATE.m2 
3.549213 3.549213

Nie wiem, jak bezpośrednio zmienić kodowanie kontrastu w powyższym modelu 2, więc aby zilustrować, w jaki sposób można bezpośrednio użyć funkcji liniowej terminów interakcji, a także jak uzyskać błąd standardowy, użyłem pakietu multcomp:

sim.data$tp <- interaction(sim.data$trt, sim.data$period)
fit.model3 <- lmer(y~tp+ (1|id), data=sim.data)
library(multcomp)
# w= tp.1.1 + (tp.2.1-tp.2.0)+(tp.3.1-tp.3.0)+(tp.4.1-tp.4.0)
# tp.x.y=interaction effect of period x and treatment y
w <- matrix(c(0, 1,-1,1,-1,1,-1,1)/time,nrow=1)
names(w)<- names(getME(fit.model3,"fixef"))
xx <- glht(fit.model3, linfct=w)
summary(xx)

A oto wynik:

 Simultaneous Tests for General Linear Hypotheses
Fit: lmer(formula = y ~ tp + (1 | id), data = sim.data)
Linear Hypotheses:
       Estimate Std. Error z value Pr(>|z|)    
1 == 0  3.54921    0.05589   63.51   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Myślę, że standardowy błąd jest uzyskiwany przez $\sqrt{w \hat{V} w^T}$ z $w$ będący powyższą formą kombinacji liniowej i $V$ oszacowana macierz wariancji-kowariancji współczynników z modelu 3.

Kodowanie odchyleń

Kolejny sposób na zrobienie $\tilde{\beta}$ mający bezpośrednią interpretację $\mathrm{ATE}$ polega na użyciu kodowania odchylenia , aby później reprezentowały zmienne towarzyszące $\mathrm{ATE}_j - \mathrm{ATE}$ porównanie:

sim.data$p2vsmean <- 0
sim.data$p3vsmean <- 0
sim.data$p4vsmean <- 0
sim.data$p2vsmean[sim.data$period==2 & sim.data$trt==1] <- 1
sim.data$p3vsmean[sim.data$period==3 & sim.data$trt==1] <- 1
sim.data$p4vsmean[sim.data$period==4 & sim.data$trt==1] <- 1
sim.data$p2vsmean[sim.data$period==1 & sim.data$trt==1] <- -1
sim.data$p3vsmean[sim.data$period==1 & sim.data$trt==1] <- -1
sim.data$p4vsmean[sim.data$period==1 & sim.data$trt==1] <- -1


fit.model4 <- lmer(y~trt+p2vsmean+p3vsmean+p4vsmean+ (1|id), data=sim.data)

Wynik:

Fixed effects:
            Estimate Std. Error t value
(Intercept)  3.48308    0.03952   88.14
trt1         3.54921    0.05589   63.51
p2vsmean    -1.14774    0.04720  -24.32
p3vsmean     1.11729    0.04720   23.67
p4vsmean     3.01025    0.04720   63.77

— jujae
źródło

Dobrze - ale jak uzyskać standardowy szacunek błędu? I czy nie powinno być możliwe zastosowanie kodowania efektów interakcji / okresu w taki sposób

\tilde{β}

$\tilde{\beta}$ (twoje beta_t) jest ATE bezpośrednio (wtedy z oszacowaniem SE)?

— tomka

@ tomka, jest możliwe, nie wiem jak bezpośrednio zmienić matrycę kontrastu terminu interakcji w modelu 2, zrobię badania i wrócę później.

— jujae

Myśląc o twojej odpowiedzi, znalazłem to. Myślę, że kodowanie odchyleń robi to, co chcę. Możesz to przetestować i dołączyć do swojej odpowiedzi. ats.ucla.edu/stat/sas/webbooks/reg/chapter5/…

— tomka

@tomka: Dokładnie o tym myślę, zobacz mój oryginalny komentarz do twojego pytania, w którym wspomniałem o kodowaniu odchyleń :) Spróbuję to zaimplementować i zaktualizuję odpowiedź później. (Mam problem z zrobieniem tego w R bez ręcznego tworzenia zmiennej zastępczej do kodowania, ale wygląda na to, że jest to jedyny sposób).

— jujae

@tomka: przepraszam za opóźnienie, zaktualizowałem część kodu odchylenia

— jujae

W przypadku pierwszego pytania rozumiem, że „fantazyjne” sposoby są potrzebne tylko wtedy, gdy nie jest od razu oczywiste, że leczenie jest niezależne od potencjalnych rezultatów. W takich przypadkach należy argumentować, że niektóre aspekty danych pozwalają na przybliżenie losowego przypisania do leczenia, co prowadzi nas do zmiennych instrumentalnych, nieciągłości regresji i tak dalej.

W twoim przypadku jednostki są losowo przydzielane do leczenia, więc wydaje się prawdopodobne, że leczenie jest niezależne od potencjalnych wyników. Następnie możemy po prostu wszystko uprościć: oszacuj model 1 ze zwykłymi najmniejszymi kwadratami, a będziesz miał spójne oszacowanie ATE. Ponieważ jednostki są losowo przydzielane do leczenia, jest to jeden z niewielu przypadków, w których wiarygodne jest założenie o efektach losowych.

— Piotr
źródło