Resampling and Logistic Regression#

Click here to run this notebook on Colab

Different ways of computing sampling distributions – and the statistics derived from them, like standard errors and confidence intervals – yield different results. None of them are right or wrong; rather, they are based on different modeling assumptions. In practice, the differences are often small compared to other sources of error, so we are free to use whichever is convenient.

This article is prompted by a recent question on Reddit:

I am trying to do parametric bootstrap for logistic regression but I don’t know what the random error terms should be. I know how to do it with linear regression since the error terms follow a normal distribution. Really appreciate any pointers to resources about this topic.

One of the responses recommends this paper: “An Application of Bootstrapping in Logistic Regression Model”, by Isaac Akpor Adjei and Rezaul Karim.

The paper suggests two ways to compute the sampling distribution of the parameters of a logistic regression by bootstrap sampling. They characterize one as parametric and the other as non-parametric, and advise:

The non-parametric bootstrap [relies] on weaker assumptions (or no assumptions at all) about the underlying distribution of the population. Statistical practitioners should use non-parametric procedures only in so far as the assumptions about the underlying distribution are seriously doubtful in their validity. … However, when assumptions are not violated, non-parametric procedures will usually have greater variance (in point estimation), less power (in hypothesis testing), wider intervals (in confidence interval estimation), lower probability of correct selection (in ranking and selection) and higher risk (in decision theory) when compared to a corresponding parametric procedure (Efron and Tibshirani, 1994 [1] ).

The premise of this advice is that the parametric and non-parametric methods are answering the same question in different ways, and that one might be preferred over the other because the results are preferable in one way or another. It is also based on the assumption that the sampling distribution of the parameters – and the statistics derived from it – are uniquely and objectively defined.

I disagree with these premises. The parametric and non-parametric methods they present are based on different modeling decisions, so they compute answers to different questions. And the models they are based on are only two of many possible models.

To explain what I mean, I will implement the methods they propose and explain the assumptions each is based on. Then I will propose a third method that is a hybrid of the two. I’ll show that all three methods yield different results, and suggest criteria for when one might be preferred over the others.

The Data#

As an example, I’ll use data from the General Social Survey (GSS). I’ll download an HDF file from the Elements of Data Science repository, which contains a subset of the GSS data that has been resampled to correct for stratified sampling.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

from os.path import basename, exists

def download(url):
    filename = basename(url)
    if not exists(filename):
        from urllib.request import urlretrieve
        local, _ = urlretrieve(url, filename)
        print('Downloaded ' + local)
download('' +
gss = pd.read_hdf('gss_eda.hdf', 'gss')
(64814, 8)

It includes 64,814 respondents and 8 variables for each respondent.

0 1972 1 23.0 16.0 2 1.0 NaN 18951.0
1 1972 2 70.0 10.0 1 1.0 NaN 24366.0
2 1972 3 48.0 12.0 2 1.0 NaN 24366.0
3 1972 4 27.0 17.0 2 1.0 NaN 30458.0
4 1972 5 61.0 12.0 2 1.0 NaN 50763.0

To demonstrate logistic regression, I’ll use on one of the questions in the General Social Survey, which asks “Do you think the use of marijuana should be made legal or not?” The responses are in a column called GRASS; here are the values.

2.0    25662
1.0    11884
Name: GRASS, dtype: int64

I’ll use StatsModels, which provides a function that does logistic regression. First we have to recode the dependent variable so 1 means “yes” and 0 means “no”. We can do that by replacing 2 with 0.

gss['GRASS'].replace(2, 0, inplace=True)

To model quadratic relationships, I’ll add columns that contain the values of AGE and EDUC squared.

gss['AGE2'] = gss['AGE']**2
gss['EDUC2'] = gss['EDUC']**2

And I’ll drop the rows that have missing values for the variables we’ll need.

data = gss.dropna(subset=['AGE', 'EDUC', 'SEX', 'GUNLAW', 'GRASS'])

Logistic Regression#

Here are the results of logistic regression with these variables.

import statsmodels.formula.api as smf

formula = 'GRASS ~ AGE + AGE2 + EDUC + EDUC2 + C(SEX)'
result_hat = smf.logit(formula, data=data).fit()
Optimization terminated successfully.
         Current function value: 0.593464
         Iterations 6
Logit Regression Results
Dep. Variable: GRASS No. Observations: 20475
Model: Logit Df Residuals: 20469
Method: MLE Df Model: 5
Date: Tue, 17 Jan 2023 Pseudo R-squ.: 0.05003
Time: 08:58:57 Log-Likelihood: -12151.
converged: True LL-Null: -12791.
Covariance Type: nonrobust LLR p-value: 1.555e-274
coef std err z P>|z| [0.025 0.975]
Intercept -1.6788 0.240 -6.988 0.000 -2.150 -1.208
C(SEX)[T.2] -0.3849 0.031 -12.394 0.000 -0.446 -0.324
AGE -0.0278 0.005 -5.399 0.000 -0.038 -0.018
AGE2 0.0001 5.28e-05 2.190 0.029 1.21e-05 0.000
EDUC 0.2000 0.031 6.412 0.000 0.139 0.261
EDUC2 -0.0029 0.001 -2.450 0.014 -0.005 -0.001

To get a sense of what the results look like, we can plot the predicted probability of saying “yes” as a function of age, for male and female respondents.

df = pd.DataFrame()
df['AGE'] = np.linspace(18, 89)
df['EDUC'] = 16

df['AGE2'] = df['AGE']**2
df['EDUC2'] = df['EDUC']**2

df['SEX'] = 1
pred1 = result_hat.predict(df)
pred1.index = df['AGE']

df['SEX'] = 2
pred2 = result_hat.predict(df)
pred2.index = df['AGE']

pred1.plot(label='Male', alpha=0.6)
pred2.plot(label='Female', alpha=0.6)

plt.title('Support for legal marijuana')

Men are more likely to support legalization than women, and younger people more likely than older people.

Boostrap sampling#

Now let’s implement the method characterized as non-parametric, which is based on bootstrap resampling. The following function samples the rows of data with replacement and runs the regression model on the resampled data.

options = dict(disp=False, start_params=result_hat.params) 

def bootstrap(i):
    bootstrapped = data.sample(n=len(data), replace=True)
    results = smf.logit(formula, data=bootstrapped).fit(**options)
    return results.params

The options sent to fit make it run faster, but don’t affect the results by much.

Each time we run this process, the result represents a single draw from the sampling distribution of the parameters.

%timeit bootstrap(0)
82.2 ms ± 1.89 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

By running it many times, we generate a sample from the sampling distribution. pqdm runs the sampling process in multiple threads.

    from pqdm.threads import pqdm
except ImportError:
    !pip install pqdm
    from pqdm.threads import pqdm
estimates = pqdm(range(101), bootstrap, n_jobs=4)
sampling_dist = pd.DataFrame(estimates)
Intercept C(SEX)[T.2] AGE AGE2 EDUC EDUC2
0 -1.747448 -0.413202 -0.027473 0.000114 0.213299 -0.003302
1 -2.055303 -0.411384 -0.029943 0.000129 0.276153 -0.005563
2 -2.050773 -0.414960 -0.023190 0.000094 0.242778 -0.004530
3 -2.031659 -0.392391 -0.029538 0.000147 0.251942 -0.004616
4 -1.847295 -0.363829 -0.031532 0.000156 0.238870 -0.004339

Here’s what the sampling distribution looks like for one of the parameters.

ci90 = sampling_dist['C(SEX)[T.2]'].quantile([0.05, 0.95])
[plt.axvline(x, ls=':') for x in ci90]
plt.title('Sampling distribution of a parameter');

The mean of the sampling distribution should be close to the parameters we estimated with the original dataset, and it is.

pd.DataFrame({"Sampling mean": sampling_dist.mean(), 
              "Estimates": result_hat.params})
Sampling mean Estimates
Intercept -1.681880 -1.678838
C(SEX)[T.2] -0.379621 -0.384919
AGE -0.027804 -0.027825
AGE2 0.000115 0.000116
EDUC 0.201137 0.200002
EDUC2 -0.002912 -0.002855

The standard deviations of the sampling distributions should be close to the standard errors computed by StatsModels, and they are. Most of them are close enough that the difference probably doesn’t matter in practice.

def standard_errors(sampling_dist, result_hat):
    df = pd.DataFrame({"Sampling std": sampling_dist.std(), 
                       "Standard error": result_hat.bse})
    num, den = df.values.T
    df['Percent diff'] = (num / den - 1) * 100
    return df
standard_errors(sampling_dist, result_hat)
Sampling std Standard error Percent diff
Intercept 0.248814 0.240243 3.567779
C(SEX)[T.2] 0.027137 0.031057 -12.622418
AGE 0.004944 0.005153 -4.066840
AGE2 0.000051 0.000053 -4.198399
EDUC 0.031756 0.031193 1.804886
EDUC2 0.001197 0.001166 2.687507

Of course, there is nothing magic about the standard errors computed by StatsModels; they are approximations based on a model of the sampling process, just like the results from resampling.

The Resampling Framework#

The purpose of the sampling distribution is to quantify the variability in an estimate due to random sampling. It is the answer to the question, “If we ran this sampling process many times, how much would we expect the results to vary?”

To answer that question, we need a model of the sampling process. We use the model to simulate the sampling process and generate counterfactual datasets that represent other samples the process could have generated.

The following figure represents this framework.

In bootstrap resampling, we treat the rows of the original dataset as if they are the population and draw a random sample from it. We use this simulated data to compute the sample statistic; in this example, it’s the parameters of the logistic regression model. We collect the results to form a sample from the sampling distribution, which we can use to compute confidence intervals and standard errors.

This way of simulating the sampling process asks, in effect, what would have happened if we had surveyed different people. But that’s not the only possible model of the sampling process. An alternative is to ask what would happen if we surveyed the same people, but they gave different answers. That’s the model that underlies the parametric method.

Parametric bootstrap#

So, how do we simulate a sampling process where we survey the same people, but they give different answers? One way is to assume, for the purposes of the simulation, that the parameters we estimated from the original data are correct. If so, we can use the regression model to compute the predicted probability for each respondent, like this:

pi_hat = result_hat.predict(data)
count    20475.000000
mean         0.317167
std          0.113654
min          0.025989
25%          0.234792
50%          0.317927
75%          0.400791
max          0.630366
dtype: float64

Now we can use these probabilities to generate a biased coin toss for each respondent.

from scipy.stats import bernoulli

simulated = bernoulli.rvs(pi_hat.values)

Then we can run the regression model with these simulated values as the dependent variable. The following function shows how.

def bootstrap2(i):
    flipped = data.assign(GRASS=bernoulli.rvs(pi_hat.values))
    results = smf.logit(formula, data=flipped).fit(**options)
    return results.params

Again, the result from a single simulation is a random value from the sampling distribution.

Intercept     -1.805556
C(SEX)[T.2]   -0.352934
AGE           -0.023101
AGE2           0.000058
EDUC           0.196758
EDUC2         -0.002599
dtype: float64

If we run it many times, we get a sample from the sampling distribution.

from pqdm.threads import pqdm

estimates = pqdm(range(101), bootstrap2, n_jobs=4)
sampling_dist2 = pd.DataFrame(estimates)
Intercept C(SEX)[T.2] AGE AGE2 EDUC EDUC2
0 -1.804015 -0.384881 -0.014708 -0.000017 0.180796 -0.002302
1 -1.432664 -0.422171 -0.023207 0.000044 0.166255 -0.001940
2 -1.415658 -0.406837 -0.033815 0.000168 0.180158 -0.001914
3 -1.840253 -0.394937 -0.028435 0.000126 0.221295 -0.003461
4 -1.642430 -0.375719 -0.034946 0.000177 0.216985 -0.003430

Here’s the sampling distribution for a single parameter.

ci90 = sampling_dist2['C(SEX)[T.2]'].quantile([0.05, 0.95])
[plt.axvline(x, ls=':') for x in ci90]
plt.title('Sampling distribution of a parameter');

Again, the standard deviations should be close to the standard errors, and they are.

standard_errors(sampling_dist2, result_hat)
Sampling std Standard error Percent diff
Intercept 0.199178 0.240243 -17.092962
C(SEX)[T.2] 0.030494 0.031057 -1.812430
AGE 0.005085 0.005153 -1.327560
AGE2 0.000054 0.000053 2.118557
EDUC 0.024899 0.031193 -20.178173
EDUC2 0.000935 0.001166 -19.819110

So we have two models of the sampling process that yield different sampling distributions. But these are not the only models.

If the first method asks, “What if we surveyed different people?” and the second asks “What if we surveyed the same people and they gave different answers?”, let’s consider a third method that asks “¿Por qué no los dos?”.

The Hybrid Model#

The following function uses bootstrap sampling to simulate surveying different people; then it uses the parametric method to simulate their responses.

def bootstrap3(i):
    bootstrapped = data.sample(n=len(data), replace=True)
    pi_hat = result_hat.predict(bootstrapped)
    flipped = bootstrapped.assign(GRASS=bernoulli.rvs(pi_hat.values))
    results = smf.logit(formula, data=flipped).fit(**options)
    return results.params
Intercept     -1.655314
C(SEX)[T.2]   -0.386624
AGE           -0.024265
AGE2           0.000083
EDUC           0.182822
EDUC2         -0.002334
dtype: float64
from pqdm.threads import pqdm

estimates = pqdm(range(101), bootstrap3, n_jobs=4)
sampling_dist3 = pd.DataFrame(estimates)
Intercept C(SEX)[T.2] AGE AGE2 EDUC EDUC2
0 -2.019658 -0.385612 -0.025974 0.000108 0.249717 -0.005012
1 -1.648528 -0.351063 -0.025192 0.000100 0.182666 -0.002312
2 -1.780499 -0.445770 -0.032862 0.000173 0.242922 -0.004569
3 -1.534347 -0.464670 -0.029799 0.000125 0.189228 -0.002218
4 -1.336757 -0.404636 -0.030217 0.000137 0.160754 -0.001453

Here’s the sampling distribution for one of the parameters.

ci90 = sampling_dist3['C(SEX)[T.2]'].quantile([0.05, 0.95])
[plt.axvline(x, ls=':') for x in ci90]
plt.title('Sampling distribution of a parameter');

Again, the standard errors are comparable to the ones computed by StatsModels.

standard_errors(sampling_dist3, result_hat)
Sampling std Standard error Percent diff
Intercept 0.236901 0.240243 -1.390859
C(SEX)[T.2] 0.035600 0.031057 14.629106
AGE 0.005769 0.005153 11.939358
AGE2 0.000061 0.000053 15.337773
EDUC 0.031469 0.031193 0.884552
EDUC2 0.001173 0.001166 0.613167

Now we have four ways to compute sampling distributions, confidence intervals, and standard errors – and they yield different results. So you might wonder…

Which One Is Right?#

None of them are right. They are based on different models of the sampling process, so they quantify different sources of uncertainty.

In some cases we might find that different methods converge on the same results, asymptotically, under certain assumptions. But that doesn’t really matter, because with finite data sets, the results are not generally the same. So the important question is whether the differences are big enough to matter in practice.

In this example, it is easy to implement multiple models and compare the results. If they were substantially different, we would need to think more carefully about the modeling assumptions they are based on and choose the one we think is the best description of the data-generating process.

But in this example, the differences are small enough that they probably don’t matter in practice. So we are free to choose whichever is easiest to implement, or fastest to compute, or convenient in some other way.

It is a common error to presume that the result of an analytic method is uniquely correct, and that results from computational methods like resampling are approximations to it. Analytic methods are often fast to compute, but they are always based on modeling assumptions and usually based on approximations, so they are no more correct than computational methods.

Copyright 2023 Allen Downey

Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)