Marginal effects calculation in R: logit model 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))) 

 A: 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).
