Motivation

This vignette illustrates how MultinomFourier can be used to identify unusual outcomes in a collection of multinomial observations. The application analyzes results from polling stations belonging to the same polling location in the 2024 Chilean municipal election. Each polling location (local de votación) hosts several polling stations (mesas), each with its own ballot box and list of registered voters.

Within a polling location, registered voters are assigned to polling stations according to alphabetical order. If voter surname is unrelated to political preference, the underlying voting distribution should be broadly similar across stations in the same polling location. This structure provides a natural multinomial goodness-of-fit setting: we can estimate a reference probability vector from the remaining stations in the location, and then test whether the result observed at a given station is statistically compatible with that null distribution.

Data

The example focuses on the Escuela Centenario polling location in the city of Iquique. This location comprises 15 polling stations, numbered 276 through 290. Most stations contain 400 registered voters, while the final two stations are smaller.

The vote matrix below records votes for four mayoral candidates, null or blank votes, and abstention. Columns correspond to individual polling stations.

library(MultinomFourier)

votes <- matrix(
  c(
    # Candidate A
     37,  34,  31,  34,  37,  37,  35,  32,  40,  37,  37,  38,  36, 22, 25,
    # Candidate B
     39,  37,  36,  36,  26,  29,  28,  39,  39,  34,  28,  32,  35, 15, 27,
    # Candidate C
     99, 107, 110, 108, 127, 107, 108, 125,  89, 100, 109, 103,  98, 74, 70,
    # Candidate D
     83,  84,   8,  87,  70, 102,  88,  74,  82,  86,  69,  86,  86, 49, 44,
    # Null and blank votes
     46,  37,  42,  41,  34,  37,  45,  26,  62,  45,  48,  49,  41, 21, 21,
    # Abstention
     96, 101, 173,  94, 106,  88,  96, 104,  88,  98, 109,  92, 104, 57, 50
  ),
  nrow = 6,
  byrow = TRUE
)

rownames(votes) <- c(
  "Candidate A",
  "Candidate B",
  "Candidate C",
  "Candidate D",
  "Null/blank",
  "Abstention"
)

colnames(votes) <- as.character(276:290)

polling_station_sizes <- colSums(votes)
polling_station_sizes
#> 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 
#> 400 400 400 400 400 400 400 400 400 400 400 400 400 238 237

The first 13 stations have exactly 400 registered voters each (\(N_m = 400\)), while stations 289 and 290 have 238 and 237 voters, respectively.

Visual Inspection

Because station sizes differ slightly, proportions provide a fairer visual comparison than raw counts. Each cell shows the percentage of votes.

prop <- sweep(votes, 2, polling_station_sizes, "/")

heat_data <- as.data.frame(as.table(prop))
names(heat_data) <- c("Option", "Station", "Proportion")

# Cell labels: percentage with one decimal, in white on dark cells
heat_data$Label <- sprintf("%.1f", 100 * heat_data$Proportion)
dark_cell <- heat_data$Proportion > min(prop) + 0.55 * diff(range(prop))
heat_data$TextColor <- ifelse(dark_cell, "white", "black")

if (requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::ggplot(
    heat_data,
    ggplot2::aes(x = Station, y = Option, fill = Proportion)
  ) +
    ggplot2::geom_tile(color = "white", linewidth = 0.25) +
    ggplot2::geom_text(
      ggplot2::aes(label = Label, color = TextColor),
      size = 2.6
    ) +
    ggplot2::scale_color_identity() +
    ggplot2::scale_fill_distiller(
      palette = "YlOrRd",
      direction = 1,
      labels = function(x) paste0(100 * x, "%")
    ) +
    ggplot2::scale_y_discrete(limits = rev(rownames(votes))) +
    ggplot2::labs(
      x = "Polling station",
      y = NULL,
      fill = "Vote proportion"
    ) +
    ggplot2::theme_minimal(base_size = 11) +
    ggplot2::theme(
      panel.grid = ggplot2::element_blank(),
      axis.text.x = ggplot2::element_text(angle = 0, hjust = 0.5)
    )
} else {
  # Rows reversed so that Candidate A appears at the top
  z <- t(prop[rev(rownames(prop)), , drop = FALSE])
  op <- par(mar = c(5, 8, 2, 2))
  image(
    seq_len(nrow(z)), seq_len(ncol(z)), z,
    col = hcl.colors(100, "YlOrRd", rev = TRUE),
    axes = FALSE,
    xlab = "Polling station",
    ylab = ""
  )
  axis(1, at = seq_len(nrow(z)), labels = rownames(z), las = 1)
  axis(2, at = seq_len(ncol(z)), labels = colnames(z), las = 1)
  lab <- matrix(sprintf("%.1f", 100 * z), nrow = nrow(z))
  dark <- z > min(z) + 0.55 * diff(range(z))
  text(
    rep(seq_len(nrow(z)), ncol(z)), rep(seq_len(ncol(z)), each = nrow(z)),
    lab, cex = 0.6, col = ifelse(dark, "white", "black")
  )
  par(op)
}

Vote proportions by polling station at Escuela Centenario.

Inspection of the heatmap highlights Station 278: Candidate D received only 8 votes (2%), compared with 70–100 votes (17%–25%) in neighboring stations, while abstention was markedly higher (173 votes).

While visually salient, graphical inspection alone does not quantify how statistically improbable this outcome is under the multinomial model.

Leave-One-Out Reference Distribution

Let \(X_m = (X_{m1}, \ldots, X_{mK})\) denote the vector of observed counts at station \(m\), with total sample size \(N_m = \sum_{k=1}^K X_{mk}\).

To evaluate whether station \(m\) is compatible with the location as a whole without allowing an anomalous station to contaminate its own reference distribution, we estimate the null probability vector using a leave-one-out scheme:

\[ \widehat p_k^{(-m)} = \frac{\sum_{j \neq m} X_{jk}}{\sum_{j \neq m} N_j}, \qquad k = 1, \ldots, K. \]

Under this formulation, station \(m\) contributes zero weight to its own reference vector \(\widehat{\mathbf{p}}^{(-m)}\).

loo_p <- sapply(
  seq_len(ncol(votes)),
  function(m) {
    other_counts <- rowSums(votes[, -m, drop = FALSE])
    other_counts / sum(other_counts)
  }
)

rownames(loo_p) <- rownames(votes)
colnames(loo_p) <- colnames(votes)

# Reference distribution for station 278
round(loo_p[, "278"], 4)
#> Candidate A Candidate B Candidate C Candidate D  Null/blank  Abstention 
#>      0.0912      0.0842      0.2700      0.2066      0.1048      0.2432

Detailed Analysis of Station 278

We model the counts at station 278 as

\[ X_{278} \sim \mathrm{Multinomial}\left(N_{278} = 400, \; \boldsymbol{\pi}\right), \]

and test whether the vote shares \(\boldsymbol{\pi}\) at this station match those of the rest of the location:

\[ \begin{aligned} H_0 &: \; \boldsymbol{\pi} = \widehat{\mathbf{p}}^{(-278)}, \\ H_1 &: \; \boldsymbol{\pi} \neq \widehat{\mathbf{p}}^{(-278)}. \end{aligned} \]

We use the log-likelihood ratio statistic, which measures how far the observed vote shares are from the reference shares:

\[ G^2(x) = 2 \sum_{k=1}^{K} x_k \ln\left(\frac{x_k}{N_{278}\, \widehat p_k^{(-278)}}\right), \]

where \(x = (x_1, \ldots, x_K)\) are the observed counts and terms with \(x_k = 0\) are taken as zero. The \(p\)-value is the probability, under \(H_0\), of a value of the statistic at least as large as the observed one:

\[ p = \mathbb{P}_{H_0}\left(G^2(X_{278}) \ge G^2(x)\right). \]

Other statistics are compared in a later section.

We first compare observed counts with expected counts under \(\widehat{\mathbf{p}}^{(-278)}\):

comparison_278 <- data.frame(
  Option = rownames(votes),
  Observed = votes[, "278"],
  Expected = round(polling_station_sizes["278"] * loo_p[, "278"], 1),
  Difference = round(votes[, "278"] - polling_station_sizes["278"] * loo_p[, "278"], 1),
  row.names = NULL
)

knitr::kable(
  comparison_278,
  format = "html",
  escape = FALSE,
  align = "l",
  caption = "Observed vs. expected counts for Polling Station 278."
)

fit_278 <- pval_fourier(
  x = votes[, "278"],
  p = loo_p[, "278"],
  rel_eps = 1e-6
)

fit_278
#> Multinomial Goodness-of-Fit Test (MultinomFourier)
#> ------------------------------------------------
#> Statistic       : Log-likelihood ratio
#> Statistic value : 160.74342
#> p-value         : 9.2015753e-33
#> log(p-value)    : -73.765933
#> 
#> Computation
#> Method          : Fourier
#> Poissonized     : Yes
#> Status          : Converged
#> Gamma           : 0.963835
#> W               : 909.592
#> Terms           : 9480
#> Gamma time      : 0.003038 s
#> Total time      : 0.4539 s

exact_fit_278 <- pval_exact(
  x = votes[, "278"],
  p = loo_p[, "278"]
)

exact_fit_278$pval
#> [1] 9.201581e-33

pval_flexible(
  x = votes[, "278"],
  p = loo_p[, "278"],
  rel_eps = 1e-6
)$method
#> [1] "fourier"

Screening All Polling Stations

We now apply pval_fourier() to all 15 stations in Escuela Centenario, each against its own leave-one-out reference vector, and check its results against the exact \(p\)-values from pval_exact(). The data are used as they are, so the comparison covers the instance actually being analyzed, including the extreme tail at station 278.

Because the \(p\)-values span more than 30 orders of magnitude, we compare the two methods through the relative error \(|\hat p_{\text{Fourier}} / p_{\text{Exact}} - 1|\) rather than the absolute difference, which would be uninformative for small \(p\)-values.

p_fourier_all <- sapply(colnames(votes), function(m) {
  pval_fourier(votes[, m], loo_p[, m], rel_eps = 1e-6)$pval
})

p_exact_all <- sapply(colnames(votes), function(m) {
  # Station 278 was already computed above
  if (m == "278") return(exact_fit_278$pval)
  pval_exact(votes[, m], loo_p[, m])$pval
})

rel_err_all <- abs(p_fourier_all / p_exact_all - 1)

screening <- data.frame(
  Station = colnames(votes),
  Voters = polling_station_sizes,
  p_Fourier = fmt_p(p_fourier_all),
  p_Exact = fmt_p(p_exact_all),
  Rel_Error = fmt_p(rel_err_all),
  row.names = NULL
)

knitr::kable(
  screening,
  format = "html",
  escape = FALSE,
  align = "l",
  caption = "Leave-one-out p-values for all stations (LLR statistic): Fourier inversion vs. exact computation."
)

The largest relative error across the 15 stations is 1.3 × 10⁻⁴, at station 285. At station 278, where \(p \approx\) 9.2 × 10⁻³³, the relative error is 6.0 × 10⁻⁷. Note that rel_eps controls the stopping rule, not the error: the observed relative errors can exceed it.

On a \(-\log_{10}(p)\) scale, where higher values indicate stronger evidence against the reference distribution, the two methods give the same profile:

screening_plot_data <- data.frame(
  Station = rep(as.integer(colnames(votes)), 2),
  Method = rep(c("Fourier (pval_fourier)", "Exact (pval_exact)"), each = ncol(votes)),
  neg_log10_p = -log10(c(p_fourier_all, p_exact_all))
)

if (requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::ggplot(
    screening_plot_data,
    ggplot2::aes(x = Station, y = neg_log10_p, color = Method, shape = Method, linetype = Method)
  ) +
    ggplot2::geom_line(linewidth = 0.7) +
    ggplot2::geom_point(size = 2.5) +
    ggplot2::scale_x_continuous(breaks = 276:290) +
    ggplot2::scale_color_manual(values = c("Fourier (pval_fourier)" = "#1f77b4", "Exact (pval_exact)" = "#ff7f0e")) +
    ggplot2::labs(
      title = "LLR p-values by polling station at Escuela Centenario",
      x = "Polling station",
      y = expression(-log[10](p)),
      color = "Method",
      shape = "Method",
      linetype = "Method"
    ) +
    ggplot2::theme_minimal(base_size = 11) +
    ggplot2::theme(legend.position = "bottom")
} else {
  matplot(
    276:290,
    -log10(cbind(p_fourier_all, p_exact_all)),
    type = "b",
    pch = c(19, 17),
    col = c("blue", "orange"),
    lty = c(1, 2),
    xlab = "Polling station",
    ylab = expression(-log[10](p)),
    xaxt = "n",
    main = "LLR p-values by polling station at Escuela Centenario"
  )
  axis(1, at = 276:290, labels = 276:290)
  legend("topright", legend = c("Fourier (pval_fourier)", "Exact (pval_exact)"),
         col = c("blue", "orange"), pch = c(19, 17), lty = c(1, 2), bty = "n")
}

Leave-one-out -log10(p) across polling stations (LLR): Fourier inversion vs. exact computation.

Station 278 stands out: its \(p\)-value is about 29 orders of magnitude smaller than that of any other station (the next smallest is station 284, with \(p \approx\) 0.003).

Comparing Multinomial Test Statistics

MultinomFourier supports several classic statistics. All of them test the same hypotheses \(H_0\) and \(H_1\) stated above for station 278; they differ in how they order the possible outcomes from least to most extreme, and therefore in which outcomes count as “at least as extreme” as the observed one.

For a station with observed counts \(x = (x_1, \ldots, x_K)\), total \(N = \sum_k x_k\) and null probabilities \(p = (p_1, \ldots, p_K)\) (here \(\widehat{\mathbf{p}}^{(-m)}\)), the multinomial sample space is

\[ \Omega_N = \left\{ y \in \{0, 1, \ldots, N\}^K : \sum_{k=1}^{K} y_k = N \right\}, \]

and the probability of an outcome \(y \in \Omega_N\) under \(H_0\) is

\[ P(y) = \frac{N!}{y_1! \cdots y_K!} \prod_{k=1}^{K} p_k^{y_k}. \]

For a statistic \(T\) in which larger values are more extreme, the \(p\)-value is

\[ p = \sum_{y \in \Omega_N} P(y) \, \mathbf{1}\left\{ T(y) \ge T(x) \right\}. \]

Probability Mass Function ("pmf"). Outcomes are ordered by their own probability under \(H_0\): an outcome is at least as extreme as \(x\) when it is at most as probable. The \(p\)-value is

\[ p_{\text{PMF}} = \sum_{y \in \Omega_N} P(y) \, \mathbf{1}\left\{ P(y) \le P(x) \right\}. \]

Pearson’s Chi-Square ("chi2"). Outcomes are ordered by the squared distance between observed and expected counts, scaled by the expected counts:

\[ \mathcal{X}^2(y) = \sum_{k=1}^{K} \frac{(y_k - N p_k)^2}{N p_k}, \qquad p_{\mathcal{X}^2} = \sum_{y \in \Omega_N} P(y) \, \mathbf{1}\left\{ \mathcal{X}^2(y) \ge \mathcal{X}^2(x) \right\}. \]

Log-Likelihood Ratio ("llr"). The statistic \(G^2\) defined above for station 278, with \(p_{G^2}\) given by the same formula with \(T = G^2\).

Power Divergence ("pd"). The family of Cressie and Read (1984), indexed by a real parameter \(\lambda\):

\[ 2 I^{\lambda}(y) = \frac{2}{\lambda(\lambda + 1)} \sum_{k=1}^{K} y_k \left[ \left( \frac{y_k}{N p_k} \right)^{\lambda} - 1 \right], \qquad \lambda \neq 0, -1, \]

with the values at \(\lambda = 0\) and \(\lambda = -1\) defined as limits. In this package, the power-divergence statistic requires \(\lambda > -1\) (values \(\lambda \le -1\) produce an error). The family contains Pearson’s statistic (\(\lambda = 1\)) and the log-likelihood ratio (\(\lambda \to 0\), available directly as stat = "llr"). We use \(\lambda = 2/3\), which lies between the two.

test_specs <- list(
  "PMF"                   = list(stat = "pmf"),
  "Pearson (Chi2)"        = list(stat = "chi2"),
  "Log-Likelihood Ratio"  = list(stat = "llr"),
  "Power Divergence (2/3)"= list(stat = "pd", lambda = 2/3)
)

# Run a p-value function on station m with the statistic described by spec
run_test <- function(fun, m, spec, ...) {
  args <- list(x = votes[, m], p = loo_p[, m], stat = spec$stat, ...)
  if (!is.null(spec$lambda)) args$lambda <- spec$lambda
  do.call(fun, args)
}

fits_by_test <- lapply(
  test_specs,
  function(spec) lapply(colnames(votes), function(m) run_test(pval_fourier, m, spec, rel_eps = 1e-6))
)

p_by_test <- sapply(fits_by_test, function(fits) sapply(fits, `[[`, "pval"))
status_by_test <- sapply(fits_by_test, function(fits) sapply(fits, `[[`, "status"))
rownames(p_by_test) <- rownames(status_by_test) <- colnames(votes)

# A dagger marks results where the series did not converge
p_display <- matrix(fmt_p(p_by_test), nrow = nrow(p_by_test), dimnames = dimnames(p_by_test))
p_display[status_by_test != 0] <- paste0(p_display[status_by_test != 0], " †")

comparison <- data.frame(
  Station = colnames(votes),
  p_display,
  check.names = FALSE,
  row.names = NULL
)

knitr::kable(
  comparison,
  format = "html",
  escape = FALSE,
  align = "l",
  caption = "Leave-one-out p-values across different test statistics (Fourier inversion). † marks results where the series did not converge (status != 0)."
)

For typical stations (e.g. 276, 277, 282, 285), all statistics yield comparable \(p\)-values. Cells marked † are those where pval_fourier() reached its limit of max_terms terms without meeting its convergence criterion (status != 0), so the value shown there is not reliable. This happens for Pearson’s chi-square at stations 278, 281, 283, 284, and for the power divergence with \(\lambda = 2/3\) at station 278; PMF and LLR converge at every station. The table below compares, for station 278, the Fourier result of each statistic with the exact \(p\)-value computed by pval_exact():

exact_278 <- sapply(
  test_specs,
  function(spec) run_test(pval_exact, "278", spec)$pval
)

fourier_278 <- lapply(fits_by_test, function(fits) fits[[which(colnames(votes) == "278")]])

comparison_278 <- data.frame(
  Statistic = names(test_specs),
  p_Fourier = fmt_p(sapply(fourier_278, `[[`, "pval")),
  Terms = sapply(fourier_278, `[[`, "n_terms"),
  Status = sapply(fourier_278, `[[`, "status"),
  p_Exact = fmt_p(exact_278),
  row.names = NULL
)

knitr::kable(
  comparison_278,
  format = "html",
  escape = FALSE,
  align = "l",
  caption = "Station 278: Fourier inversion vs. exact p-values for each statistic."
)

At this station, for PMF and LLR the Fourier series converges (7518 and 9480 terms) and agrees with the exact value. For Pearson’s chi-square and the power divergence it does not converge, and raising max_terms does not fix this: with 100,000 terms the chi-square result even becomes negative. The reason is structural. Pearson’s statistic is \(\mathcal{X}^2 = \sum_k x_k^2 / (N p_k) - N\) with integer counts \(x_k\), so its values lie close to a lattice; the terms of the Fourier series then recur instead of cancelling, and the truncated sum does not settle. The power divergence with \(\lambda = 2/3\) behaves similarly at this station. In practice, check the status of every pval_fourier() result. When it is not 0, use pval_exact() if the instance is tractable, or a statistic such as LLR or PMF, for which the series converges.

Using the exact values for the two non-converged cells, all four statistics flag station 278, with \(p\)-values between 9.2 × 10⁻³³ and 1.6 × 10⁻²⁰. At this station, the detection of the anomaly does not depend on the statistic, although the size of the \(p\)-value varies by about 12 orders of magnitude.

Methodological Scope and Limitations

When interpreting these results, three methodological aspects must be kept in mind:

1. Plug-in estimation: The leave-one-out probability vector \(\widehat{\mathbf{p}}^{(-m)}\) is estimated from data and then treated as known by pval_fourier(). The resulting numbers should strictly be interpreted as plug-in goodness-of-fit \(p\)-values.
2. Nature of the anomaly: A vanishingly small \(p\)-value indicates that station 278 is highly incompatible with the assumption that all stations share a common multinomial distribution. It does not identify the cause of the discrepancy (e.g. data-entry error, transcription error, or genuine localized demographic differences).
3. Multiple testing: When screening many stations, a multiple-testing correction (e.g. Bonferroni or Benjamini-Hochberg) should be applied over the full set of stations screened. For station 278 the conclusion does not depend on it: even a Bonferroni correction over tens of thousands of stations gives a threshold of order \(10^{-6}\), far above \(p \approx 10^{-33}\).

Summary

This example demonstrates a workflow for analyzing grouped multinomial counts:

1. Visualize proportions across groups to detect visible deviations.
2. Construct leave-one-out reference probabilities to prevent circular influence.
3. Compute multinomial \(p\)-values using pval_fourier(), checking its status and falling back to pval_exact() when the series does not converge.
4. Compare multiple statistics ("llr", "pmf", "chi2", "pd") to check that the conclusion does not depend on the choice of statistic.