margeff: Marginal effects for the multinomial logit and cumulative models

Description

Marginal effects for the multinomial logit model and cumulative logit/probit/... models: the derivative of the fitted probabilities with respect to each explanatory variable.

Usage

margeff(object, subset = NULL)

Arguments

object

A vglm multinomial or cumulative object.

subset

Numerical or logical vector, denoting the required observation(s). Recycling is used if possible. The default means all observations.

Value

A $p$ by $M+1$ by $n$ array, where $p$ is the number of explanatory variables and the (hopefully) nominal response has $M+1$ levels, and there are $n$ observations.
If is.numeric(subset) and length(subset) == 1 then a $p$ by $M+1$ matrix is returned.

Warning

Care is needed in interpretation, e.g., the change is not universally accurate for a unit change in each explanatory variable because eventually the `new' probabilities may become negative or greater than unity. Also, the `new' probabilities will not sum to one.

This function is not applicable for models with data-dependent terms such as bs and poly. Also the function should not be applied to models with any terms that have generated more than one column of the LM model matrix, such as bs and poly. For such try using numerical methods such as finite-differences. The formula in object should comprise of simple terms of the form ~ x2 + x3 + x4, etc.

Details

Computes the derivative of the fitted probabilities of a multinomial logit model or cumulative logit/probit/... model with respect to each explanatory variable.

Examples

Run this code

# Not a good example for multinomial() because the response is ordinal!!
ii <- 3; hh <- 1/100
pneumo <- transform(pneumo, let = log(exposure.time))
fit <- vglm(cbind(normal, mild, severe) ~ let, multinomial, data = pneumo)
fit <- vglm(cbind(normal, mild, severe) ~ let,
            cumulative(reverse = TRUE,  parallel = TRUE),
            data = pneumo)
fitted(fit)[ii, ]

mynewdata <- with(pneumo, data.frame(let = let[ii] + hh))
(newp <- predict(fit, newdata = mynewdata, type = "response"))

# Compare the difference. Should be the same as hh --> 0.
round(digits = 3, (newp-fitted(fit)[ii, ])/hh)  # Finite-difference approxn
round(digits = 3, margeff(fit, subset = ii)["let",])

# Other examples
round(digits = 3, margeff(fit))
round(digits = 3, margeff(fit, subset = 2)["let",])
round(digits = 3, margeff(fit, subset = c(FALSE, TRUE))["let",,])  # recycling
round(digits = 3, margeff(fit, subset = c(2, 4, 6, 8))["let",,])

Run the code above in your browser using DataLab

Get 50% off unlimited learning