Logit and probit regression in Python with marginal effects

Many outcomes in economics are yes-or-no: a person is employed or not, a firm exports or not, a borrower defaults or not. A linear regression of a 0/1 outcome, the linear probability model, is easy to read but can predict probabilities below 0 or above 1, and it forces a regressor to shift the probability by the same amount whether the baseline is 2% or 50%. Logit and probit model the probability itself, so predictions stay between 0 and 1 and effects shrink near the extremes.

What logit and probit estimate

Both models start from a latent index. An unobserved propensity y* = x′β + ε decides the outcome: y = 1 when y* > 0 and y = 0 otherwise. The two models differ only in the distribution of ε.

Pr(yᵢ = 1 | xᵢ) = F(xᵢ′β),   F = Λ (logistic) for logit, F = Φ (standard normal) for probit

The coefficients are estimated by maximum likelihood. Their scale is fixed by normalizing the error variance: the logistic error has standard deviation π/√3 ≈ 1.8 rather than 1, so logit coefficients are typically 1.6 to 1.8 times their probit counterparts. The two models usually give very similar predicted probabilities. Logit has one convenience: exp(β) is an odds ratio. Probit connects naturally to models with normal errors, such as selection and bivariate probit models.

A coefficient tells you the sign of an effect and the change in the index, not the change in the probability. The effect of a continuous regressor xₖ on the probability is f(x′β)βₖ, where f is the density of the error, so it differs from person to person. Two summaries are common:

• The average marginal effect (AME) evaluates the effect for every observation and averages. For a binary or categorical regressor it averages the discrete change in the probability from the reference level.
• The marginal effect at the means (MEM) evaluates the effect once, for a hypothetical person with average characteristics.

AMEs are the better default for reporting. Because the latent scale is normalized, coefficients from different models or samples are not directly comparable, while AMEs are less sensitive to that problem (Mood 2010). The coefficient on an interaction term is not the interaction effect on the probability either, and it can even have the opposite sign (Ai and Norton 2003).

The logit takes its name from Berkson (1944). A finite maximum-likelihood estimate does not exist when a regressor, or a combination of regressors, predicts the outcome perfectly: complete or quasi-complete separation (Albert and Anderson 1984).

In OpenEconometrics

oe.logit and oe.probit take the same arguments. The outcome must be coded 0/1. The model below corresponds to the Stata commands logit employed educ kids i.married followed by margins, dydx(*):

• oe.logit(data=df, y="employed", x=["educ", "kids", "married"], categorical=["married"]). Columns in categorical are expanded into treatment dummies named married[1], with the first level omitted. intercept=False drops the constant.
• The default covariance is "nonrobust", the inverse of the information matrix. cluster="firm" selects cluster-robust standard errors on one column; the finite-sample factor is recorded in fit.inference["correction"].
• weights= with weight_type= ("fweight", "aweight", "pweight" or "iweight") uses a weighted maximum-likelihood route that also accepts covariance="robust" and "opg". Sampling weights default to robust. Analytic weights are an OpenEconometrics extension, not a Stata logit option.
• If the regressors separate the outcome, the fit stops with an AnalysisError whose code is separation_detected, rather than returning huge coefficients.

Marginal effects come from the saved result. oe.margins(fit, ["educ", "kids", "married"], data=df) returns AMEs with delta-method standard errors and confidence limits, using the full coefficient covariance. method="mem" gives effects at the means, at={"educ": [10, 16]} evaluates on a grid, and kind="xb" reports effects on the index scale. oe.predict(fit, new_rows, kind="response", interval="mean") gives predicted probabilities with pointwise intervals. Neither function refits the model. Both need the evaluation rows explicitly: the AME is the average over the rows you pass. The prediction guide documents both functions, and the GLM guide covers other binary links, offsets and grouped data through oe.glm(family="binomial").

Supported options

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates 4,000 adults whose employment follows a true logit model. The latent index is −3.5 + 0.30·educ − 0.50·kids + 0.40·married, plus a standard logistic error. Schooling is normal with mean 13 years, the number of young children is Poisson with mean 1, and half the sample is married. The script also computes the AMEs implied by the true coefficients on the same sample, so the estimated AMEs have a target.

import numpy as np
import pandas as pd
import openecon as oe

BETA = {"Intercept": -3.5, "educ": 0.30, "kids": -0.50, "married": 0.40}


def simulate(n=4000, seed=2026):
    """Employment from a latent index with a logistic error: a true logit model."""
    rng = np.random.default_rng(seed)
    educ = rng.normal(13, 2.5, size=n)     # years of schooling
    kids = rng.poisson(1.0, size=n)        # number of young children
    married = rng.integers(0, 2, size=n)
    index = BETA["Intercept"] + BETA["educ"] * educ + BETA["kids"] * kids + BETA["married"] * married
    employed = (index + rng.logistic(size=n) > 0).astype(int)
    return pd.DataFrame({"employed": employed, "educ": educ, "kids": kids, "married": married})


def true_ame(df):
    """Average marginal effects implied by the true coefficients on this sample."""
    cdf = lambda z: 1 / (1 + np.exp(-z))
    base = BETA["Intercept"] + BETA["educ"] * df.educ + BETA["kids"] * df.kids
    index = base + BETA["married"] * df.married
    density = cdf(index) * (1 - cdf(index))
    return {"educ": np.mean(density) * BETA["educ"], "kids": np.mean(density) * BETA["kids"],
            "married[1]": np.mean(cdf(base + BETA["married"]) - cdf(base))}


if __name__ == "__main__":
    df = simulate()
    spec = dict(data=df, y="employed", x=["educ", "kids", "married"], categorical=["married"])
    logit = oe.logit(**spec)
    probit = oe.probit(**spec)
    print(logit.summary())
    print("Probit coefficients:", {c.term: round(c.estimate, 4) for c in probit.coefficients})

    effects = ["educ", "kids", "married"]
    ame_logit = oe.margins(logit, effects, data=df).set_index("variable")
    ame_probit = oe.margins(probit, effects, data=df).set_index("variable")
    table = pd.DataFrame({"true AME": pd.Series(true_ame(df)),
                          "logit AME": ame_logit["estimate"], "logit se": ame_logit["std_error"],
                          "probit AME": ame_probit["estimate"], "probit se": ame_probit["std_error"]})
    print("\nAverage marginal effects on Pr(employed = 1):")
    print(table.round(4).to_string())

With openecon 0.3.20a1 from PyPI the script prints:

OpenEconometrics LOGIT — employed
Observations: 4000  |  Covariance: nonrobust  |  Confidence: 95%

Term         Estimate  Std. error         z     P>|stat|   CI lower   CI upper
----------  ---------  ----------  --------  -----------  ---------  ---------
Intercept    -3.51802    0.202306  -17.3896  9.89224e-68   -3.91453   -3.12151
educ         0.299561   0.0151823   19.7309  1.16999e-86   0.269804   0.329318
kids        -0.485345   0.0362758  -13.3793  7.98952e-41  -0.556444  -0.414246
married[1]   0.235408   0.0687318   3.42503  0.000614741   0.100696    0.37012

pseudo_r_squared: 0.116602  |  aic: 4906.56  |  bic: 4931.73
Probit coefficients: {'Intercept': -2.1287, 'educ': 0.1811, 'kids': -0.2941, 'married[1]': 0.1449}

Average marginal effects on Pr(employed = 1):
            true AME  logit AME  logit se  probit AME  probit se
educ          0.0631     0.0636    0.0026      0.0633     0.0026
kids         -0.1052    -0.1030    0.0071     -0.1028     0.0070
married[1]    0.0845     0.0500    0.0146      0.0507     0.0146

Reading the output

• Schooling and children are recovered closely. The logit estimates 0.2996 for schooling (true 0.30) and −0.4853 for children (true −0.50), both well inside their 95% intervals.
• The marriage coefficient misses. It is 0.2354 with an interval from 0.101 to 0.370, which excludes the true 0.40; the estimate is 2.4 standard errors low. A 95% interval misses the truth in about one sample in twenty, and this is such a sample. Its AME of 0.0500 is correspondingly below the true 0.0845.
• Probit coefficients differ by a scale factor. Dividing the logit by the probit coefficients gives 1.65 for schooling, 1.65 for children and 1.62 for marriage. The models disagree about the latent scale, not about the data.
• The AMEs agree across models. One more year of schooling raises the employment probability by 6.36 percentage points in the logit and 6.33 in the probit, against a true 6.31; each additional child lowers it by about 10.3 points. Report these, not the coefficients, when the question is about probabilities.
• The pseudo R² of 0.117 is McFadden's, one minus the ratio of the log likelihoods of the model and the constant-only model. It is not a share of explained variance, so do not compare it with an OLS R².

Limitations

• Unweighted fits accept only nonrobust and cluster covariance. Heteroskedasticity-robust (robust) and outer-product (opg) standard errors require weights= and weight_type=.
• Weighted fits run in memory on the CPU, with at most 100,000 original rows and 64 expanded parameters; the out-of-core oe.Dataset route is for unweighted fits.
• The binary commands take no offset and no other link. Use oe.glm(family="binomial", ...) for offsets, grouped data and other links, or oe.cloglog.
• Separation stops the fit with separation_detected; perfectly predicted observations are not dropped automatically. Remove or recode the separating regressor yourself.
• Prediction intervals for probabilities are symmetric delta-method intervals and are not clipped to [0, 1]. Classification tables and residuals are not provided.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Ai, C. and Norton, E. C. (2003). Interaction terms in logit and probit models. Economics Letters, 80(1), 123–129. doi:10.1016/S0165-1765(03)00032-6
• Albert, A. and Anderson, J. A. (1984). On the existence of maximum likelihood estimates in logistic regression models. Biometrika, 71(1), 1–10. doi:10.1093/biomet/71.1.1
• Berkson, J. (1944). Application of the logistic function to bio-assay. Journal of the American Statistical Association, 39(227), 357–365. doi:10.1080/01621459.1944.10500699
• Mood, C. (2010). Logistic regression: Why we cannot do what we think we can do, and what we can do about it. European Sociological Review, 26(1), 67–82. doi:10.1093/esr/jcp006