Is There Position Bias in Music League Voting?

Does a song’s place in the playlist predict how many votes it receives?

1 What this report asks

Music League rounds present submitted songs to voters as a Spotify playlist. If voters work through that playlist in order, their attention and their remaining vote budget are not constant from track 1 to track n. That could produce position bias: songs in some playlist slots systematically collect more votes than equally good songs in other slots.

Three competing possibilities:

1. Primacy — early tracks do better (fresh ears, full vote budget).
2. Recency — late tracks do better (freshest in memory when voting).
3. Curvilinear — both ends beat the middle (U shape), or the middle beats both ends (inverted U, e.g. warm-up then fatigue).

All three are tested below against a null of no positional effect.

2 Analysis pipeline

The whole computation lives in this one block so that every number in the report has a single source. Unfold it to read the code; the sections that follow present its output.

Code

# ---- Load ------------------------------------------------------------------
# The export has a UTF-8 BOM, an unlabeled first column ("F"), and 19 trailing
# all-empty columns. The companion workbook's `sorted_by_vote_total` sheet has
# the same layout with that first column labeled "Total Votes", which is what
# it holds.
submissions <- read_csv(
  "raw_data/data.csv",
  col_types = cols(.default = col_character()),
  name_repair = "unique_quiet"
) |>
  select(where(~ !all(is.na(.x)))) |>
  filter(!if_all(everything(), is.na)) |>
  rename(
    votes        = 1,
    order        = `playlist_order`,
    title        = `Title`,
    album        = `Album`,
    artist       = `Artist`,
    submitter    = `Submitter`,
    round        = `Round`,
    playlist_url = `Spotify Playlist URL`
  ) |>
  mutate(
    across(c(votes, order), as.integer),
    across(c(title, album, artist, submitter, round), str_squish)
  )

n_raw <- nrow(submissions)

# ---- Round-level structure -------------------------------------------------
round_info <- submissions |>
  summarise(
    n_songs   = n(),
    n_missing = sum(is.na(votes)),
    pool      = sum(votes, na.rm = TRUE),
    mean_v    = mean(votes, na.rm = TRUE),
    sd_v      = sd(votes, na.rm = TRUE),
    .by = round
  )

# ---- Within-round measures -------------------------------------------------
dat <- submissions |>
  mutate(
    n_songs    = n(),
    vote_z     = as.numeric(scale(votes)),          # 0 = round mean, +/-1 = 1 SD
    vote_ratio = votes / mean(votes, na.rm = TRUE),  # 1 = round average
    rel_pos    = (order - 1) / (n_songs - 1),        # 0 = first, 1 = last
    .by = round
  ) |>
  mutate(round_f = factor(round), long = n_songs > 20,
         n_c = n_songs - median(round_info$n_songs))

analysis  <- dat |> drop_na(votes, rel_pos)
n_dropped <- nrow(dat) - nrow(analysis)
sd_typ    <- median(round_info$sd_v)

dupes     <- analysis |> add_count(round, order, name = "k") |> filter(k > 1)
dupe_keys <- dupes |> distinct(round, order)
contiguous_rounds <- sum(map2_lgl(round_info$round, round_info$n_songs,
  ~ identical(sort(dat$order[dat$round == .x]), 1:.y)))

# ---- Randomization check: does position depend on submitter? ---------------
sub_pos <- analysis |>
  summarise(n = n(), mean_rel = mean(rel_pos), .by = submitter) |>
  filter(n >= 10)
aov_fit <- aov(rel_pos ~ submitter,
               data = analysis |> filter(submitter %in% sub_pos$submitter))
aov_p <- summary(aov_fit)[[1]][["Pr(>F)"]][1]
aov_F <- summary(aov_fit)[[1]][["F value"]][1]

# ---- Test 1a: per-round correlations, Fisher-z meta-analysis ---------------
round_cor <- analysis |>
  summarise(
    n = n(),
    r_pearson  = cor(rel_pos, votes, method = "pearson"),
    r_spearman = cor(rel_pos, votes, method = "spearman"),
    p_round    = cor.test(rel_pos, votes)$p.value,
    .by = round
  )

n_sig     <- sum(round_cor$p_round < 0.05)
n_sig_neg <- sum(round_cor$p_round < 0.05 & round_cor$r_pearson < 0)
exp_sig   <- 0.05 * nrow(round_cor)

zt      <- round_cor |> mutate(z = atanh(r_pearson), w = n - 3)
z_bar   <- weighted.mean(zt$z, zt$w)
se_zbar <- sqrt(1 / sum(zt$w))
z_stat  <- z_bar / se_zbar
p_meta  <- 2 * pnorm(-abs(z_stat))
ci_meta <- tanh(z_bar + c(-1.96, 1.96) * se_zbar)
r_meta  <- tanh(z_bar)
t_sp    <- t.test(round_cor$r_spearman)

# ---- Cluster-robust inference (CR1), written out to avoid dependencies -----
cluster_vcov <- function(fit, cluster) {
  keep <- !is.na(coef(fit))
  X  <- model.matrix(fit)[, keep, drop = FALSE]
  u  <- residuals(fit)
  cl <- factor(cluster)
  N <- nrow(X); K <- ncol(X); G <- nlevels(cl)
  XtX_inv <- solve(crossprod(X))
  meat <- matrix(0, K, K)
  for (g in levels(cl)) {
    i <- which(cl == g)
    s <- crossprod(X[i, , drop = FALSE], u[i])
    meat <- meat + s %*% t(s)
  }
  XtX_inv %*% meat %*% XtX_inv * (G / (G - 1)) * ((N - 1) / (N - K))
}

cluster_table <- function(fit, cluster, terms) {
  V  <- cluster_vcov(fit, cluster)
  b  <- coef(fit)[!is.na(coef(fit))]
  se <- sqrt(diag(V))
  G  <- n_distinct(cluster)
  keep <- names(b) %in% terms
  tibble(
    term = names(b)[keep], estimate = b[keep], std.error = se[keep],
    statistic = (b / se)[keep],
    p.value   = 2 * pt(-abs((b / se)[keep]), df = G - 1),
    conf.low  = (b - qt(0.975, G - 1) * se)[keep],
    conf.high = (b + qt(0.975, G - 1) * se)[keep]
  )
}

# ---- Test 1b / 2: fixed-effects models ------------------------------------
m_lin  <- lm(vote_z ~ rel_pos + round_f, data = analysis)
m_quad <- lm(vote_z ~ rel_pos + I(rel_pos^2) + round_f, data = analysis)
m_raw  <- lm(votes  ~ rel_pos + round_f, data = analysis)

lin_row   <- cluster_table(m_lin,  analysis$round, "rel_pos")
raw_row   <- cluster_table(m_raw,  analysis$round, "rel_pos")
quad_rows <- cluster_table(m_quad, analysis$round, c("rel_pos", "I(rel_pos^2)"))

anova_quad <- anova(m_lin, m_quad)
bq         <- coef(m_quad)
vertex     <- -bq[["rel_pos"]] / (2 * bq[["I(rel_pos^2)"]])
shape      <- if (bq[["I(rel_pos^2)"]] > 0) "U-shaped (both ends favoured)" else
              "inverted-U (middle favoured)"

# ---- Permutation test (exact within-round null) ---------------------------
# Within-round demeaning gives the same slope as round fixed effects, cheaply.
perm_dat <- analysis |>
  mutate(
    x1 = rel_pos - mean(rel_pos),
    x2 = rel_pos^2 - mean(rel_pos^2),
    z  = vote_z - mean(vote_z),
    .by = round
  )

Xp    <- cbind(x1 = perm_dat$x1, x2 = perm_dat$x2)
A     <- solve(crossprod(Xp), t(Xp))                                  # betas = A %*% z
A_lin <- solve(crossprod(Xp[, 1, drop = FALSE]), t(Xp[, 1, drop = FALSE]))

obs_lin  <- as.numeric(A_lin %*% perm_dat$z)
obs_quad <- as.numeric(A %*% perm_dat$z)[2]

# Shuffling vote_z within a round is exactly equivalent to shuffling raw votes:
# a within-round shuffle leaves that round's mean and SD untouched, so the
# round's vote pool and vote distribution are both preserved under the null.
grp <- split(seq_len(nrow(perm_dat)), perm_dat$round)
z0  <- perm_dat$z

perm_stats <- map_dfr(seq_len(N_PERM), function(b) {
  idx <- seq_len(nrow(perm_dat))
  for (g in grp) idx[g] <- sample(g)
  zp <- z0[idx]
  tibble(lin = as.numeric(A_lin %*% zp), quad = as.numeric(A %*% zp)[2])
})

p_perm_lin  <- mean(abs(perm_stats$lin)  >= abs(obs_lin))
p_perm_quad <- mean(abs(perm_stats$quad) >= abs(obs_quad))

# ---- Shape-free smooth ----------------------------------------------------
m_gam <- gam(vote_z ~ s(rel_pos, k = 8, bs = "cs") + s(round_f, bs = "re"),
             data = analysis, method = "REML")
s_tab <- summary(m_gam)$s.table

gam_nd <- tibble(rel_pos = seq(0, 1, length.out = 200),
                 round_f = analysis$round_f[1])
gam_pr <- predict(m_gam, gam_nd, se.fit = TRUE, exclude = "s(round_f)")
gam_nd <- gam_nd |> mutate(fit = as.numeric(gam_pr$fit),
                           se  = as.numeric(gam_pr$se.fit))

# ---- Binned views --------------------------------------------------------
bins <- analysis |>
  mutate(decile = ntile(rel_pos, 10)) |>
  summarise(n = n(), m = mean(vote_z), se = sd(vote_z) / sqrt(n()), .by = decile) |>
  arrange(decile)

std_rounds <- round_info |> filter(n_songs <= 20) |> pull(round)
abs_dat    <- analysis |> filter(round %in% std_rounds)
abs_bins   <- abs_dat |>
  summarise(n = n(), m = mean(vote_z), se = sd(vote_z) / sqrt(n()), .by = order) |>
  filter(n >= 20) |> arrange(order)

edge <- abs_dat |>
  mutate(zone = case_when(order <= 3 ~ "First three",
                          order > n_songs - 3 ~ "Last three",
                          TRUE ~ "Middle"))
edge_tests <- bind_rows(
  tidy(t.test(vote_z ~ zone == "First three",
              data = edge |> filter(zone != "Last three"))) |>
    mutate(Contrast = "First three vs. middle"),
  tidy(t.test(vote_z ~ zone == "Last three",
              data = edge |> filter(zone != "First three"))) |>
    mutate(Contrast = "Last three vs. middle")
)

# ---- Does round length moderate the effect? ------------------------------
m_mod_bin  <- lm(vote_z ~ rel_pos * long + round_f, data = analysis)
m_mod_cont <- lm(vote_z ~ rel_pos * n_c  + round_f, data = analysis)
mod_bin  <- cluster_table(m_mod_bin,  analysis$round, "rel_pos:longTRUE")
mod_cont <- cluster_table(m_mod_cont, analysis$round, "rel_pos:n_c")

len_groups <- analysis |>
  mutate(grp = cut(n_songs, c(0, 16, 18, 20, 40),
                   labels = c("13-16", "17-18", "19-20", "21-37"))) |>
  summarise(Rounds = n_distinct(round), Songs = n(),
            r = cor(rel_pos, votes), .by = grp) |>
  arrange(grp)

# ---- Robustness ----------------------------------------------------------
fit_slope <- function(d, label, formula = vote_z ~ rel_pos + round_f,
                      term = "rel_pos") {
  cluster_table(lm(formula, data = d), d$round, term) |>
    transmute(Specification = label,
              Estimate = sprintf("%+.4f", estimate),
              `95% CI` = sprintf("[%+.3f, %+.3f]", conf.low, conf.high),
              p = sprintf("%.3f", p.value), Songs = nrow(d))
}

robust <- bind_rows(
  fit_slope(analysis, "Main: all songs, round FE"),
  fit_slope(anti_join(analysis, dupe_keys, by = c("round", "order")),
            "Excluding duplicated positions"),
  fit_slope(analysis |> filter(round %in% std_rounds),
            "Standard-length rounds only (13-20)"),
  fit_slope(analysis |> filter(n_songs > 20), "Long rounds only (21-37)"),
  fit_slope(analysis |> filter(!round %in% round_info$round[round_info$n_missing > 0]),
            "Excluding rounds with any missing votes"),
  fit_slope(analysis, "Adding submitter fixed effects",
            formula = vote_z ~ rel_pos + round_f + submitter),
  fit_slope(analysis, "Outcome = votes / round mean",
            formula = vote_ratio ~ rel_pos + round_f),
  fit_slope(analysis, "Outcome = raw votes", formula = votes ~ rel_pos + round_f)
)

# ---- Precision -----------------------------------------------------------
mde_sd <- 2.8 * lin_row$std.error   # ~80% power, two-sided alpha = .05

3 Summary of findings

There is a small but statistically detectable primacy effect: songs earlier in the playlist collect more votes. There is no evidence of curvature.

Code

tibble(
  Hypothesis = c(
    "Primacy (earlier slots win more votes)",
    "Recency (later slots win more votes)",
    "Curvilinear (U or inverted-U)"
  ),
  `Headline test` = c(
    sprintf("Within-round permutation p = %.3f", p_perm_lin),
    sprintf("Effect runs the other way (%s%.3f SD)",
            ifelse(lin_row$estimate < 0, "", "+"), lin_row$estimate),
    sprintf("Permutation p = %.3f; GAM edf = %.2f", p_perm_quad, s_tab[1, "edf"])
  ),
  Verdict = c(
    ifelse(p_perm_lin < 0.05 & lin_row$estimate < 0, "Supported", "Not supported"),
    ifelse(p_perm_lin < 0.05 & lin_row$estimate > 0, "Supported", "Not supported"),
    ifelse(p_perm_quad < 0.05, "Supported", "Not supported")
  )
) |>
  kable(caption = "Verdict on the three hypotheses.")

Verdict on the three hypotheses.

Hypothesis	Headline test	Verdict
Primacy (earlier slots win more votes)	Within-round permutation p = 0.009	Supported
Recency (later slots win more votes)	Effect runs the other way (-0.250 SD)	Not supported
Curvilinear (U or inverted-U)	Permutation p = 0.926; GAM edf = 1.05	Not supported

The size of the effect, stated three ways:

• Moving a song from the first slot to the last slot costs 0.25 within-round standard deviations (95% CI 0.06 to 0.44).
• On the raw vote scale that is 0.94 votes (95% CI 0.13 to 1.75), against an average song total of 8.8 votes — roughly 11% of a typical song’s haul.
• The average within-round correlation between position and votes is -0.086 (95% CI -0.147 to -0.024), and 45 of 64 rounds lean in the primacy direction.

Two important qualifications, developed in Section 11 and Section 9:

1. The estimate sits right at the edge of this dataset’s resolving power, so the direction is better established than the magnitude.
2. The effect is a gradual gradient across the whole playlist, not a bonus for the first few slots — first-three and last-three contrasts are both null.

4 Data and measures

Code

tibble(
  Quantity = c(
    "Rounds", "Songs", "Songs per round (min / median / max)",
    "Total votes cast per round (min / median / max)",
    "Votes per song per round (min / median / max)",
    "Within-round SD of votes (median)"
  ),
  Value = c(
    as.character(nrow(round_info)), as.character(n_raw),
    paste(min(round_info$n_songs), median(round_info$n_songs),
          max(round_info$n_songs), sep = " / "),
    paste(min(round_info$pool), median(round_info$pool),
          max(round_info$pool), sep = " / "),
    sprintf("%.1f / %.1f / %.1f", min(round_info$mean_v),
            median(round_info$mean_v), max(round_info$mean_v)),
    sprintf("%.2f", sd_typ)
  )
) |>
  kable(caption = "Structure of the submissions data.")

Structure of the submissions data.

Quantity	Value
Rounds	64
Songs	1185
Songs per round (min / median / max)	13 / 17 / 37
Total votes cast per round (min / median / max)	120 / 161 / 204
Votes per song per round (min / median / max)	5.2 / 10.0 / 10.0
Within-round SD of votes (median)	4.11

4.1 Why everything must be measured within round

Votes are zero-sum within a round. Each round distributes a fixed pool of voter points, so one song gaining a vote means another loses one. Round pools range from 120 to 204 points, and because round sizes differ, average votes per song differs across rounds (5.2 to 10.0).

Rounds differ in length. With 13 to 37 songs per round, “position 15” is late in a short round and mid-pack in a long one.

Comparing raw votes across rounds would therefore conflate position with round size. Every measure is defined within round:

Measure	Definition	Reading
vote_z	votes z-scored within round	0 = round average, ±1 = one within-round SD
vote_ratio	votes ÷ round mean votes	1 = round average, 1.2 = 20% above
rel_pos	\((\text{order}-1)/(n-1)\)	0 = first track, 1 = last track
order	raw playlist position	absolute slot, 1-based

5 rows are dropped for missing votes, leaving 1180 songs across 64 rounds.

4.2 Data quality

Code

tibble(
  Check = c("Rows with missing votes", "Rounds where position runs exactly 1..n",
            "Rows sharing a duplicated (round, position)", "Distinct submitters"),
  Result = c(as.character(sum(is.na(dat$votes))),
             paste0(contiguous_rounds, " of ", nrow(round_info)),
             as.character(nrow(dupes)),
             as.character(n_distinct(analysis$submitter)))
) |>
  kable(caption = "Duplicated positions are treated as ties in position, and excluded in a robustness check.")

Duplicated positions are treated as ties in position, and excluded in a robustness check.

Check	Result
Rows with missing votes	5
Rounds where position runs exactly 1..n	62 of 64
Rows sharing a duplicated (round, position)	4
Distinct submitters	53

4.3 Is playlist order plausibly random?

A position effect is only interpretable as a voting effect if position is not itself a proxy for song quality. If strong submitters reliably landed early slots, the effect would be confounded.

Code

tibble(
  Test = "ANOVA: relative position ~ submitter (submitters with >= 10 songs)",
  `F` = sprintf("%.2f", aov_F), `p` = sprintf("%.3f", aov_p)
) |>
  kable(caption = "Randomization check on playlist position.")

Randomization check on playlist position.

Test	F	p
ANOVA: relative position ~ submitter (submitters with >= 10 songs)	1.32	0.156

Code

sub_pos |>
  mutate(submitter = fct_reorder(submitter, mean_rel)) |>
  ggplot(aes(mean_rel, submitter)) +
  geom_vline(xintercept = 0.5, linetype = 2, colour = "grey40") +
  geom_point(aes(size = n), alpha = 0.7) +
  scale_size_continuous(range = c(1, 4), name = "Songs") +
  labs(title = "Mean playlist position by submitter",
       subtitle = "Reference line = 0.5 (no positional advantage)",
       x = "Mean relative position (0 = first, 1 = last)", y = NULL)

Figure 1: Mean relative playlist position by submitter. If position is assigned independently of who submitted, these scatter around 0.5 with no submitter far outside the band.

No evidence that position depends on submitter (p = 0.156), consistent with order being assigned independently of the song. The position effect can reasonably be read as a voting effect rather than a quality artifact.

5 The relationship, plotted

Code

ggplot(analysis, aes(rel_pos, vote_z)) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_point(alpha = 0.25) +
  geom_smooth(method = "lm", formula = y ~ x, se = TRUE,
              colour = "#2c7fb8", linewidth = 1) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE,
              colour = "#e6550d", linewidth = 1) +
  geom_smooth(method = "loess", span = 0.75, se = FALSE,
              colour = "grey20", linetype = 2, linewidth = 0.9) +
  labs(title = "Votes vs. playlist position",
       subtitle = "Each point is one song; votes standardized within its round",
       x = "Relative playlist position (0 = first, 1 = last)",
       y = "Votes (SD within round)")

Figure 2: Within-round standardized votes against relative playlist position. Blue = linear fit, orange = quadratic fit, dashed = loess. The fitted lines tilt down: earlier slots earn more votes.

The quadratic and loess curves track the straight line closely — a first hint that whatever is going on is monotone rather than curved.

6 Test 1: Is there a monotone (primacy or recency) effect?

6.1 Per-round correlations

The cleanest test of a consistent within-round trend: correlate position with votes separately in each of the 64 rounds, then ask whether those correlations center on zero. Spearman is reported alongside Pearson because votes are bounded counts.

Code

tibble(
  Statistic = c(
    "Mean Pearson r across rounds", "Mean Spearman rho across rounds",
    "Rounds with negative r (primacy direction)",
    "Fisher-z meta-analytic r [95% CI]", "Fisher-z test of r = 0",
    "One-sample t-test on Spearman rho"
  ),
  Value = c(
    sprintf("%.4f", mean(round_cor$r_pearson)),
    sprintf("%.4f", mean(round_cor$r_spearman)),
    sprintf("%d of %d (%.0f%%)", sum(round_cor$r_pearson < 0), nrow(round_cor),
            100 * mean(round_cor$r_pearson < 0)),
    sprintf("%.4f [%.4f, %.4f]", r_meta, ci_meta[1], ci_meta[2]),
    sprintf("z = %.2f, p = %.3f", z_stat, p_meta),
    sprintf("t(%d) = %.2f, p = %.3f", t_sp$parameter, t_sp$statistic, t_sp$p.value)
  )
) |>
  kable(caption = "Per-round correlation between playlist position and votes.")

Per-round correlation between playlist position and votes.

Statistic	Value
Mean Pearson r across rounds	-0.0665
Mean Spearman rho across rounds	-0.0694
Rounds with negative r (primacy direction)	45 of 64 (70%)
Fisher-z meta-analytic r [95% CI]	-0.0857 [-0.1471, -0.0235]
Fisher-z test of r = 0	z = -2.70, p = 0.007
One-sample t-test on Spearman rho	t(63) = -2.36, p = 0.021

Both tests reject zero. Taken one at a time, though, only 4 of 64 rounds are individually significant — barely more than the 3.2 that chance alone would deliver at α = 0.05, so the count itself proves nothing. What is mildly notable is that all 4 of them point the same way (primacy); under the null they should split evenly between directions. Either way, the effect is only reliably visible by pooling, which is exactly what a small consistent bias looks like.

Code

ggplot(round_cor, aes(r_spearman)) +
  geom_histogram(bins = 20, fill = "grey70", colour = "white") +
  geom_vline(xintercept = 0, linewidth = 1) +
  geom_vline(xintercept = mean(round_cor$r_spearman), colour = "#e6550d",
             linewidth = 1, linetype = 2) +
  labs(title = "Per-round position-votes correlation",
       subtitle = "Solid line = no effect; dashed orange = observed mean",
       x = "Spearman rho (position vs. votes) within round", y = "Rounds")

Figure 3: Distribution of the within-round position-votes correlation. The distribution straddles zero but its center of mass sits left of it, the signature of a weak primacy effect.

6.2 Pooled within-round regression

Equivalently, pool all songs with round fixed effects, absorbing every round-level difference in size and vote pool. Standard errors are cluster-robust by round (CR1), because songs in a round compete for one fixed vote pool and so are not independent observations.

Code

bind_rows(
  lin_row |> mutate(term = "Standardized votes (SD)"),
  raw_row |> mutate(term = "Raw votes")
) |>
  select(Outcome = term, Estimate = estimate, SE = std.error,
         t = statistic, p = p.value, `CI low` = conf.low, `CI high` = conf.high) |>
  kable(digits = 4, caption = "First-slot-to-last-slot effect of playlist position, round fixed effects, cluster-robust SEs by round.")

First-slot-to-last-slot effect of playlist position, round fixed effects, cluster-robust SEs by round.

Outcome	Estimate	SE	t	p	CI low	CI high
Standardized votes (SD)	-0.2504	0.0969	-2.5836	0.0121	-0.4441	-0.0567
Raw votes	-0.9369	0.4052	-2.3120	0.0241	-1.7466	-0.1271

6.3 Permutation test

The zero-sum structure makes model-based p-values worth double-checking. Shuffling votes among songs within each round generates an exact null: it preserves each round’s vote pool, its vote distribution, and its length, destroying only the link between position and votes.

Code

tibble(
  Statistic = c("Linear slope", "Quadratic coefficient"),
  Observed = sprintf("%+.4f", c(obs_lin, obs_quad)),
  `Null SD` = sprintf("%.4f", c(sd(perm_stats$lin), sd(perm_stats$quad))),
  `Two-sided p` = sprintf("%.4f", c(p_perm_lin, p_perm_quad))
) |>
  kable(caption = sprintf("Within-round permutation test, %s replicates.",
                          format(N_PERM, big.mark = ",")))

Within-round permutation test, 5,000 replicates.

Statistic	Observed	Null SD	Two-sided p
Linear slope	-0.2504	0.0953	0.0094
Quadratic coefficient	+0.0301	0.3529	0.9264

Code

perm_stats |>
  pivot_longer(everything(), names_to = "stat", values_to = "value") |>
  mutate(stat = recode(stat, lin = "Linear slope", quad = "Quadratic coefficient")) |>
  ggplot(aes(value)) +
  geom_histogram(bins = 50, fill = "grey75", colour = "white") +
  geom_vline(data = tibble(stat = c("Linear slope", "Quadratic coefficient"),
                           obs = c(obs_lin, obs_quad)),
             aes(xintercept = obs), colour = "#e6550d", linewidth = 1) +
  facet_wrap(~stat, scales = "free") +
  labs(title = "Observed effects against the within-round null",
       x = "Coefficient", y = "Permutation replicates")

Figure 4: Permutation null distributions with observed statistics marked. The linear slope sits in the left tail of its null; the quadratic coefficient sits in the middle of its own.

The permutation p of 0.0094 closely matches the cluster-robust p of 0.012, so the result does not hinge on regression assumptions.

7 Test 2: Is the relationship curvilinear?

A U or inverted-U would be invisible to the tests above. Three checks: a quadratic term, a nested model comparison, and a shape-free smooth.

Code

quad_rows |>
  mutate(term = recode(term, "rel_pos" = "Relative position (linear)",
                       "I(rel_pos^2)" = "Relative position squared")) |>
  select(Term = term, Estimate = estimate, SE = std.error, t = statistic,
         p = p.value) |>
  kable(digits = 4, caption = "Quadratic model, round fixed effects, cluster-robust SEs.")

Quadratic model, round fixed effects, cluster-robust SEs.

Term	Estimate	SE	t	p
Relative position (linear)	-0.2805	0.3969	-0.7068	0.4823
Relative position squared	0.0301	0.3907	0.0770	0.9389

Code

tibble(
  Comparison = "Adding the quadratic term to the linear model",
  `F` = sprintf("%.2f", anova_quad$F[2]),
  `p` = sprintf("%.3f", anova_quad$`Pr(>F)`[2]),
  `Permutation p` = sprintf("%.3f", p_perm_quad),
  `Adj. R2 change` = sprintf("%+.5f", summary(m_quad)$adj.r.squared -
                               summary(m_lin)$adj.r.squared)
) |>
  kable(caption = "Nested model comparison for curvature.")

Nested model comparison for curvature.

Comparison	F	p	Permutation p	Adj. R2 change
Adding the quadratic term to the linear model	0.01	0.932	0.926	-0.00094

The quadratic term is indistinguishable from zero by every route. Adding it also inflates the standard error on the linear term about fourfold (the two predictors are nearly collinear over [0, 1]), which is why the linear term loses significance here — that is multicollinearity, not evidence against the linear effect.

The fitted curvature is nominally U-shaped (both ends favoured) with a turning point at relative position 4.66. That turning point falls far outside the observed range of positions, so within the actual data the fitted curve is effectively a straight line – it is not describing a real interior peak or trough.

7.1 Shape-free check: GAM smooth

A penalized spline lets the data choose the shape, with a random intercept per round. Effective degrees of freedom near 1 means the data prefer a straight line; higher values indicate genuine curvature.

Code

tibble(
  Term = "s(relative position)",
  `Effective df` = sprintf("%.2f", s_tab[1, "edf"]),
  `F` = sprintf("%.2f", s_tab[1, "F"]),
  `p` = sprintf("%.3f", s_tab[1, "p-value"]),
  `Deviance explained` = sprintf("%.2f%%", 100 * summary(m_gam)$dev.expl)
) |>
  kable(caption = "GAM smooth for playlist position, with round random intercepts.")

GAM smooth for playlist position, with round random intercepts.

Term	Effective df	F	p	Deviance explained
s(relative position)	1.05	0.90	0.007	0.62%

The smooth uses 1.05 effective degrees of freedom — the penalty has shrunk it to essentially a straight line — while still being significant (p = 0.007). The data want a slope, not a curve. Note also that position explains only 0.6% of the deviance in votes: real, but a tiny part of what decides a round.

Code

ggplot(gam_nd, aes(rel_pos, fit)) +
  geom_hline(yintercept = 0, linetype = 2, colour = "grey40") +
  geom_ribbon(aes(ymin = fit - 1.96 * se, ymax = fit + 1.96 * se),
              fill = "#2c7fb8", alpha = 0.2) +
  geom_line(colour = "#2c7fb8", linewidth = 1) +
  labs(title = "Estimated shape of the position effect",
       subtitle = "Penalized spline with round random intercepts; ribbon = 95% CI",
       x = "Relative playlist position (0 = first, 1 = last)",
       y = "Votes (SD within round)")

Figure 5: GAM-estimated shape of the position effect with 95% confidence band. The band sits above zero early in the playlist and below it late, crossing near the middle.

7.2 Binned means

The least assumption-laden view: average standardized votes within position deciles.

Code

ggplot(bins, aes(decile, m)) +
  geom_hline(yintercept = 0, colour = "grey40", linetype = 2) +
  geom_smooth(method = "lm", se = FALSE, colour = "#e6550d", linewidth = 0.8) +
  geom_errorbar(aes(ymin = m - 1.96 * se, ymax = m + 1.96 * se), width = 0.2) +
  geom_point(size = 2.5) +
  scale_x_continuous(breaks = 1:10) +
  labs(title = "Standardized votes by position decile",
       subtitle = "Decile 1 = earliest tenth of the playlist, 10 = latest",
       x = "Decile of relative playlist position",
       y = "Mean votes (SD within round)")

Figure 6: Mean standardized votes by decile of relative playlist position, with 95% confidence intervals. The early deciles sit above zero and the late deciles below, but individual deciles are noisy.

Code

bins |>
  mutate(across(c(m, se), ~ round(.x, 3))) |>
  rename(Decile = decile, Songs = n, `Mean vote_z` = m, SE = se) |>
  kable(caption = "Position deciles. The gradient is gradual; no decile dominates.")

Position deciles. The gradient is gradual; no decile dominates.

Decile	Songs	Mean vote_z	SE
1	118	0.064	0.088
2	118	0.038	0.086
3	118	0.259	0.099
4	118	0.027	0.090
5	118	0.024	0.089
6	118	-0.009	0.093
7	118	-0.071	0.087
8	118	-0.153	0.085
9	118	-0.113	0.083
10	118	-0.066	0.093

8 Test 3: Absolute slots and playlist edges

Relative position assumes bias scales with playlist length. Fatigue might instead attach to absolute slots — the 1st track, the 20th — regardless of round length. Restricting to standard-length rounds keeps absolute slots comparable.

Code

ggplot(abs_bins, aes(order, m)) +
  geom_hline(yintercept = 0, colour = "grey40", linetype = 2) +
  geom_errorbar(aes(ymin = m - 1.96 * se, ymax = m + 1.96 * se), width = 0.2) +
  geom_point(size = 2.5) +
  geom_smooth(method = "loess", se = FALSE, colour = "#e6550d", linewidth = 0.9) +
  scale_x_continuous(breaks = seq(1, max(abs_bins$order), 2)) +
  labs(title = "Standardized votes by absolute playlist slot",
       subtitle = sprintf("Rounds of 13-20 songs (%d rounds, %d songs); slots with >= 20 songs",
                          length(std_rounds), nrow(abs_dat)),
       x = "Playlist position", y = "Mean votes (SD within round)")

Figure 7: Mean standardized votes by absolute playlist slot in rounds of 13-20 songs, with 95% confidence intervals. Slot-level means are too noisy to pin down individual slots.

Code

edge |>
  summarise(Songs = n(), `Mean vote_z` = round(mean(vote_z), 3), .by = zone) |>
  rename(Zone = zone) |>
  kable(caption = "Playlist edges versus the middle, standard-length rounds.")

Playlist edges versus the middle, standard-length rounds.

Zone	Songs	Mean vote_z
First three	165	0.067
Middle	585	0.000
Last three	165	-0.068

Code

edge_tests |>
  transmute(Contrast, `Difference (SD)` = sprintf("%+.3f", -estimate),
            `t` = sprintf("%.2f", statistic), `p` = sprintf("%.3f", p.value)) |>
  kable(caption = "Welch t-tests for edge effects.")

Welch t-tests for edge effects.

Contrast	Difference (SD)	t	p
First three vs. middle	+0.066	-0.80	0.426
Last three vs. middle	-0.068	0.82	0.415

Both edge contrasts are null. Combined with the linear result, the picture is a gradual gradient down the whole playlist, not a discrete bonus for opening tracks or a penalty for closing ones. Chopping the playlist into thirds throws away the ordering information that makes the gradient detectable.

9 Does round length moderate the effect?

Fatigue should plausibly bite harder in a 37-song playlist than an 18-song one.

Code

bind_rows(
  mod_bin  |> mutate(term = "Position x long round (>20 songs)"),
  mod_cont |> mutate(term = "Position x round length (per extra song)")
) |>
  select(Interaction = term, Estimate = estimate, SE = std.error,
         t = statistic, p = p.value) |>
  kable(digits = 4, caption = "Interaction tests for round-length moderation.")

Interaction tests for round-length moderation.

Interaction	Estimate	SE	t	p
Position x long round (>20 songs)	-0.4058	0.2163	-1.8767	0.0652
Position x round length (per extra song)	-0.0219	0.0172	-1.2707	0.2085

Code

len_groups |>
  mutate(r = round(r, 3)) |>
  rename(`Round length` = grp, `Position-votes r` = r) |>
  kable(caption = "Raw within-group correlation by round length.")

Raw within-group correlation by round length.

Round length	Rounds	Songs	Position-votes r
13-16	27	415	0.023
17-18	26	462	-0.105
19-20	2	38	-0.175
21-37	9	265	-0.151

The subgroup correlations look suggestive — near zero in the shortest rounds, clearly negative in the longest — and the robustness table below shows the effect reaching significance in long rounds but not in standard-length ones. Neither interaction test is significant, though, so this contrast is not established. Comparing two subgroup p-values is not a test of whether they differ; the standard-length estimate (-0.1620) has a confidence interval that comfortably contains the pooled estimate. The honest reading is that this dataset cannot resolve whether the effect scales with playlist length. It is a hypothesis for a larger sample, not a finding.

10 Robustness

Code

kable(robust, caption = "Robustness of the linear position effect. The first five rows are in within-round SD units; the last two change the outcome scale.")

Robustness of the linear position effect. The first five rows are in within-round SD units; the last two change the outcome scale.

Specification	Estimate	95% CI	p	Songs
Main: all songs, round FE	-0.2504	[-0.444, -0.057]	0.012	1180
Excluding duplicated positions	-0.2509	[-0.446, -0.055]	0.013	1176
Standard-length rounds only (13-20)	-0.1620	[-0.381, +0.057]	0.143	915
Long rounds only (21-37)	-0.5678	[-1.017, -0.119]	0.019	265
Excluding rounds with any missing votes	-0.2487	[-0.448, -0.050]	0.015	1151
Adding submitter fixed effects	-0.3098	[-0.507, -0.113]	0.003	1180
Outcome = votes / round mean	-0.1230	[-0.220, -0.026]	0.013	1180
Outcome = raw votes	-0.9369	[-1.747, -0.127]	0.024	1180

The effect survives excluding duplicated positions, excluding rounds with missing votes, adding submitter fixed effects (where it strengthens), and both alternative outcome scales. It is not significant in the standard-length subgroup — see Section 9 for why that should not be read as a contradiction.

Code

analysis |>
  ggplot(aes(order, votes)) +
  geom_point(alpha = 0.5, size = 0.7) +
  geom_smooth(method = "lm", se = FALSE, colour = "#e6550d", linewidth = 0.6) +
  facet_wrap(~round, scales = "free", ncol = 6) +
  labs(title = "Votes vs. playlist position, by round",
       x = "Playlist position", y = "Votes") +
  theme(strip.text = element_text(size = 6), axis.text = element_text(size = 5))

Figure 8: Every round separately: votes against playlist position with a linear fit. Most panels tilt downward, but many tilt up — the effect is a weak tendency, invisible in any single round.

11 How precisely is this pinned down?

A detected effect deserves the same scrutiny as a null one. The smallest true first-to-last effect this dataset would detect 80% of the time (two-sided, α = 0.05) is roughly 2.8 × SE.

Code

tibble(
  Quantity = c(
    "Cluster-robust SE of the first-to-last effect",
    "Minimum detectable effect (80% power)",
    "  in votes, at the median within-round SD",
    "Observed effect",
    "Ratio of observed effect to MDE"
  ),
  Value = c(
    sprintf("%.4f SD", lin_row$std.error),
    sprintf("%.3f SD", mde_sd),
    sprintf("%.2f votes", mde_sd * sd_typ),
    sprintf("%+.4f SD (%.2f votes)", lin_row$estimate,
            lin_row$estimate * sd_typ),
    sprintf("%.2f", abs(lin_row$estimate) / mde_sd)
  )
) |>
  kable(caption = "Precision of the test.")

Precision of the test.

Quantity	Value
Cluster-robust SE of the first-to-last effect	0.0969 SD
Minimum detectable effect (80% power)	0.271 SD
in votes, at the median within-round SD	1.12 votes
Observed effect	-0.2504 SD (-1.03 votes)
Ratio of observed effect to MDE	0.92

The observed effect is 0.92× the minimum detectable effect — it clears significance, but only just. Two consequences worth taking seriously:

• The direction is better established than the magnitude. The confidence interval spans 0.13 to 1.75 votes, a factor of 14 from end to end. “About one vote” is the right level of precision to quote.
• Effects estimated near the detection threshold are biased upward in magnitude. Selection on significance means the true effect is, if anything, likely smaller than 0.94 votes. Treat this as an upper-ish bound pending more rounds.

12 Conclusions

Primacy, weakly. Earlier playlist slots earn more votes. Over the full playlist the gap is about 0.9 vote (11% of an average song’s total), the average within-round correlation is -0.086, and position accounts for roughly 0.6% of the variation in votes. It is real by an exact permutation test (p = 0.0094) and survives every robustness check that keeps the full sample, but it is small enough that song quality, round theme, and taste overwhelm it in any individual round.

Not curvilinear. No U, no inverted U. The quadratic term is flat (p = 0.93), the GAM shrinks to 1.05 effective degrees of freedom, and the fitted parabola’s vertex lands outside the data. A steady decline fits.

Not an edge effect. Opening and closing tracks are not special; the decline is spread evenly across the playlist.

Practical upshot. Drawing the last slot instead of the first costs about a vote. That can flip a near-tie, and over 64 rounds it is a consistent tilt — but it will not decide a round on its own.

12.1 Caveats

• Voters may not listen in playlist order. Shuffle play, listening across several sittings, or voting from a personal shortlist all weaken the link between playlist slot and listening order. Such attenuation biases toward zero, which cuts in an interesting direction here: the true effect among voters who do listen in order is probably larger than 0.9 vote.
• Zero-sum votes. Within a round, votes must sum to a fixed pool, so mean vote_z is mechanically zero and songs within a round are negatively dependent. That is why inference rests on the within-round permutation test rather than model-based p-values alone.
• One observation per song. Vote totals are analyzed, not individual ballots. Ballot-level data would sharpen every test here and would let you ask whether specific voters drive the effect.
• Multiplicity. Many tests are reported, and the headline p is around 0.01 rather than 0.001. Two things argue against dismissing it as a multiplicity artifact: primacy was a pre-specified direction, and independent routes (per-round meta-analysis, fixed-effects regression, permutation, GAM) agree. Still, replication on future rounds is the real test.
• Position is not randomized by design. The randomization check is reassuring but observational; it cannot rule out an unmeasured association between slot and song quality.
• Round length moderation is unresolved, as Section 9 explains.

Code

sessionInfo()

R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.6.1

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] knitr_1.51      mgcv_1.9-3      nlme_3.1-168    broom_1.0.12   
 [5] lubridate_1.9.5 forcats_1.0.1   stringr_1.6.0   dplyr_1.2.0    
 [9] purrr_1.2.1     readr_2.1.6     tidyr_1.3.2     tibble_3.3.1   
[13] ggplot2_4.0.2   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] generics_0.1.4     stringi_1.8.7      lattice_0.22-7     hms_1.1.4         
 [5] digest_0.6.39      magrittr_2.0.4     evaluate_1.0.5     grid_4.5.2        
 [9] timechange_0.4.0   RColorBrewer_1.1-3 fastmap_1.2.0      jsonlite_2.0.0    
[13] Matrix_1.7-4       backports_1.5.0    scales_1.4.0       cli_3.6.6         
[17] rlang_1.2.0        crayon_1.5.3       bit64_4.6.0-1      splines_4.5.2     
[21] withr_3.0.2        yaml_2.3.12        otel_0.2.0         parallel_4.5.2    
[25] tools_4.5.2        tzdb_0.5.0         vctrs_0.7.1        R6_2.6.1          
[29] lifecycle_1.0.5    htmlwidgets_1.6.4  bit_4.6.0          vroom_1.7.0       
[33] pkgconfig_2.0.3    pillar_1.11.1      gtable_0.3.6       glue_1.8.0        
[37] xfun_0.56          tidyselect_1.2.1   farver_2.1.2       htmltools_0.5.9   
[41] labeling_0.4.3     rmarkdown_2.30     compiler_4.5.2     S7_0.2.1