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)