# modeling probability with the multinomial logit link

I am attempting to model probabilities using the multinomial logit link and I am confused about how the link works. To study the link function I have been attempting to use a deterministic system.

As an example, I am attempting to model the probability a with an intercept and a covariate (Bx.0 + Bx.1 * cov), while modeling the probability b just with an intercept. A third probability, c, is 1 - a - b.

I keep thinking that because I am modeling the probability b with just an intercept (By.0) that b should be constant, but it is not. Why is that? I suppose the obvious answer is that because Bx.0 + Bx.1 * cov is changing By.0 must change as well if b is to remain constant. Even so, it seems counter-intuitive that an intercept model (By.0) must constantly change to return a constant value (here of b). I am struggling to reconcile two ideas that seem straight-forward individually, but seem counter-intuitive when taken together.

Here is my initial example:

Bx.0 <-  0.5
Bx.1 <- -0.5
By.0 <- -0.8
cov <- -2:2

a <- (exp(Bx.0 + Bx.1 * cov) / (1 + exp(Bx.0 + Bx.1 * cov) + exp(By.0)))
b <- (exp(By.0             ) / (1 + exp(Bx.0 + Bx.1 * cov) + exp(By.0)))
c <- 1 - a - b
b
#  0.07575916 0.10781452 0.14503605 0.18344982 0.21856014


I also would like to identify the values of Bx.0 and Bx.1 that give rise to selected values of a. Here are the values of the covariate and the desired values of a I have been attempting to use to identify Bx.0 and Bx.1 algebraically. (Because I have been attempting to solve this algebraically I am unclear whether I should be posting on the Math site even though this exercise is to improve my grasp of the statistical model.)

cov = -1, a = 0.25
cov =  0, a = 0.50
cov =  1, a = 0.75


Here I define the values of a and b and obtain the values of the linear predictors:

a2   <- 0.25
b2   <- 0.20
lin.pred.a2 <- log( 1 / (((1 - a2) / a2) * (1 - b2) - b2))
lin.pred.b2 <- log( 1 / (((1 - b2) / b2) * (1 - a2) - a2))

a3   <- 0.50
b3   <- 0.20
lin.pred.a3 <- log( 1 / (((1 - a3) / a3) * (1 - b3) - b3))
lin.pred.b3 <- log( 1 / (((1 - b3) / b3) * (1 - a3) - a3))

a4   <- 0.75
b4   <- 0.20
lin.pred.a4 <- log( 1 / (((1 - a4) / a4) * (1 - b4) - b4))
lin.pred.b4 <- log( 1 / (((1 - b4) / b4) * (1 - a4) - a4))


Here are the vectors of linear predictors:

# vectors of linear predictors
a.vec   <- c(lin.pred.a2, lin.pred.a3, lin.pred.a4)
b.vec   <- c(lin.pred.b2, lin.pred.b3, lin.pred.b4)
a.vec
#  -0.7884574  0.5108256  2.7080502
b.vec
#  -1.0116009 -0.4054651  1.3862944


Here I check that my values of a and b match what I intended, which they do:

a.pred <- exp(a.vec) / (1 + exp(a.vec) + exp(b.vec))
a.pred
#  0.25 0.50 0.75

b.pred <- exp(b.vec) / (1 + exp(a.vec) + exp(b.vec))
b.pred
#  0.2 0.2 0.2


However, I remain unclear why the intercept of the linear predictor for b has three different values (-1.0116009, -0.4054651, 1.3862944) in order to return a constant value for b (0.2, 0.2, 0.2).

Furthermore, given that the intercept for the linear predictor of b varies, I suppose the intercept for the linear predictor for a probably varies as well. In which case is it possible to identify the values of Bx.0 and Bx.1 algebraically that give rise to the probabilities of a = 0.25, 0.50 and 0.75? If so, how can I do that?

In other words, I was thinking of Bx.0 + Bx.1 * cov as a linear regression, but if Bx.0 and Bx.1 take on different values for each value of a (0.25, 0.50 and 0.75) can I identify Bx.0 and Bx.1 algebraically?

Thank you for any advice or suggestions. Sorry if this post is not clear. I am familiar with multinomial logistic regression and ordinal logistic regression, but I am not certain they apply here, at least not to obtain an algebraic solution.

I have not worked out all of my concerns. However, I have worked out the connection between the multinomial logit link in my original post and multinomial logistic regression.

Suppose there are three state: a, b and c. The relationship between probabilities of states a and b to the linear predictor for state b as used in multinomial logistic regression is:

$\log(\frac{P (b)}{ P (a)}) = \beta_{0b} + \beta_{1b}X + \beta_{2b}X^2$

$\log(\frac{P (c)}{ P (a)}) = \beta_{0c} + \beta_{1c}X + \beta_{2c}X^2$

which equals:

$e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} = \frac{P (b)}{ P (a)}$

$e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2} = \frac{P (c)}{ P (a)}$

and reorganizing:

$P(a) e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} = P (b)$

$P(a) e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2} = P (c).$

Now add probabilities for all three states to obtain 1:

$P(a) + P(a) e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} + P(a) e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2} = 1.$

Now factor out $P(a)$:

$P(a) (1 + e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} + e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2}) = 1.$

and divide to isolate $P(a)$:

$P(a) = \frac{1} {1 + e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} + e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2}}.$

The above equation for $P(a)$ is the same as the one I used in my original post.

To obtain the equation for $P(b)$ used in my original post substitute the equation immediately above for $P(a)$ into the equation given a few lines up for $P(b)$:

$P(b) = \frac{1} {1 + e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} + e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2}} e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2}.$

Do the same for $P(c)$:

$P(c) = \frac{1} {1 + e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} + e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2}} e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2}.$

Lastly, recall that $P(a)$ is:

$1 - \frac{1} {1 + e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} + e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2}} e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} - \frac{1} {1 + e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} + e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2}} e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2}$

because:

$\frac{1} {1 + e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} + e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2}} + \frac{1} {1 + e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} + e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2}} e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} + \frac{1} {1 + e^{\beta_{0b} + \beta_{1b}X + \beta_{2b}X^2} + e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2}} e^{\beta_{0c} + \beta_{1c}X + \beta_{2c}X^2} = 1$

Here is an example of multinomial logistic regression in R code showing how the parameter estimates relate to the probabilities:

library(nnet)

state   cov     cov2
a       0       0
a       0       0
a       0       0
a       0       0
b       1       1
c       1       1
c       1       1
c       1       1
c       1       1
c       1       1
a       0       0
a       0       0
a       0       0
a       0       0
b       2       4
b       2       4
c       2       4
c       2       4
c       2       4
c       2       4
a       0       0
a       0       0
a       0       0
a       0       0
b       3       9
b       3       9
b       3       9
b       3       9
c       3       9
c       3       9
a       0       0
a       0       0
a       0       0
a       0       0
b       4      16
b       4      16
b       4      16
b       4      16
b       4      16
c       4      16
', header = TRUE, stringsAsFactors = FALSE)

my.data$state <- as.factor(my.data$state)
my.data$state2 <- relevel(my.data$state, ref = "a")
model1 <- multinom(state2 ~ cov + cov2, data = my.data)
summary(model1)

# Call:
# multinom(formula = state2 ~ cov + cov2, data = my.data)
#
# Coefficients:
#   (Intercept)      cov     cov2
# b  -11.644219 20.40486 1.025687
# c   -8.813602 19.25966 1.030156
#
# Std. Errors:
#   (Intercept)      cov     cov2
# b    20.11729 69.71465 69.69290
# c    19.92642 69.69671 69.69628
#
# Residual Deviance: 26.15114
# AIC: 38.15114
#

fitted <- fitted(model1)
fitted[1:10,]

# probability of b
round(log(1.557895e-01 / 8.757683e-06), 4) == round(-11.644219 + 20.40486 * 1 + 1.025687 * 1^2, 4)

# probability of c
round(log(0.8442017620 / 8.757683e-06), 4) == round(-8.813602  + 19.25966 * 1 + 1.030156 * 1^2, 4)


This post shows, I believe, that the equations and parameters used in multinomial logistic regression are the same as those I presented in my original post for the multinomial logistic link.

However, I still need to study both approaches to understand my original questions of 1). how to obtain estimates of the intercept and slope parameters algebraically (that might not be possible), 2). whether the intercept and slope estimates of b must vary as the covariate varies and 3). whether the intercept of a must vary as the covariate in b varies.

If I make progress in my understanding of those three questions I will post an update.

FIRST EDIT - Nov 28, 2015

Model estimates seem to be better if I include covariates for State a, even though I want $P(a)$ to be a constant 0.40 in this example. Note the nice standard errors.

library(nnet)

state   cov     cov2
a       1       1
a       1       1
a       1       1
a       1       1
b       1       1
c       1       1
c       1       1
c       1       1
c       1       1
c       1       1
a       2       4
a       2       4
a       2       4
a       2       4
b       2       4
b       2       4
c       2       4
c       2       4
c       2       4
c       2       4
a       3       9
a       3       9
a       3       9
a       3       9
b       3       9
b       3       9
b       3       9
b       3       9
c       3       9
c       3       9
a       4      16
a       4      16
a       4      16
a       4      16
b       4      16
b       4      16
b       4      16
b       4      16
b       4      16
c       4      16
', header = TRUE, stringsAsFactors = FALSE)

summary(model1)
# Call:
# multinom(formula = state2 ~ cov + cov2, data = my.data)
#
# Coefficients:
#   (Intercept)       cov      cov2
# b  -2.5666228 1.2589663 -0.139625
# c   0.2370914 0.1380369 -0.139825
#
# Std. Errors:
#   (Intercept)      cov      cov2
# b    2.758592 2.250422 0.4171805
# c    2.048476 1.993598 0.4172125
#
# Residual Deviance: 79.98684
# AIC: 91.98684
#

fitted <- fitted(model1)
fitted[1:10,]

# probability of b
round(log(0.09406466 / 0.3999198), 4) == round(-2.5666228 + 1.2589663 * 1 + -0.139625 * 1^2, 4)

# probability of c
round(log(0.50601555 / 0.3999198), 4) == round(0.2370914  + 0.1380369 * 1 + -0.139825 * 1^2, 4)