Creating data frame including delta method

user_17619490 0
Viewed 70

I have estimated a regression where yields are a function of nitrogen.

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)

I have used the delta method to see the marginal effects of nitrogen on yield. I can see the marginal effect of nitrogen at the average value of nitrogen =110

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

I want to make a data frame to easily see the changes in the marginal effect of nitrogen on yield at different nitrogen rates ranging from 0-110 relative to the average.

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)

I am getting an error code when I try to run B<-deltaMethod(reg, "(nitrogen*Nitrogen_rate)+(nitrogen_square*Nitrogen_rate^2)", vcov(reg)). The error says 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.

I am new to R and looking for guidance on how to get this formula in a data frame.

dataframe r
1 Answers

I do not know how to use the car package to achieve this, but the next version of the marginaleffects package will include a deltamethod() which can do this very easily. (Disclaimer: I am the author.) You can use the new function now by installing it from Github:

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

Restart R completely, then:

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

Use the tail() function to print a part of the results:

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

If you add a few parentheses and combine the two hypotheses in a one-liner, you can even estimate standard errors and confidence intervals around the difference of interest itself:

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
user_3423310
Related