提问人：Hack-R 提问时间：9/20/2014 最后编辑：Hack-R 更新时间：10/22/2014 访问量：1106
如何在 R 中对 Bootstrap 概率模型进行编码选择
How to Code Selection for Bootstrap Probit Models in R
问：
这个问题涉及如何在概率模型中对变量选择进行编码，并具有边际效应（直接或通过调用一些预先存在的包）。
作为一篇与 TLAPD 相关的博客文章，我正在对电影的免费和商业可用性对这些电影盗版水平的影响进行一些概率回归。
在 R 中运行概率的简单方法通常是通过 ，即：glm
probit <- glm(y ~ x1 + x2, data=data, family =binomial(link = "probit"))
但这对解释来说是有问题的，因为它不提供边际效应。
通常，如果我想从概率回归中获得边际效应，我会定义这个函数（我不记得原始来源，但它是一个流行的函数，经常被转发）：
mfxboot <- function(modform,dist,data,boot=500,digits=3){
x <- glm(modform, family=binomial(link=dist),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)
# start bootstrap
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 <- glm(modform, family=binomial(link=dist),samp1)
pdf1 <- ifelse(dist=="probit",
mean(dnorm(predict(x, type = "link"))),
mean(dlogis(predict(x, 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)
}
然后像这样运行回归：
mfxboot(modform = "y ~ x1 + x2",
dist = "probit",
data = piracy)
但是使用这种方法，我不知道我可以运行任何变量选择算法，如前进、后退、逐步等。
解决这个问题的最佳方法是什么？有没有更好的方法在 R 中运行概率，既能报告边际效应，又能自动选择模型？还是我应该专注于使用和执行该函数的变量选择？mfxboot
谢谢！
答：
1赞
tchakravarty
10/22/2014
#1
目前尚不清楚为什么会出现问题。给定模型的模型（变量）选择和边际效应的计算是连续的，没有理由尝试将两者结合起来。
以下是在模型（变量）选择后计算边际效应及其自举标准效应的方法：
使用您喜欢的模型选择过程（包括下面讨论的引导模型选择技术，不要与用于计算最终模型边际效应标准误差的引导程序混淆）执行变量选择。
以下是本问题中提供的数据集的示例。另请注意，这绝不是对使用逐步变量选择方法的认可。
#================================================ # read in data, and perform variable selection for # a probit model #================================================ dfE = read.csv("ENAE_Probit.csv") formE = emploi ~ genre + filiere + satisfaction + competence + anglais glmE = glm(formula = formE, family = binomial(link = "probit"), data = dfE) # perform model (variable) selection glmStepE = step(object = glmE)
- 现在，将所选模型传递给计算边际效应的函数。
#================================================ # function: compute marginal effects for logit and probit models # NOTE: this assumes that an intercept has been included by default #================================================ fnMargEffBin = function(objBinGLM) { stopifnot(objBinGLM$family$family == "binomial") vMargEff = switch(objBinGLM$family$link, probit = colMeans(outer(dnorm(predict(objBinGLM, type = "link")), coef(objBinGLM))[, -1]), logit = colMeans(outer(dlogis(predict(objBinGLM, type = "link")), coef(objBinGLM))[, -1]) ) return(vMargEff) } # test the function fnMargEffBin(glmStepE)
输出如下：
> fnMargEffBin(glmStepE) genre filiere 0.06951617 0.04571239
- 为了获得边际效应的标准误差，您可以引导边际效应，例如，使用函数中的函数，因为它提供了这样一个简洁的接口来引导从估计值中得出的统计数据。
Bootcarglm
#================================================ # compute bootstrap std. err. for the marginal effects #================================================ margEffBootE = Boot(object = glmStepE, f = fnMargEffBin, labels = names(coef(glmE))[-1], R = 100) summary(margEffBootE)
输出如下：
> summary(margEffBootE) R original bootBias bootSE bootMed genre 100 0.069516 0.0049706 0.045032 0.065125 filiere 100 0.045712 0.0013197 0.011714 0.048900
附录：
从理论上讲，有两种方法可以解释自举变量选择问题。
您可以通过使用引导模型拟合标准作为拟合度量来执行模型选择（变量选择）。Shao（1996）概述了这方面的理论，并且需要一种子抽样方法。
然后，根据上面选择的最佳模型计算边际效应及其引导标准误差。您可以对多个引导样本执行变量选择，并通过查看多个引导模型选择中保留的变量来得出一个最佳模型，或者使用模型平均估计器。Austin和Tu（2004）对此进行了讨论。
然后，根据上面选择的最佳模型计算边际效应及其引导标准误差。
评论
0赞
Hack-R
10/22/2014
非常感谢您的详细回答。我会+1，如果没有更好的答案，我会在一段时间后将其标记为解决方案。但是，出于几个原因，我实际上确实想一步到位地完成选择 + 引导，我应该更好地澄清这一点。您在最后提到的原因是其中很大一部分，另一个原因是我通常处理大数据，运行模型选择需要很长时间，引导需要更长的时间，因为我认为有某种方法可以做到这一点，我认为在一个函数/步骤中估计所有内容可能会更快。
评论