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

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

Teraz rozważam alternatywny model

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)?

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

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.