I've simulated an example in which a treatment is randomly allocated to a group of people and then some continuous measure is observed.
The treatment group sees a spike in their response and then a sharp decline. A plot of the conditional mean is show below. The control group has a flat response.
The sum of difference between treatment and control over days, plus its standard error, is of interest. I was hoping I could do this with marginaleffects
From the functional form, sum of differences should be around 3.2. However, when I use avg_comparisons
I get something closer to 0.
Here is a minimal working example
library(nimble)
#> nimble version 1.0.1 is loaded.
#> For more information on NIMBLE and a User Manual,
#> please visit https://R-nimble.org.
#>
#> Note for advanced users who have written their own MCMC samplers:
#> As of version 0.13.0, NIMBLE's protocol for handling posterior
#> predictive nodes has changed in a way that could affect user-defined
#> samplers in some situations. Please see Section 15.5.1 of the User Manual.
#>
#> Attaching package: 'nimble'
#> The following object is masked from 'package:stats':
#>
#> simulate
library(tidyverse)
library(marginaleffects)
set.seed(0)
mu <- \(x) 10*ddexp(x, location=7, scale=1) - 7*ddexp(x, location=8, scale=1)
d <- crossing(
day = 1:14,
trt = 0:1,
ix = 1:100
) %>%
mutate(
y = 10 + mu(day)*trt + rnorm(n(), 0, 1)
)
fit <- lm(y ~ factor(day)*trt, data=d)
# This gives me something close to what I want
avg_comparisons(fit,
newdata = datagrid(by = c('trt','day')),
variables = 'trt',
by='day')$estimate %>%
sum
#> [1] 3.469189
# This does not
avg_comparisons(fit,
newdata = datagrid(by = c('trt','day')),
variables = 'trt')
#>
#> Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
#> trt 1 - 0 0.248 0.0379 6.53 <0.001 33.9 0.173 0.322
#>
#> Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
#> Type: response
Created on 2023-10-27 with reprex v2.0.2
How can I estimate the effect I'm interested in using marginaleffects
?
I’m not sure exactly what you mean by
The sum of difference between treatment and control over days
But there’s a good chance you can get what you want by specifying a custom function in the comparison
argument of the comparisons()
function. Here’s an example where I compute the sum of differences between treatment and control for each day:
library(nimble)
library(tidyverse)
library(marginaleffects)
set.seed(0)
mu <- \(x) 10 * ddexp(x, location = 7, scale = 1) - 7 * ddexp(x, location = 8, scale = 1)
d <- crossing(
day = 1:14,
trt = 0:1,
ix = 1:100) %>%
mutate(y = 10 + mu(day) * trt + rnorm(n(), 0, 1))
fit <- lm(y ~ factor(day) * trt, data = d)
comparisons(fit,
variables = "trt",
by = "day",
comparison = \(hi, lo) sum(hi - lo))
#
# Term Contrast day Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
# trt 1, 0 1 -0.1180 0.284 -0.416 0.67747 0.6 -0.6743 0.438
# trt 1, 0 2 -0.0780 0.284 -0.275 0.78353 0.4 -0.6342 0.478
# trt 1, 0 3 0.0609 0.284 0.215 0.83010 0.3 -0.4953 0.617
# trt 1, 0 4 0.5200 0.284 1.832 0.06688 3.9 -0.0362 1.076
# trt 1, 0 5 1.0589 0.284 3.731 < 0.001 12.4 0.5027 1.615
# trt 1, 0 6 3.1287 0.284 11.025 < 0.001 91.5 2.5725 3.685
# trt 1, 0 7 7.1686 0.284 25.260 < 0.001 465.2 6.6123 7.725
# trt 1, 0 8 -3.2755 0.284 -11.542 < 0.001 100.0 -3.8317 -2.719
# trt 1, 0 9 -1.0755 0.284 -3.790 < 0.001 12.7 -1.6318 -0.519
# trt 1, 0 10 -0.8138 0.284 -2.867 0.00414 7.9 -1.3700 -0.258
# trt 1, 0 11 -0.1793 0.284 -0.632 0.52742 0.9 -0.7356 0.377
# trt 1, 0 12 0.1109 0.284 0.391 0.69607 0.5 -0.4454 0.667
# trt 1, 0 13 0.1139 0.284 0.401 0.68809 0.5 -0.4423 0.670
# trt 1, 0 14 0.3166 0.284 1.115 0.26466 1.9 -0.2397 0.873
#
# Columns: term, contrast, day, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
# Type: response