Как лучше всего оценить средний эффект лечения в продольном исследовании?

В продольном исследовании результаты $Y_{it}$ единиц $i$ многократно измеряются в моменты времени $t$ с общим числом фиксированных измерений $m$ (фиксированные = измерения в единицах измерения проводятся одновременно).

Единицы случайным образом назначаются либо на лечение, $G=1$ , либо на контрольную группу, $G=0$ . Я хочу оценить и проверить средний эффект лечения, то есть

A T Е знак равно Е (Y | г знак равно 1) - Е (Y | г знак равно 0),

$ATE=E(Y | G=1) - E(Y | G=0),$ где ожидания принимаются во времени и у отдельных лиц. Я рассматриваю использование для этой цели многоуровневой модели (смешанных эффектов) с фиксированной вероятностью:

Y_{я T} знак равно α + β г_{я} + U_{0 я} + е_{я T}

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

с $\alpha$ ; отсекаемый отрезок, $\beta$ ; в $ATE$ , $u$ случайного перехват через единица, а $e$ остаточные.

Сейчас я рассматриваю альтернативную модель

Y_{я T} знак равно \tilde{β} г_{я} + Σ_{J знак равно 1}^{м} κ_{J} d_{я J} + Σ_{J знак равно 1}^{м} γ_{J} d_{я J} г_{я} + {\tilde{U}}_{0 я} + {\tilde{е}}_{я 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}$

который содержит фиксированные эффекты для каждого случая где dummy если и противном случае. Кроме того, эта модель содержит взаимодействие между лечением и временем с параметрами . Таким образом, эта модель учитывает, что эффект может отличаться во времени. Это само по себе информативно, но я считаю, что это также должно повысить точность оценки параметров, поскольку учитывается неоднородность по $\kappa_j$ $t$ $d_t=1$ $j=t$ $0$ $\gamma$ $G$ $Y$

Тем не менее, в этой модели коэффициент , кажется, не равна больше. Вместо этого он представляет ATE в первом случае ( ). Таким образом, оценка может быть более эффективным , чем , но это не представляет больше. $\tilde{\beta}$ $ATE$ $t=1$ $\tilde{\beta}$ $\beta$ $ATE$

Мои вопросы :

Как лучше всего оценить эффект лечения в этом продольном дизайне исследования?
Нужно ли использовать модель 1 или есть способ использовать (возможно, более эффективную) модель 2?
Есть ли способ, чтобы интерпретировало и конкретное отклонение от случая (например, используя кодирование эффекта)? $\tilde{\beta}$ $ATE$ $\gamma$

mixed-model panel-data multilevel-analysis random-effects-model fixed-effects-model Томка
источник

В модели 2 разве ATE не равно

плюс среднее значение

\tilde{β}

$\tilde{\beta}$

γ_{j}

$\gamma_j$

jujae

Если вашей целью является исключительно оценка ATE, тогда достаточно модели 1, поскольку она будет беспристрастной. Я считаю, что добавление периода или взаимодействия в модель уменьшит дисперсию вашей оценки. И я думаю, что вы можете попытаться закодировать

как кодирование отклонения (отклонение от среднего)?

γ

$\gamma$

jujae

@jujae Основная причина для модели 2 - уменьшение дисперсии, да. Но мне интересно, как вывести ATE из модели 2. Ваш первый комментарий кажется указателем. Вы можете показать это или уточнить? Тогда это будет близко к ответу на мой вопрос!

Томка

При подборе модели 2,

имеет интерпретацию АТЕ в период 1. Коэффициенты члена взаимодействия, для рассмотрения identifiablility, будет закодирован с АТЕ на период 1 в качестве опорного уровня. Следовательно,

фактически является разницей между обработкой в периоде

и обработкой в периоде 1 из выходных данных программного обеспечения. Таким образом, в каждый период

ATE составляет

\tilde{β}

$\tilde{\beta}$

γ_{j}

$\gamma_j$

j

$j$

j

$j$

и при усреднении ATE для конкретного периода это приведет к большому среднему значению ATE, которое равно

в вашей модели 1.

\tilde{β} + γ_{j}

$\tilde{\beta} + \gamma_j$

β

$\beta$

jujae

Ответы:

Отвечая на ваш вопрос «Интересно, как вывести ATE из модели 2» в комментариях:

Прежде всего, в вашей модели 2 не все идентифицируемы, что приводит к проблеме дефицита ранга в матрице проектирования. Необходимо отбросить один уровень, например, предполагая, что для . То есть, используя контрастное кодирование и предполагая, что эффект лечения в период 1 равен 0. В R он будет кодировать термин взаимодействия с эффектом лечения в период 1 в качестве контрольного уровня, и это также является причиной, по которой интерпретирует эффекта лечения в период 1. В SAS он будет кодировать эффект лечения в период в качестве контрольного уровня, затем $\gamma_j$ $\gamma_j =0$ $j=1$ $\tilde{\beta}$ $m$ $\tilde{\beta}$ имеет интерпретацию эффекта лечения в период , а не период 1 больше. $m$

Если предположить, что контраст создается способом R, то коэффициенты, оцененные для каждого члена взаимодействия (я все равно буду обозначать это как , хотя это не совсем то, что вы определили в вашей модели), интерпретируют разницу эффекта лечения между периодами времени. и период времени 1. Обозначим ATE в каждом периоде , тогда для $\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)$ .

Я сделал простое моделирование в R, чтобы проверить это:

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)

И результаты подтверждают это:

  ATE.m1   ATE.m2 
3.549213 3.549213

Я не знаю, как напрямую изменить контрастное кодирование в модели 2 выше, поэтому, чтобы проиллюстрировать, как можно напрямую использовать линейную функцию членов взаимодействия, а также как получить стандартную ошибку, я использовал пакет 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)

И вот вывод:

 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$ предполагаемой ковариационной матрицы коэффициентов из модели 3.

Кодирование отклонения

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

Вывод:

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
источник

\tilde{β}

$\tilde{\beta}$ beta_t

@ Tomka, возможно, я не знаю, как напрямую изменить контрастную матрицу члена взаимодействия в model2, позже мы проведем некоторые исследования и вернемся.

Jujae

Размышляя о вашем ответе, я нашел это. Я думаю, что кодирование отклонений делает то, что я хочу. Вы можете проверить это и включить в свой ответ. ats.ucla.edu/stat/sas/webbooks/reg/chapter5/...

Томка

@ Tomka: Это именно то, что у меня на уме, см. мой оригинальный комментарий к вашему вопросу, где я упомянул кодирование отклонений :), я постараюсь реализовать это и обновить ответ позже. (Возникли проблемы с выполнением этого в R без ручного создания фиктивной переменной для кодирования, но, похоже, это единственный способ сделать это).

jujae

@tomka: извините за задержку, обновил часть кода отклонения

jujae

Что касается первого вопроса, я понимаю, что «причудливые» способы нужны только тогда, когда не сразу очевидно, что лечение не зависит от потенциальных результатов. В этих случаях вы должны утверждать, что некоторые аспекты данных позволяют приблизиться случайное назначение к лечению, что приводит нас к инструментальным переменным, разрыву регрессии и так далее.

В вашем случае, единицы будут рандомизированы на лечение, так что кажется правдоподобным , что лечение не зависит от возможных исходов. Тогда мы можем упростить задачу: оцените модель 1 с помощью обычных наименьших квадратов, и вы получите непротиворечивую оценку ATE. Поскольку единицы назначаются для лечения случайным образом, это один из немногих случаев, когда предположение о случайных эффектах правдоподобно.

Питер
источник