# Not a good example for multinomial() since 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((newp-fitted(fit)[ii, ]) / hh, 3) # Finite-diff approxn
round(margeff(fit, subset = ii)["let",], 3)
# Other examples
round(margeff(fit), 3)
round(margeff(fit, subset = 2)["let",], 3)
round(margeff(fit, subset = c(FALSE, TRUE))["let",,], 3) # Recycling
round(margeff(fit, subset = c(2, 4, 6, 8))["let",,], 3)
# Example 3; margeffs at a new value
mynewdata2a <- data.frame(let = 2) # New value
mynewdata2b <- data.frame(let = 2 + hh) # For finite-diff approxn
(neweta2 <- predict(fit, newdata = mynewdata2a))
fit@x[1, ] <- c(1, unlist(mynewdata2a))
fit@predictors[1, ] <- neweta2 # Needed
max(abs(margeff(fit, subset = 1)["let", ] - (
predict(fit, newdata = mynewdata2b, type = "response") -
predict(fit, newdata = mynewdata2a, type = "response")) / hh
)) # Should be 0
Run the code above in your browser using DataLab