我目前正在使用 mgcv 包在 R 中使用逻辑半参数模型。模型的输出给出了标准的对数优势系数；然而，审阅者要求边际效应（就像 Stata 中使用边距命令的那样）。我想做平均边际效应（尽管，平均边际效应——分类协变量的模式——至少是一个开始）。

我想知道是否有人在 R 中为使用 gam() 函数的模型实现了这一点。我有一个相当大的数据集（120 万个观测值，大量离散和连续协变量，以及 2,000 个人的固定效应）。考虑到几个连续协变量的非参数处理，这些估计是否容易处理？任何信息都有帮助。

这是我正在使用的内容，尽管在尝试为效果执行引导式 SE 时出现错误（来自这个有用的站点 probitlogit-marginal-effects-in-r-2/）。我不确定这如何处理非参数估计的平滑（但也许这些并不重要，因为它使用的是“预测”功能？）：

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)
}