```r
predict.mnl <- function(model, data) {
  data.model <-
    model.matrix(
      update(model$formula, 0 ~ .),
      data = data)[ , -1]
  utility <- data.model %*% model$coef
  share <- exp(utility) / sum(exp(utility))
  cbind(share, data)
}
```