R régression probit effets marginaux

Question

J'utilise R pour reproduire une étude et obtenir essentiellement les mêmes résultats les auteur a rapporté. À un moment donné, toutefois, je calcule les effets marginaux qui semblent être déraisonnablement faible. Je serais très reconnaissant si vous pouviez jeter un oeil à mon raisonnement et le code ci-dessous et voir si je me trompe à un moment ou un autre.

Mon échantillon contient 24535 observations, la variable dépendante "x028bin" est un variable binaire prenant les valeurs 0 et 1, et il y a en outre 10 variables explicatives. Neuf de ces variables indépendantes ont des niveaux numériques, la variable indépendante « f025grouped » est un facteur composé de différentes confessions religieuses.

Je voudrais lancer une régression probit, y compris pour les nuls confession religieuse et à calculer les effets marginaux. Pour ce faire, j'élimine d'abord les valeurs manquantes et utiliser les tableaux croisés entre les variables dépendantes et indépendantes pour vérifier qu'il n'y a pas de petites ou 0 cellules. Ensuite, je lance le modèle probit qui fonctionne très bien et j'obtenir également des résultats raisonnables:

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)

Cependant, lors du calcul des effets marginaux avec toutes les variables à leurs moyens à partir des coefficients de probit et un facteur d'échelle, les effets marginaux sont beaucoup plus que j'obtiens trop petit (par exemple 2.6042e-78). Le code ressemble à ceci:

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

Je suis désolé que je ne peux pas vous donner un exemple de travail que mon jeu de données est beaucoup trop grand. Tout commentaire serait grandement apprécié. Merci beaucoup.

Best,

Tobias

Answer 1

@Gavin est juste et il est préférable de demander à l'emplacement de soeur.

Dans tous les cas, voici mon truc pour interpréter les coefficients probit.

Les coefficients de régression probit sont les mêmes que les coefficients logit, jusqu'à une échelle graduée (1.6). Donc, si l'ajustement d'un modèle probit est Pr(y=1) = fi(.5 - .3*x), cela équivaut à la Pr(y=1) = invlogit(1.6(.5 - .3*x)) modèle logistique.

Et je l'utiliser pour faire un graphique, en utilisant la fonction invlogit du paquet arm. Une autre possibilité est juste de multiplier les coefficients (y compris l'interception) de 1,6, puis en appliquant la « division par 4 règle » (voir le livre de Gelman et Hill), à savoir, diviser les nouveaux coefficients de 4, et vous découvrirez une limite supérieure de la différence de prédiction correspondant à une différence d'unité en x.

Voici un exemple.

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

Answer 2

fera l'affaire pour probit ou 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)
}

Source: http://www.r-bloggers.com / / probitlogit marginaux effets-en-r