---
title: "Everyone Watches Women's Basketball: Attendance and Fan Engagement in the WNBA"
author:
  - name: Gabriela Eastwood
    corresponding: true
    email: geastwood@uchicago.edu
    orcid: 0009-0004-7969-198X
    affiliations:
      - ref: chicago
      - ref: smith
      - id: uchicago
        name: University of Chicago
        department: Harris School of Public Policy
        address: Chicago, IL
  - name: 
      given: Benjamin S. 
      family: Baumer
    orcid: 0000-0002-3279-0516
    email: bbaumer@smith.edu
    url: https://www.smith.edu/people/ben-baumer
    corresponding: false
    affiliation:
      - ref: smith
  - name: 
      given: Andrew 
      family: Zimbalist
    affiliation:
      - ref: smith
      - department: Economics
    email: azimbali@smith.edu
    url: https://www.smith.edu/people/andrew-zimbalist
    corresponding: false
affiliations:
  - id: smith
    name: Smith College
    department: Statistical & Data Sciences
    address: 44 College Lane
    city: Northampton, MA
    country: United States of America
    postal-code: 01063
    url: https://www.smith.edu/
    isni: 0000000119314868
  - id: chicago
    name: University of Chicago
    department: Harris School of Public Policy
    address: 1307 E. 60th Street
    city: Chicago, IL
    postal-code: 60637
    country: United States of America
    url: https://harris.uchicago.edu/


keywords:
  - WNBA
  - Game Attendance
abstract: |
  We examine how the relationships between attendance and three notions of fan engagement have changed over the history of the WNBA. Previous research suggests that stronger professional sports teams draw higher live game attendance. However, there is ongoing doubt about the extent to which that relationship holds in the WNBA. We improve upon previous research by modeling attendance in the broader context of the league's growth from 1997-2025, and find that the factors that drive fan demand have evolved as the league has matured. We test whether three leading theories of fan behavior---quality of play, loss aversion, and outcome uncertainty---have become more important drivers of attendance as the league has grown. We use game-level attendance data reported directly by the WNBA as well as team Elo ratings and game-by-game home win probabilities and find that in recent years fans are increasingly likely to attend games in response to quality of play and loss aversion factors. We also find strong support for the large impact of Indiana Fever guard Caitlin Clark's presence upon attendance. Our results suggest that the WNBA is entering a new era in which fans are choosing to attend games to watch good basketball rather than novelty-driven interest or broader support for women's sports. This, in turn, is consistent with a more durable popularity for the league.
  
  
key-points:
  - winning percentage and game attendance
date: last-modified
bibliography: [pkgs.bib, references.bib]
number-sections: true
---

```{r}
#| message: false
library(tidyverse)
library(knitr)
library(gt)

pkg <- function(x, ...) {
#  if (knitr::is_html_output()) {
    paste0("**", x, "** [@R-", x, "]")
#  } else {
    # isn't working:
#    https://github.com/beanumber/datascience2/issues/7
#    paste0(r"(\pkg{)", x, r"(\})")
#    paste0("**", x, "**")
#  }
}
```

```{r}
#| include: false
manual <- c("ggplot2", "tidycensus", "sandwich", "elo", "VGAM")
pkgs <- unique(c(manual, .packages()))
knitr::write_bib(pkgs, file = "pkgs.bib")
```


```{r}
wnba_gl <- read_rds(here::here("data/wnba_gl.rds")) |>
  group_by(arena) |>
  mutate(
    median_income = median_income / 1000,
    metro_pop = total_pop / 1000000,
    est_capacity = max(max(attendance), arena_capacity),
    pct_capacity = attendance / est_capacity,
    pct_capacity_safe = case_when(
      pct_capacity == 1 ~ 0.999, 
      pct_capacity == 0 ~ 0.001,
      TRUE ~ pct_capacity
    ),
    arena_size = case_when(
      arena_capacity >= 16000 ~ "arena_large",
      arena_capacity >= 11000 ~ "arena_medium",
      arena_capacity >= 5500 ~ "arena_small",
      TRUE ~ "arena_unknown"
    )
  ) |>
  ungroup() |>
  mutate(
    is_sellout = pct_capacity > 0.92,
    quality_of_play = elo_home_pregame * elo_away_pregame,
    uncertain_outcome = -sqrt(abs(elo_home_pregame - elo_away_pregame)),
    loss_aversion = elo_home_win_prob^2,
    z1_quality_of_play = (quality_of_play - mean(quality_of_play)) / sd(quality_of_play), 
    z2_uncertain_outcome = (uncertain_outcome - mean(uncertain_outcome)) / sd(uncertain_outcome),
    z3_loss_aversion = (loss_aversion - mean(loss_aversion)) / sd(loss_aversion),
    zzz_sum_all_three = z3_loss_aversion + z1_quality_of_play + z2_uncertain_outcome,
  ) |>
  ungroup()
```

```{r}
#| eval: false
ggplot(wnba_gl, aes(x = uncertain_outcome)) +
  geom_histogram() +
  scale_color_viridis_c()
ggplot(wnba_gl, aes(x = z2_uncertain_outcome)) +
  geom_histogram() +
  scale_color_viridis_c()
ggplot(wnba_gl, aes(x = elo_home_pregame, y = elo_away_pregame, color = z2_uncertain_outcome)) +
  geom_point() +
  scale_color_viridis_c()

ggplot(wnba_gl, aes(x = pct_capacity, y = pct_capacity_safe)) +
  geom_point()

wnba_gl |>
  filter(pct_capacity != pct_capacity_safe)
```



## Introduction {#sec-intro}

In the years prior to the 2020 WNBA collective bargaining agreement (CBA), NBA Commissioner Adam Silver publicly characterized the women’s league as unpopular and unprofitable [@garcia2020betting]. 
Again during the CBA renegotiations in 2025, NBA and WNBA leadership made similar negative claims about the women’s league’s long term viability. 
While leadership has been ironically pessimistic about the health of the WNBA in recent years, their claims stand in tension with recent trends: the league has enjoyed a significant increase in attendance and fan interest in the past several seasons.

While this trend is encouraging for anyone interested in the success of the WNBA, further breaking down the underlying components of this trend is necessary both to explain and sustain this growth. 
For decades, league officials and sports economists have debated the primary drivers of game attendance, identifying game-related and market-related factors impacting fan interest and game attendance. 
Much of this research, beginning with foundational work in the 1970s, assumes that fans are more likely to attend games when the home team is expected to do well. 
This relationship is widely tested and supported in men’s professional sports, but women’s leagues and the WNBA in particular are still largely underresearched.
As a relatively young league with rapid growth in viewership [@wnba2024recordseason] but limited publicly available financial data, game attendance remains the most readily available proxy for public interest in the WNBA. 
In this paper we set out to determine how the underlying causes shape game attendance and explain these more recent trends. 
Specifically, we focus on whether on-court factors are becoming more important to WNBA game attendance in recent years.

### Existing theories of fan behavior

A key line of research in sports economics links fan demand to expectations about game outcomes. 
Early work such as the uncertainty-of-outcome hypothesis suggests that attendance is higher when game results are less predictable, increasing perceived excitement [@rottenberg1956baseball; @neale1964peculiar]. @CoatesHumphreysZhou2014 use more recent behavioral economics models to suggest that people evaluate gains and losses relative to a reference point and place greater weight on avoiding expected losses than on achieving comparable gains. Rather than predicting a simple linear relationship between attendance and win probability, their model theorizes that attendance relates non-linearly to win probability. 
They argue that applying this idea of loss aversion may better explain why fans choose to go to games than the uncertainty of outcome hypothesis. 
They use MLB live-game attendance data to find patterns consistent with the hypothesis that fans are less likely to attend games in which they expect the home team to lose [@CoatesHumphreysZhou2014]. 
These frameworks together suggest that attendance decisions are sensitive not only to team quality, but to expectations about competitive performance and likely game outcomes.

In the first study in sports economics that sought to empirically identify the factors that affected attendance, @noll1974attendance considered data from MLB, the NFL, NHL, NBA and ABA (American Basketball Association) for the 1969-70 and 1970-71 seasons and found that demand was a function of population, interactions of income with star players, team win percentage, and per capita income. 
Multiple subsequent studies on men’s professional team sports continued to find a statistically significant correlation between team quality and attendance.
Differences in the literature over how to measure team quality (wins, wins lagged, number of championships in past years, win probability, star power), attendance (logs, lags, linear regression or censored), the inclusion of control variables and nonlinear forms also emerged [@wiseman2003team; @woltring2018attendance; @kim2023winning; @deschriver2002determinants; @paul2019attendance].

Studies of attendance and winning percentage most often theorize the latter as causing the former, although some work has considered a bidirectional causal relationship. 
@davis2008interaction and @horowitz2007if both found that attendance and winning percentage were related. 
Horowitz concluded a bidirectional relationship between the two, while Davis suggested that causation runs from winning percentage to attendance. More broadly, sports economists have shown that causal relationships involving team performance may operate in both directions.  
In a following study, @davis2009analyzing found winning percentage as a significant determinant of attendance for 12 National League teams. 
The study did not establish causation both ways. 
@lemke2010estimating found attendance increases as home team win probability increases. 
The bulk of the literature has found that the causal direction runs from winning percentage to attendance[^causal].

[^causal]: For a discussion of modeling causal direction for time series data in the sports context,
see @hall2002testing.

While teams and fans want to win as many games as possible, significant imbalance in winning percentages across a professional sports league will tend to decrease viewership as spectators will only watch a game if they have some uncertainty about the outcome [@PerlineNobleStoldt2018]. 
In NCAA basketball, @PerlineNobleStoldt2018 examine the standard deviation and the range of winning percentages for teams in the women’s Division I and determine that the women’s game is significantly less balanced than the men’s game. 
They hypothesize that the difference in the competitive balance of the men’s and women’s games could be one cause of lower revenue generation and attendance in the women’s game. 

### Why is the WNBA different?

@AghaBerri2023 argue that a major reason why attendance is considerably lower at WNBA games relative to NBA games is that the women’s league is still in its early years, only beginning play in 1997 compared to the NBA’s 1946. 
To account for this, they compare attendance in the tenth season of the WNBA (2006) and the NBA (1955/56) (p. 36). 
In their respective tenth years, the WNBA average attendance per home game exceeded the NBA’s by more than 2,000 fans per game. 
They do not control for the growth in population over the period 1955-2006. 
This comparison suggests that the current disparity in attendance between the women’s and men’s league should not be interpreted simply as a lack of demand, but rather a sign of the league’s maturity. 
Considering the league’s age as a factor in demand for attendance, we inquire further into whether patterns of fan behavior change as the league matures. 
@AghaBerri2023 attempt to establish whether different forces shape demand for attendance in the WNBA and the NBA, finding that team performance was not a statistically significant predictor of game attendance in the WNBA. 
As noted above, this is contrary to the virtually unanimous finding in the existing scholarship for men’s leagues. 
However, their study only employed data from 2004 through 2017 and thus cannot capture the more recent improvements in attendance shown in @fig-boxplots. 
They also measure a team's strength using the team's winning percentage from the previous year which raises issues because fans are more likely to respond to the more recent performance of the team. 
Thus, while a team’s win percentage in 2016 might have an appreciable impact on the team’s attendance at the beginning of the 2017 season, as the season wears on, this impact should diminish. 

The disparity between raw attendance numbers and attendance driving factors in the NBA and WNBA cannot be attributed to the age of the leagues alone. 
Additional important factors include: the WNBA plays during the summer months in order to promote arena use outside of the NBA season; there is less money spent on promotion; reduced media coverage; historical prejudice; poorer arena amenities and age, among others. 

Sport management literature also includes several studies that find a major motivation behind fandom for women’s sports leagues to be an area’s overall support for feminist issues, suggesting that fans in those areas may attend games out of political support rather than to see good basketball [@delia2020psychological; @guest2018fan; @leslie2022football]. 
Factors like the summer season, reduced media coverage, and the association of game attendance with feminist causes all contribute to the argument that the WNBA functions fundamentally differently to the NBA---at least in terms of attendance drivers. 
This conception of the WNBA as fundamentally different from the NBA is dangerous both to the league's reputation and future growth, motivating our more thorough examination of game attendance in the context of the league's history.

### Our contribution

It is reasonable to anticipate that as women’s sports mature and are perceived increasingly as a normal sports league (rather than a political concession or oddity) they will begin to assert similar fan demand characteristics. 
In particular, fan demand should become more responsive to team quality. 
Using game attendance and team performance data from 1997 through 2025, we test these hypotheses. 
In viewing more recent attendance figures for the WNBA in the context of the league’s history, we determine whether predictors of attendance have shifted over time. 
Specifically, we find that in recent years the strength of the effect of a game's uncertainty, the strength of both teams, and the likelihood of a home team win, on game attendance has increased, paralleling documented trends in men’s leagues. 


## Data {#sec-data}

### Game logs

[^attendance]: Note that, as in most sports leagues, there is always the possibility of finagling with game attendance figures. The WNBA appears to give attendance counting instructions to teams that define attendance as tickets out, including total ticket sales, plus complimentary tickets and promotions. There is considerable room for these factors to be manipulated, rendering the possibility of noise in the attendance figures we use.

Our analysis is based on game-level data that spans WNBA history, beginning in 1997 and continuing through the end of the 2025 season. 
Our final data set contains observations from `r format(nrow(wnba_gl), big.mark = ",")` games.

First, we collected game-level attendance data from an online repository for women's basketball data that collects attendance numbers[^attendance] directly from the league's records [@AcrossTheTimeline]. 
Attendance data was unavailable for only `r nrow(filter(wnba_gl, is.na(attendance)))` games, while (due to the COVID-19 pandemic) `r nrow(filter(wnba_gl, attendance == 0))` games (accurately) recorded attendance of 0. 
The arenas in which WNBA teams play vary widely in size.
For example, the New York Liberty alone have called both the Barclays Center (17,732 seat capacity) and the Westchester County Center (2,300 seat capacity) home in recent years. 
While teams have many reasons they might choose to play for a sold-out crowd in a smaller arena, smaller venue size will also by nature constrain attendance, making its consideration in our model important.

Second, we collected game-level performance data from the **wehoop** R package [@hutchinson_gilani_2021_wehoop], which uses the official WNBA API. 
These data include the final score of the game, among other performance metrics. 
We also used the **wehoop** package to identify whether Caitlin Clark and/or A'ja Wilson played more than 0 minutes in the game. 

@tbl-attendance shows a sample of the data. 
While many other variables are omitted from the table, we note that the attendance, the teams playing in the game, the outcome, the arena, and whether Clark or Wilson played in the game are recorded. 

```{r}
#| label: tbl-attendance
#| fig-cap: "Sample of Game-level WNBA data. The data indicates (among other things) which teams were playing, who won, the arena, the attendance, and whether Caitlin Clark and/or A'ja Wilson played in the game."

my_gt <- function(x) {
  x |>
    gt() |>
    tab_options(
      table.font.size = px(11),
      table.width = pct(100)
    )
}

wnba_gl |>
  select(
    game_date,
    home = team_abbreviation,
    away = opponent,
    arena,
    attendance,
    plus_minus,
    clark = is_cc,
    wilson = is_aw
  ) |>
  slice_sample(n = 10) |>
  mutate(game_date = as.character(game_date)) |>
  janitor::clean_names(case = "title") |>
  my_gt() |>
  tab_header(
    title = "WNBA Game Data (1997-2025)",
    subtitle = "Sample of game-level records"
  ) |>
  fmt_number(
    columns = Attendance,
    decimals = 0
  )
```

### Attendance

@fig-boxplots shows the distribution of attendance across all seasons, with a logarithmic scale on the vertical axis. 
There is no attendance data from 2020, because all games that year were played in the WNBA bubble due to the COVID-19 pandemic, and some games in 2021 were played with artificially low attendance.  

```{r}
wnba_summary <- wnba_gl |>
  group_by(season_id) |>
  summarize(
    num_games = n(),
    avg_attendance = mean(attendance)
  )

attend_hi <- max(wnba_summary$avg_attendance, na.rm = TRUE)
attend_prev_hi <- min(tail(sort(wnba_summary$avg_attendance), 2))
attend_pct_increase <- 100 * (attend_hi / attend_prev_hi - 1)
```


Several broad observations can be made from @fig-boxplots. 
First, after some initial enthusiasm during the first three years of the league, attendance largely plateaued or steadily but slightly declined from 2000 to 2017. 
(Note that @AghaBerri2023 used only data through 2017.) 
The two years preceding the pandemic (2018-2019) had notably lower attendance. 
However, since the pandemic, attendance has rebounded, exceeding its initial levels by 2024 and reaching new heights in 2025. 
In fact, average WNBA attendance in 2025 was `r format(round(attend_hi), scientific = FALSE, big.mark = ",")`, which was `r round(attend_pct_increase, 1)`% higher than the previous high in 1998. 

These last two years coincide with the arrival of Caitlin Clark, the 2024 Rookie of the Year Award winner who is widely regarded as a generational talent. 
Disentangling the "Caitlin Clark effect" from the league's overall attendance trends is one of our goals. 
To that end, we identified all games in which Clark and A'ja Wilson played. 
Wilson is a four-time WNBA MVP and three-time league champion who is currently viewed as the league's best player. 
As a superior player with less celebrity, Wilson's effect on attendance relative to Clark's provides additional benchmarking possibilities. 

```{r}
#| label: fig-boxplots
#| fig-cap: "Distribution of WNBA game attendance by season. Note that the vertical axis has a logarithmic scale. Caitlin Clark and/or A'ja Wilson's games are separated by color. "
plot_attend_eda <- ggplot(wnba_gl, aes(x = factor(season_id), y = attendance)) +
  geom_boxplot(outlier.alpha = 0.2, linewidth = 0.6, fill = "grey90") +
  geom_jitter(alpha = 0.3, width = 0.15, size = 1, aes(color = factor(is_cc | is_aw))) +
  labs(
    title = "Distribution of WNBA Game Attendance by Season",
    subtitle = "1997-2025",
    x = "Season",
    y = "Attendance"
  ) +
  scale_color_discrete(name = "Clark/Wilson") +
  theme_minimal(base_size = 10) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_attend_eda +
  scale_y_log10()
```

### Arena capacity {#sec-arenas}

While we recorded reported arena capacities, those capacities were difficult to obtain for venues not in current use, they change over time, and some games had attendance that exceeded the arena's reported capacity. 
Thus, we estimate each arena's capacity as the larger of the reported capacity, or the maximum recorded attendance at that arena in the sample. 
This ensures that our figures for the percentage of estimated arena capacity are bounded above by 1. 

We designate a sellout as a game in which the attendance exceeds 92% of the arena's estimated capacity.
Using our definition, we deem `r round(100 * nrow(filter(wnba_gl, is_sellout)) / nrow(wnba_gl), 1)`% of the games to be sellouts, with nearly half of those occurring in the past two seasons. 
`r nrow(filter(wnba_gl, is_sellout, is_cc == 1))` of Clark's `r nrow(filter(wnba_gl, is_cc == 1))` games have sold out. 
The overall average percentage of arena capacity is `r wnba_gl$pct_capacity |> mean(na.rm = TRUE) |> round(3)`.

@fig-boxplots-capacity shows the distribution of attendance as a percentage of the arena's estimated capacity. 


```{r}
#| label: fig-boxplots-capacity
#| fig-cap: "Distribution of WNBA game attendance as a percentage of estimated arena capacity by season. Clark and Wilson games are highlighted by color. "
plot_attend_eda +
  aes(y = pct_capacity) +
  scale_y_continuous("Attendance as Percentage of Estimated Arena Capacity")
```


### Estimating team strength {#sec-elo}

We used the game-level results data to model the evolution of the strength of each WNBA team over time. 
Measuring team strength is necessary to test our hypothesis that attendance is in part driven by the desire of fans to see quality basketball. 
We considered season-to-date winning percentage and a rolling 44-game (the length of a typical WNBA season) winning percentage before settling on Elo rating [@Elo1978], because Elo rating updates with each game, and thus more accurately reflects the information that fans have access to in deciding to attend games.
Originally developed for rating chess players, Elo rating is well-understood in sports analytics as a measure of team strength that updates after each game, and---unlike winning percentage and rolling winning percentage---takes into account not only the outcome but also the previous estimated team strength of both teams. 
This means that beating a team with a higher Elo rating has a greater positive effect than beating a team with a lower Elo rating, imbuing outcomes with more nuance than simple win-loss record.

Let $\theta_{i}(t)$ be the strength of team $i$ at time $t$, where Elo ratings are centered at 1500. 
Then if the probability that team $i$ beats team $j$ is modeled by a logistic curve with base 10 and scale factor 400, the win probability is given by:
$$
p_{ij}(t) = \Pr(i \text{ beats }j \text{ at time }t) = \frac{1}{1 + 10^{\frac{\theta_{j}(t) -  \theta_{i}(t)  }{400}}} \,.
$$

Each team $i$ starts the 1997 season with $\hat{\theta}_{i}(0) = 1500$, and following each game, team $i$'s rating is updated via the formula
$$
  \hat{\theta}_{i}(t+1) = \hat{\theta}_{i}(t) + K \cdot \left[ I(i \text{ beats } j)  - \hat{p}_{ij}(t) \right] 
$$
where $I(i \text{ beats } j)$ is an indicator variable that is 1 if $i$ beats $j$ and $-1$ if $j$ beats $i$, and $K = 32$ is a constant that controls how sensitive the ratings are to each game. 
For example, suppose team $i$ has an Elo rating of `r (i <- 1600)` and team $j$ `r (j <- 1450)`. 
Then team $i$  has a `r p <- 1/(1 + 10^((j-i) / 400)); round(p, 3)` probability of beating team $j$ and if they win, their rating will improve by $32 \cdot (1 - \text{`r round(p, 3)`}) = \text{`r round(32 * (1 - p), 2)`}$ points, while team $j$'s rating will decrease by the same amount. 
Note that if team $j$ upsets team $i$, the corresponding changes in rating will be `r round(32 * (1 - (1 - p)), 2)` points.
This larger amount reflects the less likely outcome. 

We fit Elo ratings using the `r pkg("elo")` package for R. 
@fig-elo shows the evolution of Elo rating $\hat{\theta}_i(t)$ by franchise for the 12 teams who were active during the 2024 season. 
The recent poor performance of the Los Angeles Sparks, the 2010 dominance of the Seattle Storm, and the recent rise of the New York Liberty and Indiana Fever are easily visible.

```{r}
#| label: fig-elo
#| fig-cap: "Evolution of team Elo ratings in the WNBA (1997-2025). Only franchises active in 2024 are shown. The recent poor performance of the Los Angeles Sparks, the 2010 dominance of the Seattle Storm, and the recent rise of the New York Liberty and Indiana Fever are easily visible. Note that all franchises begin with a 1500 rating. "
wnba_gl |>
  filter(is_active, team_abbreviation != "GSV") |>
  ggplot(aes(x = game_date, y = elo_home_pregame)) +
  geom_hline(aes(yintercept = 1500), linetype = 3) +
  geom_line() +
  facet_wrap(vars(name_current)) +
  scale_y_continuous(latex2exp::TeX(r"(Elo rating ($\hat{\theta}_i(t)$))")) +
  scale_x_date("Time (t)") +
  theme_minimal(base_size = 10) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
```

### Estimating home win probability

We also added a home court advantage equivalent to 60 points (of Elo rating) to our home win probability model, so that:
$$
  \hat{p}_{ij}(t) = \Pr(i \text{ beats }j \text{ at time }t \text{ at home}) = \frac{1}{1 + 10^{\frac{\hat{\theta}_{j}(t) - \left( \hat{\theta}_{i}(t) + 60 \right)  }{400}}} \,.
$$
We arrived at 60 points by iteratively calibrating the win probabilities to match the observed frequencies by minimizing the Brier score. 
The model fit very well for most of the games (see @tbl-calibrate). 

```{r}
#| label: tbl-calibrate
#| tbl-cap: "Calibration of home win probability model, based on Elo ratings. We note that in the most frequent situations, where the home team was a slight favorite, the observed home win probability is within 3 percentage points of the predicted home probability. However, in those rarer occasions where a home win was extremely likely or unlikely, the model was slightly less accurate. "
wnba_gl |>
  group_by(wpct_bin = cut(elo_home_win_prob, breaks = 10)) |>
  summarize(
    games = n(),
    bin_midpoint = median(elo_home_win_prob),
    obs_home_wpct = mean(win_flag, na.rm = TRUE)
  ) |>
  mutate(diff = obs_home_wpct - bin_midpoint) |>
  janitor::clean_names(case = "title") |>
  my_gt() |>
  tab_header(
    title = "WNBA Home Win Probability Model Calibration (1997-2025)",
    subtitle = "10 Equally-spaced Bins"
  ) |>
  fmt_number(
    decimals = 3
  ) |>
  fmt_number(columns = Games, decimals = 0)
```


## Methods

### Response variables

Our goal is to model the attendance ($y$) at a given WNBA game, as a function of several different measures of fan engagement detailed in @sec-fan. 
We also include several control variables to account for time and city-related effects (see @sec-control).

Since we know the arena in which the game took place and the capacity of that arena, we also developed an alternative response variable: percentage of estimated arena capacity filled ($\gamma$). 

In @fig-attendance-year, we examine the relationship between attendance and estimated home win probability. 
We add a simple quadratic regression line as a flexible approximation of the nonlinear relationship posited by @CoatesHumphreysZhou2014.
In most years, the slope appears to be positive and the curvature concave up. 
This pattern is broadly consistent with @CoatesHumphreysZhou2014's hypothesis that the expected utility of attending a game increases as the likelihood of a home-team victory rises. 
However, these observed relationships do not provide strong evidence that fans primarily seek to avoid losses. 
We do not observe a pronounced decline in attendance for games in which the home team's estimated win probability falls below 50%.

```{r}
#| label: fig-attendance-year
#| fig-cap: "WNBA attendance by season against estimated home win probability. We note that in most years, the slope of the association appears to be positive and the curvature concave up. "
plot_attend <- ggplot(wnba_gl, aes(x = elo_home_win_prob, y = attendance)) +
  geom_vline(aes(xintercept = mean(elo_home_win_prob)), linetype = 3) +
#  geom_hline(aes(yintercept = mean(attendance)), linetype = 3) +
  geom_point(aes(color = factor(is_cc | is_aw))) +
  geom_smooth(method = "lm", formula = y ~ I(x^2), se = 0, color = "black") +
  facet_wrap(vars(season_id)) +
  scale_x_continuous("Estimated Home Win Probability", limits = c(0, 1)) +
  scale_color_discrete(name = "Clark/Wilson") +
  theme_minimal(base_size = 10)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_attend +
  scale_y_log10("Attendance")
```

In @fig-capacity-year, we note similar trends in the relationship between the percentage of the estimated arena capacity filled and the estimated home win probability. 
We also note that since the percentage of arena capacity cannot exceed 1, there are edge effects created by sellouts. 
Those sellouts are much more frequent in the last two years of data, and many of the sellouts involve Caitlin Clark. 

```{r}
#| label: fig-capacity-year
#| fig-cap: "WNBA percentage of estimated arena capacity filled by season, against estimated home win probability. We note the edge effects caused by sellouts, which occur much more frequently in 2024 and 2025. "
plot_capacity <- plot_attend +
  aes(y = pct_capacity) +
  scale_y_continuous("Percentage of Estimated Arena Capacity")
plot_capacity
```

@tbl-sellout summarizes attendance figures across the 10 arenas with the most sellouts. 
We make several observations.
The Golden State Valkyries (a 2025 expansion franchise) sold out all 22 of their games at the Chase Center.
In terms of measuring fan demand, there is a profound difference between selling out a small arena (e.g., the Entertainment and Sports Arena in Washington, D.C. with a capacity of slightly more than 4000) and a large arena (e.g., Gainbridge Fieldhouse in Indianapolis, with a capacity of more than 18000). 
Teams like Washington and Connecticut, in order to generate more excitement in the arena and positive images on telecasts, may have been trying to maximize the percentage of arena capacity at the expense of drawing larger crowds. 
Or perhaps they didn't have a larger arena in which to play. 
Just two of the 27 sellouts at Gainbridge Fieldhouse occurred before Clark's arrival. 

```{r}
#| label: tbl-sellout
#| tbl-cap: "Summary statistics for attendance across the 10 WNBA arenas with the most sellouts. We note the difference between selling out a small arena (e.g., Entertainment and Sports Arena in Washington, D.C.) and a large arena (e.g., Gainbridge Fieldhouse in Indianapolis). "
wnba_gl |>
  group_by(arena) |>
  summarize(
    home_teams = paste(unique(team_abbreviation), collapse = ", "),
    num_games = n(),
    max_attendance = max(attendance),
    avg_attendance = mean(attendance, na.rm = TRUE),
    num_sellouts = sum(is_sellout, na.rm = TRUE)
  ) |>
  arrange(desc(num_sellouts)) |> 
  slice_head(n = 10) |>
  bind_rows(
    wnba_gl |>
      summarize(
        arena = paste("Total among", n_distinct(arena), "arenas"),
        home_teams = as.character(n_distinct(team_id)),
        num_games = n(),
        max_attendance = max(attendance, na.rm = TRUE),
        avg_attendance = mean(attendance, na.rm = TRUE),
        num_sellouts = sum(is_sellout, na.rm = TRUE)
      )
  ) |>
  mutate(
    pct_sellouts = num_sellouts / num_games,
    avg_capacity = avg_attendance / max_attendance
  ) |>
  janitor::clean_names(case = "title") |>
  my_gt() |>
  fmt_number(columns = `Avg Attendance`, decimals = 1) |>
  fmt_number(columns = c(`Pct Sellouts`, `Avg Capacity`), decimals = 3)
```


### Modeling fan engagement {#sec-fan}

In @sec-elo, we described our use of Elo rating to estimate the strength of the home team ($\hat{\theta}_i(t)$), the strength of the visiting team ($\hat{\theta}_j(t)$) and the home win probability ($\hat{p}_{ij}(t)$) (which is itself a function of the estimated team strengths). 

As noted in @sec-intro, there are multiple notions of what drives fan engagement relevant to the WNBA. 
Each of these notions can be captured by a simple bivariate function of the corresponding strengths of the teams $f_k(i, j, t) = f_k(\hat{\theta}_i(t), \hat{\theta}_j(t))$. 
The first and third have been rescaled so that they will have means close to 0. 

- Quality of play: fans may wish to attend a game in person if they expect the quality of play to be high. We model this as the product of the Elo rating of the two teams ($f_1(i, j, t) = \hat{\theta}_i(t) \cdot \hat{\theta}_j(t)$)
- Uncertainty of outcome: fans may wish to attend a game for which the outcome is uncertain. We model this as negative square root of the absolute difference in the two teams' Elo ratings: ($f_2(i, j, t) = -\sqrt{|\hat{\theta}_i(t) - \hat{\theta}_j(t)|}$). The negative sign ensures that closely-matched opponents will have the highest values, and the square root helps to make the spread more symmetric.[^misidentify] 
- Loss aversion: fans may wish to attend a game to see their team (the home team) win. We model this using the estimated home win probability function above ($f_3(i, j, t) = (\hat{p}_{ij}(t))^2$). Note that this term is nonlinear as in @CoatesHumphreysZhou2014.

[^misidentify]: It is possible that this formulation might misidentify how fans view outcome uncertainty. That is, fans might not seek exactly equal team strength across the two opponents, but rather would be just as satisfied as long as there is sufficient uncertainty. 

The distribution of attendance with respect to $f_1$ and $f_2$ are shown in @fig-attend-f1 and @fig-attend-f2, respectively, with the corresponding plots for percentage of estimated arena capacity shown in @fig-capacity-f1 and @fig-capacity-f2. 

#### Normalization

@fig-rescaled shows the distribution of our fan engagement variables after normalizing ($Z_k = \frac{f_k - \mu(f_k)}{\sigma(f_k)}$). 
We note that @tbl-cor shows that our three fan engagement variables are largely uncorrelated. 

```{r}
#| label: tbl-cor
#| tbl-cap: "Correlations among normalized fan engagement variables. "
wnba_gl |>
  select(matches("^z[0-9]_")) |>
  cor() |>
  janitor::clean_names(case = "title") |>
  as_tibble() |>
  my_gt() |>
  fmt_number(decimals = 3)
```

@fig-all-z illustrates how our normalized measures of fan engagement vary with respect to the estimated (pregame) strength of the home and away teams. 
The lack of correlation across these measures reported in @tbl-cor is apparent graphically. 
The 4th panel in @fig-all-z shows the simple sum of the three measures of fan engagement. 
Although this metric is not used directly in our analysis, it might be considered a composite of three competing, uncorrelated notions of fan engagement. 

```{r}
#| label: fig-all-z
#| fig-cap: "Distribution of normalized fan engagement measures as a function of home and away pregame Elo ratings. The lack of correlation across these measures is obvious. The 4th panel shows the simple sum of all three measures. "
wnba_gl |>
  select(contains("pregame"), matches("^z")) |>
  pivot_longer(cols = !contains("pregame"), names_to = "variable", values_to = "value") |>
  ggplot(aes(x = elo_home_pregame, y = elo_away_pregame, color = value)) +
  geom_point() +
  facet_wrap(vars(variable)) +
  scale_color_viridis_c(name = "Normalized\nFan Engagement\nrating") +
  scale_x_continuous(latex2exp::TeX(r"(Estimated pregame Elo rating of home team ($\hat{\theta}_i(t)$))")) +
  scale_y_continuous(latex2exp::TeX(r"(Estimated pregame Elo rating of away team ($\hat{\theta}_j(t)$))"))
```



### Control variables {#sec-control}

As time-related control variables, we include the season ($\tau_1$), the month of the year ($\tau_2$), and day of the week ($\tau_3$) as fixed effects.
To control for city-related effects, we also include the team name ($\phi_1$), as well as the metropolitan area median income ($\phi_2$, in thousands of USD) and population ($\phi_3$, in millions) of the home city[^median].
Finally, we include an indicator based on whether Caitlin Clark or A'ja Wilson played more than 0 minutes in the game ($\chi_1, \chi_2$, respectively). 

[^median]: The data for the population and the median income of cities comes from the American Community Survey, access to which was provided in R by the `r pkg("tidycensus")` package. Because the population and median income statistics were available at irregular time periods (due to the COVID-19 pandemic and other factors), impute missing values using a cubic spline for interpolated values and the last observation carried forward (or backward) for extrapolated values.  

### Interaction

Testing our primary hypotheses requires analyzing how the relationships between attendance and fan engagement have changed over time. 
To capture these changes, we include interaction terms between our measures of fan engagement ($f_k$) and season ($\tau_1$). 

### Model formulation for attendance {#sec-mod-attend}

```{r}
f_time <- "factor(season_id) + as.character(wday(game_date, label = TRUE, abbr = FALSE)) + as.character(month(game_date, label = TRUE, abbr = FALSE))"
f_place <- "factor(name_current) + median_income + metro_pop + is_cc + is_aw"
f_exp <- c("(I(z1_quality_of_play) + I(z2_uncertain_outcome) + I(z3_loss_aversion))")

wnba_df <- wnba_gl |>
  filter(!season_id %in% c(2020, 2021), is_active) |>
  filter_out(is.na(attendance))

form_basic <- paste("log(attendance + 1) ~ ", f_exp, "+", f_place, "+", f_time) |>
  as.formula()

form_interaction <- update(form_basic, paste(". ~ . + ", f_exp, "* factor(season_id)"))
```

```{r}
#| label: mod-ols
# OLS model
mod_ols_basic <- lm(formula = form_basic, data = wnba_df)
mod_ols <- lm(formula = form_interaction, data = wnba_df)
```

Like @CoatesHumphreysZhou2014, we model the natural logarithm of attendance. 

Our general regression model for attendance is:
\begin{align*}
  \ln{y} &= \beta_0 + \underbrace{\beta_1 Z_1 + \beta_2 Z_2 + \beta_3 Z_3}_{\text{fan engagement}} \\
    &+ \underbrace{\beta_4 \tau_1 + \beta_5 \tau_2 + \beta_6 \tau_3}_{\text{time-related}} + \underbrace{\beta_7 \phi_1 + \beta_8 \phi_2 + \beta_9 \phi_3}_{\text{ city-related }} \\
    &+ \underbrace{\beta_{10} \chi_1 + \beta_{11} \chi_2}_{\text{player-related}} + \underbrace{\beta_{12} Z_1 \tau_1 + \beta_{13} Z_2 \tau_2 + \beta_{14} Z_3 \tau_3}_{\text{interaction}} + \epsilon \\
    &= \symbf{X} \symbf{\beta} + \epsilon \,,
\end{align*}
where $\beta_j \in \mathbb{R}$ for $j = \{0, 1, 2, 3, 8, 9, 10, 11\}$, while the other $\beta_j$ coefficients are vectors of length equal to the number of factor levels in their corresponding categorical explanatory variable minus 1. 
Most notably $\beta_4, \beta_{12}, \beta_{13}, \beta_{14}$ have `r length(levels(factor(wnba_df$season_id))) - 1` values (corresponding to each season after 1997) and $\beta_7$ has `r length(levels(factor(wnba_df$name_current))) - 1` values (corresponding to each team). 
Only current WNBA franchises are included in the regression model. 
We also excluded all data from the 2020 and 2021 seasons, since there were no fans allowed in all or parts of both seasons. 

```{r}
#| label: mod-tobit
# Tobit model
library(VGAM)

mod_tobit_basic <- vglm(
  formula = form_basic,
  family = tobit(Upper = ifelse(wnba_df$is_sellout, log(wnba_df$attendance), Inf), imethod = 2), 
  data = wnba_df,
  na.action = "na.exclude"
)

mod_tobit <- vglm(
  formula = form_interaction,
  family = tobit(Upper = ifelse(wnba_df$is_sellout, log(wnba_df$attendance), Inf), imethod = 2), 
  data = wnba_df,
  na.action = "na.exclude"
)
```

Also like @CoatesHumphreysZhou2014, we employ a Tobit model [@tobin1958estimation; @mcdonald1980uses], because the observed attendance is censored by the arena's capacity. 
This means that what we are modeling is an unobserved, latent potential attendance, that may or may not exceed the arena's capacity. 
In this Tobit model, we have:

$$
  \mathbb{E}(\ln{y}) = \begin{cases} y^*_k & \text{if } \symbf{X} \symbf{\beta} > 0.92 \cdot y^*_k \,, \\
                         \symbf{X} \symbf{\beta} & \text{if } \symbf{X} \symbf{\beta} \leq 0.92 \cdot y^*_k
                         \end{cases} \,,
$$
where $y^*_k$ is the estimated capacity of arena $k$, as described in @sec-arenas. 

We fit the Tobit model in R using the `vglm()` and `tobit()` functions from the `r pkg("VGAM")` package. 

### Model formulation for percentage of arena capacity {#sec-mod-cap}

We fit the same model with $\gamma$ as the response, but because our response variable $\gamma$ is constrained between 0 and 1 and has significant mass at 1 exactly, we use the Papke-Wooldridge fractional response model [@papke1996econometric].
This is equivalent to fitting a generalized linear model with a logit link function and applying HC1 sandwich estimators [@mackinnon1985some] for the covariance matrix, which we accomplish using the `r pkg("sandwich")` package in R. 
$$
  \mathbb{E}{(\gamma | \symbf{X})} = G(\symbf{X} \symbf{\eta}) \,,
$$
where $G(z) = \frac{e^z}{1 + e^z}$ is the (vectorized) logistic function, $\symbf{X}$ is as above, and $\symbf{\eta}$ has the same structure as $\symbf{\beta}$ defined above. 


## Results

In keeping with recent guidance from the American Statistical Association [@wasserstein2016asa; @wasserstein2019moving], we focus our attention on the magnitude and practical impact of the fitted coefficients. 
With a sample size of over 5000 observations and dozens of coefficients, it is not so surprising that many of our statistical tests result in p-values less than 0.05. 

```{r}
#| include: false
wnba_df <- wnba_df |>
  mutate(
    y_hat_ols = exp(fitted(mod_ols)),
    y_hat_tobit = exp(fitted(mod_tobit))
  )

ggplot(data = wnba_df, aes(x = y_hat_ols, y = y_hat_tobit, color = is_sellout)) +
  geom_abline(slope = 1, linetype = 3) + 
  geom_point() +
  facet_wrap(vars(arena))

wnba_df |>
#  group_by(arena) |>
  summarize(
    cor = cor(y_hat_ols, y_hat_tobit)
  )

mod_attend <- mod_tobit
```

```{r}
#| label: tidy-vglm
library(broom)
tidy.vglm <- broom.helpers::tidy_vgam
```

### Modeling attendance {#sec-attend}

#### Baseline effects

@tbl-coefs reports estimated coefficients in the Tobit model. 
Because each measure of fan engagement interacts with season, the reported coefficients for quality of play, outcome uncertainty, and loss aversion represent their estimated effects for the reference level: a home game for the Atlanta Dream taking place on a Sunday in August 1997, in which neither Clark nor Wilson is playing.
The interaction terms, discussed in @sec-annual, capture how these relationships evolve over time. 

```{r}
#| label: sanitize
sanitize_terms <- function(x) {
  x |>
    str_replace_all("factor\\(name_current\\)", "Team: ") |>
    str_replace_all("as.character\\(wday\\(game_date\\, label = TRUE\\, abbr = FALSE\\)\\)", "Day: ") |>
    str_replace_all("as.character\\(month\\(game_date\\, label = TRUE\\, abbr = FALSE\\)\\)", "Month: ") |>
    str_replace_all("is_cc", "Player: Caitlin Clark") |>
    str_replace_all("is_aw", "Player: A'ja Wilson")
}

label_coefs <- function(x) {
  case_when(
    x == "(Intercept):1" ~ "$\\hat{\\beta}_0$",
    x == "(Intercept):2" ~ "$\\hat{\\sigma}$",
    str_detect(x, "z1") ~ "$\\hat{\\beta}_1$",
    str_detect(x, "z2") ~ "$\\hat{\\beta}_2$",
    str_detect(x, "z3") ~ "$\\hat{\\beta}_3$",
    x == "is_cc" ~ "$\\hat{\\beta}_{10}$",
    x == "is_aw" ~ "$\\hat{\\beta}_{11}$",
    x == "median_income" ~ "$\\hat{\\beta}_8$",
    x == "metro_pop" ~ "$\\hat{\\beta}_9$",
    TRUE ~ "",
  )
}
```


```{r}
#| label: tbl-coefs
#| tbl-cap: "Fitted coefficients for scalar terms in Tobit model for the natural logarithm of attendance. "
mod_attend_coefs <- mod_attend |>
  tidy(conf.int = TRUE) |>
#  map(broom::tidy, conf.int = TRUE) |>
  bind_rows(.id = "model") 

min_cc <- mod_attend_coefs |>
  filter(term == "is_cc") |>
  select(estimate) |>
  min()

beta_quality <- mod_attend_coefs |>
  filter(term == "I(z1_quality_of_play)") |>
  pull(estimate)

mod_attend_coefs |>
  filter_out(str_detect(term, "factor|name_current|character")) |>
  mutate(
    coef = label_coefs(original_term),
    term = sanitize_terms(original_term)
  ) |>
  select(coef, term, estimate, std.error, statistic, p.value, conf.low, conf.high) |>
  arrange(term) |>
  janitor::clean_names(case = "title") |>
  my_gt() |>
  fmt_markdown(columns = Coef) |>
  fmt_number(decimals = 3)
```

In the reference season, the evidence for loss aversion ($\hat{\beta}_3$) and uncertainty of outcome ($\hat{\beta}_2$) is weak, with those coefficients being of little significance, either practically or statistically. 
The coefficient for quality of play ($\hat{\beta}_1$, which *is* statistically significant) suggests that in 1997, a 1 standard deviation increase in the quality of play (as measured by the product of the Elo ratings) is associated with a modest `r round(100 * (exp(beta_quality/10) - 1), 2)`% increase in the expected latent attendance, after controlling for our other factors. 
The support for quality of play (relative to uncertainty of outcome and loss aversion) as a driver of attendance is a key finding. 

Games featuring Caitlin Clark were associated with substantially higher attendance ($\hat{\beta}_{10}$), while A'ja Wilson's presence was associated with a smaller but still nontrivial increase in attendance ($\hat{\beta}_{11}$).
The expected latent attendance more than doubles (i.e., is multiplied by $\exp(\text{`r round(min_cc, 3)`})$ = `r min_cc |> exp() |> round(2)`) when Clark plays in the game, and this is *after* controlling for the Indiana Fever team effect, the levels of fan engagement, and the season effects! 
The corresponding change in attendance for A'ja Wilson---by most accounts a better player than Clark---is multiplication by just `r mod_attend_coefs |> filter(term == "is_aw") |> pull(estimate) |> exp() |> round(2)`.

Among the control variables, median income ($\hat{\beta}_8$) and metro area population ($\hat{\beta}_9$) had small negative effects after controlling for other variables. 

@fig-attend-coefs presents the results of @tbl-coefs and some additional coefficients graphically.
Confidence intervals for the coefficient associated with the Golden State Valkyries (who have sold out all 22 of their home games) is omitted because it was so large that it renders all of the other effects indistinguishable. 
The New York Liberty have the largest team effect ($\hat{\beta}_7$), and all teams are positive relative to the Dream. 
Saturday games have the highest attendance ($\hat{\beta}_6$), as do games later in the season ($\hat{\beta}_5$). 

```{r}
#| label: fig-attend-coefs
#| fig-cap: "Relative magnitudes of fitted coefficients in attendance model, with 95% confidence intervals. The coefficient for the Golden State Valkyries is omitted. Coefficients of interest are highlighted in green. Note the relatively large size of the coefficient associated with Caitlin Clark. "
plot_attend_coefs <- mod_attend_coefs |>
  filter_out(str_detect(term, "(arena)|:|Intercept|median|metro|season|Golden")) |>
  mutate(
    hilite = str_detect(term, "is_cc|2025|quality"),
    term = sanitize_terms(term)
  ) |>
  ggplot(aes(x = exp(estimate), y = term, color = hilite)) +
  geom_vline(xintercept = 1, linetype = 3) +
  geom_point() +
  geom_errorbar(aes(xmin = exp(conf.low), xmax = exp(conf.high)), width = 0.1) +
  scale_x_continuous(latex2exp::TeX(r"(Exponentiated Fitted Coefficient ($\exp{(\hat{\beta}_i)}$))")) +
  scale_y_discrete("Term or factor level") +
  scale_color_discrete(guide = "none")
plot_attend_coefs
```


```{r}
#| label: glance-vglm
glance.vglm <- function(x, ...) {
  tibble::tibble(
    nobs        = VGAM::nobs(x),
    df.residual = VGAM::df.residual(x),
    logLik      = as.numeric(stats::logLik(x)),
    AIC         = AIC(x),
    BIC         = BIC(x),
#    deviance    = if (length(x@criterion$deviance) > 0) x@criterion$deviance else NA_real_
  )
}
```

@tbl-glance characterizes the model's fit and @tbl-anova-tobit shows the ANOVA table. 
The Tobit coefficient of `r tobit_coef <- mod_attend_coefs |> filter(original_term == "(Intercept):2") |> pull(estimate); tobit_coef |> round(2)` approximates the standard deviation of the latent attendance (on the log scale) when exponentiated (i.e., `r round(exp(tobit_coef), 3)`).

#### Annual effects {#sec-annual}

@tbl-lrt shows a likelihood ratio test between our model and a reduced model that does not include the interaction terms. 
The interactions introduce 78 additional terms corresponding to season specific slopes to our model. 
The likelihood ratio test rejects the null hypothesis that all season-specific changes in slope are equal to zero, indicating that the effects of our fan engagement measures differ significantly across seasons. 

The ANOVA table shown in @tbl-anova-tobit further support the inclusion of all three sets of interaction terms. 

```{r}
#| label: tbl-lrt
#| tbl-cap: "Model comparison: constant vs season-varying effects of fan engagement on attendance."
anova(mod_tobit_basic, mod_tobit, type = "I") |>
  my_gt() |>
  fmt_number(columns = `Pr(>Chi)`, decimals = 3)
```

@fig-season-coefs isolates the change in intercept associated with each WNBA season (the various values of the vector $\hat{\beta}_4$). 
These reveal that even after controlling for our other variables---including the Caitlin Clark effect---expected attendance has exceeded its 1997 baseline in the most recent seasons. 
Thus, even after controlling for the influence of Clark and Wilson, league expected attendance has fully recovered from the COVID-19 pandemic and the malaise of the early 2000s and 2010s. 
This is a second key finding. 

```{r}
#| label: fig-season-coefs
#| fig-cap: "Estimated season-level fixed effects for Tobit attendance model ($\\hat{\\beta}_7$). "
plot_attend_seasons <- mod_attend_coefs |>
  filter(str_detect(term, "season_id")) |>
  filter_out(str_detect(term, ":")) |>
  mutate(season_id = parse_number(str_sub(term, -4))) |>
  bind_rows(tibble(season_id = 1997, estimate = 0)) |>
  ggplot(aes(x = season_id, y = exp(estimate))) +
  geom_point() +
  geom_hline(yintercept = 1, linetype = 3) +
  geom_line() +
  geom_errorbar(aes(ymin = exp(conf.low), ymax = exp(conf.high), width = 0.2)) +
  scale_x_continuous("Season") +
  scale_y_continuous("Exponentiated Fitted Season Fixed Effect") +
  theme_minimal(base_size = 10)
plot_attend_seasons
```

@fig-season-slopes shows the exponentiated season-by-season combined slopes for our three measures of fan engagement. 
The initial values ($\hat{\beta}_1, \hat{\beta}_2, \hat{\beta}_3$, corresponding to 1997), match those in @tbl-coefs, and the subsequent interaction terms ($\hat{\beta}_{12}, \hat{\beta}_{13}, \hat{\beta}_{14}$) modify the relationships between attendance and fan engagement. 
In @fig-season-slopes, we can see that there is little to no association with attendance from the uncertainty of the outcome.[^uoh] 
The slope associated with loss aversion ($\hat{\beta}_3 + \hat{\beta}_{14}$) is always positive and generally increasing, suggesting that the attendance response to home win probability has strengthened over the league's history. 
While @AghaBerri2023 did not find a statistically significant relationship between team strength and demand between 2004 and 2017, our analysis finds that estimated home win probability is a statistically significant predictor of attendance. 
We use game-level expected probabilities and incorporate a more recent period of league growth, finding that estimated home win probability is a statistically significant predictor of attendance and that its effect has strengthened over time.
This suggests that fan behavior in the WNBA may be becoming more similar to the behavior observed in Major League Baseball by @CoatesHumphreysZhou2014.

[^uoh]: See the previous note about possible misspecification of the uncertainty of outcome variable. 

Similarly, the slope associated with quality of play  ($\hat{\beta}_1 + \hat{\beta}_{12}$) is always positive and generally increasing, following a steep fall after the league's first season. 
One possible explanation for the anomalously large fitted slope in 1997 is that a few very good and/or very bad teams stood out in the first year in the league. 
Another plausible explanation is that since all teams started with an initial Elo rating of 1500, it took a while for the relationship between attendance and team strength to settle down. 

In the case of quality of play and loss aversion, the exponentiated fitted coefficient value of approximately 0.1 in recent years means that a one standard deviation increase in these measures is associated with an increase in attendance that is about 10% higher than it would have been in the league's early years. 
This is another key finding. 

```{r}
#| label: combine-slopes
combine_slopes <- function(data) {
  slopes <- data |>
    filter(str_detect(term, "z")) |>
    mutate(
      type = ifelse(str_detect(term, "\\)$"), "main", "interaction"),
      season_id = parse_number(str_sub(term, -4)),
      variable = str_extract(term, "z[0-9][a-z,_]+")
    )
  
  main_slopes <- slopes |>
    filter(type == "main") |>
    mutate(season_id = 1997) |>
    select(season_id, variable, estimate, conf.low, conf.high)
  
  changes_in_slopes <- slopes |>
    filter(type == "interaction") |>
    select(season_id, variable, estimate, conf.low, conf.high) |>
    bind_rows(main_slopes)
  
  changes_in_slopes |>
    left_join(
      main_slopes |> select(variable, main_slope = estimate), 
      by = join_by(variable)
    ) |>
    mutate(combined_slope = main_slope + estimate)
}
```


```{r}
#| label: fig-season-slopes
#| fig-cap: "Exponentiated fitted combined (main effect + interaction effect) slope coefficients for Tobit model, by season. "
plot_attend_slopes <- mod_attend_coefs |>
  combine_slopes() |>
  ggplot(aes(x = season_id, y = exp(combined_slope), color = variable)) +
  geom_hline(yintercept = 1, linetype = 3) +
  geom_line() +
#  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.2)) +
  geom_smooth(se = 0) +
  scale_x_continuous("Season") +
  scale_y_continuous("Exponentiated Fitted Combined Slope Coefficient") +
  theme_minimal(base_size = 10)
plot_attend_slopes
```


### Modeling percentage of arena capacity

```{r}
#| label: fit-tobit-pct
# Arena capacity
mod_tobit_pct <- vglm(
  formula = update(mod_ols, pct_capacity_safe ~ .),
  family = tobit(Upper = 1, Lower = 0, imethod = 2), 
  data = wnba_df,
  na.action = "na.exclude"
)
```


```{r}
#| label: fit-glm
library(lmtest)
library(sandwich)
mod_pct <- glm(
  formula = formula(mod_tobit_pct),
  data = wnba_df, 
  family = binomial(link = "logit")
)

vcov_pct <- vcovHC(mod_pct, type = "HC1")
```

#### Main effects

@tbl-coefs2 reports the coefficients from the percentage of estimated arena capacity model, while @fig-pct-coefs graphically illustrates the relatively magnitude of the coefficients. 
The results largely mirror those of the Tobit model (compare to @tbl-coefs). 

```{r}
#| label: tbl-coefs2
#| tbl-cap: "Table of coefficients from percentage of estimated arena capacity model. "
mod_pct_coefs <- mod_pct |>
  lmtest::coeftest(vcov_pct) |>
  broom::tidy(conf.int = TRUE)

mod_pct_coefs |>
  filter_out(str_detect(term, "factor|name|character")) |>  
  mutate(
    coef = ifelse(term == "(Intercept)", "$\\hat{\\beta}_0$", label_coefs(term)),
    term = sanitize_terms(term)
  ) |>
  select(coef, everything()) |>
  arrange(term) |>
  janitor::clean_names(case = "title") |>
  my_gt() |>
  fmt_markdown(Coef) |>
  fmt_number(decimals = 3)
```

In this model, the effect associated with Caitlin Clark's appearance ($\hat{\beta}_{10}$) is even more dramatic. 
After controlling for our other factors, Clark's presence is associated with a `r cc_pct <- mod_pct_coefs |> filter(term == "is_cc") |> pull(estimate); cc_pct |> exp() |> round(2)`-fold increase in the ratio of the expected percentage of the arena filled relative to the percentage that is empty. 
In other words, for an arena that would otherwise be expected to be half-full, Clark's marginal effect corresponds to a jump in the expected percentage of arena capacity to nearly `r (0.5 + 0.5^2 * cc_pct) |> round(2)`\%. 
Clark has played all but one of her games in arena's that were more than 85% full (recall that the league average is about 50%). 

```{r}
#| label: fig-pct-coefs
#| fig-cap: "Relative magnitudes of coefficients in attendance model. Coefficients of interest are highlighted. Note the relatively large size of the coefficient associated with Caitlin Clark. "
plot_attend_coefs +
  mod_pct_coefs |>
  filter_out(str_detect(term, "(arena)|:|Intercept|median|Golden|season")) |>
  mutate(
    hilite = str_detect(term, "is_cc|2025|quality"),
    term = sanitize_terms(term)
  )
```

@tbl-glance2 and @tbl-anova2 show model fit metrics and analysis of variance for the percentage of estimated arena capacity models, respectively. 

#### Annual effects

@fig-season-coefs2 displays the estimated effects of our fan engagement measures on attendance relative to the 1997 reference season. 
Similar to the attendance trends shown in @fig-season-coefs, these effects generally remained below their 1997 levels for much of the league's history and only returned to comparable levels recently. 
In the most recent years of measurement, however, we see the league exceeding its previous highs.

```{r}
#| label: fig-season-coefs2
#| fig-cap: "Estimated season-level fixed effects for percentage of arena capacity models  ($\\hat{\\beta}_7$). "
plot_attend_seasons +
  mod_pct_coefs |>
  filter(str_detect(term, "season_id")) |>
  filter_out(str_detect(term, ":")) |>
  mutate(season_id = parse_number(str_sub(term, -4))) |>
  bind_rows(tibble(season_id = 1997, estimate = 0))
```

@fig-season-slopes2 shows that the changes in slopes associated with our measures of fan engagement follow similar patterns to those shown in @fig-season-slopes. 

```{r}
#| label: fig-season-slopes2
#| fig-cap: "Estimated season-level changes in slopes for percentage of arena capacity models. "
plot_attend_slopes +
  mod_pct_coefs |>
  combine_slopes()
```

## Conclusion

### Key results

After considering attendance across the entire history of the WNBA, we arrive at several conclusions. 

First, the effect of Caitlin Clark's presence on attendance is substantial.
As discussed above, Clark's presence in a game more than doubles the expected attendance, even after controlling for other factors. 
Although interest in a single star player is not by itself evidence of sustained development, it is notable that earlier work did not find star players to be significant drivers of WNBA game attendance [@AghaBerri2023].
 
Second, as has been widely reported, attendance has reached new heights in recent years even after controlling for Clark's presence and other game, team, and market characteristics. 
This means that the league's recent success cannot be attributed only to the "Caitlin Clark effect" and instead indicates a broader increase in demand.
The Golden State Valkyries filled all 18,064 seats in each one of their 22 home games in 2025, only one of which involved Clark. 

Finally, we find evidence that the relationships between attendance and our measurements of both quality of play and loss aversion have strengthened over time. 
The effect of quality of play appears to be more robust than the effect of our measurements of either loss aversion or uncertainty of outcome. 

These findings suggest that the WNBA is entering a more mature stage of development, in which fans increasingly attend games for the product on the court rather than primarily out of novelty or support for women's sports.

### Limitations and future work

We considered fitting a period-based model that would allow us to isolate specific periods across the league's history and craft a more narrative explanation of fan engagement trends. 
We ultimately decided against this approach because defining periods by hand would require subjective decisions about the timing of structural breaks and would risk introducing researcher bias into the analysis. 

There are reasonable objections to the way that we defined our measures of fan engagement. 
We sought definitions that were supported in the literature and could be compared directly, but other approaches could also be valid. 
In particular, our definition of the uncertainty of outcome is maximized when the two teams are equally strong, but one could argue that this should not be where the maximum occurs or that fans seek a range of uncertainty rather than a single point. 
@CoatesHumphreysZhou2014 use a reference based definition of uncertainty and find that fans only respond to uncertainty when other conditions are met, incorporating a reference based definition could be useful in future research in the WNBA context.

There is no theoretical explanation for why the slopes in @fig-season-slopes2 should converge to the same value. 

With just two seasons in our data set, while Clark's impact on attendance is undeniable, only time will tell if this impact is persistent.[^nba] 
It may very well be the case that as Clark's novelty wears off, her effect on attendance wanes. 
A follow-up study could investigate this possibility. 

[^nba]: The excitement that Clark has brought to the WNBA may be reminiscent of the infusion of interest in the NBA that occurred in the early- to mid-1980s with the arrival of Magic Johnson, Larry Bird, and Michael Jordan. 

Finally, while we chose to focus on live, in-person game attendance, fan engagement with games on television and streaming also provide important insight into league health and development. 
Future research might investigate what drives fans' engagement with games from their own homes, where the lower cost of engagement might influence fan behavior.



## References {.unnumbered}

:::{#refs}

:::

## Appendix {.appendix}

### Distribution of fan engagement


```{r}
#| label: fig-rescaled
#| fig-cap: "Distribution of fan engagement variables after rescaling. "
wnba_gl |>
  select(matches("^z")) |>
  pivot_longer(cols = everything(), values_to = "value", names_to = "fan_engagement") |>
  ggplot(aes(x = value)) +
  geom_histogram() +
  facet_wrap(vars(fan_engagement), ncol = 1) + 
  theme_minimal(base_size = 10) +
  scale_x_continuous("Standard deviations above the mean")
```




### Relationship between attendance and fan engagement

```{r}
#| label: fig-attend-f1
#| fig-cap: "WNBA game attendance by season, against normalized quality of play. "
plot_attend +
  aes(x = z1_quality_of_play) +
  scale_x_continuous("Normalized Quality of Play") +
  scale_y_log10("Attendance")
```

```{r}
#| label: fig-attend-f2
#| fig-cap: "WNBA game attendance by season, against normalized uncertainty of outcome. "
plot_attend +
  aes(x = z2_uncertain_outcome) +
  scale_x_continuous("Normalized Uncertainty of Outcome") + 
  scale_y_log10("Attendance")
```

```{r}
#| label: fig-capacity-f1
#| fig-cap: "asd"
plot_capacity +
  aes(x = z1_quality_of_play) + 
  scale_x_continuous("Normalized Quality of Play")
```

```{r}
#| label: fig-capacity-f2
#| fig-cap: "asd"
plot_capacity +
  aes(x = z2_uncertain_outcome) +
  scale_x_continuous("Normalized Uncertainty of Outcome")
```

### Model fitting metrics for attendance model {#sec-fit-tobit}

```{r}
#| label: tbl-glance
#| tbl-cap: "Model fit metrics for the Tobit model for attendance. "
mod_attend |>
  broom::glance() |>
  bind_rows(.id = "model") |>
  knitr::kable(digits = 3)
```


```{r}
#| label: tbl-anova-tobit
#| tbl-cap: "Analysis of variance for terms in models for the natural logarithm of attendance. "
mod_attend |>
  anova() |> as_tibble(rownames = "term") |>
  janitor::clean_names(case = "title") |>
  my_gt() |>
  fmt_number(columns = c(`X2 log Lik Diff`, `Log Lik`, `Pr Chi`), decimals = 3)
```

### Model fitting metrics for percentage of arena capacity model {#sec-fit-papke}

```{r}
#| label: tbl-glance2
#| tbl-cap: "Model fit metrics for the percentage of arena capacity model."
mod_pct |>
  broom::glance() |>
  knitr::kable()
```

```{r}
#| label: tbl-anova2
#| tbl-cap: "Analysis of variance for terms in models for percentage of arena capacity. "
mod_pct |>
  anova() |>
  broom::tidy() |>
  filter_out(str_detect(term, "NULL")) |>
  janitor::clean_names(case = "title") |>
  my_gt() |>
  fmt_number(columns = c(`Deviance`, `Residual Deviance`, `P Value`), decimals = 3)
```

