R пробит регрессии предельные эффекты

Question

Я использую R, чтобы повторить исследование и получить в основном те же результаты, о которых сообщил автор.Однако в какой-то момент я рассчитываю предельные эффекты, которые кажутся нереально малыми.Я был бы очень признателен, если бы вы взглянули на мои рассуждения и приведенный ниже код и выяснили, ошибаюсь ли я в тот или иной момент.

В моей выборке содержится 24535 наблюдений, зависимая переменная "x028bin" является двоичнойпеременная принимает значения 0 и 1, и, кроме того, есть 10 объясняющих переменных.Девять из этих независимых переменных имеют числовые уровни, независимая переменная "f025grouped" является фактором, состоящим из различных религиозных конфессий.

Я хотел бы провести пробитную регрессию, включая манекены для религиозной конфессии, а затем вычислить предельные эффекты.Для этого я сначала устраняю пропущенные значения и использую кросс-таблицы между зависимыми и независимыми переменными, чтобы убедиться, что нет маленьких или нулевых ячеек.Затем я запускаю пробную модель, которая работает нормально, и я также получаю разумные результаты:

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)

Однако при расчете предельных эффектов со всеми переменными по их средним значениям из пробитных коэффициентов и масштабного коэффициента, предельные эффекты Iполучить слишком мало (например, 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

Я прошу прощения, что не могу предоставить вам рабочий пример, так как мой набор данных слишком большой.Любой комментарий будет принята с благодарностью.Большое спасибо.

Лучший,

Тобиас

Answer 1

Это сделает трюк для probit или 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-в-г /

Answer 2

@ Гэвин прав, лучше спросить на родственном сайте.

В любом случае, вот мой трюк для интерпретации пробит-коэффициентов.

Пробит-коэффициенты регрессии такие же, как логит-коэффициенты, вплоть до шкалы (1.6).Итак, если подгонка пробитной модели Pr(y=1) = fi(.5 - .3*x), это эквивалентно логистической модели Pr(y=1) = invlogit(1.6(.5 - .3*x)).

. И я использую это для создания графики, используя функцию invlogit пакета arm.Другая возможность - просто умножить все коэффициенты (включая перехват) на 1,6, а затем применить «правило деления на 4» (см. Книгу Гельмана и Хилла), то есть разделить новые коэффициенты на 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