Calculating the standard error of marginal effects using the delta method, in R

Here's a little piece of code I used to calculate marginal effects on probabilities in a non-linear (bivariate probit) model, using the Delta Method.

Warning: I am not a proper statistician. Use at your own risk!

# to create the model:

library(VGAM)
mbvp.1 <- vglm(cbind(decision_matrices.1-1, decision_matrices.6-1) ~ defect_other, binom2.rho, data=ds, trace=TRUE)

# needed to compute the gradient of my function numerically

library(numDeriv)
# needed for a multivariate normal distribution
library(mvtnorm)
# bv is a vector of betas
# this calculates the change in predicted probabilities
f <- function (bv) {
    rho <- rhobit(bv[3], inverse=T)
    p <- array(dim=c(2,2,2))
    for (c1 in 0:1) {
        for (c2 in 0:1) {
            for (d_o in 0:1) {
                sgn1 <- 2*c1 - 1
                sgn2 <- 2*c2 - 1

                ############################################

                # the probability of making choice (c1,c2) 

                # when the independent variable

                # defect_other=d_o:
                ############################################
                p[c1+1,c2+1,d_o+1] <- pmvnorm(upper=c(sgn1*(bv[1]+bv[4]*d_o),sgn2*(bv[2]+bv[5]*d_o)), corr=matrix(c(1,sgn1*sgn2*rho,sgn1*sgn2*rho,1), nrow=2))
            }
        }
    }
    # I use the global firstdec so that the function conforms to what grad() expects
    # the total probability of choosing 1 for either the first
    # or the second decision:
    if (firstdec) p[2,2,2]+p[2,1,2]-p[2,2,1]-p[2,1,1] else p[2,2,2]+p[1,2,2]-p[2,2,1]-p[1,2,1]
}


firstdec <- TRUE
pred[1] <- f(coef(mbvp.1))
gradient <- grad(f, coef(mbvp.1))
# ===== this is the delta method ====
predse[1] <- sqrt(gradient %*% vcov(mbvp.1) %*% gradient)

firstdec <- FALSE
pred[2] <- f(coef(mbvp.1))
gradient <- grad(f, coef(mbvp.1))
predse[2] <- sqrt(gradient %*% vcov(mbvp.1) %*% gradient)

margprobs <- cbind(pred, se=predse, Z=pred/predse, p=2*pnorm(-abs(pred/predse)))
print(margprobs, digits=3)