1
$\begingroup$

How do you calculate marginal effects of parameters of logit model in R uging package {glm}?

Are following codes correct?

#### preparation ####
# dependent variable
yseed <- rnorm(100)
y <- ifelse(yseed > 0, 1, 0)

# independent variables
x1 <- rnorm(100, mean=100, sd=20)
x2 <- rnorm(100, mean=50, sd=20)
X <- cbind(1, x1, x2) 
Xmean <- apply(X, 2, mean)

#### analysis ####
# logit model
res <- glm(y ~ x1 + x2, family=binomial(link="logit"))
summary(res)

# Marginal effects (ME) calculation
LAMBDA <- function(x) { 1 / (1 + exp(-x))} # cdf of standard logistic distribution

# ME of (intercept)
ME_1  <- coef[1] * LAMBDA(Xmean %*% coef(res)) * (1 - LAMBDA(Xmean %*% coef(res)))
# ME of x1   
ME_x1 <- coef[2] * LAMBDA(Xmean %*% coef(res)) * (1 - LAMBDA(Xmean %*% coef(res))) 
# ME of x2
ME_x2 <- coef[3] * LAMBDA(Xmean %*% coef(res)) * (1 - LAMBDA(Xmean %*% coef(res))) 
$\endgroup$
1
  • 2
    $\begingroup$ This question appears to be off-topic because it is abouthow to use R without a reproducible example. $\endgroup$ Jan 9 '15 at 13:09
2
$\begingroup$

From a quick search, I think you must first define a function that calculates the average of the sample marginal effects for the logit model. You can have a look at this: http://www.r-bloggers.com/probitlogit-marginal-effects-in-r/

mfx <- function(x,sims=1000){
set.seed(1984)
pdf <- ifelse(as.character(x$call)[3]=="binomial(link = \"probit\")",
    mean(dnorm(predict(x, type = "link"))),
    mean(dlogis(predict(x, type = "link"))))
    pdfsd <- ifelse(as.character(x$call)[3]=="binomial(link = \"probit\")",
sd(dnorm(predict(x, type = "link"))),
sd(dlogis(predict(x, type = "link"))))
marginal.effects <- pdf*coef(x)
sim <- matrix(rep(NA,sims*length(coef(x))), nrow=sims)
for(i in 1:length(coef(x))){
  sim[,i] <- rnorm(sims,coef(x)[i],diag(vcov(x)^0.5)[i])
}
pdfsim <- rnorm(sims,pdf,pdfsd)
sim.se <- pdfsim*sim
res <- cbind(marginal.effects,sd(sim.se))
colnames(res)[2] <- "standard.error"
ifelse(names(x$coefficients[1])=="(Intercept)",
return(res[2:nrow(res),]),return(res))
}

you can decide to change the number of simulations from which std errors are calculated. Once you define your glm model (glm.model <- glm(y ~ x1+x2+x3,data=data.frame, family = binomial(link = "logit")), you can calculate the marginal effects using mfx(glm.model).

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.