¿Cuál es la mejor manera de estimar el efecto promedio del tratamiento en un estudio longitudinal?

En un estudio longitudinal, los resultados de las unidades miden repetidamente en los puntos de tiempo con un total de ocasiones de medición fija (fija = las mediciones en unidades se toman al mismo tiempo). $Y_{it}$ $i$ $t$ $m$

Las unidades se asignan aleatoriamente a un tratamiento, , o a un grupo de control, . Quiero estimar y probar el efecto promedio del tratamiento, es decir, donde las expectativas se toman a través del tiempo y de las personas. Considero usar un modelo multinivel de ocasión fija (efectos mixtos) para este propósito: $G=1$ $G=0$

A T E = E (Y | G = 1) - E (Y | G = 0),

$ATE=E(Y | G=1) - E(Y | G=0),$

Y_{i t} = α + β G_{i} + u_{0 i} + e_{i t}

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

con la intersección, la , una intersección aleatoria a través de las unidades, y el residual. $\alpha$ $\beta$ $ATE$ $u$ $e$

Ahora estoy considerando un modelo alternativo.

Y_{i t} = \tilde{β} G_{i} + \sum_{j = 1}^{m} κ_{j} d_{i j} + \sum_{j = 1}^{m} γ_{j} d_{i j} G_{i} + {\tilde{u}}_{0 i} + {\tilde{e}}_{i 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}$

$\kappa_j$ $t$ $d_t=1$ $j=t$ $0$ $\gamma$ $G$ $Y$

$\tilde{\beta}$ $ATE$ $t=1$ $\tilde{\beta}$ $\beta$ $ATE$

Mis preguntas son :

¿Cuál es la mejor manera de estimar el efecto del tratamiento en este diseño de estudio longitudinal?
¿Tengo que usar el modelo 1 o hay alguna forma de usar el modelo 2 (quizás más eficiente)?
$\tilde{\beta}$ $ATE$ $\gamma$

mixed-model panel-data multilevel-analysis random-effects-model fixed-effects-model tomka
fuente

\tilde{β}

$\tilde{\beta}$

γ_{j}

$\gamma_j$

γ

$\gamma$

@jujae La razón principal para el modelo 2 es la reducción de la varianza, sí. Pero me pregunto cómo sacar el ATE del modelo 2. Su primer comentario parece ser un puntero. ¿Puedes mostrar esto o elaborar? ¡Entonces esto sería una respuesta a mi pregunta!

tomka

\tilde{β}

$\tilde{\beta}$

γ_{j}

$\gamma_j$

j

$j$

j

$j$

\tilde{β} + γ_{j}

$\tilde{\beta} + \gamma_j$

β

$\beta$

Respuestas:

Al responder a su pregunta "Me pregunto cómo sacar el ATE del modelo 2" en los comentarios:

$\gamma_j$ $\gamma_j =0$ $j=1$ $\tilde{\beta}$ $m$ $\tilde{\beta}$ $m$

$\gamma_j$ $j$ $\mathrm{ATE}_j$ $\gamma_j= \mathrm{ATE}_j - \mathrm{ATE}_1$ $j=2,\dots, m$ $\mathrm{ATE}_j$ $\tilde{\beta} + \gamma_j$ $\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)$

Hice una simulación simple en R para verificar esto:

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)

Y los resultados lo verifican:

  ATE.m1   ATE.m2 
3.549213 3.549213

No sé cómo cambiar directamente la codificación de contraste en el modelo 2 anterior, así que para ilustrar cómo se puede usar directamente una función lineal de los términos de interacción, así como cómo obtener el error estándar, utilicé el paquete 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)

Y aquí está el resultado:

 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)

$\sqrt{w \hat{V} w^T}$ $w$ $V$

Codificación de desviación

$\tilde{\beta}$ $\mathrm{ATE}$ $\mathrm{ATE}_j - \mathrm{ATE}$

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)

Salida:

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
fuente

\tilde{β}

$\tilde{\beta}$ beta_t

@tomka, es posible, no sé cómo cambiar directamente la matriz de contraste del término de interacción en el modelo 2, investigaremos y volveremos más tarde.

jujae

Pensando en tu respuesta, encontré esto. Creo que la codificación de desviación hace lo que quiero. Puede probarlo e incluirlo en su respuesta. ats.ucla.edu/stat/sas/webbooks/reg/chapter5/…

tomka

@tomka: Eso es exactamente lo que tengo en mente, vea mi comentario original a su pregunta donde mencioné la codificación de desviación :), intentaré implementar esto y actualizar la respuesta más tarde. (Tener algunos problemas para hacerlo en R sin crear manualmente una variable ficticia para la codificación, pero parece que es la única forma de hacerlo).

jujae

@tomka: perdón por el retraso, actualicé la parte del código de desviación

jujae

Para la primera pregunta, entiendo que las formas "elegantes" solo son necesarias cuando no es inmediatamente obvio que el tratamiento es independiente de los resultados potenciales. En estos casos, debe argumentar que algún aspecto de los datos permite una aproximación de la asignación aleatoria al tratamiento, lo que nos lleva a variables instrumentales, discontinuidad de regresión, etc.

En su caso, las unidades se asignan aleatoriamente al tratamiento, por lo que parece creíble que el tratamiento es independiente de los posibles resultados. Entonces podemos mantener las cosas simples: estimar el modelo 1 con mínimos cuadrados ordinarios, y usted tiene una estimación consistente del ATE. Dado que las unidades se asignan aleatoriamente al tratamiento, este es uno de los pocos casos en los que una suposición de efectos aleatorios es creíble.

Peter
fuente