Marginal effects with logistic generalized additive model in R

Question

I am currently working with a logistic semi-parametric model in R using the mgcv package. The output from the model gives the standard log-odds coefficients; however, reviewers have requested marginal effects (like the ones in Stata using the margins command). I would like to do average marginal effects (though, marginal effects at the means--modes for categorical covariates--would at least be a start).

I was wondering if anyone has an implementation for this in R for models that use the gam() function. I have a rather large data set (1.2 million observations, a large number of discrete and continuous covariates, and fixed effects for 2,000 individuals). Are these estimates even tractable, given the non-parametric treatment of a couple of the continuous covariates? Any information would be helpful.

Here is what I am working with, though I get errors when attempting to do the bootstrapped SEs for the effects (comes from this helpful site probitlogit-marginal-effects-in-r-2/). I am not sure how this treats the non-parametrically estimated smooths (but maybe these don't matter, since it is using the "predict" function?):

mfxboot <- function(modform,dist,data,boot=1000,digits=3){ #dist is the distribution choice of logit or probit
      require(mgcv)

x <- gam(modform, family=binomial(link=dist),method="GCV.Cp",data)
      # get marginal effects

pdf <- ifelse(dist=="probit",               
 mean(dnorm(predict(x, type = "link")))            
 mean(dlogis(predict(x, type = "link")))
 marginal.effects <- pdf*coef(x)


bootvals <- matrix(rep(NA,boot*length(coef(x))), nrow=boot)  
set.seed(1111)  
for(i in 1:boot){    
samp1 <- data[sample(1:dim(data)[1],replace=T,dim(data)[1]),]    
x1 <- gam(modform, family=binomial(link=dist),method="GCV.Cp",samp1)    
pdf1 <- ifelse(dist=="probit",                   
mean(dnorm(predict(x1, type = "link"))),                   
mean(dlogis(predict(x1, type = "link"))))       
bootvals[i,] <- pdf1*coef(x1)     
}

res <- cbind(marginal.effects,apply(bootvals,2,sd),marginal.effects/apply(bootvals,2,sd))     
if(names(x$coefficients[1])=="(Intercept)"){        
res1 <- res[2:nrow(res),]    
res2 <- matrix(as.numeric(sprintf(paste("%.",paste(digits,"f",sep=""),sep=""),res1)),nrow=dim(res1)[1])         
rownames(res2) <- rownames(res1)        
} else {    
res2 <- matrix(as.numeric(sprintf(paste("%.",paste(digits,"f",sep=""),sep="")),nrow=dim(res)[1]))       
rownames(res2) <- rownames(res)    
}     
colnames(res2) <- c("marginal.effect","standard.error","z.ratio") 
      return(res2)
}