Solutions Tutorial 6: Binary Outcomes (Economics Track)

Retrieval practice

As the predicted probability that $y_i=1$; it can fall outside $[0,1]$, and its errors are inherently heteroskedastic because $\operatorname{Var}(\epsilon_i\mid X_i)=p_i(1-p_i)$.
Probit assumes $\epsilon_i$ is standard normal and logit assumes it is standard logistic; using a CDF gives an S-shaped curve that keeps the predicted probability inside $[0,1]$.
A coefficient changes the index; a continuous regressor’s probability effect is $f(X\beta)\beta_k$, while a binary regressor requires a discrete change between two predicted probabilities.
Compute the derivative for a continuous regressor, or the discrete change for a binary regressor, at each observation’s covariates, then average over the sample.
Maximum likelihood, which chooses the parameters that make the observed data most probable.
Because the effect of a regressor on the probability is scaled by the density $f(\cdot)$, which is evaluated at a different index value in each group, so the probability interaction depends on all coefficients and covariates, not only on the interaction coefficient.
When callbacks are rare, predicting no callback correctly classifies most resumes but detects none of the actual callbacks.
A threshold is absorbed into the intercept, and multiplying the latent coefficients and error by the same positive constant leaves the observed binary outcomes unchanged.

Solutions

1. Interpreting binary outcome models

a) Eligibility for a 401(k)

The coefficient on e401k, eligibility for a 401(k), is positive in all three models. It is approximately 0.0192 in the LPM ($p=0.034$), 0.1548 in logit ($p=0.0051$), and 0.0895 in probit ($p=0.0052$). All three reject a zero coefficient at 5%, but only logit and probit do so at 1%. The direction of the association agrees; significance at 1% does not.

b) Income

The income coefficient is positive and significant at 1% in all three models. The estimates are approximately 0.00623 in the LPM, 0.03368 in logit, and 0.01986 in probit. Income is measured in thousands, so the LPM estimate corresponds to about 0.623 percentage points per additional $1,000, holding the other regressors fixed. The logit and probit coefficients are changes in their indices.

c) Interpreting the LPM coefficient

The LPM coefficient on e401k is approximately 0.0192. Eligible individuals have an estimated 1.92 percentage point higher probability of IRA participation, holding the included controls fixed. This is a conditional association. The regression alone does not establish a causal effect of eligibility or whether workplace retirement saving crowds out IRA saving.

d) Comparing fit measures

The LPM’s $R^2$ measures the fraction of sample outcome variance explained by the fitted values. The logit and probit entries are McFadden pseudo-$R^2$ values, $1-\ell_{\text{full}}/\ell_{\text{null}}$, comparing each model with its intercept-only fit. They do not measure explained variance and should not be compared numerically with the LPM’s $R^2$. Compare likelihood-based fit measures for the same outcome on the same sample; they do not by themselves establish predictive performance on new data.

2. Estimating binary outcome models

a) Linear probability model

The following code estimates the three models with heteroskedasticity-robust standard errors and reports them in one table.

White-sounding names are the reference category. Every name coefficient and discrete change below refers to a black-sounding name relative to a white-sounding name, reversing the lecture’s white indicator.

library(haven)
library(fixest)
library(modelsummary)

url <- paste0(
  "https://empirical-economics.netlify.app/tutorials/datafiles/",
  "lakisha_aer.dta"
)
lakisha <- read_dta(url)
lakisha$race <- factor(lakisha$race, levels = c("w", "b"))

lpm_race <- feols(call ~ race, data = lakisha, vcov = "hc1")
logit_race <- feglm(
  call ~ race, data = lakisha,
  family = binomial(link = "logit"), vcov = "hc1"
)
probit_race <- feglm(
  call ~ race, data = lakisha,
  family = binomial(link = "probit"), vcov = "hc1"
)

modelsummary(
  list("LPM" = lpm_race, "Logit" = logit_race, "Probit" = probit_race),
  stars = c("*" = .1, "**" = .05, "***" = .01),
  coef_map = c("raceb" = "Black-sounding name"),
  gof_map = c("nobs", "r.squared")
)

	LPM	Logit	Probit
* p < 0.1, p < 0.05, * p < 0.01
Black-sounding name	-0.032***	-0.438***	-0.217***
	(0.008)	(0.107)	(0.053)
Num.Obs.	4870	4870	4870
R2	0.003	0.006	0.006

import pandas as pd
import statsmodels.formula.api as smf

url = (
    "https://empirical-economics.netlify.app/tutorials/datafiles/"
    "lakisha_aer.dta"
)
df = pd.read_stata(url)

df["race"] = pd.Categorical(
    df["race"], categories=["w", "b"], ordered=True
)

lpm = smf.ols("call ~ race", data=df).fit(cov_type="HC1")
logit = smf.logit("call ~ race", data=df).fit(disp=0, cov_type="HC1")
probit = smf.probit("call ~ race", data=df).fit(disp=0, cov_type="HC1")

term = "race[T.b]"
print(pd.DataFrame({
    "Model": ["LPM", "Logit", "Probit"],
    "Coefficient": [m.params[term] for m in (lpm, logit, probit)],
    "p-value": [m.pvalues[term] for m in (lpm, logit, probit)],
}))

    Model  Coefficient   p-value
0     LPM    -0.032033  0.000039
1   Logit    -0.438180  0.000044
2  Probit    -0.216543  0.000041

The LPM coefficient on a black-sounding name is $-0.0320$, with $p<0.001$. It is the estimated difference in callback probability relative to a white-sounding name.

b) Logit model

The logit index coefficient is $-0.4382$, with $p<0.001$. Its sign shows that a black-sounding name lowers the fitted latent index, but its magnitude is not a percentage-point effect.

c) Probit model and comparison

The probit index coefficient is $-0.2165$, with $p<0.001$. All three models therefore agree on direction and statistical significance: applicants with a black-sounding name have a lower callback probability. The three coefficient magnitudes should not be compared directly because the models use different scales.

3. Marginal effects versus coefficients

a) Calculating the average discrete changes

library(marginaleffects)
avg_comparisons(logit_race, variables = "race")


 Estimate Std. Error     z Pr(>|z|)    S   2.5 %  97.5 %
   -0.032    0.00778 -4.11   <0.001 14.7 -0.0473 -0.0168

Term: race
Type: response
Comparison: b - w

avg_comparisons(probit_race, variables = "race")


 Estimate Std. Error     z Pr(>|z|)    S   2.5 %  97.5 %
   -0.032    0.00778 -4.11   <0.001 14.7 -0.0473 -0.0168

Term: race
Type: response
Comparison: b - w

from marginaleffects import avg_comparisons

print(avg_comparisons(logit, variables="race"))

shape: (1, 7)
┌───────────┬───────────┬──────────┬──────────┬──────────┬───────────┬───────────┐
│ Estimate  ┆ Std.Error ┆ z        ┆ P(>|z|)  ┆ S        ┆ 2.5%      ┆ 97.5%     │
│ ---       ┆ ---       ┆ ---      ┆ ---      ┆ ---      ┆ ---       ┆ ---       │
│ f64       ┆ f64       ┆ f64      ┆ f64      ┆ f64      ┆ f64       ┆ f64       │
╞═══════════╪═══════════╪══════════╪══════════╪══════════╪═══════════╪═══════════╡
│ -0.032033 ┆ 0.007783  ┆ -4.11555 ┆ 0.000039 ┆ 14.66008 ┆ -0.047288 ┆ -0.016778 │
└───────────┴───────────┴──────────┴──────────┴──────────┴───────────┴───────────┘
Term: race
Contrast: b - w

print(avg_comparisons(probit, variables="race"))

shape: (1, 7)
┌───────────┬───────────┬──────────┬──────────┬──────────┬───────────┬───────────┐
│ Estimate  ┆ Std.Error ┆ z        ┆ P(>|z|)  ┆ S        ┆ 2.5%      ┆ 97.5%     │
│ ---       ┆ ---       ┆ ---      ┆ ---      ┆ ---      ┆ ---       ┆ ---       │
│ f64       ┆ f64       ┆ f64      ┆ f64      ┆ f64      ┆ f64       ┆ f64       │
╞═══════════╪═══════════╪══════════╪══════════╪══════════╪═══════════╪═══════════╡
│ -0.032033 ┆ 0.007783  ┆ -4.11555 ┆ 0.000039 ┆ 14.66008 ┆ -0.047288 ┆ -0.016778 │
└───────────┴───────────┴──────────┴──────────┴──────────┴───────────┴───────────┘
Term: race
Contrast: b - w

b) Interpreting the probit discrete change

The average discrete change is approximately $-0.0320$. Switching from a white-sounding to a black-sounding name reduces the fitted callback probability by 3.20 percentage points. The comparison is between the two name categories, rather than a derivative with respect to a dummy.

c) Comparing the three probability effects

The two discrete changes equal the LPM coefficient, approximately $-0.0320$. With an intercept and one dummy, each model can fit the two group probabilities separately. Logit and probit maximize one binary log-likelihood per group at that group’s observed callback rate; OLS also reproduces the two group means. All three therefore give exactly the same probability difference, up to numerical precision. Equality follows from this specification, rather than from small effects or approximately linear probability curves.

4. Interactions in a nonlinear model

a) The difference between the two name gaps

The name gap is $p_{10}-p_{00}$ for men and $p_{11}-p_{01}$ for women. The interaction on the probability scale is therefore

\[ \Delta=(p_{11}-p_{01})-(p_{10}-p_{00}). \]

Both regressors are binary, so this is a difference between discrete changes. In this model there are four fitted cell probabilities and one such interaction contrast.

b) Probability and index interactions

Under probit,

\[ \Delta=\Phi(\beta_0+\beta_1+\beta_2+\beta_3) -\Phi(\beta_0+\beta_2) -\Phi(\beta_0+\beta_1) +\Phi(\beta_0). \]

This expression depends on all four coefficients. It can be non-zero when $\beta_3=0$, and its sign can differ from the sign of $\beta_3$. Testing $\beta_3=0$ tests an interaction in the index. Testing $\Delta=0$ instead requires the uncertainty of this probability contrast, which software can calculate using the coefficient covariance matrix.

For an LPM with the same regressors, substituting its four cell probabilities into the contrast gives $\Delta=\beta_3$. Its interaction coefficient and associated test therefore refer directly to this difference in probability changes.

5. The name gap for men and women

a) The LPM interaction

library(marginaleffects)

lakisha$female <- as.integer(lakisha$sex == "f")

lpm_int <- feols(call ~ race * female, data = lakisha, vcov = "hc1")
logit_int <- glm(call ~ race * female, data = lakisha, family = binomial("logit"))

modelsummary(
  list("LPM" = lpm_int, "Logit" = logit_int),
  stars = c("*" = .1, "**" = .05, "***" = .01),
  coef_map = c("raceb" = "Black-sounding name", "female" = "Female",
               "raceb:female" = "Black-sounding name x female"),
  gof_map = c("nobs")
)

	LPM	Logit
* p < 0.1, p < 0.05, * p < 0.01
Black-sounding name	-0.030*	-0.453*
	(0.016)	(0.234)
Female	0.010	0.120
	(0.014)	(0.166)
Black-sounding name x female	-0.002	0.017
	(0.018)	(0.263)
Num.Obs.	4870	4870

import numpy as np
import statsmodels.formula.api as smf
from marginaleffects import avg_comparisons, datagrid

df["female"] = (df["sex"] == "f").astype(int)

lpm_int = smf.ols("call ~ race * female", data=df).fit(cov_type="HC1")
logit_int = smf.logit("call ~ race * female", data=df).fit(disp=0)

print(lpm_int.summary().tables[1])

====================================================================================
                       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------------
Intercept            0.0887      0.012      7.478      0.000       0.065       0.112
race[T.b]           -0.0304      0.016     -1.960      0.050      -0.061     3.3e-06
female               0.0102      0.014      0.745      0.456      -0.017       0.037
race[T.b]:female    -0.0022      0.018     -0.125      0.901      -0.037       0.033
====================================================================================

print(logit_int.summary().tables[1])

====================================================================================
                       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------------
Intercept           -2.3297      0.147    -15.882      0.000      -2.617      -2.042
race[T.b]           -0.4526      0.234     -1.935      0.053      -0.911       0.006
female               0.1204      0.166      0.726      0.468      -0.205       0.446
race[T.b]:female     0.0165      0.263      0.063      0.950      -0.499       0.533
====================================================================================

The LPM interaction coefficient is -0.0022. For men, a black-sounding name lowers the callback probability by 3.0 percentage points; for women the change is 0.22 percentage points larger in absolute value. With $p=0.90$, there is no evidence that the name gap differs between men and women.

b) The logit interaction

The logit interaction coefficient is 0.0165, with $p=0.95$. It is positive, whereas the LPM coefficient is negative. Both are far from significant, so the substantive conclusion is the same, but the sign alone already shows that the two coefficients do not measure the same thing.

c) The interaction on the probability scale

avg_comparisons(logit_int, variables = "race", by = "female")


 female Estimate Std. Error     z Pr(>|z|)    S   2.5 %    97.5 %
      0  -0.0304    0.01551 -1.96   0.0499  4.3 -0.0608 -9.36e-06
      1  -0.0326    0.00899 -3.63   <0.001 11.8 -0.0503 -1.50e-02

Term: race
Type: response
Comparison: b - w

avg_comparisons(logit_int, variables = "race", by = "female",
                hypothesis = ~ pairwise)


 Hypothesis Estimate Std. Error      z Pr(>|z|)   S   2.5 % 97.5 %
  (1) - (0) -0.00224     0.0179 -0.125    0.901 0.2 -0.0374 0.0329

Type: response

print(avg_comparisons(logit_int, variables="race", by="female"))

shape: (2, 8)
┌────────┬───────────┬───────────┬───────────┬──────────┬───────────┬───────────┬───────────┐
│ female ┆ Estimate  ┆ Std.Error ┆ z         ┆ P(>|z|)  ┆ S         ┆ 2.5%      ┆ 97.5%     │
│ ---    ┆ ---       ┆ ---       ┆ ---       ┆ ---      ┆ ---       ┆ ---       ┆ ---       │
│ i64    ┆ f64       ┆ f64       ┆ f64       ┆ f64      ┆ f64       ┆ f64       ┆ f64       │
╞════════╪═══════════╪═══════════╪═══════════╪══════════╪═══════════╪═══════════╪═══════════╡
│ 0      ┆ -0.030408 ┆ 0.01551   ┆ -1.960557 ┆ 0.049931 ┆ 4.323927  ┆ -0.060807 ┆ -0.000009 │
│ 1      ┆ -0.032647 ┆ 0.008985  ┆ -3.633336 ┆ 0.00028  ┆ 11.803418 ┆ -0.050258 ┆ -0.015036 │
└────────┴───────────┴───────────┴───────────┴──────────┴───────────┴───────────┴───────────┘
Term: race
Contrast: b - w

print(avg_comparisons(logit_int, variables="race", by="female",
                      hypothesis="b1 - b0 = 0"))

shape: (1, 7)
┌───────────┬───────────┬───────────┬──────────┬──────────┬───────────┬──────────┐
│ Estimate  ┆ Std.Error ┆ z         ┆ P(>|z|)  ┆ S        ┆ 2.5%      ┆ 97.5%    │
│ ---       ┆ ---       ┆ ---       ┆ ---      ┆ ---      ┆ ---       ┆ ---      │
│ f64       ┆ f64       ┆ f64       ┆ f64      ┆ f64      ┆ f64       ┆ f64      │
╞═══════════╪═══════════╪═══════════╪══════════╪══════════╪═══════════╪══════════╡
│ -0.002239 ┆ 0.017925  ┆ -0.124914 ┆ 0.900591 ┆ 0.151055 ┆ -0.037371 ┆ 0.032893 │
└───────────┴───────────┴───────────┴──────────┴──────────┴───────────┴──────────┘
Term: b1-b0=0

A black-sounding name lowers the callback probability by 3.04 percentage points for men and by 3.26 for women. The difference, -0.0022, is the interaction effect on the probability. It has the sign of the LPM coefficient, not of the logit coefficient, and is equal to the LPM coefficient.

d) Why the logit and LPM agree here

race and female are both dummies, and the model includes their interaction, so it has four parameters for four cells (white men, white women, black men, black women). Such a model is saturated. As with the single dummy in the lecture, the log-likelihood splits into one coin log-likelihood per cell, and maximum likelihood sets each fitted probability equal to that cell’s callback rate. OLS does the same, because the LPM is a regression on a full set of cell dummies. Both models therefore fit the same four probabilities, and every difference between them, including the difference in the name gap between men and women, is identical. The logit interaction coefficient is different because it is measured in log-odds: it is the change in the log odds ratio, not in the probability difference.

e) An interaction without an interaction term

logit_exp <- glm(call ~ race + yearsexp, data = lakisha,
                 family = binomial("logit"))
avg_comparisons(
  logit_exp, variables = "race",
  newdata = datagrid(yearsexp = c(2, 20), grid_type = "counterfactual"),
  by = "yearsexp"
)


 yearsexp Estimate Std. Error     z Pr(>|z|)    S   2.5 %  97.5 %
        2  -0.0260    0.00648 -4.01   <0.001 14.0 -0.0387 -0.0133
       20  -0.0462    0.01180 -3.92   <0.001 13.4 -0.0694 -0.0231

Term: race
Type: response
Comparison: b - w

logit_exp = smf.logit("call ~ race + yearsexp", data=df).fit(disp=0)
print(avg_comparisons(
    logit_exp, variables="race",
    newdata=datagrid(model=logit_exp, yearsexp=[2, 20],
                     grid_type="counterfactual"),
    by="yearsexp",
))

shape: (2, 8)
┌──────────┬───────────┬───────────┬───────────┬──────────┬───────────┬───────────┬───────────┐
│ yearsexp ┆ Estimate  ┆ Std.Error ┆ z         ┆ P(>|z|)  ┆ S         ┆ 2.5%      ┆ 97.5%     │
│ ---      ┆ ---       ┆ ---       ┆ ---       ┆ ---      ┆ ---       ┆ ---       ┆ ---       │
│ i64      ┆ f64       ┆ f64       ┆ f64       ┆ f64      ┆ f64       ┆ f64       ┆ f64       │
╞══════════╪═══════════╪═══════════╪═══════════╪══════════╪═══════════╪═══════════╪═══════════╡
│ 2        ┆ -0.025968 ┆ 0.006476  ┆ -4.009773 ┆ 0.000061 ┆ 14.006111 ┆ -0.038661 ┆ -0.013275 │
│ 20       ┆ -0.046237 ┆ 0.011803  ┆ -3.917363 ┆ 0.00009  ┆ 13.447382 ┆ -0.069371 ┆ -0.023103 │
└──────────┴───────────┴───────────┴───────────┴──────────┴───────────┴───────────┴───────────┘
Term: race
Contrast: b - w

With 2 years of experience the name gap is 2.6 percentage points; with 20 years it is 4.6 percentage points. The index coefficient on race is the same for both, but more experience raises the index and moves the resumes up the lower part of the S-curve, where the curve is steeper. The same index shift then produces a larger change in probability. Thus a regressor’s probability effect can depend on another regressor even without an interaction term in the index.

6. Where the logit estimates come from

a) The intercept-only logit

logit_null <- glm(call ~ 1, data = lakisha, family = binomial("logit"))
alpha_hat <- coef(logit_null)[[1]]
c(alpha_hat = alpha_hat, fitted = plogis(alpha_hat),
  sample_mean = mean(lakisha$call))

  alpha_hat      fitted sample_mean 
-2.43566996  0.08049281  0.08049281

logit_null = smf.logit("call ~ 1", data=df).fit(disp=0)
alpha_hat = logit_null.params.iloc[0]
print(alpha_hat, 1 / (1 + np.exp(-alpha_hat)), df["call"].mean())

-2.4356699573682383 0.08049281314168372 0.08049282

The estimate is $\hat\alpha=-2.4357$ and $\Lambda(\hat\alpha)=0.0805$, which is the share of resumes with a callback, $392/4870$. With only an intercept, every resume has the same callback probability, so each resume is a flip of the same coin. The coin example showed that the likelihood is then maximized at the sample share of successes. The logit only expresses that share on the log-odds scale: $\hat\alpha=\log\{\hat p/(1-\hat p)\}$.

b) Evaluating the log-likelihood

With $n_1$ callbacks and $n_0$ resumes without a callback,

\[ \ell(p)=n_1\log p+n_0\log(1-p). \]

At $\hat p=0.0805$ this gives $392\log(0.0805)+4478\log(0.9195)=-1363.46$, which equals the value reported by logLik(logit_null) (R) or logit_null.llf (Python): -1363.46.

c) One probability per group

logit_race_ml <- glm(call ~ race, data = lakisha, family = binomial("logit"))
b <- coef(logit_race_ml)
rates <- tapply(lakisha$call, lakisha$race, mean)

c(fitted_white = plogis(b[["(Intercept)"]]),
  fitted_black = plogis(b[["(Intercept)"]] + b[["raceb"]]))

fitted_white fitted_black 
  0.09650924   0.06447639

rates

         w          b 
0.09650924 0.06447639

# Coefficient on race versus the log of the ratio of the two odds
c(coefficient = b[["raceb"]],
  log_odds_ratio = log((rates[["b"]] / (1 - rates[["b"]])) /
                         (rates[["w"]] / (1 - rates[["w"]]))))

   coefficient log_odds_ratio 
    -0.4381802     -0.4381802

logit_race_ml = smf.logit("call ~ race", data=df).fit(disp=0)
a, g = logit_race_ml.params.iloc[0], logit_race_ml.params.iloc[1]
rates = df.groupby("race", observed=True)["call"].mean()
odds = rates / (1 - rates)

print(1 / (1 + np.exp(-a)), 1 / (1 + np.exp(-(a + g))))

0.09650924024640664 0.06447638603696082

print(rates)

race
w    0.096509
b    0.064476
Name: call, dtype: float32

print(g, np.log(odds["b"] / odds["w"]))

-0.43818021345657443 -0.43818015

The fitted probabilities are 0.0965 for white-sounding and 0.0645 for black-sounding names, exactly the two callback rates. The log-likelihood splits into one coin log-likelihood for each group, and each is maximized at that group’s share of callbacks.

Because $\log\{p/(1-p)\}=\alpha+\gamma\,\text{race}_i$, the coefficient $\gamma$ is the difference between the two groups’ log odds, which is the log of their odds ratio: both numbers equal -0.4382.

d) The likelihood-ratio test

\[ LR=2(\ell_U-\ell_R)=2\{-1354.97-(-1363.46)\}=16.98. \]

The restricted model sets the coefficient on race to zero, so there is one restriction and the statistic is compared with a $\chi^2_1$ distribution, whose 5% critical value is 3.84. Since 16.98 > 3.84 ($p<0.001$), we reject that the name has no effect on the callback probability.

The squared $z$-statistic on race is $(-4.08)^2=16.67$. The two statistics test the same hypothesis and are asymptotically equivalent; in a sample of this size they are close but not identical, because the $z$-test uses the curvature at the unrestricted maximum and the LR test compares the heights of the two maxima.

Self-study

7. Average effects and effects at the mean

a) Observation-specific effects of experience

The experience coefficient changes the logit index by a constant amount. Its probability effect is $\hat p_i(1-\hat p_i)\hat\beta_{\text{experience}}$, which varies because the fitted probabilities differ across resumes.

logit_full <- glm(
  call ~ race + female + yearsexp + factor(occupbroad),
  data = lakisha, family = binomial("logit")
)
p_full <- predict(logit_full, type = "response")
beta_exp <- coef(logit_full)[["yearsexp"]]
me_exp <- p_full * (1 - p_full) * beta_exp

knitr::kable(data.frame(
  Quantity = c("Minimum", "Maximum"),
  `Percentage points per year` = 100 * range(me_exp),
  check.names = FALSE
), digits = 3)

Quantity	Percentage points per year
Minimum	0.160
Maximum	0.878

logit_full = smf.logit(
    "call ~ race + female + yearsexp + C(occupbroad)", data=df
).fit(disp=0)
p_full = logit_full.predict(df)
beta_exp = logit_full.params["yearsexp"]
me_exp = p_full * (1 - p_full) * beta_exp

print(pd.DataFrame({
    "Quantity": ["Minimum", "Maximum"],
    "Percentage points per year": 100 * np.array([me_exp.min(), me_exp.max()]),
}).round(3))

  Quantity  Percentage points per year
0  Minimum                       0.160
1  Maximum                       0.878

The effects range from 0.160 to 0.878 percentage points per year. These are local slopes, so they approximate the change from one additional year rather than calculating an exact finite change.

b) The average marginal effect

Average the slopes calculated at the observed resumes:

\[ AME_{\text{experience}}=\frac{1}{n}\sum_i \hat p_i(1-\hat p_i)\hat\beta_{\text{experience}}. \]

ame_exp <- mean(me_exp)
100 * ame_exp

[1] 0.3484451

ame_exp = me_exp.mean()
print(100 * ame_exp)

0.3484451357086889

The AME is 0.348 percentage points per year: the average local change in fitted callback probability per additional year of experience, holding each resume’s other characteristics fixed.

c) The marginal effect at the mean

The MEM evaluates the derivative at the mean of the model’s regressor vector, including its dummy columns:

\[ \bar p=\Lambda(\bar x'\hat\beta),\qquad MEM_{\text{experience}}=\bar p(1-\bar p)\hat\beta_{\text{experience}}. \]

x_bar <- colMeans(model.matrix(logit_full))
p_at_mean <- plogis(sum(x_bar * coef(logit_full)))
mem_exp <- p_at_mean * (1 - p_at_mean) * beta_exp

knitr::kable(data.frame(
  Effect = c("AME", "MEM"),
  `Percentage points per year` = 100 * c(ame_exp, mem_exp),
  check.names = FALSE
), digits = 3)

Effect	Percentage points per year
AME	0.348
MEM	0.335

x_bar = logit_full.model.exog.mean(axis=0)
p_at_mean = 1 / (1 + np.exp(-(x_bar @ logit_full.params.to_numpy())))
mem_exp = p_at_mean * (1 - p_at_mean) * beta_exp

print(pd.DataFrame({
    "Effect": ["AME", "MEM"],
    "Percentage points per year": 100 * np.array([ame_exp, mem_exp]),
}).round(3))

  Effect  Percentage points per year
0    AME                       0.348
1    MEM                       0.335

The MEM is 0.335 percentage points per year. It need not equal the AME because evaluating a non-linear function at the average covariates need not equal averaging its values. The dummy-column means are sample proportions: an actual resume has one name category, one sex category, and one occupation category. The mean regressor vector combines fractional indicators and need not represent an observed resume.

8. The latent-variable framework

a) From a latent propensity to an observed outcome

The latent variable $y_i^*=X_i'\beta+\epsilon_i$ is an unobserved continuous propensity. We observe $y_i=1$ when it exceeds a threshold and $y_i=0$ otherwise. For example, the lecture represents the employer’s valuation of a resume as latent and the callback decision as observed. With the threshold normalized to zero, $y_i=1\{y_i^*>0\}$.

b) What distinguishes probit from logit

The standard normal and standard logistic distributions are symmetric around zero. For either choice, continuity and symmetry give

\[ \begin{aligned} P(y_i=1\mid X_i) &=P(\epsilon_i>-X_i'\beta\mid X_i)\\ &=1-F(-X_i'\beta)=F(X_i'\beta). \end{aligned} \]

Probit uses $F=\Phi$, the standard normal CDF; logit uses $F=\Lambda$, the standard logistic CDF. The symmetry step applies to these distributions, rather than to probability distributions in general.

c) Threshold and scale normalizations

If the original threshold is $a$, the event $X_i'\beta+\epsilon_i>a$ can be written with a zero threshold by subtracting $a$ from the intercept. This leaves the observed outcomes unchanged.

The scale is also unidentified. For any $c>0$,

\[ 1\{X_i'\beta+\epsilon_i>0\} =1\{cX_i'\beta+c\epsilon_i>0\}. \]

We therefore fix the error variance to define the latent scale. Probit assumes variance 1; standard logit assumes variance $\pi^2/3$. Their raw coefficients use different scales, so compare probability effects when interpreting magnitudes across these models.

9. Evaluating callback predictions

a) Classification at a threshold of 0.5

The fitted probabilities come from the full logit in Question 7. The following helper counts all four cells, including cells with zero observations.

classify <- function(actual, probability, threshold) {
  predicted <- as.integer(probability > threshold)
  confusion <- table(
    Actual = factor(actual, levels = 0:1),
    Predicted = factor(predicted, levels = 0:1)
  )
  tn <- confusion[1, 1]
  fp <- confusion[1, 2]
  fn <- confusion[2, 1]
  tp <- confusion[2, 2]
  metrics <- c(
    Accuracy = (tn + tp) / sum(confusion),
    Sensitivity = tp / (tp + fn),
    Specificity = tn / (tn + fp),
    Precision = if (tp + fp > 0) tp / (tp + fp) else NA_real_
  )
  list(confusion = confusion, metrics = metrics)
}

result_05 <- classify(lakisha$call, p_full, 0.5)
max(p_full)

[1] 0.2437187

knitr::kable(result_05$confusion, caption = "Threshold 0.5: rows actual, columns predicted")

Threshold 0.5: rows actual, columns predicted
	0	1
0	4478	0
1	392	0

def classify(actual, probability, threshold):
    predicted = (probability > threshold).astype(int)
    confusion = pd.crosstab(
        pd.Categorical(actual, categories=[0, 1]),
        pd.Categorical(predicted, categories=[0, 1]),
        rownames=["Actual"], colnames=["Predicted"], dropna=False,
    )
    tn, fp, fn, tp = confusion.to_numpy().ravel()
    metrics = pd.Series({
        "Accuracy": (tn + tp) / confusion.to_numpy().sum(),
        "Sensitivity": tp / (tp + fn),
        "Specificity": tn / (tn + fp),
        "Precision": tp / (tp + fp) if tp + fp > 0 else np.nan,
    })
    return confusion, metrics

confusion_05, metrics_05 = classify(df["call"], p_full, 0.5)
print("Largest fitted probability:", p_full.max())

Largest fitted probability: 0.24371868025986707

print(confusion_05)

Predicted     0  1
Actual            
0          4478  0
1           392  0

The largest fitted probability is 0.244, below 0.5. The classifier predicts no callbacks: there are 4478 true negatives and 392 false negatives, with zero true positives and false positives.

b) Accuracy and the majority-class baseline

baseline <- mean(lakisha$call == 0)
knitr::kable(data.frame(
  Metric = names(result_05$metrics),
  Value = unname(result_05$metrics)
), digits = 4)

Metric	Value
Accuracy	0.9195
Sensitivity	0.0000
Specificity	1.0000
Precision	NA

baseline

[1] 0.9195072

baseline = (df["call"] == 0).mean()
print(metrics_05.round(4))

Accuracy       0.9195
Sensitivity    0.0000
Specificity    1.0000
Precision         NaN
dtype: float64

print("Majority-class accuracy:", baseline)

Majority-class accuracy: 0.9195071868583162

Accuracy is 91.95%, exactly the majority-class baseline. Sensitivity, $TP/(TP+FN)$, is zero; specificity, $TN/(TN+FP)$, is one. Precision, $TP/(TP+FP)$, is undefined because no positive predictions are made. High accuracy here reflects the rarity of callbacks; none of the actual callbacks are detected.

c) Classification at a threshold of 0.1

result_01 <- classify(lakisha$call, p_full, 0.1)
knitr::kable(result_01$confusion, caption = "Threshold 0.1: rows actual, columns predicted")

Threshold 0.1: rows actual, columns predicted
	0	1
0	3522	956
1	250	142

knitr::kable(data.frame(
  Metric = names(result_01$metrics),
  Value = unname(result_01$metrics)
), digits = 4)

Metric	Value
Accuracy	0.7524
Sensitivity	0.3622
Specificity	0.7865
Precision	0.1293

confusion_01, metrics_01 = classify(df["call"], p_full, 0.1)
print(confusion_01)

Predicted     0    1
Actual              
0          3522  956
1           250  142

print(metrics_01.round(4))

Accuracy       0.7524
Sensitivity    0.3622
Specificity    0.7865
Precision      0.1293
dtype: float64

At 0.1, the classifier detects 142 callbacks and makes 956 false positive predictions. Sensitivity rises to 36.2%, while specificity falls to 78.7%. Accuracy is 75.2%, and precision is 12.9%.

The lower threshold detects some callbacks at the cost of more false positives. Choosing a threshold requires the relative costs of those errors. Accuracy alone weights correct positive and negative classifications equally and, with this rare outcome, favours predicting no callbacks.

d) Evaluation on new resumes

These are training-sample results: the observations were used to estimate the coefficients. Reserve test data before fitting or selecting models. Choose the specification and threshold using training data, with a validation set or cross-validation within those data if needed, then evaluate the chosen procedure on the untouched test set. Keep related observations, such as resumes sent to the same vacancy, together when constructing splits.

A larger McFadden pseudo-$R^2$ indicates better likelihood-based fit relative to the intercept-only model on the same estimation sample. It does not establish better performance on new observations. Evaluate probability predictions as well as threshold-based classifications on the test data.

Takeaways

What should you be able to do now?

Interpret a linear-probability-model coefficient as a change in probability in percentage points, and name the two costs of that simplicity.
Estimate LPM, logit, and probit models in R or Python, and compare their coefficients, significance, and reported fit measures.
Compute and interpret average discrete changes for binary regressors and distinguish a continuous regressor’s AME from its MEM.
Explain the latent-variable motivation, distributional assumptions, and threshold and scale normalizations of probit and logit.
Compute an interaction effect on the probability scale and explain why it can differ from the interaction coefficient of a logit or probit.
Explain how maximum likelihood produces the logit estimates in simple cases, evaluate a log-likelihood, and carry out a likelihood-ratio test.
Evaluate classification against the majority-class baseline, interpret a confusion matrix, and distinguish training-sample fit from performance on new observations.