What are likelihood functions, again?

Likelihood functions are essentially probability models of our data. They provide a mathematical model of the process that generates our data. Remember, the goal of statistical modeling is to build models that can make accurate predictions and describe data! We do this by choosing likelihood functions that describe observed data.

What if the normal distribution doesn’t describe our data?

As an example, we’re going to use data on flights in the US:

head(flights)

## # A tibble: 6 × 23
##    year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##   <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
## 1  2013    12     3     1809           1815        -6     2051           2127
## 2  2013     9    27     1516           1520        -4     1739           1740
## 3  2013     5     2     1350           1330        20     1645           1622
## 4  2013     9    21     1434           1420        14     1555           1553
## 5  2013     8    23     1939           1830        69     2238           2126
## 6  2013    10     4     1900           1900         0     2248           2235
## # … with 15 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>, flight_delayed <chr>,
## #   delayed <dbl>, distance_c <dbl>, distance_2sd <dbl>

Let’s try to see if the distance of a flight can predict whether the flight is delayed or not. What probability model would describe flight delays? The normal distribution?

ggplot(flights, aes(x = delayed)) +
  geom_density()

The normal distribution doesn’t look right to me, not to mention the fact that the normal distribution describes a data generating process that produces continuous values from \(-\infty\) to \(\infty\). Let’s run with it anyway

normal_fit <- glm(delayed ~ distance_2sd,
                  family = gaussian(),
                  data = flights)

library(bayesplot)
pp_check(normal_fit)

The normal distribution doesn’t model the data very well. How can we tell? Well, in the plot about we see the distribution of the outcome in the darker line labeled \(y\) and simulated draws from the model in the light blue lines labeled \(y_{rep}\). If our model was capturing the data generating process, the simulated draws \(y_{rep}\) would match the actual data \(y\).

What we need is a likelihood function that describes a data generating process that produces discrete values of either 0 or 1. That is exactly what the binomial and Bernoulli distributions describe. Let’s write out that model:

\[ \begin{aligned} y_i &\sim \text{Binomial}(1, p_i) \\ \text{logit}(p_i) &= \alpha + \beta \times x_i \end{aligned} \]

This formulation shows the likelihood and link functions.

Wait, what is a link function?

Link functions are used to map the range of a linear equation to the range of the outcome variable. That means that even if we model the outcome as a linear combination \(\mu = \alpha + \beta * x\), it result of this equation will always be restricted by the link function that we use. For example, the logit link function restricts the range of the linear equation to fall between 0 and 1 and the log link function restricts the range of the linear equation to be positive, continuous, and greater than 0.

Logit Link Function

logit_fit1 <- glm(Survived ~ Class,
                  family = binomial(link = "logit"),
                  data = titanic_data)
tidy(logit_fit1, conf.int = TRUE)

## # A tibble: 4 × 7
##   term         estimate std.error statistic p.value conf.low conf.high
##   <chr>           <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)  6.62e-16     0.707  9.36e-16    1.00    -1.44      1.44
## 2 Class2nd    -5.80e-16     1     -5.80e-16    1       -2.00      2.00
## 3 Class3rd    -8.81e-16     1     -8.81e-16    1.00    -2.00      2.00
## 4 ClassCrew   -7.21e-16     1     -7.21e-16    1.00    -2.00      2.00

The model ran, but the results look suspicious. Z-scores of 0? P-values of 1? Let’s inspect the data to see if we can figure out the problem.

glimpse(titanic_data)

## Rows: 32
## Columns: 5
## $ Class    <fct> 1st, 2nd, 3rd, Crew, 1st, 2nd, 3rd, Crew, 1st, 2nd, 3rd, Crew…
## $ Sex      <fct> Male, Male, Male, Male, Female, Female, Female, Female, Male,…
## $ Age      <fct> Child, Child, Child, Child, Child, Child, Child, Child, Adult…
## $ Survived <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, N…
## $ Freq     <dbl> 0, 0, 35, 0, 0, 0, 17, 0, 118, 154, 387, 670, 4, 13, 89, 3, 5…

This data is aggregated! Remember how we wrote of the equation for the binomial regression?

\[ \begin{aligned} y_i &\sim \text{Binomial}(1, p_i) \\ \text{logit}(p_i) &= \alpha + \beta \times x \end{aligned} \]

The \(1\) governs the number of trials for each event. Usually, we just set this value equal to \(1\) because each row of our data is one event (e.g. the person voted or they didn’t). If we wanted to model binomial data where each event only has \(1\) trial it would be better to use the Bernoulli distribution which assumes that the number of trials is \(1\). It would be written like this:

\[ \begin{aligned} y_i &\sim \text{Bernoulli}(p_i) \\ \text{logit}(p_i) &= \alpha + \beta \times x \end{aligned} \]

A regression with a Bernoulli likelihood can be implemented in a package like brms. For the aggregated data, we need to set the number of trials for each event and we do this with the Freq variable. The easiest way to do this with our code is to set the weight = argument equal to Freq:

logit_fit2 <- glm(Survived ~ Class,
                  family = binomial(link = "logit"),
                  data = titanic_data,
                  weights = Freq)
tidy(logit_fit2, conf.int = TRUE)

## # A tibble: 4 × 7
##   term        estimate std.error statistic  p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)    0.509     0.115      4.44 8.79e- 6    0.287     0.736
## 2 Class2nd      -0.856     0.166     -5.16 2.51e- 7   -1.18     -0.533
## 3 Class3rd      -1.60      0.144    -11.1  1.07e-28   -1.88     -1.32 
## 4 ClassCrew     -1.66      0.139    -12.0  4.97e-33   -1.94     -1.39

Now that looks a lot better.

Binomial Exercises

Now you practice. Run the code block below to disaggregate the Titanic data:

titanic_data_2 <- 
  titanic_data %>% 
  mutate(New_Response = map(Freq, ~ rep_len(1, .x))) %>% 
  unnest(cols = c(New_Response)) %>% 
  select(-Freq, -New_Response)

With the data disaggregated, try answering the following questions:

1. Run a standard logistic regression model with the disaggregated data (binomial likelihood function, logit link function).

2. Run a probit regression with the disaggregated data (binomial likelihood function, probit link function).

Doesn’t everyone use Negative-Binomial instead of Poisson?

Usually, where we see someone using a count model, they use a negative binomial regression instead of a Poisson regression. Well, it turns out that the negative Binomial model doesn’t involve a Binomial distribution. It’s much more accurate to call it a Gamma-Poisson model. Why Gamma-Poisson and why do we use it? One limitation of the Poisson distribution is that it assumes that the mean is equal to the variance. When this isn’t the case, we have either over, or under, dispersion. This is where the Gamma-Poisson model helps and we can write it like this:

\[ \begin{aligned} y_i &\sim \text{Gamma-Poisson}(\lambda_i, \phi) \\ \text{log}(\lambda_i) &= \alpha + \beta \times x \\ \phi &\sim \text{Gamma}(0.01, 0.01) \end{aligned} \]

Instead of assuming the mean and variance are equal, the Gamma-Poisson model has an additional parameter (\(\phi\)) that allows each observation to have its own rate (it lets the mean differ from the variance). And, the rate parameter follows a gamma distribution because the gamma distribution is a positive continuous distribution defined between (0 and \(\infty\)). Let’s run it.

#install.packages("MASS")
library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

gp_fit1 <- glm.nb(lax_tbit_main_checkpoint ~ lockdown,
                  data = tsa_data)

tidy(gp_fit1, conf.int = TRUE)

## # A tibble: 2 × 7
##   term        estimate std.error statistic     p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>       <dbl>    <dbl>     <dbl>
## 1 (Intercept)    6.44    0.00685    941.   0              6.43      6.45 
## 2 lockdown      -0.220   0.0428      -5.13 0.000000288   -0.302    -0.135

Because the Gamma-Poisson model involves two likelihood functions, it is considered a mixture model. Mixture models are models that combine two likelihood functions in order to describe a data generating process.

So, how do we interpret this model?

(1 - exp(coef(gp_fit1)[[2]])) * 100

## [1] 19.72062

Air travel decreased by about 20 in March of 2020.