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))
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.
What are the possible reasons for this slight inconsistency?
Given the inconsistency, can I continue using regression model to calculate probability and confidence interval?