marginalEffects <- function(model, ...){
    UseMethod("marginalEffects")
    }

marginalEffects.glm <- function(model, method=c("averaged", "atmeans"), ...){
    method <- match.arg(method)
    link <- model$family$link
    if (!(link %in% c("logit", "probit"))) stop ("must use logit or probit link")
    f <- if (link == "logit") dlogis else dnorm
    if (method == "averaged"){
        Xb <- predict(model, type="link")
        colMeans(outer(f(Xb), coef(model)))
        }
    else {
        b <- coef(model)
        xbar <- colMeans(model.matrix(model))
        xbarb <- sum(xbar * b)
        f(xbarb)*b
        }
    }

SLID <- read.table(        
"http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/SLID-women.txt", 
    header=TRUE)
    
SLID$region <- factor(SLID$region,       # re-order levels
    levels=c("Atlantic", "Quebec", "Ontario", "Prairies", "BC"))

mod.slid.logit <- glm(working ~ region + kids0004 + kids0509 + familyIncome 
    + education, family=binomial, data=SLID)

mod.slid.probit <- glm(working ~ region + kids0004 + kids0509 + familyIncome 
    + education, family=binomial(link=probit), data=SLID)
    
marginalEffects(mod.slid.logit)
marginalEffects(mod.slid.logit, method="atmeans")

marginalEffects(mod.slid.probit)
marginalEffects(mod.slid.probit, method="atmeans")