I have a dataset with this structure:
# 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
I run a LS model interacting time and treatment:
z <- zls$new()
z$zelig(out ~ time * treatment, data = dat)
summary(z)
Here a trimmed output...
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
....
I would like to estimate the first differences (treatment = 1, treatment = 0) for each time so that I can plot the effects by time.
Any ideas? Thank you in advance
Here a solution using a loop.
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