10

Let me use UCLA example on multinominal logit as a running example---

library(nnet)
library(foreign)

ml <- read.dta("http://www.ats.ucla.edu/stat/data/hsbdemo.dta")
ml$prog2 <- relevel(ml$prog, ref = "academic")
test <- multinom(prog2 ~ ses + write, data = ml)

dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
predict(test, newdata = dses, "probs")

I wonder how can I get 95% confidence interval?

Johan Larsson
  • 3,496
  • 18
  • 34
user2966726
  • 101
  • 1
  • 1
  • 3

2 Answers2

7

This can be accomplished with the effects package, which I showcased for another question at Cross Validated here.

Let's look at your example.

library(nnet)
library(foreign)

ml <- read.dta("http://www.ats.ucla.edu/stat/data/hsbdemo.dta")
ml$prog2 <- relevel(ml$prog, ref = "academic")
test <- multinom(prog2 ~ ses + write, data = ml)

Instead of using the predict() from base, we use Effect() from effects

require(effects)

fit.eff <- Effect("ses", test, given.values = c("write" = mean(ml$write)))

data.frame(fit.eff$prob, fit.eff$lower.prob, fit.eff$upper.prob)

  prob.academic prob.general prob.vocation L.prob.academic L.prob.general L.prob.vocation U.prob.academic
1     0.4396845    0.3581917     0.2021238       0.2967292     0.23102295      0.10891758       0.5933996
2     0.4777488    0.2283353     0.2939159       0.3721163     0.15192359      0.20553211       0.5854098
3     0.7009007    0.1784939     0.1206054       0.5576661     0.09543391      0.05495437       0.8132831
  U.prob.general U.prob.vocation
1      0.5090244       0.3442749
2      0.3283014       0.4011175
3      0.3091388       0.2444031

If we want to, we can also plot the predicted probabilities with their respective confidence intervals using the facilities in effects.

plot(fit.eff)

Imgur

Community
  • 1
  • 1
Johan Larsson
  • 3,496
  • 18
  • 34
  • Johan, is there a way to use the effects function when you are just trying to get confidence intervals for the intercept only model (e.g., the formula was: ses ~ 1). It doesn't seem like writing just '1' or 'Intercept' in the effect function works. – Jamie Mar 04 '23 at 09:43
4

Simply use the confint function on your model object.

ci <- confint(test, level=0.95)

Note that confint is a generic function and a specific version is run for multinom, as you can see by running 

> methods(confint)
[1] confint.default   confint.glm*      confint.lm*       confint.multinom*
[5] confint.nls*

EDIT:

as for the matter of calculating confidence interval for the predicted probabilities, I quote from: https://stat.ethz.ch/pipermail/r-help/2004-April/048917.html

Is there any possibility to estimate confidence intervalls for the probabilties with the multinom function?

No, as confidence intervals (sic) apply to single parameters not probabilities (sic). The prediction is a probability distribution, so the uncertainty would have to be some region in Kd space, not an interval. Why do you want uncertainty statements about predictions (often called tolerance intervals/regions)? In this case you have an event which happens or not and the meaningful uncertainty is the probability distribution. If you really have need of a confidence region, you could simulate from the uncertainty in the fitted parameters, predict and summarize somehow the resulting empirical distribution.

nico
  • 50,859
  • 17
  • 87
  • 112
  • What I need are not CI of multinominal coefficients. Rather, I need CI for predicted probabilities--- predict(test, newdata = dses, "probs"). – user2966726 Nov 07 '13 at 22:06
  • 5
    Don't understand why Prof Ripley (author of quote above) is fussing about this. Without thinking *too* deeply, seems to me that his objections would apply to any "confidence intervals" (or whatever you want to call them), but this is a fairly standard practice. I dug into `nnet:::predict.multinom`, `nnet:::predict.nnet`, and `nnet:::vcov.multinom`. *If* you can construct model matrix `X` such that the predicted probabilities are `X %*% coef(model)`, then the variances (on the logit scale are `X %*% vcov(model) %*% t(X)`, and you would then construct CIs on the logit scale and back-transform – Ben Bolker Nov 07 '13 at 23:41
  • However, the aforementioned functions are not quite transparent enough for me to think I can hack them in a short period of time. If you can justify an ordinal response, `ordinal::predict.clm()` does allow an `se.fit` option. – Ben Bolker Nov 07 '13 at 23:42
  • @BenBolker The emmeans and development version of the marginaleffects package support nnet::multinom and mclogit::mblogit multinomial models and can calculate 95% confidence intervals on the response scale via the Delta method. Backtransforming from the logit or latent (centered logit) scale via the softmax transform in principle should be better (the inverse link function for multinomial is softmax, not inverse logit), but still implementing this in the marginaleffects package... – Tom Wenseleers Sep 04 '22 at 20:20