Quelle est la meilleure façon d'estimer l'effet moyen du traitement dans une étude longitudinale?

Dans une étude longitudinale, les résultats des unités sont mesurés à plusieurs reprises aux points temporels avec un total de occasions de mesure fixes (fixe = les mesures sur les unités sont prises en même temps). $Y_{it}$ $i$ $t$ $m$

Les unités sont assignées au hasard soit à un traitement, , soit à un groupe témoin, . Je veux estimer et tester l'effet moyen du traitement, c'est-à-dire où les attentes sont prises à travers le temps et les individus. J'envisage d'utiliser à cet effet un modèle multiniveau à effets fixes (effets mixtes): $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}$

avec l'ordonnée à l'origine, l' , une interception aléatoire entre les unités, et le résiduel. $\alpha$ $\beta$ $ATE$ $u$ $e$

J'envisage maintenant un modèle alternatif

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}$

qui contient les effets fixes pour chaque occasion où factice si et sinon. De plus, ce modèle contient une interaction entre le traitement et le temps avec les paramètres . Ce modèle prend donc en compte le fait que l'effet de peut différer dans le temps. C'est informatif en soi, mais je pense que cela devrait aussi augmenter la précision d'estimation des paramètres, car l'hétérogénéité de est prise en compte. $\kappa_j$ $t$ $d_t=1$ $j=t$ $0$ $\gamma$ $G$ $Y$

Cependant, dans ce modèle, le coefficient ne semble plus égaler l' . Au lieu de cela, il représente l'ATE à la première occasion ( ). Donc, l'estimation de peut être plus efficace que mais elle ne représente plus l' . $\tilde{\beta}$ $ATE$ $t=1$ $\tilde{\beta}$ $\beta$ $ATE$

Mes questions sont :

Quelle est la meilleure façon d'estimer l'effet du traitement dans ce plan d'étude longitudinale?
Dois-je utiliser le modèle 1 ou existe-t-il un moyen d'utiliser (peut-être plus efficace) le modèle 2?
Existe-t-il un moyen pour que ait l'interprétation de l' et l'écart spécifique à l'occasion (par exemple en utilisant le codage d'effet)? $\tilde{\beta}$ $ATE$ $\gamma$

mixed-model panel-data multilevel-analysis random-effects-model fixed-effects-model tomka
la source

Dans le modèle 2, l'ATE n'est-il pas égal à plus la moyenne de ?

\tilde{β}

$\tilde{\beta}$

γ_{j}

$\gamma_j$

jujae

Si votre objectif est d'estimer exclusivement l'ETA, le modèle 1 suffira, car il ne sera pas biaisé. L'ajout de période ou d'interaction dans le modèle réduira la variance de votre estimation, je crois. Et je pense que vous voudrez peut-être essayer de coder comme codage d' écart (écart par rapport à la moyenne)?

γ

$\gamma$

jujae

@jujae La raison principale du modèle 2 est la réduction de la variance, oui. Mais je me demande comment sortir l'ATE du modèle 2. Votre premier commentaire semble être un pointeur. Pouvez-vous montrer cela ou élaborer? Ce serait alors proche d'une réponse à ma question!

tomka

Lorsque vous ajustez le modèle 2, a l'interprétation d'ATE dans la période 1. Les coefficients du terme d'interaction, pour la considération de l'identifiablilité, seront codés avec ATE à la période 1 comme niveau de référence. Par conséquent, est en fait la différence entre le traitement à la période et le traitement à la période 1 à partir de la sortie du logiciel. Donc, à chaque période , l'ATE est et lorsque la moyenne de l'ATE spécifique à la période, cela conduira à l'ATE moyen, qui est dans votre modèle 1.

\tilde{β}

$\tilde{\beta}$

γ_{j}

$\gamma_j$

j

$j$

j

$j$

\tilde{β} + γ_{j}

$\tilde{\beta} + \gamma_j$

β

$\beta$

jujae

Réponses:

Répondant à votre question "Je me demande comment sortir l'ATE du modèle 2" dans les commentaires:

Tout d'abord, dans votre modèle 2, tous les sont pas identifiables, ce qui pose le problème du manque de rang dans la matrice de conception. Il est nécessaire de supprimer un niveau, par exemple en supposant que pour . Autrement dit, en utilisant le codage de contraste et en supposant que l'effet du traitement à la période 1 est 0. Dans R, il codera le terme d'interaction avec l'effet du traitement à la période 1 comme niveau de référence, et c'est également la raison pour laquelle a l'interprétation de l'effet du traitement à la période 1. Dans SAS, il codera l'effet du traitement à la période comme niveau de référence, puis a l'interprétation de l'effet du traitement à la période $\gamma_j$ $\gamma_j =0$ $j=1$ $\tilde{\beta}$ $m$ $\tilde{\beta}$ $m$ , plus la période 1.

En supposant que le contraste est créé de la manière R, les coefficients estimés pour chaque terme d'interaction (je désignerai toujours ceci par , bien que ce ne soit pas précisément ce que vous avez défini dans votre modèle) ont l'interprétation de la différence d'effet de traitement entre les périodes et période 1. Noter ATE à chaque période , puis pour . Par conséquent, un estimateur pour $\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)$

J'ai fait une simulation simple en R pour vérifier cela:

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)

Et les résultats le vérifient:

  ATE.m1   ATE.m2 
3.549213 3.549213

Je ne sais pas comment changer directement le codage de contraste dans le modèle 2 ci-dessus, donc pour illustrer comment on peut directement utiliser une fonction linéaire des termes d'interaction, ainsi que comment obtenir l'erreur standard, j'ai utilisé le package 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)

Et voici la sortie:

 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$

Codage des écarts

$\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)

Production:

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
la source

\tilde{β}

$\tilde{\beta}$ beta_t

@tomka, il est possible, je ne sais pas comment changer directement la matrice de contraste du terme d'interaction dans model2, fera des recherches et reviendra plus tard.

jujae

En réfléchissant à votre réponse, j'ai trouvé cela. Je pense que le codage de déviation fait ce que je veux. Vous pouvez le tester et l'inclure dans votre réponse. ats.ucla.edu/stat/sas/webbooks/reg/chapter5/…

tomka

@tomka: C'est exactement ce que je pense, voir mon commentaire d'origine à votre question où j'ai mentionné le codage de l'écart :), je vais essayer de l'implémenter et de mettre à jour la réponse plus tard. (Avoir du mal à le faire dans R sans créer manuellement une variable fictive pour le codage, mais il semble que c'est la seule façon de le faire).

jujae

@tomka: désolé pour le retard, mise à jour de la partie code de déviation

jujae

Pour la première question, je crois comprendre que des moyens «fantaisistes» ne sont nécessaires que lorsqu'il n'est pas immédiatement évident que le traitement est indépendant des résultats potentiels. Dans ces cas, vous devez faire valoir que certains aspects des données permettent une approximation de l'assignation aléatoire au traitement, ce qui nous amène aux variables instrumentales, à la discontinuité de la régression, etc.

Dans votre cas, les unités sont assignées au hasard au traitement, il semble donc crédible que le traitement soit indépendant des résultats potentiels. Ensuite, nous pouvons simplement garder les choses simples: estimez le modèle 1 avec les moindres carrés ordinaires, et vous avez une estimation cohérente de l'ATE. Étant donné que les unités sont assignées au hasard au traitement, il s'agit de l'un des rares cas où une hypothèse d'effets aléatoires est crédible.

Peter
la source