如何计算日级边际效应的总和？-解网

问：

我模拟了一个例子，在这个例子中，将治疗随机分配给一组人，然后观察到一些连续的测量。

治疗组的反应激增，然后急剧下降。条件均值图如下所示。对照组的反应是平坦的。

治疗和对照组在天数内的差异之和，加上其标准误差，是有意义的。我希望我能用marginaleffects

从函数形式来看，差异之和应在 3.2 左右。但是，当我使用时，我得到的东西接近 0。avg_comparisons

这是一个最小的工作示例

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

^{创建于 2023-10-27 使用 reprex v2.0.2}

我怎样才能估计我有兴趣使用的效果？marginaleffects

r 边际效应

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

如何计算日级边际效应的总和？

How do I compute the sum of day level marginal effects?

评论

评论