# 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, 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(dig=3, (newp-fitted(fit)[ii,])/hh) # Finite-difference approximation
round(dig=3, margeff(fit, subset=ii)["let",])
# Other examples
round(dig=3, margeff(fit))
round(dig=3, margeff(fit, subset=2)["let",])
round(dig=3, margeff(fit, subset=c(FALSE,TRUE))["let",,]) # recycling
round(dig=3, margeff(fit, subset=c(2,4,6,8))["let",,])
Run the code above in your browser using DataLab