3

数据：数据

代码：

#function that calculates ‘the average of the sample marginal effects’.
mfxboot <- function(modform,dist,data,boot=1000,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)
}

## Regression
probit_enae = glm(emploi ~ genre + filiere + satisfaction + competence + anglais, family=binomial(link="probit"),
              data=ENAE_Probit.df)
summary(probit_enae) #Summary output of the regression
confint(probit_enae) #Gives the 95% confidence interval for the estimated coefficients

## Running the mfxboot for Marginal effects
mfx_enae = mfxboot(emploi ~ genre + filiere + satisfaction + competence + anglais,"probit",ENAE_Probit.df)

问题：

当我运行 mfxboot 函数时，我收到以下错误消息：

bootvals[i, ] <- pdf1 * coef(x1) 中的错误：要替换的项目数不是替换长度的倍数

知道为什么会这样吗？以及如何解决这个问题的任何建议？

谢谢。

4

1 回答 1

6

我无法重现您的错误。也许您应该将sessionInfo()输出添加到您的问题中。尽管如此，建议的功能增强mfxboot如下。

我的建议是将mfxboot函数重构为两个函数——一个返回给定glm对象的边际效应，第二个引导它。

您可以使用包Boot中的函数轻松完成此操作，因为这是一个很好的引导对象car的前端。glm

这是一些演示此过程的代码，阅读起来更清晰：

第 1 步：估计概率模型

library(car)

#================================================
# read in data, and estimate 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)

第 2 步：编写一个返回边际效应的函数

以下函数将族glm对象作为输入，binomial并计算 logit 和 probit 链接的适当边际效应。

#================================================
# 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(glmE)

第 3 步：引导边际效应

以下代码使用包中的Boot函数car来引导边际效应。请注意 的接口如何针对从和对象Boot派生的统计信息进行优化。lmglm

#================================================
# compute bootstrap std. err. for the marginal effects
#================================================
margEffBootE = Boot(object = glmE, f = fnMargEffBin, 
     labels = names(coef(glmE))[-1], R = 100)
summary(margEffBootE)

这是输出：

> summary(margEffBootE)
               R  original    bootBias   bootSE   bootMed
genre        100  0.070733  0.00654606 0.042162  0.074563
filiere      100  0.043173  0.00060356 0.014064  0.043486
satisfaction 100  0.050773 -0.00110501 0.037737  0.048310
competence   100 -0.020144  0.00407027 0.034194 -0.014987
anglais      100 -0.018906 -0.00170887 0.033522 -0.019164
于 2014-10-21T20:48:51.367 回答