mfxboot函数对概率回归的边际影响? [英] mfxboot function for marginal effects for probit regressions?

问题描述

数据:数据

代码:

#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函数时，出现以下错误消息:

When I run the mfxboot function, I get the following error message:

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

Error in bootvals[i, ] <- pdf1 * coef(x1) : number of items to replace is not a multiple of replacement length

关于为什么会发生的任何想法?还有关于如何解决此问题的建议吗?

Any idea as to why that happens? And Any suggestion of how to go around this issue?

谢谢.

推荐答案

我无法重现您的错误.也许您应该在问题中添加sessionInfo()输出.尽管如此，下面还是建议对mfxboot函数进行增强.

I am not able to reproduce your error. Perhaps you should add sessionInfo() output to your question. A suggested enhancement to the mfxboot function follows below nonetheless.

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

My suggestion would be to refactor the mfxboot function into two functions -- one that returns the marginal effects given a glm object, and the second which bootstraps it.

您可以使用 Boot 轻松完成此操作在car包中起作用，因为 这是引导glm对象的不错的前端.

You can do this easily using the Boot function in the car package since that is a nice front-end for bootstrapping glm objects.

这里有一些代码演示了此过程，阅读起来更加简洁:

Here is some code that demonstrates this process, which is much cleaner to read:

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步:编写一个返回边际效应的函数

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

Step 2: Write a function that returns the marginal effects

The following function takes as input a glm object of the binomial family and computes appropriate marginal effects for logit and probit links.

#================================================
# 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步:引导边际效应

以下代码使用car包中的Boot函数来引导边缘效应.请注意Boot的接口如何针对从lm和glm对象得出的统计信息进行优化.

Step 3: Bootstrap the marginal effects

The following code uses the Boot function from the car package to bootstrap the marginal effects. Note how the interface of Boot is optimized for statistics derived from lm and glm objects.

#================================================
# 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

这篇关于mfxboot函数对概率回归的边际影响?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！