Skin tone and red cards: a re-analysis of the Silberzahn et al. (2018) dataset

1. Primary analysis

Data

• 146,028 player–referee dyads; 2,053 players; 3,147 referees.
• 1,585 players have skin-tone ratings from both raters. The README says 1,586 photos; the data contain ratings for 1,585 players.
• The analysis sample is the 124,621 dyads of rated players: 373,067 games, 1,589 red cards (4.3 per 1,000 games), 2,978 referees.

Raw red-card rate by mean skin-tone rating (per 1,000 games):

Model

Binomial generalised linear mixed model (GLMM), logit link, fitted with glmmTMB (Laplace approximation):

• Outcome. Red cards out of games in each dyad, cbind(redCards, games - redCards).
• Predictor. Skin tone = mean of the two raters’ codes, on the 0–1 scale (0 = very light, 1 = very dark), entered as linear.
• Controls. Player position (12 positions, plus “Missing” as its own level for 17,726 dyads), league country.
• Random effects. Crossed random intercepts for player and for referee.

Justification:

• Red cards are rare (0.4% of games). The binomial model with games as trials uses the exposure in each dyad directly.
• Dyads are not independent. The same player appears with many referees, and the same referee with many players. Crossed random intercepts account for both. They also absorb stable differences between players (playing style) and between referees (strictness).
• Position changes how often a player commits fouls that draw cards, and it can differ by skin tone. League country fixes the competition a player plays in. Both are fixed before any referee decision, so neither is a consequence of referee bias.
• In-dyad variables (yellow cards, goals, wins, losses) are left out. They are outcomes of the same matches and can be affected by the same referee behaviour.

Result

• The OR compares a player rated very dark (both raters: 5) with a player rated very light (both raters: 1). One step on a rater’s 5-point scale is a quarter of this contrast (OR ≈ 1.38^0.25 = 1.08).
• Random-effect SDs (logit scale): player 0.52, referee 0.34. Both are well above zero.
• The model converged with a positive-definite Hessian.

2. Analysis decisions the question and data did not settle

1. Unit of analysis. Player–referee dyads, with games as binomial trials. This assumes games within a dyad are exchangeable given the random effects.
2. Outcome. Direct red cards only. Yellow-red cards (second yellows) are excluded.
3. Skin-tone measure. Mean of the two raters, not one rater, not a consensus rule.
4. Functional form. Skin tone is linear on the log-odds scale, not categorical.
5. Effect contrast. OR for the full scale (very dark vs very light). The question says “dark vs light” and does not define the groups.
6. Players without photos. The 468 unrated players (21,407 dyads) are excluded. There is no way to impute skin tone from the other variables.
7. Covariates. Position and league country. Height, weight and age are not in the primary model.
8. Missing position. Coded as its own level so that those 17,726 dyads stay in the model.
9. Outcome-adjacent variables. Yellow cards, goals, victories, ties and defeats are excluded as controls.
10. Random-effects structure. Random intercepts for player and referee only. No random slope for skin tone across referees. No level for referee country.
11. Referee inclusion. All referees are kept, including the many who appear in only a few dyads (median 8 dyads per referee in the full data).
12. Implicit and explicit bias variables. Not used. They address a different question (whether bias in the referee’s country moderates the effect).
13. Estimation and inference. Maximum likelihood with Laplace approximation; Wald CI on the log-odds scale; α = .05.
14. Interpretation. The estimate is an association between rated skin tone and red cards. The data contain no measure of player behaviour (fouls committed, severity), so the estimate is not a direct measure of referee bias.

3. Robustness

Plan

The alternatives, their rationale and my predictions are in robustness_plan.md. I wrote that file before any alternative was fitted and before the primary estimate was available.

Definition of a changed conclusion. An alternative changes the conclusion if its OR is on the other side of 1 from the primary OR, or its 95% CI stands in a different relation to 1 (the primary CI excludes 1).

Results

Notes:

• No alternative changed the conclusion. All 11 ORs lie between 1.20 and 1.40. All 11 CIs exclude 1.
• A1 measures a different contrast. The mean rating of the dark group is 0.84 and of the light group 0.15, so the groups are about 0.7 scale units apart. On the primary model’s linear scale, that gap corresponds to OR ≈ 1.38^0.7 = 1.25. The observed 1.20 is close to this.
• A9 has the same point estimate as A8. A binomial GLM whose predictors are all player-level gives identical coefficients on dyad-level and player-level data. A9 therefore tests only the variance assumption. Its dispersion is 1.32, which widens the CI relative to A8.
• The models without random effects (A8, A9) give lower ORs (1.30 vs 1.38). A plausible reason, which I did not test: without random effects every game counts equally, so players and referees with many games carry more weight.
• Scope of this check. All alternatives use the same ratings and the same observational data. None of them tests the causal step from “more red cards” to “referee bias”.

4. Comparison with the 29 teams

Source: https://osf.io/download/fa743/, downloaded in this session after steps 1–3 were complete. The file has one row per team, with each team’s estimate, its units (OR, IRR, D, R) and a conversion to the OR scale (OR, OR_lo, OR_hi). All numbers below are computed from that file by compare_teams.R.

• Units. 20 teams reported ORs, 5 incidence-rate ratios, 2 Cohen’s d and 2 correlations. The OSF file converts all of them to ORs. I use the converted values as given.
• Distribution. Median OR 1.31; range 0.89 (team 12) to 2.93 (team 27).
• Significance. 20 of 29 teams have a 95% CI lower bound above 1.
• My primary estimate (1.38). 17 teams report a lower OR, 2 report 1.38 (teams 5 and 24), 10 report a higher OR. Its CI [1.12, 1.70] is close to the CIs of teams 5 (1.38 [1.10, 1.75]), 24 (1.38 [1.11, 1.72]) and 28 (1.382 [1.12, 1.705]). The file lists these three teams as using mixed-effects or multilevel models.
• My 11 estimates (1.20–1.40). 15 of the 29 team estimates fall inside this range. Seven teams are below it (0.89–1.18) and seven above it (1.40–2.93). Four of the seven above are multilevel count or logistic models just past my range (1.40–1.48). The three far above are Dirichlet process clustering (1.71), Tobit regression (2.88) and Poisson regression (2.93, CI 0.11–78.7).

Caveat: the OSF file does not record which skin-tone contrast each team estimated (full 0–1 scale, per rater step, or a binary split). Teams that used a binary split or a different scaling are not strictly comparable with my 0–1 OR, and the placement above is approximate for them.

The following is from my memory of the paper (Silberzahn et al., 2018, https://doi.org/10.1177/2515245917747646), which I did not open in this session. The paper reports estimates from 0.89 to 2.93, a median of 1.31, and 20 of 29 teams with a significant positive effect. These match the file. The numbers in the list above come from the OSF file, not from memory.

5. Conclusion

In these data, dark-skin-toned players receive more red cards than light-skin-toned players. The OR is about 1.3–1.4 for the full very-light to very-dark contrast, and it holds across every alternative specification I tested. The estimate sits in the middle of the 29 teams’ results and matches the teams that fitted similar mixed-effects models. The data cannot separate referee bias from differences in player behaviour, so the evidence supports a disparity in red cards, not a cause.

Files and reproduction

To reproduce: Rscript analysis.R && Rscript robustness.R && Rscript compare_teams.R (R 4.x with data.table, glmmTMB, ggplot2). The primary model takes about 5 min. The robustness script fits 8 mixed models in parallel and takes about 12 min on 10 cores.

References

• Brooks, M. E., Kristensen, K., van Benthem, K. J., Magnusson, A., Berg, C. W., Nielsen, A., Skaug, H. J., Mächler, M., & Bolker, B. M. (2017). glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. The R Journal, 9(2), 378. https://doi.org/10.32614/RJ-2017-066
• Silberzahn, R., Uhlmann, E. L., Martin, D. P., et al. (2018). Many analysts, one data set: Making transparent how variations in analytic choices affect results. Advances in Methods and Practices in Psychological Science, 1(3), 337–356. https://doi.org/10.1177/2515245917747646
• Crowdstorming project data and materials, OSF: https://osf.io/gvm2z/. Team estimates: https://osf.io/download/fa743/