ANOVA Models with Logistic Regression

The UC Berkeley data set is built into R already in the UCBAdmission data set!

library(tidyverse); library(broom); library(car)
theme_set(theme_bw())

# Creating the Berkeley Data set:
UCB <- 
  UCBAdmissions |> 
  data.frame() |> 
  pivot_wider(
    names_from = "Admit",
    values_from = "Freq"
  ) |> 
  mutate(
    applied = Admitted + Rejected,
    accept_rate = Admitted/applied
  ) |> 
  janitor::clean_names() |> 
  arrange(dept, gender) |> 
  # Have the results match the slides
  mutate(gender = (gender == 'Female')*1)

# Just looking at the three departments in the slides:
UCB3 <- UCB |> filter(dept %in% c("A", "B", "C")) 

Let’s look at the UCB data set:

UCB3

## # A tibble: 6 × 6
##   gender dept  admitted rejected applied accept_rate
##    <dbl> <fct>    <dbl>    <dbl>   <dbl>       <dbl>
## 1      0 A          512      313     825       0.621
## 2      1 A           89       19     108       0.824
## 3      0 B          353      207     560       0.630
## 4      1 B           17        8      25       0.68 
## 5      0 C          120      205     325       0.369
## 6      1 C          202      391     593       0.341

So we want to see how gender (\(X\)) and department (\(Z\)) change the probability someone is accepted (\(Y\))

We have two main logistic models we can use that includes all 3:

1. Additive Model: \(\log \left ( \frac{\pi_{ik}}{1-\pi_{ik}} \right) = \alpha + \beta^X_i +\beta^Z_k\)

2. Interaction Model: \(\log \left ( \frac{\pi_{ik}}{1-\pi_i} \right) = \alpha + \beta^X_i +\beta^Z_k + \beta^{XZ}_{ik}\)

   • This is our Saturated Model if \(X\) and \(Z\) are factors

Comparing the additive vs interaction model

Let’s start by fitting the two models:

# Additive Model
berk_add <- 
  glm(
    formula = accept_rate ~ gender + dept,
    weight = applied,
    family = binomial,
    data = UCB3
  )

tidy(berk_add, exponentiate = F)

## # A tibble: 4 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   0.571     0.0696     8.21  2.26e-16
## 2 gender        0.197     0.117      1.68  9.23e- 2
## 3 deptB        -0.0366    0.110     -0.333 7.40e- 1
## 4 deptC        -1.32      0.117    -11.3   1.86e-29

# Interaction model
berk_int <- 
  glm(
    formula = accept_rate ~ gender * dept,
    weight = applied,
    family = binomial,
    data = UCB3
  )

tidy(berk_int,
     exponentiate = F)

## # A tibble: 6 × 5
##   term         estimate std.error statistic  p.value
##   <chr>           <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    0.492     0.0717     6.86  6.94e-12
## 2 gender         1.05      0.263      4.00  6.21e- 5
## 3 deptB          0.0416    0.113      0.368 7.13e- 1
## 4 deptC         -1.03      0.135     -7.58  3.34e-14
## 5 gender:deptB  -0.832     0.510     -1.63  1.03e- 1
## 6 gender:deptC  -1.18      0.300     -3.93  8.53e- 5

If there is only 1 interaction term (aka both \(X\) and \(Z\) are binary) we can just conduct a hypothesis test to see if that term is 0 using tidy().

However, if at least one of \(X\) or \(Z\) are not binary, we’ll have more than 1 model term, so we can’t conduct a simple z-test :(

What we can do is compare the deviances of the two models:

fit_stats <- 
  bind_rows(
    "add" = glance(berk_add),
    "int" = glance(berk_int),
    .id = "model"
  ) 

fit_stats |> 
  dplyr::select(model, deviance, df.residual, nobs) |> 
  mutate(deviance = round(deviance, 3))

## # A tibble: 2 × 4
##   model deviance df.residual  nobs
##   <chr>    <dbl>       <int> <int>
## 1 add       17.2           2     6
## 2 int        0             0     6

A difference in the deviances will follow a \(\chi^2\) model with a difference in the number of model parameters

The difference in parameters between an additive and interaction model is the number of interaction terms that we have: df = \((I - 1)(K - 1)\)

For our example, \(I = 2\), \(K = 3\) then df = 2, which is also the difference in the df.redisual from glance()

fit_stats |> 
  summarize(
    comparison = paste(model,  collapse = ' vs '),
    G = deviance |> diff() |> abs(),
    df = df.residual |> diff() |> abs()
  ) |> 
  mutate(
    p_val = pchisq(q = G,df = df, lower = F)
  )

## # A tibble: 1 × 4
##   comparison     G    df    p_val
##   <chr>      <dbl> <int>    <dbl>
## 1 add vs int  17.2     2 0.000184

So we have strong evidence that the interaction term is needed :(

Instead of finding the test stat and p-value by hand, we can use Anova() in the car package with type = 3

car::Anova(mod = berk_int, type = 3)

## Analysis of Deviance Table (Type III tests)
## 
## Response: accept_rate
##             LR Chisq Df Pr(>Chisq)    
## gender        19.054  1  1.271e-05 ***
## dept          69.791  2  6.999e-16 ***
## gender:dept   17.197  2  0.0001844 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The interaction row of the results match exactly what we have!

If you have two or more nested models that you want to compare directly, we can use the anova() function in base R.

Let’s fit a model that just has gender and another model that just has department:

berk_gen <- 
  glm(
    formula = accept_rate ~ gender,
    weight = applied,
    family = binomial,
    data = UCB3
  )


berk_dept <- 
  glm(
    formula = accept_rate ~ dept,
    weight = applied,
    family = binomial,
    data = UCB3
  )

Now we can use anova() by giving it multiple models from simplest to the most complicated:

options(width = 400)
anova(berk_gen, berk_add, berk_int, 
      test = "LRT") |> tidy()

## # A tibble: 3 × 6
##   term                        df.residual residual.deviance    df deviance   p.value
##   <chr>                             <dbl>             <dbl> <dbl>    <dbl>     <dbl>
## 1 accept_rate ~ gender                  4          1.67e+ 2    NA     NA   NA       
## 2 accept_rate ~ gender + dept           2          1.72e+ 1     2    150.   3.08e-33
## 3 accept_rate ~ gender * dept           0          1.79e-13     2     17.2  1.84e- 4

# Can we drop Gender entirely?
# That is, are Admission and Gender independent accounting for Department?
anova(berk_gen, berk_int, 
      test = "LRT")|> tidy()

## # A tibble: 2 × 6
##   term                        df.residual residual.deviance    df deviance   p.value
##   <chr>                             <dbl>             <dbl> <dbl>    <dbl>     <dbl>
## 1 accept_rate ~ gender                  4          1.67e+ 2    NA      NA  NA       
## 2 accept_rate ~ gender * dept           0          1.79e-13     4     167.  4.79e-35

# Can we just drop the interaction term?
anova(berk_dept, berk_int, 
      test = "LRT")|> tidy()

## # A tibble: 2 × 6
##   term                        df.residual residual.deviance    df deviance   p.value
##   <chr>                             <dbl>             <dbl> <dbl>    <dbl>     <dbl>
## 1 accept_rate ~ dept                    3          2.01e+ 1    NA     NA   NA       
## 2 accept_rate ~ gender * dept           0          1.79e-13     3     20.1  0.000165

anova() compares the model on a row with the in the previous row

So the test results for row 2 compares gender only vs the additive model and row 3 compares the additive model vs the interaction model

Testing Departments B - F

We saw in chapter 2 that department A was the problematic department with the odds a man is accepted compared to a woman being around 1/3.

The code chunk below removes department A. Test to see if gender is needed when department is in the model. That is, perform a test to see if Gender and Admission are conditional independent on Department!

UCB5 <- 
  UCB |> 
  filter(dept != 'A')

tibble(UCB5)

## # A tibble: 10 × 6
##    gender dept  admitted rejected applied accept_rate
##     <dbl> <fct>    <dbl>    <dbl>   <dbl>       <dbl>
##  1      0 B          353      207     560      0.630 
##  2      1 B           17        8      25      0.68  
##  3      0 C          120      205     325      0.369 
##  4      1 C          202      391     593      0.341 
##  5      0 D          138      279     417      0.331 
##  6      1 D          131      244     375      0.349 
##  7      0 E           53      138     191      0.277 
##  8      1 E           94      299     393      0.239 
##  9      0 F           22      351     373      0.0590
## 10      1 F           24      317     341      0.0704