Multinomial Logistic Regression - Independent Regressions

Question

Given the standard multinomial logistic regression model with reference class K:

$P(Y=j|X^{(i)}) = \frac{\exp(\theta_j^TX^{(i)})}{1+ \sum_{m=1}^{K-1}\exp(\theta_m^T X^{(i)})}$

Would I run into any problems if I estimate K-1 binary regressions independently (R: glm(,family = "binomial") vs. nnet::multinom()) on filtered data sets where the outcomes are K and j? This implies there will be different sets of predictors and possibly different transformations (e.g. floors, caps, splines, etc.).

Answer 1

You will get close, except that multinomial logit will enforce the constraint that all predicted probabilities add up to one, while different logistic regressions cannot enforce that constraint.

Answer 2

Just for the record - apparently my original question referred to a methodology first described in an article by Begg and Gray "Calculation of polychotomous logistic regression parameters using individualized regressions" from 1984. Hosmer mentions it on page 277 of Applied Logistic Regression...

Answer 3

I think you made a small typo. The sum in the denominator should be running from $1$ to $k-1$, not $k$. Edit: nevermind you fixed it :)

Fitting $k-1$ independent regressions, you get estimates for $\{\theta_j\}$ with $j=1, \ldots, k-1$. The $k-1$ models are parametrized as: $$ P(Y=j|X^{(i)}) = \frac{\exp(\theta_j^TX^{(i)})}{1+ \exp(\theta_j^T X^{(i)})}. $$ Looking at each row, the dependent variable will either be a succes ($=j$) or not $(\neq j)$. You're basically ignoring all the other outcomes here.

On the other hand if you fit a multinomial logistic regression, it will enforce, as @MaartenBuis says, the constraint that all the probabilities sum to $1$. You're also assuming, probably more correctly, that each dependent variable can be any value from $1$ to $k$, and it will follow a multinomial distribution with the parameter vector being dependent on the $X^{(i)}$ data. In other words, the model assumed is $$ P(Y=j|X^{(i)}) = \frac{\exp(\theta_j^TX^{(i)})}{1+ \sum_{m=1}^{k-1}\exp(\theta_m^T X^{(i)})} $$ with $j=1,\ldots,K-1$, and you still get the same amount of coefficient estimates. But you can see that $$ \sum_{i=1}^k P(Y=i|X^{(i)}) = \sum_{i=1}^{k-1} \frac{\exp(\theta_j^TX^{(i)})}{1+ \sum_{m=1}^{k-1}\exp(\theta_m^T X^{(i)})} + \frac{1}{1+\sum_{m=1}^{k}\exp(\theta_m^T X^{(i)})} = 1. $$

I can't prove when the estimates between these two procedures will be close or not, but I can tell you that these two procedures are assuming totally different models to be true.