Where Are We At?

Binary Response Models

Binary Outcomes

• Outcomes restricted to FALSE vs TRUE, or 0 vs 1.

• We'd have $y\in \{0,1\}$.

• In those situations we are primarily interested in estimating the response probability or the probability of success:

$$p\left(x\right)=Pr(y=1|x)$$

• how does $p\left(x\right)$ change as we change $x$?

• we ask

  If we increase $x$ by one unit, how would the probability of $y=1$ change?

Remembering Bernoulli Fun

The Linear Probability Model (LPM)

• The simplest option. Model the response probability as

$$Pr(y=1|x)=p\left(x\right)={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_K{x}_K$$

• Interpretation: a 1 unit change in $x_1$, say, results in a change of $p\left(x\right)$ of ${\beta }_1$.

Example: Mroz (1987)

• Female labor market participation

• How does inlf (in labor force) status depend on non-wife household income, her education, age and number of small children?

Mroz 1987

data(mroz, package = "wooldridge")
plot(factor(inlf) ~ age, data = mroz, 
     ylevels = 2:1,
     ylab = "in labor force?")

Running the LPM

LPM = lm(inlf ~ nwifeinc + educ + exper 
         + I(exper^2) + age +I(age^2) + kidslt6, mroz)
broom::tidy(LPM)

## # A tibble: 8 × 5
##   term         estimate std.error statistic  p.value
##   <chr>           <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  0.322     0.486       0.662  5.08e- 1
## 2 nwifeinc    -0.00343   0.00145    -2.36   1.86e- 2
## 3 educ         0.0375    0.00735     5.10   4.33e- 7
## 4 exper        0.0383    0.00577     6.63   6.44e-11
## 5 I(exper^2)  -0.000565  0.000189   -2.98   2.96e- 3
## 6 age         -0.00112   0.0225     -0.0497 9.60e- 1
## 7 I(age^2)    -0.000182  0.000258   -0.706  4.80e- 1
## 8 kidslt6     -0.260     0.0341     -7.64   6.72e-14

LPM: Predicting negative probabilities?!

LPM in Saturated Model: No Problem!

library(dplyr)
mroz %<>% 
  # classify age into 3 and huswage into 2 classes
  mutate(age_fct = cut(age,breaks = 3,labels = FALSE),
         huswage_fct = cut(huswage, breaks = 2,labels = FALSE)) %>%
  mutate(classes = paste0("age_",age_fct,"_hus_",huswage_fct))
LPM_saturated = mroz %>%
  lm(inlf ~ classes, data = .)
broom::tidy(LPM_saturated)

## # A tibble: 6 × 5
##   term               estimate std.error statistic  p.value
##   <chr>                 <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)          0.611     0.0277    22.0   2.98e-83
## 2 classesage_1_hus_2  -0.611     0.350     -1.75  8.11e- 2
## 3 classesage_2_hus_1  -0.0257    0.0404    -0.635 5.25e- 1
## 4 classesage_2_hus_2  -0.277     0.203     -1.37  1.72e- 1
## 5 classesage_3_hus_1  -0.149     0.0494    -3.01  2.72e- 3
## 6 classesage_3_hus_2  -0.111     0.350     -0.317 7.51e- 1

LPM in Saturated Model: No Problem!

Nonlinear Binary Response Models

In this class of models we change the way we model the response probability $p\left(x\right)$. Instead of the simple linear structure from above, we write

$$Pr(y=1|x)=p\left(x\right)=G({\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_K{x}_K)$$

• almost identical to LPM!

• except the linear index ${\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_K{x}_K$ is now inside some function $G(\cdot )$.

• Main property of $G$: transforms any $z\in R$ into a number in the interval $(0,1)$.

• This immediately solves our problem of getting weird predictions for probabilities.

$G$: probit and logit

Running probit and logit in R: the glm function

• We use the glm function to run a generalized linear model

• This generalizes our standard linear model. We have to specify a family and a link:

probit <- glm(inlf ~ age, 
                    data = mroz, 
                    family = binomial(link = "probit"))
logit <- glm(inlf ~ age, 
                    data = mroz, 
                    family = binomial(link = "logit"))

Interpretation

modelsummary::modelsummary(list("probit" = probit,"logit" = logit))

Interpretation

The model is

$$Pr(y=1|\text{age})=G\left(x\beta \right)=G({\beta }_0+{\beta }_1\text{age})$$ and the marginal effect of age on the response probability is

$$\frac{\partial Pr(y=1|\text{age})}{\partial \text{age}}=g({\beta }_0+{\beta }_1\text{age}){\beta }_1$$

• function $g$ is defined as $g\left(z\right)=\frac{dG}{dz}\left(z\right)$ - the first derivative function of $G$ (i.e. the slope of $G$).

• given $G$ that is nonlinear, this means that $g$ will be non-constant. You are able to try this out yourself using this app here:

ScPoApps::launchApp("marginal_effects_of_logit_probit")

or online

Interpretation

So you can see that there is not one single marginal effect in those models, as that depends on where we evaluate the previous expression. In practice, there are two common approaches:

1. report effect at the average values of $x$: $$g\left(\overline{x}\beta \right){\beta }_j$$

2. report the sample average of all marginal effects: $$\frac{1}{n}\sum _{i=1}^{N}g\left(x_i\beta \right){\beta }_j$$

Thankfully there are packages available that help us to compute those marginal effects fairly easily. One of them is called mfx, and we would use it as follows:

Interpretation

f <- "inlf ~ age + kidslt6 + nwifeinc" # setup a formula
glms <- list()
glms$probit <- glm(formula = f, 
                    data = mroz, 
                    family = binomial(link = "probit"))
glms$logit <- glm(formula = f, 
                    data = mroz, 
                    family = binomial(link = "logit"))
# now the marginal effects versions
glms$probitMean <- mfx::probitmfx(formula = f, 
                    data = mroz, atmean = TRUE)
glms$probitAvg <- mfx::probitmfx(formula = f, 
                    data = mroz, atmean = FALSE)
glms$logitMean <- mfx::logitmfx(formula = f, 
                    data = mroz, atmean = TRUE)
glms$logitAvg <- mfx::logitmfx(formula = f, 
                    data = mroz, atmean = FALSE)

Interpretation

GOF in Binary Models

• There is no universally accepted $R^2$ for binary models.

• We can think of a pseudo $R^2$ which compares our model to one without any regressors:

glms$probit0 <- update(glms$probit, formula = . ~ 1)  # intercept model only
1 - as.vector(logLik(glms$probit)/logLik(glms$probit0))

## [1] 0.07084972

• But that's not super informative (unlike the standard $R^2$). Changes in likelihood value are highly non-linear, so that's not great.

• Let's check accuracy - what's the proportion correctly predicted! round(fitted(x)) assigns 1 if the predicted prob $>0.5$.

prop.table(table(true = mroz$inlf, pred = round(fitted(glms$probit))))

##     pred
## true         0         1
##    0 0.1699867 0.2616202
##    1 0.1221780 0.4462151

GOF in Binary Models: ROC Curves

• The 0.5 cutoff is arbitrary. What if all predicted probs are $>0.5$ but in the data there are about 50% of zeros?

• Let's choose an arbitrary cutoff $c\in (0,1)$ and check accuracy for each value. This gives a better overview.

• Also, we can confront the true positives rate (TPR) with the false positives rate (FPR).

  1. TPR: number of women correctly predicted to work divided by num of working women.
  2. FPR: number of women incorrectly predicted to work divided by num of non-working women.

• Plotting FPR vs TPR for each $c$ defines the ROC (Receiver Operating Characteristics) Curve.

• A good model has a ROC curve in the upper left corner: FPR = 0, TPR = 1.

GOF in Binary Models: ROC Curves

library(ROCR)
pred <- prediction(fitted(glms$probit), mroz$inlf)
par(mfrow = c(1,2), mar = lowtop)
plot(performance(pred,"acc"))
plot(performance(pred,"tpr","fpr"))
abline(0,1,lty = 2, col = "red")