3

我正在使用 R 来复制一项研究，并获得与作者报告的几乎相同的结果。然而，在某一时刻，我计算的边际效应似乎小得不切实际。如果您能看看我的推理和下面的代码，看看我是否在某个地方或另一个地方弄错了，我将不胜感激。

我的样本包含 24535 个观测值，因变量“x028bin”是取值 0 和 1 的二元变量，此外还有 10 个解释变量。其中九个自变量具有数字级别，自变量“f025grouped”是由不同宗教教派组成的因素。

我想运行一个概率回归，包括宗教教派的假人，然后计算边际效应。为此，我首先消除缺失值并使用因变量和自变量之间的交叉表来验证没有小单元格或 0 单元格。然后我运行运行良好的概率模型，我也获得了合理的结果：

probit4AKIE <- glm(x028bin ~ x003 + x003squ + x025secv2 + x025terv2 + x007bin + x04chief + x011rec + a009bin + x045mod + c001bin + f025grouped, family=binomial(link="probit"), data=wvshm5red2delna, na.action=na.pass)

summary(probit4AKIE)

但是，当从概率系数和比例因子计算所有变量的均值时，我获得的边际效应太小（例如 2.6042e-78）。代码如下所示：

ttt <- cbind(wvshm5red2delna$x003,
wvshm5red2delna$x003squ,
wvshm5red2delna$x025secv2,
wvshm5red2delna$x025terv2,
wvshm5red2delna$x007bin,
wvshm5red2delna$x04chief,
wvshm5red2delna$x011rec,
wvshm5red2delna$a009bin,
wvshm5red2delna$x045mod,
wvshm5red2delna$c001bin,
wvshm5red2delna$f025grouped,
wvshm5red2delna$f025grouped,
wvshm5red2delna$f025grouped,
wvshm5red2delna$f025grouped,
wvshm5red2delna$f025grouped,
wvshm5red2delna$f025grouped,
wvshm5red2delna$f025grouped,
wvshm5red2delna$f025grouped,
wvshm5red2delna$f025grouped) #I put variable "f025grouped" 9 times because this variable consists of 9 levels

ttt <- as.data.frame(ttt)

xbar <- as.matrix(mean(cbind(1,ttt[1:19]))) #1:19 position of variables in dataframe ttt

betaprobit4AKIE <- probit4AKIE$coefficients

zxbar <- t(xbar) %*% betaprobit4AKIE

scalefactor <- dnorm(zxbar)

marginprobit4AKIE <- scalefactor * betaprobit4AKIE[2:20] #2:20 are the positions of variables in the output of the probit model 'probit4AKIE' (variables need to be in the same ordering as in data.frame ttt), the constant in the model occupies the first position

marginprobit4AKIE #in this step I obtain values that are much too small

很抱歉，由于我的数据集太大，我无法为您提供工作示例。任何评论将不胜感激。非常感谢。

最好的，

托拜厄斯

4

2 回答 2

1

这将为probit or logit解决问题：

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

资料来源：http ://www.r-bloggers.com/probitlogit-marginal-effects-in-r/

于 2014-09-10T19:31:13.957 回答
1

@Gavin 是对的，最好在姐妹网站上询问。

无论如何，这是我解释概率系数的技巧。

概率回归系数与 logit 系数相同，最高为 (1.6)。因此，如果概率模型的拟合为Pr(y=1) = fi(.5 - .3*x)，则相当于逻辑模型Pr(y=1) = invlogit(1.6(.5 - .3*x))。

我用它来制作图形，使用 package 的函数 invlogit arm。另一种可能性是将所有系数（包括截距）乘以 1.6，然后应用“除以 4 规则”（参见 Gelman 和 Hill 的书），即将新系数除以 4，您会发现对应于 x 的单位差的预测差的上限。

这是一个例子。

x1 = rbinom(100,1,.5)
x2 = rbinom(100,1,.3)
x3 = rbinom(100,1,.9)
ystar = -.5  + x1 + x2 - x3 + rnorm(100)
y = ifelse(ystar>0,1,0)
probit = glm(y~x1 + x2 + x3, family=binomial(link='probit'))
xbar <- as.matrix(mean(cbind(1,ttt[1:3])))

# now the graphic, i.e., the marginal effect of x1, x2 and x3
library(arm)
curve(invlogit(1.6*(probit$coef[1] + probit$coef[2]*x + probit$coef[3]*xbar[3] + probit$coef[4]*xbar[4]))) #x1
curve(invlogit(1.6*(probit$coef[1] + probit$coef[2]*xbar[2] + probit$coef[3]*x + probit$coef[4]*xbar[4]))) #x2
curve(invlogit(1.6*(probit$coef[1] + probit$coef[2]*xbar[2] + probit$coef[3]*xbar[3] + probit$coef[4]*x))) #x3
于 2011-04-27T15:29:49.003 回答