Slight inconsistency between probability estimated from logistic regression and manual calculation from cross table

Question

I want to estimate (using Titanic data as an example) the probability of being in each of the classes (class 1, class 2, class 3) by sex (male, female) and port of embarkation (C, Q, S). I used multinomial logistic regression model to estimate the probability of each class among each group defined by the interaction of sex and embarkation. The R scripts are as below:

library(titanic)
library(nnet)

d1 <- titanic_train
d1 <- subset(d1, d1$Embarked != "")
d1$Pclass <- factor(d1$Pclass, levels=c("1", "2", "3"))
d1$Sex <- factor(d1$Sex, levels=c("female", "male"))
d1$Embarked <- factor(d1$Embarked)

mod <- multinom(Pclass ~ Sex + Embarked + Sex*Embarked, data=d1)
summary(mod)

The output is:

I calculated, for example, the probability of being in class 2 among EmbarkedQ and male:

exp(-1.8151727 + 0.3801844 + 2.507507 - 1.0735799)/(1 + exp(-1.8151727  
    + 0.3801844 + 2.507507 - 1.0735799) + exp(-0.6257148 + 0.6492546 + 
      4.121453 - 0.4819176))

= 0.02437628

On the other hand, I manually calculated the probability for the same group:

table(d1$Pclass, d1$Sex, d1$Embarked)

1/(1 + 1 + 39) = 0.02439024

The two approaches yielded very close results and practically negligible, but are not exactly the same. In theory, to my understanding, the probability from the manual calculation should be the correct one. I am concerned about the result from the regression model.

Questions:

1. What are the possible reasons for this slight inconsistency?

2. Given the inconsistency, can I continue using regression model to calculate probability and confidence interval?

Answer 1

As Henry mentioned, this is because multinom() uses a somewhat relaxed convergence threshold so that the algorithm stops before reaching to true maximim likelihood. If you set reltol to a lower value, e.g., to 1e-15, the model predictions are equal to the observed proportions:

mod <- multinom(Pclass ~ Sex + Embarked + Sex*Embarked, data=d1,
                reltol = 1e-15)

predict(mod, newdata = data.frame(Sex = "male", Embarked = "Q"),
        type = "probs")[["2"]]
#> [1] 0.02439024

In practice, this difference is meaningless, but one should get in the habit of adjusting the tolerance on solvers. Defaults often lead to fast but approximate performance, and it can be beneficial to examine lower tolerances when results are suspect.

Answer 2

Logistic regression uses maximum likelihood methods, and presumably the algorithm stops when the optimisation is "close enough". Your $0.02437628 \approx \frac{1}{41.023}$, which is close to your desired $\frac1{41}$.

You would find similar results for the other observations. If you converted these into (log-)likelihoods, the actual proportions would give a marginally higher overall likelihood but without a substantial difference.

A more relevant issue is the small number of 1st and 2nd class passengers embarking at Queenstown as your table shows, which if this is regarded as a random process means there is considerable uncertainty in the estimated probabilities, far larger than any optimisation errors. The large standard errors reported for some of the coefficients hint at this.

(Incidentally, my great-grandmother, great-aunt and great-uncle and others embarked on the Titanic at Southampton and disembarked at Cherbourg - whether they should be included as survivors is another methodological issue.)

Answer 3

There isn't any inconsistency. The results of a model aren't going to exactly match the data. The interesting thing here is how close they come. So, the key word is "slight".

One reason it is so small is that you have a problem with the data -- of 6 cells, 4 have either 1 or 2 people. That's a recipe for overfitting. I think R should have given some kind of warning.

As for your question 2: I'd say "no" but not because the results don't match perfectly, but because the model is overfit.

One way to look for overfitting is to run your model on titanic_test.

testpred <- predict(mod, newdata=titanic_test)
table(testpred,titanic_test$Pclass)

and it doesn't do nearly as well as it did on train.

testpred   1   2   3
       1  28   4   8
       2   0   0   0
       3  79  89 210