0

次の構造のデータセットがあります。

# libraries
library(Zelig) # 5.0-12
library(datatable)

# create data 
time <- factor(rep(-12:12, 50))
treatment <- rbinom(length(time), 1, .75)
outcome <- rnorm(length(time), 1, 3) + 3 * treatment

dat <- data.table(outcome, time, treatment)
dat

            outcome time treatment
   1: 5.2656458  -12         0
   2: 4.8888805  -11         1
   3: 2.6322592  -10         1
   4: 8.2449092   -9         1
   5: 0.5752739   -8         0
  ---                         
1246: 2.1865924    8         0
1247: 1.6028838    9         1
1248: 2.4056725   10         1
1249: 2.0257008   11         1
1250: 6.1503307   12         1

時間と治療を相互作用させる LS モデルを実行します。

z <- zls$new()
z$zelig(out ~ time * treatment, data = dat)
summary(z)

ここでトリミングされた出力...

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)
(Intercept)        2.40264    0.71552   3.358  0.00081
time-11           -1.61292    1.08177  -1.491  0.13622
time-10           -1.03283    0.99850  -1.034  0.30116
time-9            -1.47934    1.02667  -1.441  0.14987
time-8            -0.35614    1.02667  -0.347  0.72874
time-7            -1.05803    1.04304  -1.014  0.31061
time-6            -2.25316    1.16178  -1.939  0.05269
.... 
treatment          1.28097    0.89440   1.432  0.15234
time-11:treatment  2.86965    1.30927   2.192  0.02859
time-10:treatment  1.69479    1.25788   1.347  0.17813
time-9:treatment   1.78684    1.27330   1.403  0.16078
time-8:treatment   0.82332    1.27330   0.647  0.51801
time-7:treatment   1.62808    1.28334   1.269  0.20482
time-6:treatment   2.64653    1.36895   1.933  0.05344
time-5:treatment   3.08572    1.36895   2.254  0.02437
....

時間ごとに効果をプロットできるように、各時間の最初の差 (治療 = 1、治療 = 0) を推定したいと思います。

何か案は？前もって感謝します

4

1 に答える 1

1

ここでは、ループを使用したソリューションです。

m <- zelig(outcome ~ time * treatment, model = "ls", data = dat)

output <- NULL

for (i in unique(dat$time)) {

t0 <-  setx(m, treatment = 0, time = i)
t1 <-  setx(m, treatment = 1, time = i)

ss <- sim(m, x = t0, x1 = t1, num = 10000)
fd <- unlist(ss$sim.out[["x1"]][["fd"]])

r <- data.table(time = i, mean = mean(fd), low = quantile(fd, .025), high = quantile(fd, 0.975))
output <- rbind(output, r)
}

output
    time     mean         low     high
 1:  -12 1.506365 -0.30605416 3.347631
 2:  -11 1.013915 -0.83479749 2.817791
 3:  -10 2.673004  0.72371241 4.645537
 4:   -9 1.291547 -0.62162353 3.183365
 5:   -8 2.985348  0.59834003 5.351312
 6:   -7 3.911258  1.95825840 5.878157
 7:   -6 4.222870  1.86773822 6.567400
 8:   -5 3.152967  0.81620039 5.483884
 9:   -4 3.893867  1.77629999 6.003647
10:   -3 2.319123  0.35445149 4.278032
11:   -2 1.942848  0.03771276 3.844245
12:   -1 3.879313  1.92915419 5.852765
13:    0 1.388601 -0.93881332 3.703387
14:    1 3.576107  1.54679622 5.567298
15:    2 2.413652  0.58863014 4.225094
16:    3 2.160988  0.03251586 4.266438
17:    4 2.203825  0.28985053 4.080361
18:    5 4.445642  2.40569051 6.510071
19:    6 1.504513 -0.27797349 3.251175
20:    7 2.542558  0.77794333 4.269277
21:    8 2.682681  0.93322199 4.449863
22:    9 4.271228  2.39189897 6.137469
23:   10 2.540004  0.66875643 4.454354
24:   11 3.454584  1.54938921 5.340096
25:   12 3.682521  1.85539403 5.501669
    time     mean         low     high
于 2016-04-11T21:29:24.630 に答える