# Details of binary logistic regression, estimating $P(Y=1|X)$

I am trying to understand logistic regression, but most sources I have found tend to leave the actual computational step as sort of a "black box", like using r glm(y~x,...) which obscures the underlying compuation. As an exercise, I want to write my own routine for doing (binary) logistic regression which requires me knowing the details.

For binary outcomes $$Y$$, and some predictor data $$X$$, we attempt to model the conditional probability

$$p(X) = P(Y=1|X) = \frac{1}{1+ e^-{\beta X}}$$

and the corresponding linear model

$$\mathcal{l} = \text{log}\bigg( \frac{p}{1-p} \bigg) = \beta X$$

How do we now proceed in practice to compute $$\beta$$? In my understanding we have something like a linear system, so that for $$n$$ measurements & $$k$$ predictive variables

$$l_1 = \beta_0 + \beta_1x_{11} + ... +\beta_k x_{1k}$$ $$l_2 = \beta_0 + \beta_1x_{21} + ... +\beta_k x_{21}$$  $$l_n = \beta_0 + \beta_1x_{n1} + ... +\beta_k x_{nk}$$

which could then be solved using standard methods (such as SGD for an appropriate cost function). But how do we compute/estimate the values in $$l_i$$? The observed value of $$P(Y=1|X)$$ would simply be the scalar $$p* = \frac{\sum_{i=1}^n y_i}{n}$$ for each row, meaning that the system to be solved would be

$$\text{log}(\frac{p*}{1-p*})\times(1,1,....,1)^T = \beta X$$

Is this correct? I have been testing with the following r code:

set.seed(1234)
x <- rnorm(1000)
y <- rbinom(1000, 1, exp(-2.3 + 0.1*x)/(1+exp(-2.3 + 0.1*x)))

fit = glm(y ~ x, family=binomial(link='logit'))
fit$coefficients p = sum(y)/length(y) z = rep(p, 1000) z_star = log(z/(1-z)) fit2 = lm(z_star ~ x) fit2$coefficients


which gives the coefficients estimates for $$\beta$$

> fit$coefficients (Intercept) x -2.2261215 0.1651474 > fit2$coefficients
(Intercept)             x
-2.219647e+00 -5.338566e-16


which differ noticably in the estimates of the $$x$$ coefficent. I thought this may be becuse of the different optimization methods used in the glm() vs lm() methods, but is my understand of the modelling procedure correct?

• You can have a read of en.wikipedia.org/wiki/Logistic_regression#Model_fitting to see how we fit the model (obtain estimates for $\beta$). Essentially, we model $y_i\mid x_i, \beta$ as being $\mathsf{Bernoulli}(p_i)$, where $p_i = \frac{1}{1+e^{-\beta^T x_i}}$. From this, we form the log-likelihood of the data and maximise this with respect to $\beta$ (i.e. we perform maximum likelihood estimation). The equations for maximising the log-likelihood are solved numerically (there is no closed-form solution for them in this case). Jun 18 '19 at 12:47
• Looks like there is a mistake in your code. P= sum(y)/length(y) assumes that p is the same for all observations. This is why the relationship with x in the linear model disappears.
– jsk
Jun 18 '19 at 16:11
• @jsk That is part of my question, is this not the case? The scalar value computed above is the same for each row, since it is just the sum of 1s in Y divided by n Jun 18 '19 at 18:09
• No, there's problems with the notation used that is creating the confusion. Every observation has a different $p_i$ that changes based on $x_i$.
– jsk
Jun 18 '19 at 20:04

The fit from the linear model is very different because it's a linear probability model, i.e. you're directly estimating

$$P(Y=1|X=x) = \beta_0 + \beta_1 x$$

whereas in logistic regression you're modeling the log odds as a linear model

$$\log \left( \frac{P(Y=1|X=x)}{P(Y=0|X=x)} \right) = \beta_0 + \beta_1 x$$

which is very different, i.e., the coefficient $$\beta_1$$ has a totally different meaning.

Also, to clarify: you estimate the coefficients in a logistic regression model by maximum likelihood--i.e. by maximizing the binomial (actually Bernoulli) likelihood, not by solving a system of equations.

• +1 Nicely and simply explained. Welcome to the site; hope to see more contributions from you.
– EdM
Jun 18 '19 at 15:33
• In the setup/code above you'll see that I find a linear model for the log-odds, not the probability itself. What is unclear is the part in bold, how do we compute the values for log(p/(1-p)) to be used in the model for the log-odds? Jun 18 '19 at 18:12

Generally, we need to get the maximum likelihood estimate with following steps:

1) Write the log-likelihood function $$L(\beta)$$.

2) Get $$\partial {L(\beta)}\over \partial {\beta}$$.

3) Get $$\partial^2 {L(\beta)}\over \partial {\beta^2}$$.

4) Use Newton–Raphson method to get the maximum likelihood estimate (MLE) of $$\beta$$.

The $$l_i$$ do not correspond to anything in your dataset. Logistic regression is fit via Maximum Likelihood which seeks to assign $$p^*$$ to observations such that the binomial log likelihood is maximized.

That is, given covariates $$x$$ and a binary outcome, we can assign a probability $$p^*$$ that future observations with covariates $$x$$ will yield an outcome of 1.

So we don't fit the $$l_i$$, we fit the $$p^*$$.