r语言 - 创建包含delta方法的数据帧

Question

我估计了一个回归，其中产量是氮的函数。

yields<-rnorm(50, mean=2000, sd=10)
nitrogen<-rnorm(50, 110, sd=5)
nitrogen_square <- nitrogen^2
reg<-lm(yields~nitrogen+nitrogen_square)
summary(reg)

我用delta法来观察氮对产量的边际效应。我可以看到氮的边际效应在氮的平均值=110

deltaMethod(reg, "(nitrogen*110)+(nitrogen_square*110^2)", vcov(reg))

我想做一个数据框架，以方便地看到不同施氮率下氮对产量的边际效应的变化，相对于平均值，范围从0-110。

Nitrogen_avg<-110
Nitrogen_rate<-0:110
A<-deltaMethod(reg, "(nitrogen*Nitrogen_avg)+(nitrogen_square*Nitrogen_avg^2)", vcov(reg))
B<-deltaMethod(reg, "(nitrogen*Nitrogen_rate)+(nitrogen_square*Nitrogen_rate^2)", vcov(reg))
Marg_yield_diff<- A-B
df<- data.frame(Nitrogen_avg, Nitrogen_rate, A, B, Marg_yield_diff)

当我尝试运行B<-deltaMethod(reg, "(nitrogen*Nitrogen_rate)+(nitrogen_square*Nitrogen_rate^2)", vcov(reg))时，我得到一个错误代码。错误显示为Error in row.names < 1L : comparison is not allowed for expressions In addition: Warning messages: 1: In gd[i] <- eval(reg(g., para.names[i]), envir) : number of items to replace is not a multiple of replacement length 2: In gd[i] <- eval(reg(g., para.names[i]), envir) : number of items to replace is not a multiple of replacement length。

我是R的新手，正在寻找如何在数据框架中获得此公式的指导。

Answer 1

我不知道如何使用car包来实现这一点，但是marginaleffects包的下一个版本将包括一个deltamethod()，它可以很容易地做到这一点。(免责声明:我是作者。)你现在可以通过从Github安装它来使用这个新功能:

library(remotes)
install_github("vincentarelbundock/marginaleffects")

完全重启R，然后:

library(marginaleffects)
yields <- rnorm(50, mean = 2000, sd = 10)
nitrogen <- rnorm(50, 110, sd = 5)
nitrogen_square <- nitrogen^2
Nitrogen_avg <- 110
Nitrogen_rate <- 0:110
reg <- lm(yields ~ nitrogen + nitrogen_square)
A <- deltamethod(
model = reg,
hypothesis = "(nitrogen*Nitrogen_avg)+(nitrogen_square*Nitrogen_avg^2) = 0")
B <- deltamethod(
model = reg,
hypothesis = "(nitrogen*Nitrogen_rate)+(nitrogen_square*Nitrogen_rate^2) = 0")
B$Marg_yield_diff <- A$estimate - B$estimate

使用tail()函数打印部分结果:

tail(B)
#>                                                               term
#> 106 (nitrogen*Nitrogen_rate)+(nitrogen_square*Nitrogen_rate^2) = 0
#> 107 (nitrogen*Nitrogen_rate)+(nitrogen_square*Nitrogen_rate^2) = 0
#> 108 (nitrogen*Nitrogen_rate)+(nitrogen_square*Nitrogen_rate^2) = 0
#> 109 (nitrogen*Nitrogen_rate)+(nitrogen_square*Nitrogen_rate^2) = 0
#> 110 (nitrogen*Nitrogen_rate)+(nitrogen_square*Nitrogen_rate^2) = 0
#> 111 (nitrogen*Nitrogen_rate)+(nitrogen_square*Nitrogen_rate^2) = 0
#>      estimate std.error  statistic   p.value  conf.low conf.high
#> 106 -156.6524  574.6930 -0.2725846 0.7851726 -1283.030  969.7252
#> 107 -156.7363  575.1396 -0.2725187 0.7852232 -1283.989  970.5166
#> 108 -156.7935  575.4914 -0.2724516 0.7852748 -1284.736  971.1488
#> 109 -156.8242  575.7485 -0.2723832 0.7853274 -1285.270  971.6220
#> 110 -156.8283  575.9109 -0.2723136 0.7853809 -1285.593  971.9362
#> 111 -156.8059  575.9785 -0.2722426 0.7854355 -1285.703  972.0913
#>     Marg_yield_diff
#> 106     -0.15343896
#> 107     -0.06961516
#> 108     -0.01235936
#> 109      0.01832843
#> 110      0.02244822
#> 111      0.00000000

如果你加上一些括号，并将两个假设合并成一行，你甚至可以估计兴趣差本身的标准误差和置信区间:

hyp <- "((nitrogen*Nitrogen_avg)+(nitrogen_square*Nitrogen_avg^2)) - ((nitrogen*Nitrogen_rate)+(nitrogen_square*Nitrogen_rate^2)) = 0"
deltamethod(
model = reg,
hypothesis = hyp) |>
head() |>
subset(select = 2:7)
#>    estimate std.error  statistic   p.value  conf.low conf.high
#> 1 -156.8059  575.9785 -0.2722426 0.7854355 -1285.703  972.0913
#> 2 -153.9324  565.5740 -0.2721702 0.7854911 -1262.437  954.5723
#> 3 -151.0855  555.2645 -0.2720965 0.7855478 -1239.384  937.2129
#> 4 -148.2652  545.0500 -0.2720213 0.7856056 -1216.543  920.0131
#> 5 -145.4714  534.9304 -0.2719446 0.7856646 -1193.916  902.9728
#> 6 -142.7042  524.9058 -0.2718664 0.7857248 -1171.501  886.0922