%run functions.ipynb

Functions loaded.

13. Logistic regression

Linear is additive, doesnt work for binary outcomes (0 or 1). To model binary we need to add two features to base model: $y = a + bx$: a non-linear transformation to bound output between 0 and 1 and a model that trats output as a probability and maps into random binary outcomes.

13.1 Logistic regression with a single predictor

Logit function:

$logit(x) = log(\frac{x}{1-x})$ - maps the range (0,1) to the range (-inf, inf) and is useful for modeling probabilities. The inverse logit function is: $logit^{-1}(x) = \frac{1}{1 + e^{-x}}$ - maps back to the unit range.

Example: modelling political preference given income

Conservatives traditionally more support with voters with higher income. National election data from 1992. For each respondent i, we label $y_i = 1$ if they voted for Bush and $y_i = 0$ if they voted for Clinton and exclude any others. Predict preferences given respondents income on a 5 point scale.

# load nes5200 .dta and replicate nes_setup.R preprocessing
brdata = pd.read_stata("../ros_data/nes5200_processed_voters_realideo.dta")
brdata = brdata.dropna(subset=["black", "female", "educ1", "age", "income", "state"])
nes = brdata[brdata["year"].between(1952, 2000)].copy()

# extract numeric codes from categorical labels (e.g. "1. democrat" -> 1)
nes["presvote_code"] = nes["presvote"].astype(str).str[0].astype(float)
nes["income"] = nes["income"].astype(str).str[0].astype(int)

# correct rvote: 0=Clinton (dem), 1=Bush (rep), NaN=other
nes["rvote"] = nes["presvote_code"].map({1.0: 0, 2.0: 1})

nes92 = nes[(nes["year"] == 1992) & nes["rvote"].notna()].copy()
nes92.head()

	year	resid	weight1	weight2	weight3	age	gender	race	educ1	urban	...	perfin	presadm	age_10	age_sq_10	newfathe	newmoth	parent_party	white	presvote_code	rvote
32092	1992.0	1.0	1.0	1.0	1.0	41.0	2. female	4. native american	2. high school (12 grades or fewer, incl	2. suburban areas	...	1.0	1.0	4.1	16.809999	1.0	1.0	2.0	0.0	2.0	1.0
32093	1992.0	2.0	1.0	1.0	1.0	83.0	2. female	1. white	3. some college(13 grades or more,but no	3. rural, small towns, outlying and adja	...	2.0	1.0	8.3	68.889999	-1.0	-1.0	-2.0	1.0	2.0	1.0
32095	1992.0	4.0	1.0	1.0	1.0	67.0	2. female	2. black	2. high school (12 grades or fewer, incl	1. central cities	...	2.0	1.0	6.7	44.889996	-1.0	-1.0	-2.0	0.0	1.0	0.0
32096	1992.0	6.0	1.0	1.0	1.0	28.0	1. male	1. white	4. college or advanced degree (no cases	2. suburban areas	...	2.0	1.0	2.8	7.840000	-1.0	-1.0	-2.0	1.0	2.0	1.0
32097	1992.0	8.0	1.0	1.0	1.0	71.0	2. female	4. native american	2. high school (12 grades or fewer, incl	2. suburban areas	...	2.0	1.0	7.1	50.410000	NaN	NaN	NaN	0.0	1.0	0.0

5 rows × 64 columns

sns.lmplot(data=nes92, x='income', y='rvote', logistic=True, scatter_kws={'s': 1,'alpha': 0.3}, x_jitter=0.2, y_jitter=0.05)

<seaborn.axisgrid.FacetGrid at 0x1376f7fe0>

# drop nans and reformat output and predictor
df = nes92[["income", "rvote"]].dropna()
x = df["income"].values.astype(int)
y = df["rvote"].values.astype(int)

# fit logistic regression model
with pm.Model() as fit_1:
        
    # priors
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)
    income = pm.Normal("income", mu=0, sigma=2.5)
    
    # logistic likelihood
    p = pm.math.invlogit(intercept + income * x)
    pm.Bernoulli("rvote", p=p, observed=y)
    
    trace = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

summary = pm.summary(trace)
print(summary)
# az.summary(trace, var_names=['intercept', 'income'])

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, income]

Output()

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 1 seconds.

            mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \
intercept -1.402  0.195  -1.751   -1.015      0.007    0.005     784.0   
income     0.326  0.059   0.204    0.426      0.002    0.002     795.0   

           ess_tail  r_hat  
intercept     773.0   1.01  
income        803.0   1.01  

We write this as $Pr(y=1) = logit^{-1}(-1.438 + 0.262 \cdot income)$, which means that the log odds of voting for Bush increases by 0.262 for each unit increase in income. There is no sigma term - the uncertainty is captured in its probabilistic prediction of binary outcomes.

from scipy.special import expit  # equivalent to invlogit

# get all posterior samples
alpha_samples = trace.posterior['intercept'].values.flatten()
beta_samples  = trace.posterior['income'].values.flatten()

x_range = np.linspace(df['income'].min(), df['income'].max(), 200)

# posterior means for the coefficients
alpha_hat_mean = trace.posterior['intercept'].mean().item()
beta_hat_mean = trace.posterior['income'].mean().item()

# scatter plot
x_jitter = np.random.uniform(-0.2, 0.2, size=len(df))
y_jitter = np.random.uniform(-0.02, 0.02, size=len(df))
plt.scatter(df['income'] + x_jitter, df['rvote'] + y_jitter, alpha=0.3, s=5)

# plot all posterior curves
for a, b in zip(alpha_samples, beta_samples):
    plt.plot(x_range, expit(a + b * x_range), color='black', alpha=0.01, linewidth=0.5)

# logistic curve
x_range = np.linspace(df['income'].min(), df['income'].max(), 200)
plt.plot(x_range, expit(alpha_hat_mean + beta_hat_mean * x_range), color='red')

plt.xlabel('income')
plt.ylabel('rvote')
plt.show()

The logistic regression model

We want to model an outcome bounded between 0 and 1, not continuous, as we could get negative probabilities or probabilities above 1. We also loose data if we pretend that binary is a continuous variable.

Instead of modelling y, we model probability that y=1:

$Pr(y=1) = logit^{-1}(X_i \beta)$

The $logit^{-1}$ (inverse logit or logistic function) transforms $X\beta$ (any number between -∞ and ∞) to between 0 and 1. So output is always a valid probability. $X\beta$ is the linear predictor.

Two ways to write the same thing:

$Pr(y=1) = logit^{-1}(X_i \beta)$ - go from linear predictor to probability $logit(p) = X_i \beta$ - go from probability to linear predictor

Key intuition: the curve is steepest at the middle (p=0.5). Because logistic function is S-shaped, the same step on logit scale does not always correspond to the same step on probability scale. It moves a lot in the middle, but very little at the tails.

Near middle: +0.4 on logit scale goes from 50% to 60% Near top: +0.4 on logit scale goes from 90% to 93% Near bottom: +0.4 on logit scale goes from 7% to 10%

Curve is steepest at the middle and flattens out at the tails. So a change in logit scale has a bigger effect on probability when we are near the middle than when we are near the tails.

Important: $\beta$ is a constant effect on the log-odds scale / logit scale, but not on the probability scale. Its effect on probability depends on where we are on the curve - we have to have the base probability to know how much a change in logit scale will change the probability. So we can say that a change in logit scale has a bigger effect on probability when we are near the middle than when we are near the tails.

n_sims = len(alpha_samples)
for j in np.random.choice(n_sims, size=20,replace=False):
    plt.plot(x_range, expit(alpha_samples[j] + beta_samples[j] * x_range), color='gray', linewidth=0.5)

13.2 Interpreting logistic regression coefficients and the divide-by-4 rule

Because of the nonlinearity, logistic regression coefficients are not easy to interpret.

Interpretation is explained using the model: Pr(Bush support)=logit−1(−1.40 + 0.33 ∗income)

Evaluation at and near the mean of the data

Due to nonlinearity we need to pick where to evaluate changes if we want to interpret on probability scale. We can evaluate at the mean of the data, which is a common choice.

1. Intercept: as with linear, we interpret the intercept as the predicted value when all predictors are 0. But with Bush example, there is no 0, so we can interpret for highest income group (income=5) and lowest income group (income=1). For lowest income group, logit is -1.40 + 0.331 = -1.07, which corresponds to a probability of 0.25. For highest income group, logit is -1.40 + 0.335 = 0.25, which corresponds to a probability of 0.56. We can also evaluate at the mean of the data, which is 3.1 for income. At mean, logit is -1.40 + 0.33*3.1 = -0.37, which corresponds to a probability of 0.41.
2. A difference of 1 in income corresponds to a change of 0.33 in logit probability of supporting Bush. We can evaluate this at the mean of the data, which is 3.1. At mean, logit is -1.40 + 0.333.1 = -0.37, which corresponds to a probability of 0.41. If we increase income by 1, logit becomes -1.40 + 0.334.1 = -0.04, which corresponds to a probability of 0.49. So a change of 1 in income corresponds to a change of 0.08 in probability of supporting Bush at the mean of the data.

The divide-by-4 rule

Logistic curve is steepest at the middle - the point $\alpha + \beta x = 0$ (where p=0.5). The slope of the curve - derivative of the logistic function is $\beta e^0 / (1 + e^0)^2 = \beta / 4$. So the maximum effect of a change in logit scale on probability scale is $\beta / 4$. So we can divide $\beta$ by 4 to get an upper bound on how much a change in logit scale will change probability. This is a useful rule of thumb for interpreting logistic regression coefficients. For the Bush example, $\beta$ is 0.33, so the maximum effect on probability is 0.33 / 4 = 0.08, which is consistent with our evaluation at the mean of the data. So we can say that a change of 1 in income corresponds to a change of at most 0.08 in probability of supporting Bush.

Interpretation of coefficients as odds ratios

We can also interpret logistic regression coefficients as odds ratios. If two outcomes have the probabilities (p, 1-p) then $p/(1-p)$ is called the odds. Odds of 1 is the same as a probability of 0.5 - equal likelihood of both outcomes. Odds of 0.5 or 2 represent the probabilities (1/3 and 2/3) - one outcome is twice as likely as the other. Dividing two odds gives us an odds ratio $\frac{p_1/(1-p_1)}{p_2/(1-p_2)}$. A odds ratio of 2 corresponds to a change from p=0.33 to p=0.5 or from p=0.5 to p=0.67.

An advantage of odds ratios (instead of probabilities) is that they can be scaled indefinitely without hitting the boundary of 0 or 1. For example going from odds 2 to 4 increases probability from 2/3 to 4/5 and doubling the odds again to 8 increases probability from 4/5 to 8/9.

Exponentiated logistic regression coefficients can be interpreted as odds ratios.

Odds and odds ratios can be difficult to interpret. Authors prefer to use original scale. For example, adding 0.2 on logit scale corresponds to a change in probability from $logit^{-1}(0)$ to $logit^{-1}(0.2)$.

Coefficient estimates and standard errors

Coefficients in classical logistic regression are estimated using maximum likelihood estimation (MLE) - works well for models with few predictors and large datasets. As with linear, standard errors represent estimation uncertainty.

We can roughly say that coefficient estimates within 2 SE of $\hat{\beta}$ are consistent with the data. But this is a rough rule of thumb and should not be taken too literally. We can also use confidence intervals to get a more precise estimate of the range of values that are consistent with the data.

Statistical significance

As with linear, an estimate is conventionally labeled as statistically significant if it is more than 2 standard errors away from 0 or more formally if the 95% confidence interval does not include 0.

Authors do not recommend using statistical significance about whether to include predictors or interpret coefficients as real. Statistical significance could be from chance alone and non-significant results could be from a real effect. Standard errors convey uncertainty in estimation without trying to make a one-to-one mapping between statistical significance and real effects.

When considering multiple predictors, we should follow the same approach as with linear regression.

Displaying the results of several logistic regressions

As with linear, we can display estimates from a series of logistic regressions in a single plot.

13.3 Predictions and comparisons

Logistic regression predictions are probabilistic, so for each unobserved data point, we get a predicted probability rather than a point prediction $p_i^new = Pr(y_i^new = 1) = logit^{-1}(X_i^new \beta)$.

Point prediction

# Predicting a new data point
new_income = 5  # example income value
new_p_samples = expit(alpha_samples + beta_samples * new_income)
print(f"Predicted probabilities of voting for Bush at income={new_income}: {new_p_samples[:5]}...")  # show first 5 samples
new_p_mean = new_p_samples.mean()
print(f"Predicted probability of voting for Bush at income={new_income}: {new_p_mean:.3f}")
print(f"95% CI: [{np.percentile(new_p_samples, 2.5):.3f}, {np.percentile(new_p_samples, 97.5):.3f}]")

Predicted probabilities of voting for Bush at income=5: [0.55124018 0.54722858 0.57218149 0.57138985 0.53081649]...
Predicted probability of voting for Bush at income=5: 0.557
95% CI: [0.500, 0.618]

Linear predictor and expected outcome with uncertainties

# Linear predictor (logit scale) for a new observation
new_income = 5
linpred = alpha_samples + beta_samples * new_income  # shape: (n_draws,)

print(linpred[:5])  # show first 5 linear predictor samples
print(f"Linear predictor (logit) for income={new_income}: mean={linpred.mean():.3f}, 95% CI=[{np.percentile(linpred, 2.5):.3f}, {np.percentile(linpred, 97.5):.3f}]")

p_pred = expit(linpred)  # shape: (n_draws,) — equivalent to posterior_epred in R

print(f"Predicted probability of voting for Bush at income={new_income}: mean={p_pred.mean():.3f}, S.D={p_pred.std():.3f}, 95% CI=[{np.percentile(p_pred, 2.5):.3f}, {np.percentile(p_pred, 97.5):.3f}]")

[0.2056828  0.1894792  0.29075719 0.28752395 0.12342238]
Linear predictor (logit) for income=5: mean=0.229, 95% CI=[-0.002, 0.479]
Predicted probability of voting for Bush at income=5: mean=0.557, S.D=0.030, 95% CI=[0.500, 0.618]

Logistic regression with just an intercept

Linear regression with just an intercept is the average and with single binary predictor is the difference in means. Similarly, with logistic it is the estimate of a proportion.

Simple example: random sample of 50 people tested and 10 have a disease. Proportion is 10/50 = 0.2 with standard error sqrt(0.2*0.8/50) = 0.056.

Model below

y = np.repeat([0, 1], [40, 10])
print(y)

# fit logistic regression model (intercept only)
with pm.Model() as fit_1:

    # prior
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)

    # logistic likelihood
    p = pm.math.invlogit(intercept)
    pm.Bernoulli("y", p=p, observed=y)

    trace = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

summary = pm.summary(trace)
print(summary)

Initializing NUTS using jitter+adapt_diag...

[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 1 1 1 1 1 1 1 1 1 1]

Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept]

Output()

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 0 seconds.

            mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \
intercept -1.386  0.352  -2.058   -0.775      0.009    0.006    1384.0   

           ess_tail  r_hat  
intercept    1388.0    1.0  

# get all posterior samples
alpha_samples = trace.posterior['intercept'].values.flatten()

print(alpha_samples[:5])  # show first 5 samples

mean_alpha = alpha_samples.mean()

p_pred = expit(alpha_samples)  # shape: (n_draws,) — equivalent to posterior_epred in R

print(f"Predicted probability={new_income}: mean={p_pred.mean():.3f}, S.D={p_pred.std():.3f}, 95% CI=[{np.percentile(p_pred, 2.5):.3f}, {np.percentile(p_pred, 97.5):.3f}]")

[-1.63977661 -2.45984146 -2.29669381 -2.04400165 -1.5051266 ]
Predicted probability=5: mean=0.206, S.D=0.056, 95% CI=[0.109, 0.321]

intercept -1.386 0.352

We can transform the prediction to probability scale: $logit^{-1}(-1.386) = 0.2$, which is the same as the proportion of people with the disease in the sample. We can also transform the standard error on logit scale to probability scale using the delta method, which gives us a standard error of 0.056, which is the same as the standard error of the proportion.

Logistic regression with a single binary predictor

Logistic regression on an indicator variable is equivalent to comparison of proportions.

Simple example: Test for disease on two populations. Population A 10 out of 50 have the disease and population B 20 out of 60 have the disease. Proportion in A is 10/50 = 0.2 with SE sqrt(0.20.8/50) = 0.056 and proportion in B is 20/60 = 0.33 with SE sqrt(0.330.67/60) = 0.061. The difference in proportions is 0.33 - 0.2 = 0.13 with SE sqrt(0.056^2 + 0.061^2) = 0.083.

x = np.repeat([0, 1], [50, 60])
print(x)
y = np.repeat([0, 1, 0, 1], [40, 10, 40, 20])
print(y)

# fit logistic regression model (intercept only)
with pm.Model() as fit_1:
    # prior
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)
    x_co = pm.Normal("x", mu=0, sigma=2.5)

    # logistic likelihood
    p = pm.math.invlogit(intercept + x_co * x)
    pm.Bernoulli("y", p=p, observed=y)

    trace = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

summary = pm.summary(trace)
print(summary)

# get all posterior samples
alpha_samples = trace.posterior['intercept'].values.flatten()
beta_samples  = trace.posterior['x'].values.flatten()

print(alpha_samples[:5])  # show first 5 samples
print(beta_samples[:5])  # show first 5 samples

mean_alpha = alpha_samples.mean()
mean_beta = beta_samples.mean()

print(f"Posterior mean of intercept: {mean_alpha:.3f}")
print(f"Posterior mean of x coefficient: {mean_beta:.3f}")

p_pred_1 = expit(alpha_samples + beta_samples * 1)  # predicted logit for x=1
p_pred_0 = expit(alpha_samples + beta_samples * 0)  # predicted logit for x=0

print(f"Predicted probability for x=1: mean={p_pred_1.mean():.3f}, S.D={p_pred_1.std():.3f}, 95% CI=[{np.percentile(p_pred_1, 2.5):.3f}, {np.percentile(p_pred_1, 97.5):.3f}]")
print(f"Predicted probability for x=0: mean={p_pred_0.mean():.3f}, S.D={p_pred_0.std():.3f}, 95% CI=[{np.percentile(p_pred_0, 2.5):.3f}, {np.percentile(p_pred_0, 97.5):.3f}]")

diff = p_pred_1 - p_pred_0
print(f"Difference in predicted probabilities (x=1 - x=0): mean={diff.mean():.3f}, S.D={diff.std():.3f}, 95% CI=[{np.percentile(diff, 2.5):.3f}, {np.percentile(diff, 97.5):.3f}]")

Initializing NUTS using jitter+adapt_diag...

[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]

Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, x]

Output()

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 1 seconds.

            mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \
intercept -1.373  0.355  -2.035   -0.724      0.010    0.007    1262.0   
x          0.654  0.449  -0.210    1.507      0.013    0.009    1213.0   

           ess_tail  r_hat  
intercept    1417.0    1.0  
x            1348.0    1.0  
[-1.60103958 -1.87830767 -1.51071735 -1.2356058  -0.95301136]
[ 1.26476646  1.19283626  0.84703384  0.80792874 -0.34978388]
Posterior mean of intercept: -1.373
Posterior mean of x coefficient: 0.654
Predicted probability for x=1: mean=0.330, S.D=0.060, 95% CI=[0.219, 0.452]
Predicted probability for x=0: mean=0.208, S.D=0.057, 95% CI=[0.109, 0.329]
Difference in predicted probabilities (x=1 - x=0): mean=0.122, S.D=0.083, 95% CI=[-0.046, 0.286]

Comparing two proportions and running a logistic regression on a 0/1 group indicator are the same analysis, just expressed differently.

13.4 Latent-data formulation

We can interpret logistic regression directly - nonlinear model for probability of binary (yes or success), but also indirectly using unobserved or latent variables.

Below formulation $y_i = 1$ if $z_i^* > 0$ and $y_i = 0$ if $z_i^* < 0$, $z_i^* = X_i \beta + \epsilon_i$ (linear predictor) with independent errors $\epsilon_i$ that follow a logistic distribution. So we can think of logistic regression as a linear regression on an unobserved latent variable $z_i^*$, which is then transformed into a binary outcome $y_i$ using a threshold at 0. This is an alternative way to understand logistic regression and can be useful for certain applications, such as when we want to model the underlying process that generates the binary outcome.

the latent variable formulation reproduces the standard logistic regression probability model.

Example: someone in the lowest income group (income=1). Their linear predictor is -1.40 + 0.33*1 = -1.07, with logistic noise scattered around it (bell curve around -1.07). Most of its mass is below 0, but tail pokes above 0. $Pr (y_i = 1) = logit^{-1}(-1.07) = 0.25$. Center is below 0, so probability of y=1 is less than 0.5.

Practical payoff it recasts nonlinear probability as a linear model on hidden scale.

Interpretation of latent variables

E.g., in pre-election survey, $y_i = 1$ for Bush supporters and $y_i = 0$ for Clinton supporters. Unobserved continuous variable $z_i^*$ could be interpreted as a respondents "utility" or preference for Bush relative to Clinton: the sign tells us who they support and the magnitude tells us how strongly they support them.

Only the sign, not the magnitude can be determined from binary data. But we can ask respondents to rate their preference on a scale, which gives us more information about the magnitude of their preference.

Nonidentifiability of the latent scale parameter

Identification is where the data do not supply precise information about parameters in the model. Identification arises with scaling of the latent model in logistic regression.

Behind each yes/no is a continuous $z_i = X_i \beta + \epsilon_i$, we observe a 1 when $z_i$ crosses 0. Error has a logistic distribution, but the logistic curve is nearly identical to a normal curve with mean 0 and standard deviation 1.6. So we can write the model with normal noise: $z_i = X_i \beta + \epsilon_i$ with $\epsilon_i \sim N(0, 1.6^2)$. We fix 1.6 because its not identifiable (data cannot distinguish between parameters - multiple settings produce same predictions). The issue is that multiplying every parameter by the same constant does not change the predictions.

$z_i = -1.4 + 0.33*x_i + \epsilon_i$ \sigma = 1.6 is equivalent to $z_i = -14 + 3.3*x_i + \epsilon_i$ \sigma = 16, etc. Scaling doesnt matter because we are observing the binary outcome, not the latent variable. The fix is to set the scale $\sigma$ to a constant value of 1.6. Once fixed, the other parameters are identifiable.

13.5 Maximum likelihood and Bayesian inference for logistic regression

Classical linear regression can be done using purely algorithmic methods (as least squares) or maximum likelihood estimation (MLE). For logistic regression, usual place to start is maximum likelihood, followed by Bayesian inference and other regularization methods.

Maximum likelihood using iteratively weighted least squares

Binary logistic model predicts probability of y=1 and the probability is $logit^{-1}(X_i \beta)$ - number between 0 and 1 due to inverse logit (sigmoid) transformation. $X_i \beta$ is linear combination of predictors and coefficients.

Likelihood asks given some candidate coefficients $\beta$, how likely is the observed data? For each data point: if $y_i = 1$, probability is $logit^{-1}(X_i \beta)$ and if $y_i = 0$, probability is $1 - logit^{-1}(X_i \beta)$. Multiply all these per-point probabilities together to get the overall probability.

Finding best $\beta$. We want to find the $\beta$ that maximizes the likelihood of the observed data. Standard calculus move is take derivative, set to 0 and solve for $\beta$. In practice, take log of the likelihood, as it turns the product into a sum, easier to differentiate. The $\beta$ that maximizes the log-likelihood also maximises the likelihood, since log is monotonic.

Unlike linear, there is no closed-form solution for the MLE of $\beta$ in logistic regression. Instead, we use numerical methods, iteratively adjusting $\beta$ to find where derivative is 0. Sometimes no maximum exists - separation (predictor perfectly separates the data) or collinearity (predictor is a linear combination of other predictors).

Bayesian inference with a uniform prior distribution

The posterior (what we believe about $\beta$ after seeing the data) is proportional to the likelihood (how likely is the observed data given $\beta$) times the prior (what we believed about $\beta$ before seeing the data). With a uniform prior, each value of $\beta$ is equally likely a priori, multiplying likelihood by a constant does not change the shape of the posterior, so the posterior is proportional to the likelihood.

Bayesian inference gives us a full distribution over $\beta$, which allows us to quantify uncertainty in our estimates and make probabilistic statements about $\beta$.

However, a flat prior offers no guidance and so if data alone cannot pin down, we get no estimate. This is separation / non-identifiability. In cases where a predictor perfectly separates the data, even a mildly informative prior can help to regularize the estimates and provide more reasonable estimates of $\beta$.

Default prior in stan_glm

A logistic regression with several predictors:

$Pr(y_i = 1) = logit^{-1}(\alpha + \beta_1 x_{i1} + \beta_2 x_{i2} + ... + \beta_k x_{ik})$

a is intercept and bx is a slope coefficient for predictor x\dots

Each bx gets a normal prior with mean 0 and standard deviation 2.5. Mean 0 - before seeing data, we expect no relationship between predictors and outcome. SD 2.5 - makes prior scale-aware.

Intercept is different. Prior is not directly on intercept, but on the linear predictor when all predictors are at their mean. The 'centered intercept' gets a normal prior with mean 0 and SD 2.5. Intercept is value when all predictors are 0, the centered intercept is value when all predictors are at their mean. This makes the prior more interpretable.

Takeaway: zero-centered (shrink to zero) and scale-aware (unit independent), with intercept parameterized to something interpretable.

Bayesian inference with some prior information

If we know what to expect, we build that into the prior to get tighter more precise estimates. E.g., one predictor logistic regression, where from domain knowledge we expect slope b to be between 0 and 1. Soft vs hard constraints: hard would forbid values outside of 0 and 1, soft would allow them but with low probability - normal with mean 0.5 and SD 0.5. Normal distribution puts 68% of its mass between 0 and 1, and 95% between -0.5 and 1.5, so it allows for some uncertainty.

Comparing maximum likelihood and Bayesian inference using a simulation study

Consider examples, assuming true parameter a=-2 and b=0.8 gathered from experiment where x are from uniform distribution between 0 and 1.

import statsmodels.formula.api as smf

def bayes_sim(n, a=-2, b=0.8):
    rng = np.random.default_rng()
    x = rng.uniform(-1, 1, n)
    z = rng.logistic(loc=a + b * x, scale=1)   # rlogis(n, a + b*x, 1)
    y = (z > 0).astype(int)
    fake = pd.DataFrame({'x': x, 'y': y})

    # glm(y ~ x, family=binomial(link="logit"))
    glm_fit = smf.logit('y ~ x', data=fake).fit(disp=False)
    coef_df = pd.DataFrame({'coef': glm_fit.params.round(1), 'se': glm_fit.bse.round(1)})
    print("GLM (MLE):\n", coef_df, "\n")

    # stan_glm(y ~ x, family=binomial, prior=normal(0.5, 0.5))
    with pm.Model() as model:
        intercept = pm.Normal("intercept", mu=0.5, sigma=0.5)
        x_coef    = pm.Normal("x",         mu=0.5, sigma=0.5)
        p = pm.math.invlogit(intercept + x_coef * x)
        pm.Bernoulli("y", p=p, observed=y)
        trace = pm.sample(1000, tune=1000, target_accept=0.9, progressbar=False, random_seed=42)

    print("Bayesian (stan_glm):\n", pm.summary(trace)[['mean', 'sd']].round(1))
    return glm_fit, trace

print("Simulating with n=10:")
bayes_sim(10, a=-2, b=0.8)

print("\nSimulating with n=100:")
bayes_sim(100, a=-2, b=0.8)

print("\nSimulating with n=1000:")
bayes_sim(1000, a=-2, b=0.8)

Initializing NUTS using jitter+adapt_diag...

Simulating with n=10:
GLM (MLE):
            coef   se
Intercept  -2.4  1.3
x          -2.3  2.1 

Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, x]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 0 seconds.
Initializing NUTS using jitter+adapt_diag...

Bayesian (stan_glm):
            mean   sd
intercept  -0.3  0.4
x           0.2  0.5

Simulating with n=100:
GLM (MLE):
            coef   se
Intercept  -2.0  0.4
x           1.6  0.6 

Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, x]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 0 seconds.
Initializing NUTS using jitter+adapt_diag...

Bayesian (stan_glm):
            mean   sd
intercept  -1.3  0.2
x           0.8  0.3

Simulating with n=1000:
GLM (MLE):
            coef   se
Intercept  -2.1  0.1
x           0.9  0.2 

Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, x]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 1 seconds.

Bayesian (stan_glm):
            mean   sd
intercept  -2.0  0.1
x           0.8  0.2

(<statsmodels.discrete.discrete_model.BinaryResultsWrapper at 0x150e58260>,
 Inference data with groups:
 	> posterior
 	> sample_stats
 	> observed_data)

For small n GLM prediction is outside what we expect (and large SE), but Bayesian is closer due to priors. With n=100, GLM is closer. With n=1000, both are close and prior doesnt make a big difference.

13.6 Cross validation and log score for logistic regression

Model predicts probabilities (e.g., 70% yes), but if we only look at whether the rounded prediction is correct or not, we lose information. A confident (95% yes) prediction that is wrong is worse than a less confident (55% yes) prediction that is wrong and should be penalized more.

Brier score is to take the probability, subtract the actual outcome (0 or 1) and square it and take the average across all predictions. This is squared residuals, but on probability scale.

Bayesian R2 - in linear regression R2 uses models residual standard deviation to measure how much of the variation in the outcome is explained by the model. In logistic regression, we can use a similar approach, but its less interpretable because they dont add up to 1.

Log score and deviance - more natural fit is reward model by log of probability it assigned to what actually happened. Confident and right = high score, confident and wrong = low score. Deviance is this multiplied by -2.

Cross validation - if we judge on the same data we fit the model on, we will overfit and get an overly optimistic estimate of how well the model will perform on new data. So we need to evaluate the model by leaving out a point, fitting the model on the rest of the data and predicting the left out point. Repeat for all points and average the log score across all points. This is leave-one-out log score.

Understanding the log score for discrete predictions

Setup: fit logistic model to predict probabilities of new cases, then see what actually happened. For each new case take the log of the probability the model gave to the actual outcome - if 1, take log p, if 0, take log (1-p). Sum those up - this is the out-of-sample log score. Its always negative, closer to 0 is better. If model is confident and right, log score is close to 0. If model is confident and wrong, log score is very negative.

Calibration point 1 - no information baseline. Useless model says 50% for every case, this gives log(0.5) = -0.693 for each case. So a coin flip scores -0.693 * n. This is reference for knows nothing.

Calibration point 2 - slightly informed model. Model says 0.8 for half and 0.6 for other half, and crucially is calibrated, meaning that when it says 0.8, 80% of the time it is right and when it says 0.6, 60% of the time it is right. To get the expected score per case we average over what actually happened: for 0.8, 80% of the time it is right, 20% of time its wrong, contributing 0.8 * ln(0.8) + 0.2 * ln(0.2) = -0.5. For 0.6 it is 0.6 * ln(0.6) + 0.4 * ln(0.4) = -0.673. So the average score is 0.5 * -0.5 + 0.5 * -0.673 = -0.586. Going from coin flip (-0.693) to slightly informed (-0.586) is an improvement.

Log score for logistic regression

Null model gives 0.5 probability for each case, for 1179 survey respondents, log score is -0.693 * 1179 = -817.

Using the fact that 59.5% survey respondents supported Clinton and 40.5% supported Bush, we can assign 40.5% probability to Pr(y=1). This improves the log score to 477log(0.405) + 702log(0.595) = -787, which is better than the null model.

# Intercept-only model

# drop nans and reformat output and predictor
df = nes92[["income", "rvote"]].dropna()
y = df["rvote"].values.astype(int)

# fit logistic regression model
with pm.Model() as fit_1:
        
    # priors
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)
    
    # logistic likelihood
    p = pm.math.invlogit(intercept)
    pm.Bernoulli("rvote", p=p, observed=y)
    
    trace = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

summary = pm.summary(trace)
print(summary)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept]

Output()

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 0 seconds.

            mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \
intercept -0.386  0.059  -0.492   -0.275      0.002    0.001    1482.0   

           ess_tail  r_hat  
intercept    2103.0    1.0  

alpha_samples = trace.posterior['intercept'].values.flatten()

# intercept-only model: p is the same scalar for every observation,
# so we average expit(alpha) over all posterior draws
predp_bayes = np.mean(expit(alpha_samples))

logscore_bayes = np.sum(y * np.log(predp_bayes) + (1 - y) * np.log(1 - predp_bayes))
print(f"Log score: {logscore_bayes:.1f}")

Log score: -795.6

# With income as a predictor

# drop nans and reformat output and predictor
df = nes92[["income", "rvote"]].dropna()
x = df["income"].values.astype(int)
y = df["rvote"].values.astype(int)

# fit logistic regression model
with pm.Model() as fit_1:
        
    # priors
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)
    income = pm.Normal("income", mu=0, sigma=2.5)
    
    # logistic likelihood
    p = pm.math.invlogit(intercept + income * x)
    pm.Bernoulli("rvote", p=p, observed=y)
    
    trace = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

summary = pm.summary(trace)
print(summary)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, income]

Output()

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 1 seconds.

            mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \
intercept -1.402  0.195  -1.751   -1.015      0.007    0.005     784.0   
income     0.326  0.059   0.204    0.426      0.002    0.002     795.0   

           ess_tail  r_hat  
intercept     773.0   1.01  
income        803.0   1.01  

alpha_samples = trace.posterior['intercept'].values.flatten()
beta_samples  = trace.posterior['income'].values.flatten()

# income model: p varies per observation, average expit(a + b*x) over all posterior draws
# shape: (n_draws, n_obs) → mean over draws → (n_obs,)
predp_bayes = np.mean(
    [expit(a + b * x) for a, b in zip(alpha_samples, beta_samples)],
    axis=0
)

logscore_bayes = np.sum(y * np.log(predp_bayes) + (1 - y) * np.log(1 - predp_bayes))
print(f"Log score: {logscore_bayes:.1f}")

Log score: -778.5

Including income as a predictor improves the log score per observation from -796 to -778 = (796-778)/1179 = 0.015 per observation.

with fit_1:
    pm.compute_log_likelihood(trace)

loo_result = az.loo(trace)
print(loo_result)

Output()

Computed from 4000 posterior samples and 1179 observations log-likelihood matrix.

         Estimate       SE
elpd_loo  -780.49     8.53
p_loo        2.03        -
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)     1179  100.0%
   (0.70, 1]   (bad)         0    0.0%
   (1, Inf)   (very bad)    0    0.0%

13.7 Building a logistic regression model: wells in Bangladesh

Modelling decisions of households in Bangladesh about whether to change their source of drinking water.

Researchers measured all the wells and labelled with arsenic level as safe or unsafe. People with unsafe were encouraged to switch to a safe well. A few years later they surveyed the households again to see if they had switched. We carry out logistic regression to understand factors that predict whether a household switched to a safe well.

Outcome variable is 1 if household i switched to a safe well and 0 if they did not. We have the following predictors:

• distance to the nearest safe well
• arsenic level of the well they were using
• whether members of household were in community organisations
• education level of the household head

# load wells.csv
wells = pd.read_csv('../ros_data/wells.csv')

wells.head()

	switch	arsenic	dist	dist100	assoc	educ	educ4
0	1	2.36	16.826000	0.16826	0	0	0.0
1	1	0.71	47.321999	0.47322	0	0	0.0
2	0	2.07	20.966999	0.20967	0	10	2.5
3	1	1.15	21.486000	0.21486	0	12	3.0
4	1	1.10	40.874001	0.40874	1	14	3.5

# drop nans and reformat output and predictor
df = wells[["switch", "dist100"]].dropna()
x = df["dist100"].values.astype(float)
y = df["switch"].values.astype(int)

# fit logistic regression model
with pm.Model() as fit_1:
        
    # priors
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)
    dist = pm.Normal("dist", mu=0, sigma=2.5)
    
    # logistic likelihood
    p = pm.math.invlogit(intercept + dist * x)
    pm.Bernoulli("switch", p=p, observed=y)
    
    model_dist = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

summary = pm.summary(model_dist, hdi_prob=0.95)
print(summary)

with fit_1:
    pm.compute_log_likelihood(model_dist)

loo_result = az.loo(model_dist)
print(loo_result)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, dist]

Output()

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.

Output()

            mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  ess_bulk  \
intercept  0.606  0.058     0.490      0.719      0.002    0.001    1040.0   
dist      -0.622  0.095    -0.807     -0.434      0.003    0.002    1156.0   

           ess_tail  r_hat  
intercept    1383.0   1.01  
dist         1404.0   1.00  

Computed from 4000 posterior samples and 3020 observations log-likelihood matrix.

         Estimate       SE
elpd_loo -2040.01    10.41
p_loo        1.89        -
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)     3020  100.0%
   (0.70, 1]   (bad)         0    0.0%
   (1, Inf)   (very bad)    0    0.0%

Pr(switching) = logit⁻¹(0.6 − 0.62 × dist100)

Interpretation:

1. Constant when distance is 0 (next door), logit is 0.6, which corresponds to a probability of 0.65. 65% probability of switching when next door.
2. When a well is at the mean distance of 48m, distance is 0.48, so logit is 0.6 - 0.62*0.48 = 0.31. Slope at this point is -0.62e(0.31)/(1 + e(0.31))^2 = -0.14, which corresponds to a change in probability of -0.14. So a well at the mean distance of 48m has a 14% lower probability of switching than a well next door. Thus adding 100m to the distance to the nearest safe well decreases the probability of switching by 14%.
3. The divide by 4 rule gives us -0.62/4 = -0.155, which is consistent with our evaluation at the mean of the data. So we can say that adding 100m to the distance to the nearest safe well decreases the probability of switching by at most 15.5%.

The SD of 0.095 is small compared with the slope of -0.62, so we can be confident that the effect of distance is negative.

Graphing the fitted model

# jitter 0s toward 0 and 1s toward 1 so overlapping points are visible
switch_j = np.where(y == 0, np.random.uniform(0, 0.05, len(y)), np.random.uniform(0.95, 1, len(y)))

# x values for the smooth logistic curve
x_range = np.linspace(x.min(), x.max(), 200)
# print(x_range)
# posterior mean coefficients
a, b = trace.posterior['intercept'].values.mean(), trace.posterior['dist'].values.mean()

# raw data with jittered y
plt.scatter(x, switch_j, s=2, alpha=0.3)
# fitted logistic curve using posterior mean
plt.plot(x_range, expit(a + b * x_range))
plt.xlabel('dist100')
plt.ylabel('switch')
plt.show()

Adding a second input variable

Add arsenic level

# drop nans and reformat output and predictor
df = wells[["switch", "dist100", "arsenic"]].dropna()
dist = df["dist100"].values.astype(float)
arsenic = df["arsenic"].values.astype(float)
y = df["switch"].values.astype(int)

# fit logistic regression model
with pm.Model() as fit_2:
        
    # priors
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)
    dist_prior = pm.Normal("dist", mu=0, sigma=2.5)
    arsenic_prior = pm.Normal("arsenic", mu=0, sigma=2.5)
    
    # logistic likelihood
    p = pm.math.invlogit(intercept + dist_prior * dist + arsenic_prior * arsenic)
    pm.Bernoulli("switch", p=p, observed=y)
    
    model_arsenic = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

summary = pm.summary(model_arsenic, hdi_prob=0.95)
print(summary)

with fit_2:
    pm.compute_log_likelihood(model_arsenic)

loo_result = az.loo(model_arsenic)
print(loo_result)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, dist, arsenic]

Output()

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 3 seconds.

Output()

            mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  ess_bulk  \
intercept  0.006  0.078    -0.148      0.152      0.002    0.002    2002.0   
dist      -0.901  0.106    -1.107     -0.690      0.002    0.002    1974.0   
arsenic    0.460  0.041     0.381      0.544      0.001    0.001    2271.0   

           ess_tail  r_hat  
intercept    1897.0    1.0  
dist         1996.0    1.0  
arsenic      1913.0    1.0  

Computed from 4000 posterior samples and 3020 observations log-likelihood matrix.

         Estimate       SE
elpd_loo -1968.52    15.65
p_loo        3.29        -
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)     3020  100.0%
   (0.70, 1]   (bad)         0    0.0%
   (1, Inf)   (very bad)    0    0.0%

az.compare({"dist": model_dist, "dist+arsenic": model_arsenic})

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
dist+arsenic	0	-1968.521601	3.291553	0.000000	0.939109	15.649363	0.000000	False	log
dist	1	-2040.007534	1.885423	71.485933	0.060891	10.412671	12.122895	False	log

Arsenic below 0.5 is safe

To convert logit to probability = 1 / (1 + exp(-logit))

Pr(switching) = logit⁻¹(0.006 − 0.9 × dist100 + 0.46 × arsenic)

If dist is 0 and arsenic is 1, logit is 0.006 + 0.46 = 0.466, which corresponds to a probability of 0.61. So a well with arsenic level of 1 has a 61% probability of switching when next door.

If dist is 1 and arsenic is 0, logit is 0.006 - 0.9 = -0.894, which corresponds to a probability of 0.29. So a well with distance of 100m and arsenic level of 0 has a 29% probability of switching.

If dist is 1 and arsenic is 1, logit is 0.006 - 0.9 + 0.46 = -0.434, which corresponds to a probability of 0.39. So a well with distance of 100m and arsenic level of 1 has a 39% probability of switching.

If dist is 1 and arsenic is 2, logit is 0.006 - 0.9 + 0.92 = 0.026, which corresponds to a probability of 0.51. So a well with distance of 100m and arsenic level of 2 has a 51% probability of switching.

From formula: comparing two wells with same arsenic level, but one is 100m further away, the logit difference is -0.9, which corresponds to a probability difference of 0.39 - 0.29 = 0.10. So a well that is 100m further away has a 10% lower probability of switching. Similarly, comparing two wells with same distance, but one has arsenic level 1 higher, the logit difference is 0.46, which corresponds to a probability difference of 0.51 - 0.39 = 0.12. So a well that has arsenic level 1 higher has a 12% higher probability of switching.

For quick interpretation, we can use the divide by 4 rule. For distance, -0.9/4 = -0.225, which means that a well that is 100m further away has at most a 22.5% lower probability of switching. For arsenic level, 0.46/4 = 0.115, which means that a well that has arsenic level 1 higher has at most a 11.5% higher probability of switching.

Comparing the coefficient estimates when adding a predictor

Distance coefficient changes from -0.62 to -0.9 when we add arsenic level as a predictor. Because wells that are far from the nearest safe well are likely to be unsafe.

#####

13.8 Bibliographic note

13.9 Exercises