# Marginal effects with logistic generalized additive model in R [closed]

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),replace=T,dim(data)),]
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)=="(Intercept)"){
res1 <- res[2:nrow(res),]
res2 <- matrix(as.numeric(sprintf(paste("%.",paste(digits,"f",sep=""),sep=""),res1)),nrow=dim(res1))
rownames(res2) <- rownames(res1)
} else {
res2 <- matrix(as.numeric(sprintf(paste("%.",paste(digits,"f",sep=""),sep="")),nrow=dim(res)))
rownames(res2) <- rownames(res)
}
colnames(res2) <- c("marginal.effect","standard.error","z.ratio")
return(res2)
}


## closed as off-topic by mkt, Peter Flom♦Jul 17 at 11:57

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – mkt, Peter Flom
If this question can be reworded to fit the rules in the help center, please edit the question.