14. Working with logistic regression

%run functions.ipynb

14.1 Graphing logistic regression and binary data

It can be useful to plot the discrete data points along with the logistic regression curve to visualize how well the model fits the data. However, plotting the discrete data points is not as visually informative as when plotting continuous data. We can use jittering to spread out the discrete data points.

14.2 Logistic regression with interactions

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

wells.head()

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

# fit logistic regression model
with pm.Model() as fit_1:
        
    # priors
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)
    dist_coef = pm.Normal("dist", mu=0, sigma=2.5)
    arsenic_coef = pm.Normal("arsenic_coef", mu=0, sigma=2.5)
    interaction_coef = pm.Normal("interaction_coef", mu=0, sigma=2.5)
    # logistic likelihood
    p = pm.math.invlogit(intercept + dist_coef * dist100 + arsenic_coef * arsenic + interaction_coef * dist100_arsenic)
    pm.Bernoulli("switch", p=p, observed=outcome)
    
    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)

To understand the model, we can:

1. Evaluate predictions and interactions at the mean of the data (0.48 for dist100 and 1.66 for arsenic). That is 48 meters to nearest safe well and a mean of 1.66 amoung unsafe wells.
2. Dividing by 4 to get approximate differences on probability scale

switching probability = $logit^{-1}(-0.14 + (-0.594) * dist100 + 0.553 * arsenic + (-0.172) * dist100 * arsenic)$

We interpret each in turn:

1. Constant term: $logit^{-1}(-0.14) = (1/(1 + exp(0.14))) = 0.47 - probability of switching when distance to nearest safe well and arsenic level in current well are at 0. This is impossible, since all arsenic levels are above 0. We dont interpret this term in isolation. Instead at the average distance to nearest safe well and average arsenic level, the probability of switching is $logit^{-1}(-0.14 + (-0.594) * 0.48 + 0.553 * 1.66 + (-0.172) * 0.48 * 1.66) = logit^{-1}(0.355) = 0.59$. This means that at the average distance to nearest safe well and average arsenic level, the probability of switching is 59%.
2. Coefficient for dist100: -0.594. This corresponds to comparing two wells that differ by 1 unit in distance when arsenic levels are 0. Again, we should not interpret this term in isolation. Instead we consider the effect of distance at the average arsenic level. At the average arsenic level, the effect of distance is -0.594 + (-0.172) * 1.66 = -0.88 on logit scale. Dividing by 4 gives us an approximate difference of -0.22 on the probability scale. This means that at the average arsenic level, a 1 unit increase in distance to nearest safe well is associated with a 22% decrease in the probability of switching.
3. Coefficient for arsenic: 0.553. This corresponds to comparing two wells that differ by 1 unit in arsenic level when distance to nearest safe well is 0. Again, we should not interpret this term in isolation. Instead we consider the effect of arsenic at the average distance to nearest safe well. At the average distance to nearest safe well, the effect of arsenic is 0.553 + (-0.172) * 0.48 = 0.47 on logit scale. Dividing by 4 gives us an approximate difference of 0.12 on the probability scale. This means that at the average distance to nearest safe well, a 1 unit increase in arsenic level is associated with a 12% increase in the probability of switching.
4. Coefficient for interaction: -0.172. Can be interpreted in two ways. For each additional unit of arsenic, the -0.172 is added to the coefficient for distance. We have already seen that the coefficient for distance is -0.88 at the average arsenic level, and so we can interpret this to say that the importance of distance as a predictor increases for households with higher arsenic. Alternatively, for each additional unit of distance, the -0.172 is added to the coefficient for arsenic. We have already seen that the coefficient for arsenic is 0.47 at the average distance to nearest safe well, and so we can interpret this to say that the importance of arsenic as a predictor decreases for households with higher distance to nearest safe well.

Centering the inputs

Before fitting interactions, makes sense to center predictors, this makes them more interpretable. We dont standardise them, since we want to interpret the coefficients in terms of the original units. We can center the predictors by subtracting the mean from each value.

wells["dist100_centered"] = wells["dist100"] - wells["dist100"].mean()
wells["arsenic_centered"] = wells["arsenic"] - wells["arsenic"].mean()
wells["dist100_arsenic_centered"] = wells["dist100_centered"] * wells["arsenic_centered"]

Refitting the model with centered inputs

We model using centered inputs, not predictors - interaction is not centered, but the interaction between the centered predictors.

dist100_centered = wells["dist100_centered"].values.astype(float)
arsenic_centered = wells["arsenic_centered"].values.astype(float)
dist100_arsenic_centered = wells["dist100_arsenic_centered"].values.astype(float)

# fit logistic regression model
with pm.Model() as fit_2:
        
    # priors
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)
    dist_coef = pm.Normal("dist", mu=0, sigma=2.5)
    arsenic_coef = pm.Normal("arsenic_coef", mu=0, sigma=2.5)
    interaction_coef = pm.Normal("interaction_coef", mu=0, sigma=2.5)
    # logistic likelihood
    p = pm.math.invlogit(intercept + dist_coef * dist100_centered + arsenic_coef * arsenic_centered + interaction_coef * dist100_arsenic_centered)
    pm.Bernoulli("switch", p=p, observed=outcome)
    
    model_interact = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

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

with fit_2:
    pm.compute_log_likelihood(model_interact)

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

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

Formula for the model is:

probability of switching = logit^{-1}(intercept + dist100_{centered} * dist100_{coefficient} + arsenic_{centered} * arsenic_{coefficient} + dist100_{arsenic}{centered} * interaction{coefficient})

probability of switching = logit^{-1}(0.352 + (-0.875) * dist100_{centered} + (0.471) * arsenic_{centered} + (-0.182) * dist100_{arsenic}_{centered})

Interpretation:

1. Constant term: logit^{-1}(0.352) = 0.59 - probability of switching when distance to nearest safe well and arsenic level in current well are at their mean values.
2. Coefficient for dist100_{centered}: -0.875. Coefficient on logit scale for distance if arsenic is at its mean. Quick interpretation divide by 4 to get approximate difference on probability scale. -0.875 / 4 = -0.22. This means that at the average arsenic level, a 1 unit increase in distance to nearest safe well is associated with a 22% decrease in the probability of switching.
3. Coefficient for arsenic_{centered}: 0.471. Coefficient on logit scale for arsenic if distance is at its mean. Quick interpretation divide by 4 to get approximate difference on probability scale. 0.471 / 4 = 0.12. This means that at the average distance to nearest safe well, a 1 unit increase in arsenic level is associated with a 12% increase in the probability of switching.
4. Coefficient for interaction unchanged by centering and has same interpretation as before.

Statistical significance of the interaction

Interaction term interaction_coef -0.182 standard error 0.102, not quite 2 standard errors away from 0, so not statistically significant at the 5% level. However, we this does not mean that the interaction is not important and there is no effect. Rather accept uncertainty in context that there is a plausable effect. LOO log is -1968 and elpd diff is very small, so the interaction does not improve the model fit.

from scipy.special import expit

post = model_dist.posterior.mean(dim=["chain", "draw"])
b0, b1, b2, b3 = float(post["intercept"]), float(post["dist"]), float(post["arsenic_coef"]), float(post["interaction_coef"])

rng = np.random.default_rng(42)
switch_jitter = wells["switch"] + rng.uniform(-0.05, 0.05, size=len(wells))

x = np.linspace(0, wells["dist"].max(), 300)
x100 = x / 100

fig, ax = plt.subplots()
ax.scatter(wells["dist"], switch_jitter, alpha=0.1, s=5, color="gray")
ax.plot(x, expit(b0 + b1*x100 + b2*0.5 + b3*0.5*x100), label="arsenic = 0.5")
ax.plot(x, expit(b0 + b1*x100 + b2*1.0 + b3*1.0*x100), label="arsenic = 1.0")
ax.set_xlim(0, wells["dist"].max())
ax.set_xlabel("dist")
ax.set_ylabel("P(switch)")
ax.legend()
plt.show()

Adding social predictors

Are well users more likely to switch if they have more education or community connections?

educ = years of education = educ4 = education/4 assoc = 1 is member of a community association, 0 otherwise

dist100 = df["dist100"].values.astype(float)
outcome = df["switch"].values.astype(int)
arsenic = df["arsenic"].values.astype(float)
educ4 = wells["educ4"].values.astype(float)
assoc = wells["assoc"].values.astype(int)

# fit logistic regression model
with pm.Model() as fit_1:
        
    # priors
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)
    dist_coef = pm.Normal("dist", mu=0, sigma=2.5)
    arsenic_coef = pm.Normal("arsenic_coef", mu=0, sigma=2.5)
    educ4_coef = pm.Normal("educ4_coef", mu=0, sigma=2.5)
    assoc_coef = pm.Normal("assoc_coef", mu=0, sigma=2.5)

    # logistic likelihood
    p = pm.math.invlogit(intercept + dist_coef * dist100 + arsenic_coef * arsenic + educ4_coef * educ4 + assoc_coef * assoc)
    pm.Bernoulli("switch", p=p, observed=outcome)
    
    model_2 = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

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

with fit_1:
    pm.compute_log_likelihood(model_2)

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

Respondents with higher education are more likely to switch logit-1(0.170), divide by 4 to get approximate difference on probability scale, 0.170 / 4 = 0.043, so a 1 unit increase in education (4 years) is associated with a 4.3% increase in the probability of switching.

Assoc has a negative coefficient, but not estimated precisely, so we remove from the model.

dist100 = df["dist100"].values.astype(float)
outcome = df["switch"].values.astype(int)
arsenic = df["arsenic"].values.astype(float)
educ4 = wells["educ4"].values.astype(float)

# fit logistic regression model
with pm.Model() as fit_1:
        
    # priors
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)
    dist_coef = pm.Normal("dist", mu=0, sigma=2.5)
    arsenic_coef = pm.Normal("arsenic_coef", mu=0, sigma=2.5)
    educ4_coef = pm.Normal("educ4_coef", mu=0, sigma=2.5)

    # logistic likelihood
    p = pm.math.invlogit(intercept + dist_coef * dist100 + arsenic_coef * arsenic + educ4_coef * educ4)
    pm.Bernoulli("switch", p=p, observed=outcome)
    
    model_3 = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

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

with fit_1:
    pm.compute_log_likelihood(model_3)

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

az.compare({"model_dist": model_dist, "model_3": model_3})

Adding education improves predictive log score, but with large uncertainty

Adding further interactions

When inputs have large main effects, general practice is to include interactions with these inputs.

# centering
wells["dist100_centered"] = wells["dist100"] - wells["dist100"].mean()
wells["arsenic_centered"] = wells["arsenic"] - wells["arsenic"].mean()
wells["educ4_centered"] = wells["educ4"] - wells["educ4"].mean()

#interactions
wells["arsenic_educ4_centered"] = wells["arsenic_centered"] * wells["educ4_centered"]
wells["dist100_educ4_centered"] = wells["dist100_centered"] * wells["educ4_centered"]

# reformat output and predictor
dist100_centered = wells["dist100_centered"].values.astype(float)
arsenic_centered = wells["arsenic_centered"].values.astype(float)
educ4_centered = wells["educ4_centered"].values.astype(float)
educ4 = wells["educ4"].values.astype(float)
arsenic_educ4_centered = wells["arsenic_educ4_centered"].values.astype(float)
dist100_educ4_centered = wells["dist100_educ4_centered"].values.astype(float)

# fit logistic regression model
with pm.Model() as fit_1:
        
    # priors
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)
    dist_coef = pm.Normal("dist", mu=0, sigma=2.5)
    arsenic_coef = pm.Normal("arsenic_coef", mu=0, sigma=2.5)
    educ4_coef = pm.Normal("educ4_coef", mu=0, sigma=2.5)
    arsenic_educ4_coef = pm.Normal("arsenic_educ4_coef", mu=0, sigma=2.5)
    dist100_educ4_coef = pm.Normal("dist100_educ4_coef", mu=0, sigma=2.5)

    # logistic likelihood
    p = pm.math.invlogit(intercept + dist_coef * dist100_centered + arsenic_coef * arsenic_centered + educ4_coef * educ4 + arsenic_educ4_coef * arsenic_educ4_centered + dist100_educ4_coef * dist100_educ4_centered)
    pm.Bernoulli("switch", p=p, observed=outcome)
    
    model_4 = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

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

with fit_1:
    pm.compute_log_likelihood(model_4)

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

Formula for the model is:

probability of switching = logit^{-1}(intercept + dist100_centered * dist100_coef + arsenic_centered * arsenic_coef + educ4 * educ4_coef + arsenic_educ4_centered * arsenic_educ4_coef + dist100_educ4_centered * dist100_educ4_coef)

probability of switching = logit^{-1}(0.119 + (-0.918) * dist100_centered + (0.492) * arsenic_centered + (0.188) * educ4 + (0.077) * arsenic_educ4_centered + (0.330) * dist100_educ4_centered)

1. Interaction between distance and education: 0.330. A difference of 4 years of education is associated with a 0.330 increase in the coefficient for distance. dist100 had a negative coefficient, thus positive changes in education reduce distances negative effect on switching. Makes sense, as people with more education might have other resources, so walking the extra distance to a safe well is less of a burden.
2. Interaction between arsenic and education: 0.077. A difference of 4 years of education is associated with a 0.077 increase in the coefficient for arsenic. Arsenic had a positive coefficient, thus positive changes in education increase arsenics positive effect on switching. Makes sense, as people with more education might be more aware of the dangers of arsenic, so they are more likely to switch to a safe well when their current well has high arsenic levels.

az.compare({"model_4": model_4, "model_3": model_3})

Stanadrdising the inputs

We should seriously consider standardising the inputs as a default in models with interactions.

14.3 Predictive simulation

# fit logistic regression model
with pm.Model() as fit_1:
        
    # priors
    intercept = pm.Normal("intercept", mu=0, sigma=2.5)
    dist_coef = pm.Normal("dist", mu=0, sigma=2.5)

    # logistic likelihood
    p = pm.math.invlogit(intercept + dist_coef * dist100)
    pm.Bernoulli("switch", p=p, observed=outcome)
    
    simple_model = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

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

with fit_1:
    pm.compute_log_likelihood(simple_model)

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

b0_sims = simple_model.posterior["intercept"].values.flatten()
b1_sims = simple_model.posterior["dist"].values.flatten()

rng = np.random.default_rng(42)
switch_jitter = wells["switch"] + rng.uniform(-0.05, 0.05, size=len(wells))
x = np.linspace(0, wells["dist100"].max(), 200)
# print(x)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# beta_0 vs beta_1
ax1.scatter(b0_sims, b1_sims, alpha=0.1, s=5)
ax1.set_xlabel(r"$\beta_0$")
ax1.set_ylabel(r"$\beta_1$")

# data + uncertainty curves
ax2.scatter(wells["dist100"], switch_jitter, alpha=0.1, s=5, color="gray")
for i in rng.choice(len(b0_sims), 20, replace=False):
    ax2.plot(x, expit(b0_sims[i] + b1_sims[i] * x), color="gray", lw=0.5)
ax2.plot(x, expit(b0_sims.mean() + b1_sims.mean() * x), color="black")
ax2.set_xlabel("dist100")
ax2.set_ylabel("P(switch)")

plt.tight_layout()
plt.show()

Predictive simulation using the binomial distribution

Now predict for a new home.

for x in np.arange(0, 1.1, 0.1):
    print("Home", x * 100, "m:", expit(b0_sims.mean() + b1_sims.mean() * x))

# X_new: design matrix [intercept, dist100]
X_new = np.column_stack([np.ones(len(wells)), wells["dist100"].values])
# print(X_new)
sims_mat = np.column_stack([b0_sims, b1_sims])  # (n_sims, 2)

n_new = X_new.shape[0]
# print(n_new)  # number of new observations
n_sims = len(b0_sims)
rng = np.random.default_rng(42)

y_new = np.zeros((n_sims, n_new), dtype=int)
for s in range(n_sims):
    p_new = expit(X_new @ sims_mat[s])   # X_new %*% sims[s,]
    y_new[s] = rng.binomial(1, p_new)    # rbinom(n_new, 1, p_new)

print(y_new.shape)          # (n_sims, n_new)
print(y_new.mean(axis=0)[:5])  # posterior predicted P(switch) for first 5 rows

Predictive simulation using the latent logistic distribution

Here we can add noise to the linear predictor and then transform to the probability scale.

# latent variable formulation: y=1 if z>0, where z = X@beta + epsilon, epsilon ~ Logistic(0,1)
y_new2 = np.zeros((n_sims, n_new), dtype=int)
for s in range(n_sims):
    u = rng.uniform(0, 1, size=n_new)
    epsilon_new = np.log(u / (1 - u))           # logit(runif(n_new, 0, 1))
    z_new = X_new @ sims_mat[s] + epsilon_new   # X_new %*% t(sims[s,]) + epsilon_new
    y_new2[s] = (z_new > 0).astype(int)         # ifelse(z_new > 0, 1, 0)

print(y_new2.shape)
print(y_new2.mean(axis=0)[:5])  # should match y_new closely

14.4 Average predictive comparisons on the probability scale

Logistic regression is non-linear on the probability scale and linear on the logit scale. The parameters live on the linear scale, but nonlinear in the relation of inputs to outputs. Logistic regression coefficients cannot be interpreted in original scale of data.

Where we have many predictors a simple graph of x vs y is not possible. Instead we use summary comparable to linear regression coefficients, which gives the expected, or average, difference in Pr(y=1) for a 1 unit change in x.

Focus on one input at a time, the one we care about is 'u' all others are 'v', so full input vector is (x = u, v). Pick two values of u (e.g., high and low) and ask what changes moving from one to the other while holding v constant. In well example, moving from household 100 meters from safe well versus 0 meters, while holding all other inputs constant.

The predictive difference is model probability at high u minus model probability at low u. Predictive difference = E(y | u_high, v) - E(y | u_low, v). E(y | ...) means expected output given the inputs - for binary that is Py(y=1 | ...). We can compute the predictive difference for each household in the data, because they are all different, and then average over all households to get the average predictive difference. That says, for the example, moving from high to low is associated with an average change of X percentage points change in the outcome.

We can also do this as a rate - change in probability divided by change in input. This is the average predictive comparison, which is the average change in probability for a 1 unit change in input. Predictive comparison = [E(y | u_high, v) - E(y | u_low, v)] / (u_high - u_low). This is the average change in probability for a 1 unit change in input - slope. However, depends on u, v, and \theta, and different at every point in the input space. So, to report a single number, either fix for representative values of x and \theta or average over the data.

Problems with evaluating predictive comparisons at a central value

Often a rough approach is to evaluate E(y|u,v,\theta) at a central value of v, such as the mean or median. However, this is problematic when inputs are highly varied with no single central value, or binary or bimodal inputs. This approach cannot be automated as it requires choice about how to setup the range for each input.

The main issue is that the average of v is not the average of the effects of v. This is because averages can fall in regions of the input space where there are no data, or can undershoot. The only reliable way is to compute predictive comparison for each real data point and then average over all data points. This is the average predictive comparison.

14.5 Residuals for discrete-data regression

14.6 Identification and separation

14.7 Bibliographic note

14.8 Exercises