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

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?

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"]]
#>  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.

• I compared 1/41=0.024390243902439 to the estimate (0.0243902438057923) from the regression model with reltol = 1e-100. It seems that there is still slight inconsistency. I understand this inconsistency in practice can be ignorable. However, given this example, it seems that the model estimate is in general unable to PERFECTLY predict the actual data (which is used to fit the model). Given this inconsistency, I am still curious if I can use regression model to estimate the probability of the outcome and its confidence interval, instead of manual calculation (like using 1/41). Thank you! Aug 10 at 7:59
• "I am still curious if I can use regression model to estimate the probability of the outcome and its confidence interval" You can. If you think an error in the 10th decimal point invalidates your inferences, you should quit statistics today.
– Noah
Aug 10 at 22:40

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.)

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

• Thank you for your helpful answer. Because my goal is to estimate the probability in the titanic_train dataset, not in the titanic_test dataset. I am wondering if I can still use the overfitted model to estimate the probability in the dataset where the model was fitted. If not, what could be the reasons not to do so? Aug 8 at 11:22
• Well ... I guess if that's your goal, OK, but I don't think it's a sensible goal (sorry). First, the usual reasons to do a regression are 1) Explanation and 2) Inference (onto some other data). If you want to estimate the probability in the actual data, why model? You have the data! Both 1) and 2) are compromised by overfitting. Running the model on test data shows that inference is problematic. But the parameters will also be off. Second, one big purpose of dividing into train and test is precisely to avoid this. Aug 8 at 11:31
• You wrote: "In the real dataset I am using," PLEASE don't ask questions about one data set when you want to analyze another. If you can't tell us all the details, fine, but tell us what you can tell us. And tell us in the question1 Your real data set has missing data; AFAIK, the Titanic data set does not. And, the problems that arose from your analysis of Titanic data were specific to that data set. I think you should ask a new question, with your real data set. Aug 8 at 14:39