这是用mgcv包配合写的，你看下能不能用
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)
}
扫码加好友，拉您进群
全部版块
我的主页
收藏
