创建包含增量方法的数据框架

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)

用δ法考察了氮素对产量的边际效应。在平均氮值为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包来实现这一点，但是下一个版本的the marginaleffects package将包括一个非常容易做到这一点的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