Bachelor_Thesis_Analysis

1. Setup

2. Data Import

# list.files(here("Raw Data needed"))

# Manually collected special-effects data
effects_raw <- read_csv2(
  here("Raw Data needed", "performance_effects.csv"),
  show_col_types = FALSE
)

# Song and musical-feature data
songs_raw <- read_csv2(
  here("Raw Data needed", "complete_songs.csv"),
  show_col_types = FALSE
)

# Performance-order data
performance_order_raw <- read_csv2(
  here("Raw Data needed", "performance_order.csv"),
  show_col_types = FALSE
)

# Historical voting data through 2024
votes_old_raw <- read_csv(
  here("Raw Data needed", "Votes_1957_2024.csv"),
  show_col_types = FALSE
)

# Voting extension for 2025
votes_2025_raw <- read_csv2(
  here("Raw Data needed", "Votes_2025.csv"),
  show_col_types = FALSE
)

dim(effects_raw)

## [1] 1486   13

dim(songs_raw)

## [1] 1754   37

dim(performance_order_raw)

## [1] 1756    7

dim(votes_old_raw)

## [1] 51452     7

dim(votes_2025_raw)

## [1] 962   8

nrow(votes_2025_raw)

## [1] 962

names(votes_2025_raw)

## [1] "year"            "round"           "from_country"    "to_country"     
## [5] "pairwise_points" "tele_points"     "jury_points"     "...8"

3. Raw Data Inspection

glimpse(effects_raw)

## Rows: 1,486
## Columns: 13
## $ year                  <dbl> 2025, 2025, 2025, 2025, 2025, 2025, 2025, 2025, …
## $ to_country            <chr> "Norway", "Luxembourg", "Estonia", "Israel", "Li…
## $ song                  <chr> "Lighter", "La poupée monte le son", "Espresso M…
## $ artist_name           <chr> "Kyle Alessandro", "Laura Thorn", "Tommy Cash", …
## $ round                 <chr> "finale", "finale", "finale", "finale", "finale"…
## $ performance_order     <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1…
## $ fire                  <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, …
## $ smoke                 <dbl> 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, …
## $ wind_machine          <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, …
## $ confetti              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ special_effects_count <dbl> 2, 0, 0, 0, 2, 0, 0, 0, 1, 1, 1, 0, 1, 0, 2, 0, …
## $ coding_complete       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ Notes_Data_collection <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …

glimpse(songs_raw)

## Rows: 1,754
## Columns: 37
## $ year                   <dbl> 2025, 2025, 2025, 2025, 2025, 2025, 2025, 2025,…
## $ to_country             <chr> "Albania", "Armenia", "Australia", "Austria", "…
## $ song                   <chr> "Zjerm", "SURVIVOR", "Go-Jo", "Wasted Love", "R…
## $ artist_name            <chr> "Shkodra Elektronike", "PARG", "Milkshake Man",…
## $ artist_size            <dbl> 2, 1, 1, 1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1,…
## $ gender                 <chr> "Mixed", "Male", "Male", "Male", "All-male", "M…
## $ main_singers           <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ backing_singers        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,…
## $ backing_dancers        <dbl> 0, 0, 2, 0, 3, 3, 4, 4, 4, 4, 5, 0, 0, 4, 3, 1,…
## $ backing_instruments    <dbl> 1, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ number_of_performers   <dbl> 2, 1, 3, 1, 6, 4, 5, 5, 5, 5, 6, 1, 1, 5, 6, 2,…
## $ country_debut          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ returning_artist       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ artist_nationality     <chr> "Albania", "Armenia", "Australia", "Austria", "…
## $ nationality_match      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1,…
## $ LGBTQ_category         <chr> NA, NA, NA, "Queer", NA, "Queer", "Gay", NA, "Q…
## $ LGBTQ_out_time         <chr> NA, NA, NA, "Openly", NA, "Openly", "Openly", N…
## $ performance_notes      <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ genre                  <chr> "Pop", "Rock", "Pop", "Pop", "Pop", "Dance", "P…
## $ returning              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ host                   <chr> "Switzerland", "Switzerland", "Switzerland", "S…
## $ language               <chr> "Albanian", "Armenian, English", "English", "En…
## $ dialect_sign_imaginary <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,…
## $ language_count         <dbl> 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 2, 1, 1,…
## $ has_english            <dbl> 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0,…
## $ tempo                  <dbl> 9493670654, 1304347839, 1401869202, 625, 141509…
## $ key                    <dbl> 0, 2, 11, 9, 2, 6, 7, 10, 1, 0, 4, 4, 7, 5, 3, …
## $ mode                   <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0,…
## $ danceability           <chr> "0.783782766", "0.746351126", "0.620062281", "0…
## $ energy                 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ loudness               <chr> "-11.70", "-8.81", "-3.61", "-3.82", "-4.57", "…
## $ speechiness            <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ acousticness           <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ instrumentalness       <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ liveness               <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ valence                <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ time_signature         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…

glimpse(performance_order_raw)

## Rows: 1,756
## Columns: 7
## $ year                       <dbl> 1956, 1956, 1956, 1956, 1956, 1956, 1956, 1…
## $ to_country                 <chr> "Switzerland", "Netherlands", "Belgium", "G…
## $ performer                  <chr> "Lys Assia", "Jetty Paerl", "Fud Leclerc", …
## $ song                       <chr> "Refrain", "De Vogels Van Holland", "Messie…
## $ semifinal                  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ PerformanceOrder_final     <dbl> 2, 1, 3, 4, 5, 13, 7, 8, 10, 11, 12, 14, 6,…
## $ PerformanceOrder_semifinal <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…

glimpse(votes_old_raw)

## Rows: 51,452
## Columns: 7
## $ year            <dbl> 1957, 1957, 1957, 1957, 1957, 1957, 1957, 1957, 1957, …
## $ round           <chr> "final", "final", "final", "final", "final", "final", …
## $ from_country    <chr> "Austria", "Austria", "Austria", "Austria", "Austria",…
## $ to_country      <chr> "Netherlands", "France", "Denmark", "Luxembourg", "Ger…
## $ pairwise_points <dbl> 6, 0, 0, 3, 0, 0, 1, 0, 0, 5, 2, 0, 0, 1, 1, 1, 0, 0, …
## $ tele_points     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ jury_points     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…

glimpse(votes_2025_raw)

## Rows: 962
## Columns: 8
## $ year            <dbl> 2025, 2025, 2025, 2025, 2025, 2025, 2025, 2025, 2025, …
## $ round           <chr> "final", "final", "final", "final", "final", "final", …
## $ from_country    <chr> "Luxembourg", "Estonia", "Israel", "Lithuania", "Spain…
## $ to_country      <chr> "Norway", "Norway", "Norway", "Norway", "Norway", "Nor…
## $ pairwise_points <dbl> 0, 0, 0, 2, 2, 10, 0, 0, 14, 1, 0, 0, 1, 5, 1, 4, 4, 0…
## $ tele_points     <dbl> 0, 0, 0, 2, 2, 10, 0, 0, 8, 1, 0, 0, 1, 3, 0, 4, 4, 0,…
## $ jury_points     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 2, 1, 0, 0, 0, …
## $ ...8            <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…

# Number of rows with and without a year
effects_raw |>
  summarise(
    total_rows = n(),
    rows_with_year = sum(!is.na(year)),
    rows_without_year = sum(is.na(year))
  )

## # A tibble: 1 × 3
##   total_rows rows_with_year rows_without_year
##        <int>          <int>             <int>
## 1       1486            413              1073

songs_raw |>
  summarise(
    total_rows = n(),
    rows_with_year = sum(!is.na(year)),
    rows_without_year = sum(is.na(year))
  )

## # A tibble: 1 × 3
##   total_rows rows_with_year rows_without_year
##        <int>          <int>             <int>
## 1       1754           1754                 0

performance_order_raw |>
  summarise(
    total_rows = n(),
    rows_with_year = sum(!is.na(year)),
    rows_without_year = sum(is.na(year))
  )

## # A tibble: 1 × 3
##   total_rows rows_with_year rows_without_year
##        <int>          <int>             <int>
## 1       1756           1756                 0

votes_old_raw |>
  summarise(
    total_rows = n(),
    rows_with_year = sum(!is.na(year)),
    rows_without_year = sum(is.na(year))
  )

## # A tibble: 1 × 3
##   total_rows rows_with_year rows_without_year
##        <int>          <int>             <int>
## 1      51452          51452                 0

votes_2025_raw |>
  summarise(
    total_rows = n(),
    rows_with_year = sum(!is.na(year)),
    rows_without_year = sum(is.na(year))
  )

## # A tibble: 1 × 3
##   total_rows rows_with_year rows_without_year
##        <int>          <int>             <int>
## 1        962            962                 0

songs_raw |>
  select(year, to_country, song, tempo, energy, danceability, valence) |>
  slice_head(n = 20)

## # A tibble: 20 × 7
##     year to_country song                       tempo energy danceability valence
##    <dbl> <chr>      <chr>                      <dbl> <chr>  <chr>        <chr>  
##  1  2025 Albania    Zjerm                 9493670654 <NA>   0.783782766  <NA>   
##  2  2025 Armenia    SURVIVOR              1304347839 <NA>   0.746351126  <NA>   
##  3  2025 Australia  Go-Jo                 1401869202 <NA>   0.620062281  <NA>   
##  4  2025 Austria    Wasted Love                  625 <NA>   0.660020048  <NA>   
##  5  2025 Azerbaijan Run With U            1415094299 <NA>   0.767386435  <NA>   
##  6  2025 Belgium    Strobe Lights         8522727203 <NA>   0.677794755  <NA>   
##  7  2025 Croatia    Poison Cake            728155365 <NA>   0.899225528  <NA>   
##  8  2025 Cyprus     Shh                    116279068 <NA>   0.97480672   <NA>   
##  9  2025 Czechia    Kiss Kiss Goodbye      681818161 <NA>   0.351895209  <NA>   
## 10  2025 Denmark    Hallucination         7731958771 <NA>   0.980092473  <NA>   
## 11  2025 Estonia    Espresso Macchiato           120 <NA>   0.886701067  <NA>   
## 12  2025 Finland    ICH KOMME             7894736481 <NA>   0.448452318  <NA>   
## 13  2025 France     maman                 5357143021 <NA>   0.266809981  <NA>   
## 14  2025 Georgia    Freedom                131578949 <NA>   0.757988785  <NA>   
## 15  2025 Germany    Baller                7009346008 <NA>   0.994010366  <NA>   
## 16  2025 Greece     Asteromáta            1219512177 <NA>   0.661153399  <NA>   
## 17  2025 Iceland    RÓA                   1376146851 <NA>   0.990902424  <NA>   
## 18  2025 Ireland    Laika Party           5357143021 <NA>   0.53104055   <NA>   
## 19  2025 Israel     New Day Will Rise     8823529053 <NA>   0.588710481  <NA>   
## 20  2025 Italy      Volevo Essere Un Duro 1485148468 <NA>   0.548016483  <NA>

songs_raw |>
  summarise(
    tempo_missing = sum(is.na(tempo)),
    energy_missing = sum(is.na(energy)),
    danceability_missing = sum(is.na(danceability)),
    valence_missing = sum(is.na(valence))
  )

## # A tibble: 1 × 4
##   tempo_missing energy_missing danceability_missing valence_missing
##           <int>          <int>                <int>           <int>
## 1             0            211                    0             211

songs_raw |>
  summarise(
    tempo_min = min(tempo, na.rm = TRUE),
    tempo_max = max(tempo, na.rm = TRUE)
  )

## # A tibble: 1 × 2
##   tempo_min  tempo_max
##       <dbl>      <dbl>
## 1        75 9740259552

songs_raw |>
  filter(year >= 2016, year <= 2025) |>
  group_by(year) |>
  summarise(
    n = n(),
    tempo_missing = sum(is.na(tempo)),
    energy_missing = sum(is.na(energy)),
    danceability_missing = sum(is.na(danceability)),
    valence_missing = sum(is.na(valence)),
    .groups = "drop"
  )

## # A tibble: 9 × 6
##    year     n tempo_missing energy_missing danceability_missing valence_missing
##   <dbl> <int>         <int>          <int>                <int>           <int>
## 1  2016    42             0              0                    0               0
## 2  2017    42             0              0                    0               0
## 3  2018    43             0              0                    0               0
## 4  2019    41             0              1                    0               1
## 5  2021    39             0              0                    0               0
## 6  2022    40             0              0                    0               0
## 7  2023    37             0              0                    0               0
## 8  2024    37             0              0                    0               0
## 9  2025    37             0             37                    0              37

songs_raw |>
  filter(year >= 2016, year <= 2025) |>
  group_by(year) |>
  summarise(
    tempo_min = min(tempo, na.rm = TRUE),
    tempo_max = max(tempo, na.rm = TRUE),
    .groups = "drop"
  )

## # A tibble: 9 × 3
##    year tempo_min  tempo_max
##   <dbl>     <dbl>      <dbl>
## 1  2016     11497     205008
## 2  2017      8488     181758
## 3  2018     10908     189743
## 4  2019      1009 9615384674
## 5  2021      8056     175719
## 6  2022     76905     169985
## 7  2023     13796     160071
## 8  2024      1351     160068
## 9  2025       120 9493670654

# Check duplicate year-country combinations in performance-level datasets
effects_raw |>
  filter(!is.na(year)) |>
  count(year, to_country) |>
  filter(n > 1)

## # A tibble: 0 × 3
## # ℹ 3 variables: year <dbl>, to_country <chr>, n <int>

songs_raw |>
  count(year, to_country) |>
  filter(n > 1)

## # A tibble: 7 × 3
##    year to_country      n
##   <dbl> <chr>       <int>
## 1  1956 Belgium         2
## 2  1956 France          2
## 3  1956 Germany         2
## 4  1956 Italy           2
## 5  1956 Luxembourg      2
## 6  1956 Netherlands     2
## 7  1956 Switzerland     2

performance_order_raw |>
  count(year, to_country) |>
  filter(n > 1)

## # A tibble: 5 × 3
##    year to_country      n
##   <dbl> <chr>       <int>
## 1  1956 Belgium         2
## 2  1956 France          2
## 3  1956 Germany         2
## 4  1956 Italy           2
## 5  1956 Netherlands     2

# Check duplicate voting observations
votes_old_raw |>
  count(year, round, from_country, to_country) |>
  filter(n > 1)

## # A tibble: 0 × 5
## # ℹ 5 variables: year <dbl>, round <chr>, from_country <chr>, to_country <chr>,
## #   n <int>

votes_2025_raw |>
  count(year, round, from_country, to_country) |>
  filter(n > 1)

## # A tibble: 0 × 5
## # ℹ 5 variables: year <dbl>, round <chr>, from_country <chr>, to_country <chr>,
## #   n <int>

# Inspect round labels
unique(effects_raw$round)

## [1] "finale" NA

unique(votes_old_raw$round)

## [1] "final"        "semi-final"   "semi-final-1" "semi-final-2"

unique(votes_2025_raw$round)

## [1] "final"

4. Data Cleaning

4.1 Special Effects Data

effects_clean <- effects_raw |>
  # Standardize column names
  clean_names() |>
  
  # Remove empty rows created during the Excel-to-CSV export
  filter(!is.na(year)) |>
  
  # Standardize the Grand Final label
  mutate(
    round = if_else(round == "finale", "final", round)
  )

dim(effects_clean)

## [1] 413  13

unique(effects_clean$round)

## [1] "final"

effects_clean |>
  count(year) |>
  arrange(desc(year))

## # A tibble: 16 × 2
##     year     n
##    <dbl> <int>
##  1  2025    26
##  2  2024    26
##  3  2023    26
##  4  2022    25
##  5  2021    26
##  6  2019    26
##  7  2018    26
##  8  2017    26
##  9  2016    26
## 10  2015    27
## 11  2014    26
## 12  2013    26
## 13  2012    26
## 14  2011    25
## 15  2010    25
## 16  2009    25

# Check allowed values for each binary effect variable
unique(effects_clean$fire)

## [1]  1  0 NA

unique(effects_clean$smoke)

## [1]  1  0 NA

unique(effects_clean$wind_machine)

## [1]  0  1 NA

unique(effects_clean$confetti)

## [1]  0 NA  1

# Check whether the calculated total matches the four component variables
effects_clean |>
  mutate(
    effects_sum_check = fire + smoke + wind_machine + confetti
  ) |>
  filter(
    !is.na(effects_sum_check),
    special_effects_count != effects_sum_check
  ) |>
  select(
    year,
    to_country,
    song,
    fire,
    smoke,
    wind_machine,
    confetti,
    special_effects_count,
    effects_sum_check
  )

## # A tibble: 0 × 9
## # ℹ 9 variables: year <dbl>, to_country <chr>, song <chr>, fire <dbl>,
## #   smoke <dbl>, wind_machine <dbl>, confetti <dbl>,
## #   special_effects_count <dbl>, effects_sum_check <dbl>

effects_clean |>
  filter(
    is.na(fire) |
    is.na(smoke) |
    is.na(wind_machine) |
    is.na(confetti) |
    is.na(special_effects_count)
  ) |>
  select(
    year,
    to_country,
    song,
    fire,
    smoke,
    wind_machine,
    confetti,
    special_effects_count,
    coding_complete,
    notes_data_collection
  )

## # A tibble: 1 × 10
##    year to_country song   fire smoke wind_machine confetti special_effects_count
##   <dbl> <chr>      <chr> <dbl> <dbl>        <dbl>    <dbl>                 <dbl>
## 1  2024 Netherlan… Euro…    NA    NA           NA       NA                    NA
## # ℹ 2 more variables: coding_complete <dbl>, notes_data_collection <chr>

4.2 Voting Data

# Clean historical voting data
votes_old_clean <- votes_old_raw |>
  clean_names()

# Clean 2025 voting extension and remove the empty export column
votes_2025_clean <- votes_2025_raw |>
  clean_names() |>
  select(-x8)

# Compare the column names of both datasets
names(votes_old_clean)

## [1] "year"            "round"           "from_country"    "to_country"     
## [5] "pairwise_points" "tele_points"     "jury_points"

names(votes_2025_clean)

## [1] "year"            "round"           "from_country"    "to_country"     
## [5] "pairwise_points" "tele_points"     "jury_points"

# Compare their dimensions
dim(votes_old_clean)

## [1] 51452     7

dim(votes_2025_clean)

## [1] 962   7

# Check 1: pairwise points must equal televote + jury points
votes_2025_clean |>
  filter(pairwise_points != tele_points + jury_points)

## # A tibble: 0 × 7
## # ℹ 7 variables: year <dbl>, round <chr>, from_country <chr>, to_country <chr>,
## #   pairwise_points <dbl>, tele_points <dbl>, jury_points <dbl>

# Check 2: every finalist should have exactly 37 voting-source rows
votes_2025_clean |>
  count(to_country, name = "number_of_voting_rows") |>
  filter(number_of_voting_rows != 37)

## # A tibble: 0 × 2
## # ℹ 2 variables: to_country <chr>, number_of_voting_rows <int>

# Check 3: calculate total points received by each finalist
votes_2025_clean |>
  group_by(to_country) |>
  summarise(
    tele_total = sum(tele_points),
    jury_total = sum(jury_points),
    total_points = sum(pairwise_points),
    .groups = "drop"
  ) |>
  arrange(desc(total_points))

## # A tibble: 26 × 4
##    to_country  tele_total jury_total total_points
##    <chr>            <dbl>      <dbl>        <dbl>
##  1 Austria            178        258          436
##  2 Israel             297         60          357
##  3 Estonia            258         98          356
##  4 Sweden             195        126          321
##  5 Italy               97        159          256
##  6 Greece             126        105          231
##  7 France              50        180          230
##  8 Albania            173         45          218
##  9 Ukraine            158         60          218
## 10 Switzerland          0        214          214
## # ℹ 16 more rows

# Combine historical voting data with the 2025 extension
votes_clean <- bind_rows(
  votes_old_clean,
  votes_2025_clean
)

# Check the dimensions of the combined dataset
dim(votes_clean)

## [1] 52414     7

# Check the year range
range(votes_clean$year, na.rm = TRUE)

## [1] 1957 2025

# Check that 2025 contains exactly 962 observations
votes_clean |>
  filter(year == 2025) |>
  summarise(
    number_of_rows = n()
  )

## # A tibble: 1 × 1
##   number_of_rows
##            <int>
## 1            962

# Check the available round labels
unique(votes_clean$round)

## [1] "final"        "semi-final"   "semi-final-1" "semi-final-2"

votes_clean |>
  filter(
    round == "final",
    year >= 2009,
    year <= 2025
  ) |>
  group_by(year) |>
  summarise(
    rows = n(),
    tele_missing = sum(is.na(tele_points)),
    jury_missing = sum(is.na(jury_points)),
    tele_available = sum(!is.na(tele_points)),
    jury_available = sum(!is.na(jury_points)),
    .groups = "drop"
  )

## # A tibble: 16 × 6
##     year  rows tele_missing jury_missing tele_available jury_available
##    <dbl> <int>        <int>        <int>          <int>          <int>
##  1  2009  1025         1025         1025              0              0
##  2  2010   950          950          950              0              0
##  3  2011  1050         1050         1050              0              0
##  4  2012  1066         1066         1066              0              0
##  5  2013   988          988          988              0              0
##  6  2014   936          936          936              0              0
##  7  2015  1053         1053         1053              0              0
##  8  2016  1066            0            0           1066           1066
##  9  2017  1066            0            0           1066           1066
## 10  2018  1092            0            0           1092           1092
## 11  2019  1040            0            0           1040           1040
## 12  2021   988            0            0            988            988
## 13  2022   975            0            0            975            975
## 14  2023   962            0           26            962            936
## 15  2024   925            0            0            925            925
## 16  2025   962            0            0            962            962

4.3 Musical Features Data

songs_clean <- songs_raw |>
  clean_names() |>
  mutate(
    energy = as.numeric(energy),
    danceability = as.numeric(danceability),
    valence = as.numeric(valence)
  )

# Check the cleaned variable types
songs_clean |>
  select(tempo, energy, danceability, valence) |>
  glimpse()

## Rows: 1,754
## Columns: 4
## $ tempo        <dbl> 9493670654, 1304347839, 1401869202, 625, 1415094299, 8522…
## $ energy       <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ danceability <dbl> 0.7837828, 0.7463511, 0.6200623, 0.6600200, 0.7673864, 0.…
## $ valence      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…

tempo_check <- songs_clean |>
  mutate(
    tempo_proposed = case_when(
      tempo <= 220 ~ tempo,
      tempo / 10 <= 220 ~ tempo / 10,
      tempo / 100 <= 220 ~ tempo / 100,
      tempo / 1000 <= 220 ~ tempo / 1000,
      tempo / 10000 <= 220 ~ tempo / 10000,
      tempo / 100000 <= 220 ~ tempo / 100000,
      tempo / 1000000 <= 220 ~ tempo / 1000000,
      tempo / 10000000 <= 220 ~ tempo / 10000000,
      tempo / 100000000 <= 220 ~ tempo / 100000000,
      TRUE ~ NA_real_
    )
  )

# Inspect the first 20 original and proposed tempo values
tempo_check |>
  select(year, to_country, song, tempo, tempo_proposed) |>
  slice_head(n = 20)

## # A tibble: 20 × 5
##     year to_country song                       tempo tempo_proposed
##    <dbl> <chr>      <chr>                      <dbl>          <dbl>
##  1  2025 Albania    Zjerm                 9493670654           94.9
##  2  2025 Armenia    SURVIVOR              1304347839          130. 
##  3  2025 Australia  Go-Jo                 1401869202          140. 
##  4  2025 Austria    Wasted Love                  625           62.5
##  5  2025 Azerbaijan Run With U            1415094299          142. 
##  6  2025 Belgium    Strobe Lights         8522727203           85.2
##  7  2025 Croatia    Poison Cake            728155365           72.8
##  8  2025 Cyprus     Shh                    116279068          116. 
##  9  2025 Czechia    Kiss Kiss Goodbye      681818161           68.2
## 10  2025 Denmark    Hallucination         7731958771           77.3
## 11  2025 Estonia    Espresso Macchiato           120          120  
## 12  2025 Finland    ICH KOMME             7894736481           78.9
## 13  2025 France     maman                 5357143021           53.6
## 14  2025 Georgia    Freedom                131578949          132. 
## 15  2025 Germany    Baller                7009346008           70.1
## 16  2025 Greece     Asteromáta            1219512177          122. 
## 17  2025 Iceland    RÓA                   1376146851          138. 
## 18  2025 Ireland    Laika Party           5357143021           53.6
## 19  2025 Israel     New Day Will Rise     8823529053           88.2
## 20  2025 Italy      Volevo Essere Un Duro 1485148468          149.

# Check the resulting range
tempo_check |>
  summarise(
    proposed_min = min(tempo_proposed, na.rm = TRUE),
    proposed_max = max(tempo_proposed, na.rm = TRUE),
    proposed_missing = sum(is.na(tempo_proposed))
  )

## # A tibble: 1 × 3
##   proposed_min proposed_max proposed_missing
##          <dbl>        <dbl>            <int>
## 1         48.5         206.                0

tempo_check |>
  filter(year >= 2016, year <= 2025) |>
  group_by(year) |>
  summarise(
    n = n(),
    tempo_min = min(tempo_proposed, na.rm = TRUE),
    tempo_mean = mean(tempo_proposed, na.rm = TRUE),
    tempo_max = max(tempo_proposed, na.rm = TRUE),
    .groups = "drop"
  )

## # A tibble: 9 × 5
##    year     n tempo_min tempo_mean tempo_max
##   <dbl> <int>     <dbl>      <dbl>     <dbl>
## 1  2016    42      67.9      118.       205.
## 2  2017    42      75.0      123.       182.
## 3  2018    43      74.9      115.       190.
## 4  2019    41      72.0      119.       184.
## 5  2021    39      66.4      114.       176.
## 6  2022    40      76.9      113.       170.
## 7  2023    37      75.0      124.       160.
## 8  2024    37      82.0      124.       160.
## 9  2025    37      53.6       99.9      149.

songs_clean <- songs_clean |>
  mutate(
    tempo = case_when(
      tempo <= 220 ~ tempo,
      tempo / 10 <= 220 ~ tempo / 10,
      tempo / 100 <= 220 ~ tempo / 100,
      tempo / 1000 <= 220 ~ tempo / 1000,
      tempo / 10000 <= 220 ~ tempo / 10000,
      tempo / 100000 <= 220 ~ tempo / 100000,
      tempo / 1000000 <= 220 ~ tempo / 1000000,
      tempo / 10000000 <= 220 ~ tempo / 10000000,
      tempo / 100000000 <= 220 ~ tempo / 100000000,
      TRUE ~ NA_real_
    )
  )

# Check the final musical-feature variable types
songs_clean |>
  select(tempo, energy, danceability, valence) |>
  glimpse()

## Rows: 1,754
## Columns: 4
## $ tempo        <dbl> 94.93671, 130.43478, 140.18692, 62.50000, 141.50943, 85.2…
## $ energy       <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ danceability <dbl> 0.7837828, 0.7463511, 0.6200623, 0.6600200, 0.7673864, 0.…
## $ valence      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…

# Check final tempo range
songs_clean |>
  summarise(
    tempo_min = min(tempo, na.rm = TRUE),
    tempo_max = max(tempo, na.rm = TRUE),
    tempo_missing = sum(is.na(tempo))
  )

## # A tibble: 1 × 3
##   tempo_min tempo_max tempo_missing
##       <dbl>     <dbl>         <int>
## 1      48.5      206.             0

# Overall ranges and missing values of the four H2 variables
songs_clean |>
  summarise(
    tempo_min = min(tempo, na.rm = TRUE),
    tempo_max = max(tempo, na.rm = TRUE),
    tempo_missing = sum(is.na(tempo)),
    
    energy_min = min(energy, na.rm = TRUE),
    energy_max = max(energy, na.rm = TRUE),
    energy_missing = sum(is.na(energy)),
    
    danceability_min = min(danceability, na.rm = TRUE),
    danceability_max = max(danceability, na.rm = TRUE),
    danceability_missing = sum(is.na(danceability)),
    
    valence_min = min(valence, na.rm = TRUE),
    valence_max = max(valence, na.rm = TRUE),
    valence_missing = sum(is.na(valence))
  )

## # A tibble: 1 × 12
##   tempo_min tempo_max tempo_missing energy_min energy_max energy_missing
##       <dbl>     <dbl>         <int>      <dbl>      <dbl>          <int>
## 1      48.5      206.             0     0.0314      0.999            211
## # ℹ 6 more variables: danceability_min <dbl>, danceability_max <dbl>,
## #   danceability_missing <int>, valence_min <dbl>, valence_max <dbl>,
## #   valence_missing <int>

# Missing values by year for the period relevant to our analysis
songs_clean |>
  filter(year >= 2016, year <= 2025) |>
  group_by(year) |>
  summarise(
    n = n(),
    tempo_missing = sum(is.na(tempo)),
    energy_missing = sum(is.na(energy)),
    danceability_missing = sum(is.na(danceability)),
    valence_missing = sum(is.na(valence)),
    .groups = "drop"
  )

## # A tibble: 9 × 6
##    year     n tempo_missing energy_missing danceability_missing valence_missing
##   <dbl> <int>         <int>          <int>                <int>           <int>
## 1  2016    42             0              0                    0               0
## 2  2017    42             0              0                    0               0
## 3  2018    43             0              0                    0               0
## 4  2019    41             0              1                    0               1
## 5  2021    39             0              0                    0               0
## 6  2022    40             0              0                    0               0
## 7  2023    37             0              0                    0               0
## 8  2024    37             0              0                    0               0
## 9  2025    37             0             37                    0              37

4.4 Performance Order Data

performance_order_clean <- performance_order_raw |>
  clean_names() |>
  select(
    year,
    to_country,
    performer,
    song,
    performance_order_final
  ) |>
  rename(
    performance_order = performance_order_final
  )

# Inspect the cleaned structure
glimpse(performance_order_clean)

## Rows: 1,756
## Columns: 5
## $ year              <dbl> 1956, 1956, 1956, 1956, 1956, 1956, 1956, 1956, 1956…
## $ to_country        <chr> "Switzerland", "Netherlands", "Belgium", "Germany", …
## $ performer         <chr> "Lys Assia", "Jetty Paerl", "Fud Leclerc", "Walter A…
## $ song              <chr> "Refrain", "De Vogels Van Holland", "Messieurs Les N…
## $ performance_order <dbl> 2, 1, 3, 4, 5, 13, 7, 8, 10, 11, 12, 14, 6, 8, 9, 2,…

performance_order_check <- effects_clean |>
  select(
    year,
    to_country,
    song,
    performance_order_effects = performance_order
  ) |>
  left_join(
    performance_order_clean |>
      select(
        year,
        to_country,
        performance_order_department = performance_order
      ),
    by = c("year", "to_country")
  )

# Check whether any Grand Final entries are missing department performance-order data
performance_order_check |>
  filter(is.na(performance_order_department))

## # A tibble: 31 × 5
##     year to_country     song       performance_order_ef…¹ performance_order_de…²
##    <dbl> <chr>          <chr>                       <dbl>                  <dbl>
##  1  2025 Norway         Lighter                         1                     NA
##  2  2025 Luxembourg     La poupée…                      2                     NA
##  3  2025 Estonia        Espresso …                      3                     NA
##  4  2025 Israel         New Day W…                      4                     NA
##  5  2025 Lithuania      Tavo akys                       5                     NA
##  6  2025 Spain          Esa diva                        6                     NA
##  7  2025 Ukraine        Bird of P…                      7                     NA
##  8  2025 United Kingdom What the …                      8                     NA
##  9  2025 Austria        Wasted Lo…                      9                     NA
## 10  2025 Iceland        Róa                            10                     NA
## # ℹ 21 more rows
## # ℹ abbreviated names: ¹​performance_order_effects,
## #   ²​performance_order_department

# Check whether the two performance-order sources disagree
performance_order_check |>
  filter(
    !is.na(performance_order_department),
    performance_order_effects != performance_order_department
  )

## # A tibble: 0 × 5
## # ℹ 5 variables: year <dbl>, to_country <chr>, song <chr>,
## #   performance_order_effects <dbl>, performance_order_department <dbl>

performance_order_check |>
  filter(is.na(performance_order_department)) |>
  count(year, name = "missing_department_orders") |>
  arrange(desc(year))

## # A tibble: 6 × 2
##    year missing_department_orders
##   <dbl>                     <int>
## 1  2025                        26
## 2  2024                         1
## 3  2022                         1
## 4  2019                         1
## 5  2018                         1
## 6  2016                         1

performance_order_check |>
  filter(
    is.na(performance_order_department),
    year < 2025
  )

## # A tibble: 5 × 5
##    year to_country     song        performance_order_ef…¹ performance_order_de…²
##   <dbl> <chr>          <chr>                        <dbl>                  <dbl>
## 1  2024 Netherlands    Europapa                         5                     NA
## 2  2022 Czech Republic Lights Off                       1                     NA
## 3  2019 Czech Republic Friend of …                      3                     NA
## 4  2018 Czech Republic Lie to Me                       14                     NA
## 5  2016 Czech Republic I Stand                          2                     NA
## # ℹ abbreviated names: ¹​performance_order_effects,
## #   ²​performance_order_department

5. Variable Construction

5.1 Voting Outcomes

votes_songlevel <- votes_clean |>
  filter(
    round == "final",
    year >= 2016,
    year <= 2025
  ) |>
  group_by(year, to_country) |>
  summarise(
    tele_points = sum(tele_points, na.rm = TRUE),
    jury_points = sum(jury_points, na.rm = TRUE),
    .groups = "drop"
  )

# Inspect the resulting dataset
glimpse(votes_songlevel)

## Rows: 232
## Columns: 4
## $ year        <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016…
## $ to_country  <chr> "Armenia", "Australia", "Austria", "Azerbaijan", "Belgium"…
## $ tele_points <dbl> 134, 191, 120, 73, 51, 180, 33, 53, 0, 109, 24, 10, 56, 11…
## $ jury_points <dbl> 115, 320, 31, 44, 130, 127, 40, 43, 41, 148, 80, 1, 52, 12…

# Count the number of finalists represented in each year
votes_songlevel |>
  count(year, name = "number_of_finalists") |>
  arrange(year)

## # A tibble: 9 × 2
##    year number_of_finalists
##   <dbl>               <int>
## 1  2016                  26
## 2  2017                  26
## 3  2018                  26
## 4  2019                  26
## 5  2021                  26
## 6  2022                  25
## 7  2023                  26
## 8  2024                  25
## 9  2025                  26

voting_entities <- votes_clean |>
  filter(
    round == "final",
    year >= 2016,
    year <= 2025
  ) |>
  group_by(year) |>
  summarise(
    jury_voting_entities = n_distinct(
      from_country[
        !is.na(jury_points) &
        from_country != "Rest of the World"
      ]
    ),
    tele_voting_entities = n_distinct(
      from_country[!is.na(tele_points)]
    ),
    .groups = "drop"
  )

voting_entities

## # A tibble: 9 × 3
##    year jury_voting_entities tele_voting_entities
##   <dbl>                <int>                <int>
## 1  2016                   42                   42
## 2  2017                   42                   42
## 3  2018                   43                   43
## 4  2019                   41                   41
## 5  2021                   39                   39
## 6  2022                   40                   40
## 7  2023                   37                   38
## 8  2024                   37                   38
## 9  2025                   37                   38

5.2 Normalized Scores

# Calculate the number of eligible jury and televote voting entities
# separately for every finalist
voting_opportunities <- votes_clean |>
  filter(
    round == "final",
    year >= 2016,
    year <= 2025
  ) |>
  group_by(year, to_country) |>
  summarise(
    jury_opportunities = n_distinct(
      from_country[
        !is.na(jury_points) &
        from_country != "Rest of the World"
      ]
    ),
    tele_opportunities = n_distinct(
      from_country[!is.na(tele_points)]
    ),
    .groups = "drop"
  )

# Add the voting opportunities to the song-level voting outcomes
# and construct normalized scores
votes_normalized <- votes_songlevel |>
  left_join(
    voting_opportunities,
    by = c("year", "to_country")
  ) |>
  mutate(
    max_jury_points = jury_opportunities * 12,
    max_tele_points = tele_opportunities * 12,
    
    jury_score_norm = jury_points / max_jury_points,
    tele_score_norm = tele_points / max_tele_points
  )

# Inspect the resulting structure
glimpse(votes_normalized)

## Rows: 232
## Columns: 10
## $ year               <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 201…
## $ to_country         <chr> "Armenia", "Australia", "Austria", "Azerbaijan", "B…
## $ tele_points        <dbl> 134, 191, 120, 73, 51, 180, 33, 53, 0, 109, 24, 10,…
## $ jury_points        <dbl> 115, 320, 31, 44, 130, 127, 40, 43, 41, 148, 80, 1,…
## $ jury_opportunities <int> 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41,…
## $ tele_opportunities <int> 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41,…
## $ max_jury_points    <dbl> 492, 492, 492, 492, 492, 492, 492, 492, 492, 492, 4…
## $ max_tele_points    <dbl> 492, 492, 492, 492, 492, 492, 492, 492, 492, 492, 4…
## $ jury_score_norm    <dbl> 0.23373984, 0.65040650, 0.06300813, 0.08943089, 0.2…
## $ tele_score_norm    <dbl> 0.27235772, 0.38821138, 0.24390244, 0.14837398, 0.1…

# Check for missing or impossible denominators
votes_normalized |>
  filter(
    is.na(max_jury_points) |
    is.na(max_tele_points) |
    max_jury_points <= 0 |
    max_tele_points <= 0
  )

## # A tibble: 0 × 10
## # ℹ 10 variables: year <dbl>, to_country <chr>, tele_points <dbl>,
## #   jury_points <dbl>, jury_opportunities <int>, tele_opportunities <int>,
## #   max_jury_points <dbl>, max_tele_points <dbl>, jury_score_norm <dbl>,
## #   tele_score_norm <dbl>

# Check whether normalized scores fall outside the valid 0-1 range
votes_normalized |>
  filter(
    jury_score_norm < 0 |
    jury_score_norm > 1 |
    tele_score_norm < 0 |
    tele_score_norm > 1
  )

## # A tibble: 0 × 10
## # ℹ 10 variables: year <dbl>, to_country <chr>, tele_points <dbl>,
## #   jury_points <dbl>, jury_opportunities <int>, tele_opportunities <int>,
## #   max_jury_points <dbl>, max_tele_points <dbl>, jury_score_norm <dbl>,
## #   tele_score_norm <dbl>

# Inspect voting opportunities by year
votes_normalized |>
  group_by(year) |>
  summarise(
    jury_opportunities_min = min(jury_opportunities),
    jury_opportunities_max = max(jury_opportunities),
    tele_opportunities_min = min(tele_opportunities),
    tele_opportunities_max = max(tele_opportunities),
    .groups = "drop"
  )

## # A tibble: 9 × 5
##    year jury_opportunities_min jury_opportunities_max tele_opportunities_min
##   <dbl>                  <int>                  <int>                  <int>
## 1  2016                     41                     41                     41
## 2  2017                     41                     41                     41
## 3  2018                     42                     42                     42
## 4  2019                     40                     40                     40
## 5  2021                     38                     38                     38
## 6  2022                     39                     39                     39
## 7  2023                     36                     36                     37
## 8  2024                     36                     36                     37
## 9  2025                     36                     36                     37
## # ℹ 1 more variable: tele_opportunities_max <int>

# Inspect selected 2025 examples
votes_normalized |>
  filter(
    year == 2025,
    to_country %in% c(
      "Austria",
      "Israel",
      "United Kingdom",
      "Switzerland"
    )
  ) |>
  select(
    year,
    to_country,
    tele_points,
    jury_points,
    tele_opportunities,
    jury_opportunities,
    tele_score_norm,
    jury_score_norm
  )

## # A tibble: 4 × 8
##    year to_country tele_points jury_points tele_opportunities jury_opportunities
##   <dbl> <chr>            <dbl>       <dbl>              <int>              <int>
## 1  2025 Austria            178         258                 37                 36
## 2  2025 Israel             297          60                 37                 36
## 3  2025 Switzerla…           0         214                 37                 36
## 4  2025 United Ki…           0          88                 37                 36
## # ℹ 2 more variables: tele_score_norm <dbl>, jury_score_norm <dbl>

5.3 Jury–Televote Gap

votes_normalized <- votes_normalized |>
  mutate(
    jury_televote_gap = tele_score_norm - jury_score_norm
  )

# Inspect the overall range of the gap
votes_normalized |>
  summarise(
    gap_min = min(jury_televote_gap, na.rm = TRUE),
    gap_mean = mean(jury_televote_gap, na.rm = TRUE),
    gap_max = max(jury_televote_gap, na.rm = TRUE),
    gap_missing = sum(is.na(jury_televote_gap))
  )

## # A tibble: 1 × 4
##   gap_min   gap_mean gap_max gap_missing
##     <dbl>      <dbl>   <dbl>       <int>
## 1  -0.495 -0.0000469   0.607           0

# Inspect selected 2025 examples
votes_normalized |>
  filter(
    year == 2025,
    to_country %in% c(
      "Austria",
      "Israel",
      "United Kingdom",
      "Switzerland"
    )
  ) |>
  select(
    year,
    to_country,
    tele_score_norm,
    jury_score_norm,
    jury_televote_gap
  )

## # A tibble: 4 × 5
##    year to_country     tele_score_norm jury_score_norm jury_televote_gap
##   <dbl> <chr>                    <dbl>           <dbl>             <dbl>
## 1  2025 Austria                  0.401           0.597            -0.196
## 2  2025 Israel                   0.669           0.139             0.530
## 3  2025 Switzerland              0               0.495            -0.495
## 4  2025 United Kingdom           0               0.204            -0.204

6. Dataset Merging

6.1 Identifier Harmonization

# Effects observations from 2016 onward that do not match voting data
effects_clean |>
  filter(year >= 2016) |>
  anti_join(
    votes_normalized,
    by = c("year", "to_country")
  ) |>
  select(year, to_country, song) |>
  arrange(desc(year), to_country)

## # A tibble: 5 × 3
##    year to_country     song              
##   <dbl> <chr>          <chr>             
## 1  2024 Netherlands    Europapa          
## 2  2022 Czech Republic Lights Off        
## 3  2019 Czech Republic Friend of a Friend
## 4  2018 Czech Republic Lie to Me         
## 5  2016 Czech Republic I Stand

# Effects observations that do not match musical-features data
effects_clean |>
  anti_join(
    songs_clean,
    by = c("year", "to_country")
  ) |>
  select(year, to_country, song) |>
  arrange(desc(year), to_country)

## # A tibble: 4 × 3
##    year to_country     song              
##   <dbl> <chr>          <chr>             
## 1  2022 Czech Republic Lights Off        
## 2  2019 Czech Republic Friend of a Friend
## 3  2018 Czech Republic Lie to Me         
## 4  2016 Czech Republic I Stand

effects_clean <- effects_clean |>
  mutate(
    to_country = if_else(
      to_country == "Czech Republic",
      "Czechia",
      to_country
    )
  )

# Recheck unmatched voting identifiers
effects_clean |>
  filter(year >= 2016) |>
  anti_join(
    votes_normalized,
    by = c("year", "to_country")
  ) |>
  select(year, to_country, song) |>
  arrange(desc(year), to_country)

## # A tibble: 1 × 3
##    year to_country  song    
##   <dbl> <chr>       <chr>   
## 1  2024 Netherlands Europapa

# Recheck unmatched musical-feature identifiers
effects_clean |>
  anti_join(
    songs_clean,
    by = c("year", "to_country")
  ) |>
  select(year, to_country, song) |>
  arrange(desc(year), to_country)

## # A tibble: 0 × 3
## # ℹ 3 variables: year <dbl>, to_country <chr>, song <chr>

6.2 Merge Effects and Voting Data

analysis_data <- effects_clean |>
  left_join(
    votes_normalized |>
      select(
        year,
        to_country,
        tele_points,
        jury_points,
        jury_opportunities,
        tele_opportunities,
        jury_score_norm,
        tele_score_norm,
        jury_televote_gap
      ),
    by = c("year", "to_country")
  )

# Check dimensions after the merge
dim(analysis_data)

## [1] 413  20

# Inspect the resulting structure
glimpse(analysis_data)

## Rows: 413
## Columns: 20
## $ year                  <dbl> 2025, 2025, 2025, 2025, 2025, 2025, 2025, 2025, …
## $ to_country            <chr> "Norway", "Luxembourg", "Estonia", "Israel", "Li…
## $ song                  <chr> "Lighter", "La poupée monte le son", "Espresso M…
## $ artist_name           <chr> "Kyle Alessandro", "Laura Thorn", "Tommy Cash", …
## $ round                 <chr> "final", "final", "final", "final", "final", "fi…
## $ performance_order     <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1…
## $ fire                  <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, …
## $ smoke                 <dbl> 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, …
## $ wind_machine          <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, …
## $ confetti              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ special_effects_count <dbl> 2, 0, 0, 0, 2, 0, 0, 0, 1, 1, 1, 0, 1, 0, 2, 0, …
## $ coding_complete       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ notes_data_collection <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ tele_points           <dbl> 67, 24, 258, 297, 62, 10, 158, 0, 178, 33, 42, 4…
## $ jury_points           <dbl> 22, 23, 98, 60, 34, 27, 60, 88, 258, 0, 116, 133…
## $ jury_opportunities    <int> 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, …
## $ tele_opportunities    <int> 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, …
## $ jury_score_norm       <dbl> 0.05092593, 0.05324074, 0.22685185, 0.13888889, …
## $ tele_score_norm       <dbl> 0.150900901, 0.054054054, 0.581081081, 0.6689189…
## $ jury_televote_gap     <dbl> 0.0999749750, 0.0008133133, 0.3542292292, 0.5300…

# Count missing normalized voting outcomes by year
analysis_data |>
  group_by(year) |>
  summarise(
    n = n(),
    tele_score_missing = sum(is.na(tele_score_norm)),
    jury_score_missing = sum(is.na(jury_score_norm)),
    gap_missing = sum(is.na(jury_televote_gap)),
    .groups = "drop"
  ) |>
  arrange(year)

## # A tibble: 16 × 5
##     year     n tele_score_missing jury_score_missing gap_missing
##    <dbl> <int>              <int>              <int>       <int>
##  1  2009    25                 25                 25          25
##  2  2010    25                 25                 25          25
##  3  2011    25                 25                 25          25
##  4  2012    26                 26                 26          26
##  5  2013    26                 26                 26          26
##  6  2014    26                 26                 26          26
##  7  2015    27                 27                 27          27
##  8  2016    26                  0                  0           0
##  9  2017    26                  0                  0           0
## 10  2018    26                  0                  0           0
## 11  2019    26                  0                  0           0
## 12  2021    26                  0                  0           0
## 13  2022    25                  0                  0           0
## 14  2023    26                  0                  0           0
## 15  2024    26                  1                  1           1
## 16  2025    26                  0                  0           0

# Inspect post-2015 observations with missing voting outcomes
analysis_data |>
  filter(
    year >= 2016,
    is.na(tele_score_norm) | is.na(jury_score_norm)
  ) |>
  select(
    year,
    to_country,
    song,
    tele_score_norm,
    jury_score_norm,
    jury_televote_gap
  )

## # A tibble: 1 × 6
##    year to_country  song     tele_score_norm jury_score_norm jury_televote_gap
##   <dbl> <chr>       <chr>              <dbl>           <dbl>             <dbl>
## 1  2024 Netherlands Europapa              NA              NA                NA

6.3 Merge Musical Features

analysis_data <- analysis_data |>
  left_join(
    songs_clean |>
      select(
        year,
        to_country,
        tempo,
        energy,
        danceability,
        valence
      ),
    by = c("year", "to_country")
  )

# Check dimensions after the merge
dim(analysis_data)

## [1] 413  24

# Inspect the newly added musical-feature variables
analysis_data |>
  select(
    year,
    to_country,
    tempo,
    energy,
    danceability,
    valence
  ) |>
  glimpse()

## Rows: 413
## Columns: 6
## $ year         <dbl> 2025, 2025, 2025, 2025, 2025, 2025, 2025, 2025, 2025, 202…
## $ to_country   <chr> "Norway", "Luxembourg", "Estonia", "Israel", "Lithuania",…
## $ tempo        <dbl> 120.96774, 116.27907, 120.00000, 88.23529, 75.75758, 75.7…
## $ energy       <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ danceability <dbl> 0.8429474, 0.8348589, 0.8867011, 0.5887105, 0.1796271, 0.…
## $ valence      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…

analysis_data |>
  group_by(year) |>
  summarise(
    n = n(),
    tempo_missing = sum(is.na(tempo)),
    energy_missing = sum(is.na(energy)),
    danceability_missing = sum(is.na(danceability)),
    valence_missing = sum(is.na(valence)),
    .groups = "drop"
  ) |>
  arrange(year)

## # A tibble: 16 × 6
##     year     n tempo_missing energy_missing danceability_missing valence_missing
##    <dbl> <int>         <int>          <int>                <int>           <int>
##  1  2009    25             0              0                    0               0
##  2  2010    25             0              0                    0               0
##  3  2011    25             0              0                    0               0
##  4  2012    26             0              0                    0               0
##  5  2013    26             0              0                    0               0
##  6  2014    26             0              0                    0               0
##  7  2015    27             0              0                    0               0
##  8  2016    26             0              0                    0               0
##  9  2017    26             0              0                    0               0
## 10  2018    26             0              0                    0               0
## 11  2019    26             0              0                    0               0
## 12  2021    26             0              0                    0               0
## 13  2022    25             0              0                    0               0
## 14  2023    26             0              0                    0               0
## 15  2024    26             0              0                    0               0
## 16  2025    26             0             26                    0              26

6.4 Merge Validation

# Check final dimensions
dim(analysis_data)

## [1] 413  24

# Check for duplicate year-country identifiers
analysis_data |>
  count(year, to_country) |>
  filter(n > 1)

## # A tibble: 0 × 3
## # ℹ 3 variables: year <dbl>, to_country <chr>, n <int>

# Check for missing essential identifiers
analysis_data |>
  summarise(
    year_missing = sum(is.na(year)),
    country_missing = sum(is.na(to_country)),
    song_missing = sum(is.na(song)),
    performance_order_missing = sum(is.na(performance_order))
  )

## # A tibble: 1 × 4
##   year_missing country_missing song_missing performance_order_missing
##          <int>           <int>        <int>                     <int>
## 1            0               0            0                         0

# Check final year coverage and observations per year
analysis_data |>
  count(year, name = "number_of_entries") |>
  arrange(year)

## # A tibble: 16 × 2
##     year number_of_entries
##    <dbl>             <int>
##  1  2009                25
##  2  2010                25
##  3  2011                25
##  4  2012                26
##  5  2013                26
##  6  2014                26
##  7  2015                27
##  8  2016                26
##  9  2017                26
## 10  2018                26
## 11  2019                26
## 12  2021                26
## 13  2022                25
## 14  2023                26
## 15  2024                26
## 16  2025                26

7. Final Sample Construction

7.1 H1 Sample: Special Effects and Televote

sample_h1 <- analysis_data |>
  filter(year >= 2016) |>
  drop_na(
    tele_score_norm,
    special_effects_count
  )

# Inspect the H1 sample size and year coverage
sample_h1 |>
  summarise(
    observations = n(),
    first_year = min(year),
    last_year = max(year),
    number_of_contests = n_distinct(year)
  )

## # A tibble: 1 × 4
##   observations first_year last_year number_of_contests
##          <int>      <dbl>     <dbl>              <int>
## 1          232       2016      2025                  9

# Inspect the number of H1 observations per contest
sample_h1 |>
  count(year, name = "observations") |>
  arrange(year)

## # A tibble: 9 × 2
##    year observations
##   <dbl>        <int>
## 1  2016           26
## 2  2017           26
## 3  2018           26
## 4  2019           26
## 5  2021           26
## 6  2022           25
## 7  2023           26
## 8  2024           25
## 9  2025           26

# Check missingness in variables relevant to H1
sample_h1 |>
  summarise(
    tele_score_missing = sum(is.na(tele_score_norm)),
    effects_count_missing = sum(is.na(special_effects_count)),
    fire_missing = sum(is.na(fire)),
    smoke_missing = sum(is.na(smoke)),
    wind_machine_missing = sum(is.na(wind_machine)),
    confetti_missing = sum(is.na(confetti))
  )

## # A tibble: 1 × 6
##   tele_score_missing effects_count_missing fire_missing smoke_missing
##                <int>                 <int>        <int>         <int>
## 1                  0                     0            0             0
## # ℹ 2 more variables: wind_machine_missing <int>, confetti_missing <int>

# Check the range of the main H1 variables
sample_h1 |>
  summarise(
    tele_score_min = min(tele_score_norm),
    tele_score_max = max(tele_score_norm),
    effects_count_min = min(special_effects_count),
    effects_count_max = max(special_effects_count)
  )

## # A tibble: 1 × 4
##   tele_score_min tele_score_max effects_count_min effects_count_max
##            <dbl>          <dbl>             <dbl>             <dbl>
## 1              0          0.938                 0                 2

# Check distribution of the special-effects count
sample_h1 |>
  count(special_effects_count, name = "observations") |>
  arrange(special_effects_count)

## # A tibble: 3 × 2
##   special_effects_count observations
##                   <dbl>        <int>
## 1                     0          157
## 2                     1           58
## 3                     2           17

7.2 H2 Sample: Musical Characteristics and Jury

sample_h2 <- analysis_data |>
  filter(year >= 2016) |>
  drop_na(
    jury_score_norm,
    tempo,
    energy,
    danceability,
    valence
  )

# Inspect the H2 sample size and year coverage
sample_h2 |>
  summarise(
    observations = n(),
    first_year = min(year),
    last_year = max(year),
    number_of_contests = n_distinct(year)
  )

## # A tibble: 1 × 4
##   observations first_year last_year number_of_contests
##          <int>      <dbl>     <dbl>              <int>
## 1          206       2016      2024                  8

# Inspect the number of H2 observations per contest
sample_h2 |>
  count(year, name = "observations") |>
  arrange(year)

## # A tibble: 8 × 2
##    year observations
##   <dbl>        <int>
## 1  2016           26
## 2  2017           26
## 3  2018           26
## 4  2019           26
## 5  2021           26
## 6  2022           25
## 7  2023           26
## 8  2024           25

# Show observations from 2016 onward that are not included in H2
analysis_data |>
  filter(year >= 2016) |>
  mutate(
    missing_jury = is.na(jury_score_norm),
    missing_tempo = is.na(tempo),
    missing_energy = is.na(energy),
    missing_danceability = is.na(danceability),
    missing_valence = is.na(valence)
  ) |>
  filter(
    missing_jury |
    missing_tempo |
    missing_energy |
    missing_danceability |
    missing_valence
  ) |>
  select(
    year,
    to_country,
    song,
    missing_jury,
    missing_tempo,
    missing_energy,
    missing_danceability,
    missing_valence
  ) |>
  arrange(year, to_country)

## # A tibble: 27 × 8
##     year to_country  song              missing_jury missing_tempo missing_energy
##    <dbl> <chr>       <chr>             <lgl>        <lgl>         <lgl>         
##  1  2024 Netherlands Europapa          TRUE         FALSE         FALSE         
##  2  2025 Albania     Zjerm             FALSE        FALSE         TRUE          
##  3  2025 Armenia     Survivor          FALSE        FALSE         TRUE          
##  4  2025 Austria     Wasted Love       FALSE        FALSE         TRUE          
##  5  2025 Denmark     Hallucination     FALSE        FALSE         TRUE          
##  6  2025 Estonia     Espresso Macchia… FALSE        FALSE         TRUE          
##  7  2025 Finland     Ich komme         FALSE        FALSE         TRUE          
##  8  2025 France      Maman             FALSE        FALSE         TRUE          
##  9  2025 Germany     Baller            FALSE        FALSE         TRUE          
## 10  2025 Greece      Asteromáta        FALSE        FALSE         TRUE          
## # ℹ 17 more rows
## # ℹ 2 more variables: missing_danceability <lgl>, missing_valence <lgl>

7.3 H3 Sample: Performance Order and Televote

sample_h3 <- analysis_data |>
  filter(year >= 2016) |>
  drop_na(
    tele_score_norm,
    performance_order
  )

# Inspect the H3 sample size and year coverage
sample_h3 |>
  summarise(
    observations = n(),
    first_year = min(year),
    last_year = max(year),
    number_of_contests = n_distinct(year)
  )

## # A tibble: 1 × 4
##   observations first_year last_year number_of_contests
##          <int>      <dbl>     <dbl>              <int>
## 1          232       2016      2025                  9

# Inspect the number of H3 observations per contest
sample_h3 |>
  count(year, name = "observations") |>
  arrange(year)

## # A tibble: 9 × 2
##    year observations
##   <dbl>        <int>
## 1  2016           26
## 2  2017           26
## 3  2018           26
## 4  2019           26
## 5  2021           26
## 6  2022           25
## 7  2023           26
## 8  2024           25
## 9  2025           26

# Check missingness and overall range of H3 variables
sample_h3 |>
  summarise(
    tele_score_missing = sum(is.na(tele_score_norm)),
    performance_order_missing = sum(is.na(performance_order)),
    performance_order_min = min(performance_order),
    performance_order_max = max(performance_order)
  )

## # A tibble: 1 × 4
##   tele_score_missing performance_order_missing performance_order_min
##                <int>                     <int>                 <dbl>
## 1                  0                         0                     1
## # ℹ 1 more variable: performance_order_max <dbl>

# Check performance order coverage by contest
sample_h3 |>
  group_by(year) |>
  summarise(
    observations = n(),
    first_position = min(performance_order),
    last_position = max(performance_order),
    .groups = "drop"
  ) |>
  arrange(year)

## # A tibble: 9 × 4
##    year observations first_position last_position
##   <dbl>        <int>          <dbl>         <dbl>
## 1  2016           26              1            26
## 2  2017           26              1            26
## 3  2018           26              1            26
## 4  2019           26              1            26
## 5  2021           26              1            26
## 6  2022           25              1            25
## 7  2023           26              1            26
## 8  2024           25              1            26
## 9  2025           26              1            26

# Confirm the 2024 performance order positions retained in H3
sample_h3 |>
  filter(year == 2024) |>
  select(
    year,
    to_country,
    performance_order
  ) |>
  arrange(performance_order)

## # A tibble: 25 × 3
##     year to_country performance_order
##    <dbl> <chr>                  <dbl>
##  1  2024 Sweden                     1
##  2  2024 Ukraine                    2
##  3  2024 Germany                    3
##  4  2024 Luxembourg                 4
##  5  2024 Israel                     6
##  6  2024 Lithuania                  7
##  7  2024 Spain                      8
##  8  2024 Estonia                    9
##  9  2024 Ireland                   10
## 10  2024 Latvia                    11
## # ℹ 15 more rows

7.4 H4 Sample: Jury versus Televote Gap

sample_h4 <- analysis_data |>
  filter(year >= 2016) |>
  drop_na(
    jury_televote_gap,
    special_effects_count
  )

# Inspect the H4 sample size and year coverage
sample_h4 |>
  summarise(
    observations = n(),
    first_year = min(year),
    last_year = max(year),
    number_of_contests = n_distinct(year)
  )

## # A tibble: 1 × 4
##   observations first_year last_year number_of_contests
##          <int>      <dbl>     <dbl>              <int>
## 1          232       2016      2025                  9

# Inspect the number of H4 observations per contest
sample_h4 |>
  count(year, name = "observations") |>
  arrange(year)

## # A tibble: 9 × 2
##    year observations
##   <dbl>        <int>
## 1  2016           26
## 2  2017           26
## 3  2018           26
## 4  2019           26
## 5  2021           26
## 6  2022           25
## 7  2023           26
## 8  2024           25
## 9  2025           26

# Check missingness in variables required for H4
sample_h4 |>
  summarise(
    gap_missing = sum(is.na(jury_televote_gap)),
    effects_count_missing = sum(is.na(special_effects_count))
  )

## # A tibble: 1 × 2
##   gap_missing effects_count_missing
##         <int>                 <int>
## 1           0                     0

# Check the observed ranges of the H4 variables
sample_h4 |>
  summarise(
    gap_min = min(jury_televote_gap),
    gap_max = max(jury_televote_gap),
    effects_count_min = min(special_effects_count),
    effects_count_max = max(special_effects_count)
  )

## # A tibble: 1 × 4
##   gap_min gap_max effects_count_min effects_count_max
##     <dbl>   <dbl>             <dbl>             <dbl>
## 1  -0.495   0.607                 0                 2

7.5 Sample Comparison

sample_comparison <- tibble(
  hypothesis = c("H1", "H2", "H3", "H4"),
  observations = c(
    nrow(sample_h1),
    nrow(sample_h2),
    nrow(sample_h3),
    nrow(sample_h4)
  ),
  first_year = c(
    min(sample_h1$year),
    min(sample_h2$year),
    min(sample_h3$year),
    min(sample_h4$year)
  ),
  last_year = c(
    max(sample_h1$year),
    max(sample_h2$year),
    max(sample_h3$year),
    max(sample_h4$year)
  ),
  number_of_contests = c(
    n_distinct(sample_h1$year),
    n_distinct(sample_h2$year),
    n_distinct(sample_h3$year),
    n_distinct(sample_h4$year)
  )
)

sample_comparison

## # A tibble: 4 × 5
##   hypothesis observations first_year last_year number_of_contests
##   <chr>             <int>      <dbl>     <dbl>              <int>
## 1 H1                  232       2016      2025                  9
## 2 H2                  206       2016      2024                  8
## 3 H3                  232       2016      2025                  9
## 4 H4                  232       2016      2025                  9

8. Data Validation

8.1 Final Dataset Integrity

# Final dimensions of the master analytical dataset
dim(analysis_data)

## [1] 413  24

# Check year coverage
sort(unique(analysis_data$year))

##  [1] 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2021 2022 2023 2024
## [16] 2025

# Check round labels
unique(analysis_data$round)

## [1] "final"

# Check for duplicate performance identifiers
analysis_data |>
  count(year, to_country) |>
  filter(n > 1)

## # A tibble: 0 × 3
## # ℹ 3 variables: year <dbl>, to_country <chr>, n <int>

# Check essential identifier completeness
analysis_data |>
  summarise(
    year_missing = sum(is.na(year)),
    country_missing = sum(is.na(to_country)),
    song_missing = sum(is.na(song)),
    artist_missing = sum(is.na(artist_name)),
    performance_order_missing = sum(is.na(performance_order))
  )

## # A tibble: 1 × 5
##   year_missing country_missing song_missing artist_missing
##          <int>           <int>        <int>          <int>
## 1            0               0            0              0
## # ℹ 1 more variable: performance_order_missing <int>

8.2 Hypothesis Sample Validation

hypothesis_validation <- tibble(
  hypothesis = c("H1", "H2", "H3", "H4"),
  
  observations = c(
    nrow(sample_h1),
    nrow(sample_h2),
    nrow(sample_h3),
    nrow(sample_h4)
  ),
  
  outcome_missing = c(
    sum(is.na(sample_h1$tele_score_norm)),
    sum(is.na(sample_h2$jury_score_norm)),
    sum(is.na(sample_h3$tele_score_norm)),
    sum(is.na(sample_h4$jury_televote_gap))
  ),
  
  main_predictors_missing = c(
    sum(is.na(sample_h1$special_effects_count)),
    
    sum(
      is.na(sample_h1$tempo) |
      is.na(sample_h2$energy) |
      is.na(sample_h3$danceability) |
      is.na(sample_h4$valence)
    ),
    
    sum(is.na(sample_h3$performance_order)),
    
    sum(is.na(sample_h4$special_effects_count))
  )
)

hypothesis_validation

## # A tibble: 4 × 4
##   hypothesis observations outcome_missing main_predictors_missing
##   <chr>             <int>           <int>                   <int>
## 1 H1                  232               0                       0
## 2 H2                  206               0                      26
## 3 H3                  232               0                       0
## 4 H4                  232               0                       0

8.3 Variable Range and Missingness Checks

range_validation <- tibble(
  check = c(
    "Normalized televote outside 0-1",
    "Normalized jury outside 0-1",
    "Jury-televote gap outside -1 to 1",
    "Special-effects count outside 0-4",
    "Performance order below 1",
    "Tempo at or below 0",
    "Energy outside 0-1",
    "Danceability outside 0-1",
    "Valence outside 0-1"
  ),
  
  violations = c(
    sum(
      analysis_data$tele_score_norm < 0 |
      analysis_data$tele_score_norm > 1,
      na.rm = TRUE
    ),
    
    sum(
      analysis_data$jury_score_norm < 0 |
      analysis_data$jury_score_norm > 1,
      na.rm = TRUE
    ),
    
    sum(
      analysis_data$jury_televote_gap < -1 |
      analysis_data$jury_televote_gap > 1,
      na.rm = TRUE
    ),
    
    sum(
      analysis_data$special_effects_count < 0 |
      analysis_data$special_effects_count > 4,
      na.rm = TRUE
    ),
    
    sum(
      analysis_data$performance_order < 1,
      na.rm = TRUE
    ),
    
    sum(
      analysis_data$tempo <= 0,
      na.rm = TRUE
    ),
    
    sum(
      analysis_data$energy < 0 |
      analysis_data$energy > 1,
      na.rm = TRUE
    ),
    
    sum(
      analysis_data$danceability < 0 |
      analysis_data$danceability > 1,
      na.rm = TRUE
    ),
    
    sum(
      analysis_data$valence < 0 |
      analysis_data$valence > 1,
      na.rm = TRUE
    )
  )
)

range_validation

## # A tibble: 9 × 2
##   check                             violations
##   <chr>                                  <int>
## 1 Normalized televote outside 0-1            0
## 2 Normalized jury outside 0-1                0
## 3 Jury-televote gap outside -1 to 1          0
## 4 Special-effects count outside 0-4          0
## 5 Performance order below 1                  0
## 6 Tempo at or below 0                        0
## 7 Energy outside 0-1                         0
## 8 Danceability outside 0-1                   0
## 9 Valence outside 0-1                        0

9. Export Analysis Dataset

9.1 Export Master Dataset

write_csv(
  analysis_data,
  here("data output", "analysis_data_master.csv"),
  na = "NA"
)

# Confirm that the file was successfully created
file.exists(
  here("data output", "analysis_data_master.csv")
)

## [1] TRUE

9.2 Export Hypothesis Samples

# Export the four hypothesis-specific analytical samples
write_csv(
  sample_h1,
  here("data output", "sample_h1.csv"),
  na = "NA"
)

write_csv(
  sample_h2,
  here("data output", "sample_h2.csv"),
  na = "NA"
)

write_csv(
  sample_h3,
  here("data output", "sample_h3.csv"),
  na = "NA"
)

write_csv(
  sample_h4,
  here("data output", "sample_h4.csv"),
  na = "NA"
)

# Confirm that all four files were successfully created
file.exists(
  c(
    here("data output", "sample_h1.csv"),
    here("data output", "sample_h2.csv"),
    here("data output", "sample_h3.csv"),
    here("data output", "sample_h4.csv")
  )
)

## [1] TRUE TRUE TRUE TRUE

9.3 Analysis Dataset Overview

dataset_overview <- tibble(
  dataset = c(
    "analysis_data_master.csv",
    "sample_h1.csv",
    "sample_h2.csv",
    "sample_h3.csv",
    "sample_h4.csv"
  ),
  purpose = c(
    "Master merged dataset",
    "H1: Special effects and televote",
    "H2: Musical features and jury",
    "H3: Running order and televote",
    "H4: Special effects and jury-televote gap"
  ),
  observations = c(
    nrow(analysis_data),
    nrow(sample_h1),
    nrow(sample_h2),
    nrow(sample_h3),
    nrow(sample_h4)
  )
)

dataset_overview

## # A tibble: 5 × 3
##   dataset                  purpose                                  observations
##   <chr>                    <chr>                                           <int>
## 1 analysis_data_master.csv Master merged dataset                             413
## 2 sample_h1.csv            H1: Special effects and televote                  232
## 3 sample_h2.csv            H2: Musical features and jury                     206
## 4 sample_h3.csv            H3: Running order and televote                    232
## 5 sample_h4.csv            H4: Special effects and jury-televote g…          232

10. Descriptive Analysis

10.1 Descriptive Statistics

descriptive_statistics <- bind_rows(
  
  # Variables from the 2016-2025 analytical sample
  sample_h1 |>
    summarise(
      variable = "Normalized televote score",
      n = sum(!is.na(tele_score_norm)),
      mean = mean(tele_score_norm, na.rm = TRUE),
      sd = sd(tele_score_norm, na.rm = TRUE),
      min = min(tele_score_norm, na.rm = TRUE),
      max = max(tele_score_norm, na.rm = TRUE)
    ),
  
  sample_h1 |>
    summarise(
      variable = "Special-effects count",
      n = sum(!is.na(special_effects_count)),
      mean = mean(special_effects_count, na.rm = TRUE),
      sd = sd(special_effects_count, na.rm = TRUE),
      min = min(special_effects_count, na.rm = TRUE),
      max = max(special_effects_count, na.rm = TRUE)
    ),
  
  sample_h3 |>
    summarise(
      variable = "Performance order",
      n = sum(!is.na(performance_order)),
      mean = mean(performance_order, na.rm = TRUE),
      sd = sd(performance_order, na.rm = TRUE),
      min = min(performance_order, na.rm = TRUE),
      max = max(performance_order, na.rm = TRUE)
    ),
  
  sample_h4 |>
    summarise(
      variable = "Normalized jury score",
      n = sum(!is.na(jury_score_norm)),
      mean = mean(jury_score_norm, na.rm = TRUE),
      sd = sd(jury_score_norm, na.rm = TRUE),
      min = min(jury_score_norm, na.rm = TRUE),
      max = max(jury_score_norm, na.rm = TRUE)
    ),
  
  sample_h4 |>
    summarise(
      variable = "Jury-televote gap",
      n = sum(!is.na(jury_televote_gap)),
      mean = mean(jury_televote_gap, na.rm = TRUE),
      sd = sd(jury_televote_gap, na.rm = TRUE),
      min = min(jury_televote_gap, na.rm = TRUE),
      max = max(jury_televote_gap, na.rm = TRUE)
    ),
  
  # H2 musical-feature variables
  sample_h2 |>
    summarise(
      variable = "Tempo",
      n = sum(!is.na(tempo)),
      mean = mean(tempo, na.rm = TRUE),
      sd = sd(tempo, na.rm = TRUE),
      min = min(tempo, na.rm = TRUE),
      max = max(tempo, na.rm = TRUE)
    ),
  
  sample_h2 |>
    summarise(
      variable = "Energy",
      n = sum(!is.na(energy)),
      mean = mean(energy, na.rm = TRUE),
      sd = sd(energy, na.rm = TRUE),
      min = min(energy, na.rm = TRUE),
      max = max(energy, na.rm = TRUE)
    ),
  
  sample_h2 |>
    summarise(
      variable = "Danceability",
      n = sum(!is.na(danceability)),
      mean = mean(danceability, na.rm = TRUE),
      sd = sd(danceability, na.rm = TRUE),
      min = min(danceability, na.rm = TRUE),
      max = max(danceability, na.rm = TRUE)
    ),
  
  sample_h2 |>
    summarise(
      variable = "Valence",
      n = sum(!is.na(valence)),
      mean = mean(valence, na.rm = TRUE),
      sd = sd(valence, na.rm = TRUE),
      min = min(valence, na.rm = TRUE),
      max = max(valence, na.rm = TRUE)
    )
) |>
  mutate(
    across(c(mean, sd, min, max), ~ round(.x, 3))
  )

descriptive_statistics

## # A tibble: 9 × 6
##   variable                      n    mean     sd    min     max
##   <chr>                     <int>   <dbl>  <dbl>  <dbl>   <dbl>
## 1 Normalized televote score   232   0.192  0.202  0       0.938
## 2 Special-effects count       232   0.397  0.623  0       2    
## 3 Performance order           232  13.5    7.48   1      26    
## 4 Normalized jury score       232   0.192  0.166  0       0.845
## 5 Jury-televote gap           232   0      0.189 -0.495   0.607
## 6 Tempo                       206 120.    25.4   66.4   205.   
## 7 Energy                      206   0.68   0.178  0.095   0.966
## 8 Danceability                206   0.577  0.141  0.17    0.918
## 9 Valence                     206   0.439  0.189  0.08    0.901

10.2 Special Effects Frequencies

effects_frequencies <- tibble(
  effect = c(
    "Fire",
    "Smoke",
    "Wind machine",
    "Confetti",
    "Any special effect"
  ),
  
  performances = c(
    sum(sample_h1$fire == 1),
    sum(sample_h1$smoke == 1),
    sum(sample_h1$wind_machine == 1),
    sum(sample_h1$confetti == 1),
    sum(sample_h1$special_effects_count > 0)
  )
) |>
  mutate(
    percentage = round(
      performances / nrow(sample_h1) * 100,
      1
    )
  )

effects_frequencies

## # A tibble: 5 × 3
##   effect             performances percentage
##   <chr>                     <int>      <dbl>
## 1 Fire                         28       12.1
## 2 Smoke                        51       22  
## 3 Wind machine                 12        5.2
## 4 Confetti                      1        0.4
## 5 Any special effect           75       32.3

10.3 Outcome Distributions

outcome_distribution_data <- sample_h4 |>
  select(
    jury_score_norm,
    tele_score_norm
  ) |>
  pivot_longer(
    cols = everything(),
    names_to = "outcome",
    values_to = "normalized_score"
  ) |>
  mutate(
    outcome = dplyr::recode(
  outcome,
  jury_score_norm = "Jury",
  tele_score_norm = "Televote"
    )
  )

ggplot(
  outcome_distribution_data,
  aes(
    x = normalized_score,
    linetype = outcome
  )
) +
  geom_density(
    linewidth = 1
  ) +
  labs(
  title = "Distribution of Normalised Jury and Televote Scores",
  x = "Normalised score",
  y = "Density",
  linetype = "Outcome"
) +
  theme_minimal()

ggplot(
  sample_h4,
  aes(x = jury_televote_gap)
) +
  geom_histogram(
    bins = 25,
    boundary = 0
  ) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed"
  ) +
  labs(
  title = "Distribution of the Jury versus Televote Gap",
  x = "Normalised televote score minus normalised jury score",
  y = "Number of performances"
) +
  theme_minimal()

10.4 Correlation Analysis

h2_correlations <- sample_h2 |>
  select(
    jury_score_norm,
    tempo,
    energy,
    danceability,
    valence
  ) |>
  cor(
    use = "complete.obs",
    method = "pearson"
  ) |>
  round(3)

h2_correlations

##                 jury_score_norm  tempo energy danceability valence
## jury_score_norm           1.000 -0.020 -0.227       -0.122  -0.121
## tempo                    -0.020  1.000  0.058       -0.235   0.016
## energy                   -0.227  0.058  1.000        0.309   0.392
## danceability             -0.122 -0.235  0.309        1.000   0.483
## valence                  -0.121  0.016  0.392        0.483   1.000

main_correlations <- sample_h4 |>
  select(
    tele_score_norm,
    jury_score_norm,
    jury_televote_gap,
    special_effects_count,
    performance_order
  ) |>
  cor(
    use = "complete.obs",
    method = "pearson"
  ) |>
  round(3)

main_correlations

##                       tele_score_norm jury_score_norm jury_televote_gap
## tele_score_norm                 1.000           0.487             0.643
## jury_score_norm                 0.487           1.000            -0.357
## jury_televote_gap               0.643          -0.357             1.000
## special_effects_count           0.087          -0.098             0.179
## performance_order               0.175           0.076             0.120
##                       special_effects_count performance_order
## tele_score_norm                       0.087             0.175
## jury_score_norm                      -0.098             0.076
## jury_televote_gap                     0.179             0.120
## special_effects_count                 1.000             0.091
## performance_order                     0.091             1.000

11. Regression Analysis

11.1 Model Specification

regression_specifications <- tibble(
  hypothesis = c("H1", "H2", "H3", "H4"),
  outcome = c(
    "tele_score_norm",
    "jury_score_norm",
    "tele_score_norm",
    "jury_televote_gap"
  ),
  main_predictor = c(
    "special_effects_count",
    "tempo + energy + danceability + valence",
    "performance_order",
    "special_effects_count"
  ),
  control = c(
    "performance_order",
    "None",
    "special_effects_count",
    "performance_order"
  ),
  year_fixed_effects = c("Yes", "Yes", "Yes", "Yes"),
  observations = c(
    nrow(sample_h1),
    nrow(sample_h2),
    nrow(sample_h3),
    nrow(sample_h4)
  )
)

regression_specifications

## # A tibble: 4 × 6
##   hypothesis outcome      main_predictor control year_fixed_effects observations
##   <chr>      <chr>        <chr>          <chr>   <chr>                     <int>
## 1 H1         tele_score_… special_effec… perfor… Yes                         232
## 2 H2         jury_score_… tempo + energ… None    Yes                         206
## 3 H3         tele_score_… performance_o… specia… Yes                         232
## 4 H4         jury_televo… special_effec… perfor… Yes                         232

11.2 H1: Special Effects and Televote

model_h1 <- feols(
  tele_score_norm ~ special_effects_count + performance_order | year,
  data = sample_h1,
  vcov = "HC1"
)

summary(model_h1)

## OLS estimation, Dep. Var.: tele_score_norm
## Observations: 232
## Fixed-effects: year: 9
## Standard-errors: Heteroskedasticity-robust 
##                       Estimate Std. Error t value  Pr(>|t|)    
## special_effects_count 0.025604   0.022939 1.11617 0.2655626    
## performance_order     0.004544   0.001707 2.66194 0.0083399 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.198316     Adj. R2: -0.00722 
##                  Within R2:  0.036139

h1_results <- tidy(model_h1) |>
  select(
    term,
    estimate,
    std.error,
    statistic,
    p.value
  )

h1_results

## # A tibble: 2 × 5
##   term                  estimate std.error statistic p.value
##   <chr>                    <dbl>     <dbl>     <dbl>   <dbl>
## 1 special_effects_count  0.0256    0.0229       1.12 0.266  
## 2 performance_order      0.00454   0.00171      2.66 0.00834

11.3 H2: Musical Features and Jury

model_h2 <- feols(
  jury_score_norm ~ tempo + energy + danceability + valence | year,
  data = sample_h2,
  vcov = "HC1"
)

summary(model_h2)

## OLS estimation, Dep. Var.: jury_score_norm
## Observations: 206
## Fixed-effects: year: 8
## Standard-errors: Heteroskedasticity-robust 
##               Estimate Std. Error   t value Pr(>|t|)    
## tempo        -0.000169   0.000549 -0.308221 0.758245    
## energy       -0.202588   0.081343 -2.490538 0.013595 *  
## danceability -0.079366   0.097233 -0.816243 0.415362    
## valence       0.000501   0.071126  0.007043 0.994388    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.162885     Adj. R2: 0.004079
##                  Within R2: 0.057132

# Equivalent OLS specification with explicit year indicators
# for the purpose of calculating VIFs
h2_vif_model <- lm(
  jury_score_norm ~
    tempo +
    energy +
    danceability +
    valence +
    factor(year),
  data = sample_h2
)

# Variance Inflation Factors
vif(h2_vif_model)

##                  GVIF Df GVIF^(1/(2*Df))
## tempo        1.170208  1        1.081762
## energy       1.313037  1        1.145878
## danceability 1.555184  1        1.247070
## valence      1.557424  1        1.247968
## factor(year) 1.264050  7        1.016878

h2_results <- tidy(model_h2) |>
  select(
    term,
    estimate,
    std.error,
    statistic,
    p.value
  )

h2_results

## # A tibble: 4 × 5
##   term          estimate std.error statistic p.value
##   <chr>            <dbl>     <dbl>     <dbl>   <dbl>
## 1 tempo        -0.000169  0.000549  -0.308    0.758 
## 2 energy       -0.203     0.0813    -2.49     0.0136
## 3 danceability -0.0794    0.0972    -0.816    0.415 
## 4 valence       0.000501  0.0711     0.00704  0.994

# Joint test of all four musical characteristics in H2
wald(
  model_h2,
  keep = c(
    "tempo",
    "energy",
    "danceability",
    "valence"
  )
)

## Wald test, H0: joint nullity of tempo, energy, danceability and valence
##  stat = 2.43095, p-value = 0.049002, on 4 and 194 DoF, VCOV: Heteroskedasticity-robust.

11.4 H3: Running Order and Televote

model_h3 <- feols(
  tele_score_norm ~ performance_order + special_effects_count | year,
  data = sample_h3,
  vcov = "HC1"
)

summary(model_h3)

## OLS estimation, Dep. Var.: tele_score_norm
## Observations: 232
## Fixed-effects: year: 9
## Standard-errors: Heteroskedasticity-robust 
##                       Estimate Std. Error t value  Pr(>|t|)    
## performance_order     0.004544   0.001707 2.66194 0.0083399 ** 
## special_effects_count 0.025604   0.022939 1.11617 0.2655626    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.198316     Adj. R2: -0.00722 
##                  Within R2:  0.036139

h3_results <- tidy(model_h3) |>
  select(
    term,
    estimate,
    std.error,
    statistic,
    p.value
  )

h3_results

## # A tibble: 2 × 5
##   term                  estimate std.error statistic p.value
##   <chr>                    <dbl>     <dbl>     <dbl>   <dbl>
## 1 performance_order      0.00454   0.00171      2.66 0.00834
## 2 special_effects_count  0.0256    0.0229       1.12 0.266

11.5 H4: Special Effects and Jury–Televote Gap

model_h4 <- feols(
  jury_televote_gap ~ special_effects_count + performance_order | year,
  data = sample_h4,
  vcov = "HC1"
)

summary(model_h4)

## OLS estimation, Dep. Var.: jury_televote_gap
## Observations: 232
## Fixed-effects: year: 9
## Standard-errors: Heteroskedasticity-robust 
##                       Estimate Std. Error t value  Pr(>|t|)    
## special_effects_count 0.056863   0.021305 2.66900 0.0081721 ** 
## performance_order     0.002633   0.001542 1.70746 0.0891413 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.184534     Adj. R2: 0.002769
##                  Within R2: 0.045939

h4_results <- tidy(model_h4) |>
  select(
    term,
    estimate,
    std.error,
    statistic,
    p.value
  )

h4_results

## # A tibble: 2 × 5
##   term                  estimate std.error statistic p.value
##   <chr>                    <dbl>     <dbl>     <dbl>   <dbl>
## 1 special_effects_count  0.0569    0.0213       2.67 0.00817
## 2 performance_order      0.00263   0.00154      1.71 0.0891

11.6 Regression Diagnostics

# Equivalent lm() models used only for regression diagnostics

diag_h1 <- lm(
  tele_score_norm ~
    special_effects_count +
    performance_order +
    factor(year),
  data = sample_h1
)

diag_h2 <- lm(
  jury_score_norm ~
    tempo +
    energy +
    danceability +
    valence +
    factor(year),
  data = sample_h2
)

diag_h3 <- lm(
  tele_score_norm ~
    performance_order +
    special_effects_count +
    factor(year),
  data = sample_h3
)

diag_h4 <- lm(
  jury_televote_gap ~
    special_effects_count +
    performance_order +
    factor(year),
  data = sample_h4
)

# Count observations exceeding the common Cook's distance threshold 4/N
influence_summary <- tibble(
  hypothesis = c("H1", "H2", "H3", "H4"),
  observations = c(
    nrow(sample_h1),
    nrow(sample_h2),
    nrow(sample_h3),
    nrow(sample_h4)
  ),
  cooks_threshold = c(
    4 / nrow(sample_h1),
    4 / nrow(sample_h2),
    4 / nrow(sample_h3),
    4 / nrow(sample_h4)
  ),
  influential_observations = c(
    sum(cooks.distance(diag_h1) > 4 / nrow(sample_h1)),
    sum(cooks.distance(diag_h2) > 4 / nrow(sample_h2)),
    sum(cooks.distance(diag_h3) > 4 / nrow(sample_h3)),
    sum(cooks.distance(diag_h4) > 4 / nrow(sample_h4))
  ),
  maximum_cooks_distance = c(
    max(cooks.distance(diag_h1)),
    max(cooks.distance(diag_h2)),
    max(cooks.distance(diag_h3)),
    max(cooks.distance(diag_h4))
  )
)

influence_summary

## # A tibble: 4 × 5
##   hypothesis observations cooks_threshold influential_observations
##   <chr>             <int>           <dbl>                    <int>
## 1 H1                  232          0.0172                       14
## 2 H2                  206          0.0194                       10
## 3 H3                  232          0.0172                       14
## 4 H4                  232          0.0172                       16
## # ℹ 1 more variable: maximum_cooks_distance <dbl>

ggplot(
  data.frame(
    fitted = fitted(diag_h1),
    residuals = residuals(diag_h1)
  ),
  aes(
    x = fitted,
    y = residuals
  )
) +
  geom_point() +
  geom_smooth(
    method = "loess",
    se = FALSE
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dashed"
  ) +
  labs(
    title = "Residuals versus Fitted Values: H1/H3 Model",
    x = "Fitted normalised televote score",
    y = "Residual"
  ) +
  theme_minimal()

ggplot(
  data.frame(
    fitted = fitted(diag_h2),
    residuals = residuals(diag_h2)
  ),
  aes(
    x = fitted,
    y = residuals
  )
) +
  geom_point() +
  geom_smooth(
    method = "loess",
    se = FALSE
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dashed"
  ) +
  labs(
    title = "Residuals versus Fitted Values: H2 Model",
    x = "Fitted normalised jury score",
    y = "Residual"
  ) +
  theme_minimal()

ggplot(
  data.frame(
    fitted = fitted(diag_h4),
    residuals = residuals(diag_h4)
  ),
  aes(
    x = fitted,
    y = residuals
  )
) +
  geom_point() +
  geom_smooth(
    method = "loess",
    se = FALSE
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dashed"
  ) +
  labs(
    title = "Residuals versus Fitted Values: H4 Model",
    x = "Fitted jury–televote gap",
    y = "Residual"
  ) +
  theme_minimal()

11.7 Robustness Checks

11.7.1 Binary Special-Effects Measure

# Create binary special-effects variable for H1
sample_h1_binary <- sample_h1 |>
  mutate(
    any_special_effect = if_else(
      special_effects_count > 0,
      1,
      0
    )
  )

# Create the same binary variable for H4
sample_h4_binary <- sample_h4 |>
  mutate(
    any_special_effect = if_else(
      special_effects_count > 0,
      1,
      0
    )
  )

# Check the distribution
sample_h1_binary |>
  count(
    any_special_effect,
    name = "observations"
  ) |>
  mutate(
    percentage = round(
      observations / nrow(sample_h1_binary) * 100,
      1
    )
  )

## # A tibble: 2 × 3
##   any_special_effect observations percentage
##                <dbl>        <int>      <dbl>
## 1                  0          157       67.7
## 2                  1           75       32.3

# H1 robustness: any special effect -> normalized televote score
model_h1_binary <- feols(
  tele_score_norm ~ any_special_effect + performance_order | year,
  data = sample_h1_binary,
  vcov = "HC1"
)

summary(model_h1_binary)

## OLS estimation, Dep. Var.: tele_score_norm
## Observations: 232
## Fixed-effects: year: 9
## Standard-errors: Heteroskedasticity-robust 
##                    Estimate Std. Error t value  Pr(>|t|)    
## any_special_effect 0.032281   0.030486 1.05887 0.2908131    
## performance_order  0.004506   0.001709 2.63684 0.0089621 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.19838     Adj. R2: -0.007867
##                 Within R2:  0.035519

# H4 robustness: any special effect -> jury-televote gap
model_h4_binary <- feols(
  jury_televote_gap ~ any_special_effect + performance_order | year,
  data = sample_h4_binary,
  vcov = "HC1"
)

summary(model_h4_binary)

## OLS estimation, Dep. Var.: jury_televote_gap
## Observations: 232
## Fixed-effects: year: 9
## Standard-errors: Heteroskedasticity-robust 
##                    Estimate Std. Error t value Pr(>|t|)    
## any_special_effect 0.063785   0.027968 2.28062 0.023524 *  
## performance_order  0.002604   0.001543 1.68693 0.093028 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.185432     Adj. R2: -0.006969
##                  Within R2:  0.036623

11.7.2 Influence Robustness

# H3 robustness sample:
# exclude observations flagged by Cook's distance in the H1/H3 diagnostic model
sample_h3_influence <- sample_h3 |>
  mutate(
    cooks_distance = cooks.distance(diag_h3)
  ) |>
  filter(
    cooks_distance <= 4 / nrow(sample_h3)
  )

# H4 robustness sample:
# exclude observations flagged by Cook's distance in the H4 diagnostic model
sample_h4_influence <- sample_h4 |>
  mutate(
    cooks_distance = cooks.distance(diag_h4)
  ) |>
  filter(
    cooks_distance <= 4 / nrow(sample_h4)
  )

# Compare primary and influence-robustness sample sizes
tibble(
  hypothesis = c("H3", "H4"),
  primary_n = c(
    nrow(sample_h3),
    nrow(sample_h4)
  ),
  robustness_n = c(
    nrow(sample_h3_influence),
    nrow(sample_h4_influence)
  ),
  excluded = c(
    nrow(sample_h3) - nrow(sample_h3_influence),
    nrow(sample_h4) - nrow(sample_h4_influence)
  )
)

## # A tibble: 2 × 4
##   hypothesis primary_n robustness_n excluded
##   <chr>          <int>        <int>    <int>
## 1 H3               232          218       14
## 2 H4               232          216       16

# H3 influence robustness
model_h3_influence <- feols(
  tele_score_norm ~ performance_order + special_effects_count | year,
  data = sample_h3_influence,
  vcov = "HC1"
)

summary(model_h3_influence)

## OLS estimation, Dep. Var.: tele_score_norm
## Observations: 218
## Fixed-effects: year: 9
## Standard-errors: Heteroskedasticity-robust 
##                       Estimate Std. Error  t value   Pr(>|t|)    
## performance_order     0.004606   0.001342 3.430792 0.00072643 ***
## special_effects_count 0.015553   0.017355 0.896183 0.37119612    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.151199     Adj. R2: 0.016075
##                  Within R2: 0.053135

# H4 influence robustness
model_h4_influence <- feols(
  jury_televote_gap ~ special_effects_count + performance_order | year,
  data = sample_h4_influence,
  vcov = "HC1"
)

summary(model_h4_influence)

## OLS estimation, Dep. Var.: jury_televote_gap
## Observations: 216
## Fixed-effects: year: 9
## Standard-errors: Heteroskedasticity-robust 
##                       Estimate Std. Error t value   Pr(>|t|)    
## special_effects_count 0.071066   0.015451 4.59942 7.4104e-06 ***
## performance_order     0.003626   0.001233 2.93970 3.6623e-03 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.140621     Adj. R2: 0.085328
##                  Within R2: 0.117675

11.7.3 H2 Nonlinearity

model_h2_nonlinear <- feols(
  jury_score_norm ~
    tempo + I(tempo^2) +
    energy + I(energy^2) +
    danceability + I(danceability^2) +
    valence + I(valence^2) |
    year,
  data = sample_h2,
  vcov = "HC1"
)

summary(model_h2_nonlinear)

## OLS estimation, Dep. Var.: jury_score_norm
## Observations: 206
## Fixed-effects: year: 8
## Standard-errors: Heteroskedasticity-robust 
##                      Estimate Std. Error   t value Pr(>|t|) 
## tempo              0.00058653   0.004029  0.145586  0.88440 
## I(tempo^2)        -0.00000364   0.000016 -0.230127  0.81824 
## energy            -0.66368161   0.484309 -1.370369  0.17219 
## I(energy^2)        0.38018907   0.367637  1.034143  0.30238 
## danceability      -0.58716161   0.587399 -0.999596  0.31878 
## I(danceability^2)  0.46035797   0.501465  0.918026  0.35977 
## valence            0.15203006   0.294396  0.516413  0.60617 
## I(valence^2)      -0.16803508   0.290067 -0.579297  0.56307 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.161775     Adj. R2: -0.003079
##                  Within R2:  0.069936

12. Random Forest Analysis

12.1 Model Setup and Validation Strategy

rf_data <- analysis_data |>
  filter(year >= 2016, year <= 2024) |>
  drop_na(
    tele_score_norm,
    jury_score_norm,
    performance_order,
    fire,
    smoke,
    wind_machine,
    tempo,
    energy,
    danceability,
    valence
  ) |>
  select(
    year,
    to_country,
    tele_score_norm,
    jury_score_norm,
    performance_order,
    fire,
    smoke,
    wind_machine,
    tempo,
    energy,
    danceability,
    valence
  )

# Check dimensions and contest coverage
rf_data |>
  summarise(
    observations = n(),
    first_year = min(year),
    last_year = max(year),
    number_of_contests = n_distinct(year)
  )

## # A tibble: 1 × 4
##   observations first_year last_year number_of_contests
##          <int>      <dbl>     <dbl>              <int>
## 1          206       2016      2024                  8

# Check observations per contest
rf_data |>
  count(year, name = "observations") |>
  arrange(year)

## # A tibble: 8 × 2
##    year observations
##   <dbl>        <int>
## 1  2016           26
## 2  2017           26
## 3  2018           26
## 4  2019           26
## 5  2021           26
## 6  2022           25
## 7  2023           26
## 8  2024           25

# Confirm that all RF variables are complete
rf_data |>
  summarise(
    across(
      everything(),
      ~ sum(is.na(.x))
    )
  )

## # A tibble: 1 × 12
##    year to_country tele_score_norm jury_score_norm performance_order  fire smoke
##   <int>      <int>           <int>           <int>             <int> <int> <int>
## 1     0          0               0               0                 0     0     0
## # ℹ 5 more variables: wind_machine <int>, tempo <int>, energy <int>,
## #   danceability <int>, valence <int>

set.seed(123)

# Create leave-one-contest-year-out cross-validation folds
rf_folds <- group_vfold_cv(
  rf_data,
  group = year,
  v = 8
)

# Inspect the folds
rf_folds

## # Group 8-fold cross-validation 
## # A tibble: 8 × 2
##   splits           id       
##   <list>           <chr>    
## 1 <split [180/26]> Resample1
## 2 <split [180/26]> Resample2
## 3 <split [181/25]> Resample3
## 4 <split [180/26]> Resample4
## 5 <split [180/26]> Resample5
## 6 <split [180/26]> Resample6
## 7 <split [181/25]> Resample7
## 8 <split [180/26]> Resample8

# Verify which contest year is held out in each fold
rf_fold_summary <- tibble(
  fold = rf_folds$id,
  test_year = map_dbl(
    rf_folds$splits,
    ~ unique(assessment(.x)$year)
  ),
  training_observations = map_int(
    rf_folds$splits,
    ~ nrow(analysis(.x))
  ),
  test_observations = map_int(
    rf_folds$splits,
    ~ nrow(assessment(.x))
  )
) |>
  arrange(test_year)

rf_fold_summary

## # A tibble: 8 × 4
##   fold      test_year training_observations test_observations
##   <chr>         <dbl>                 <int>             <int>
## 1 Resample1      2016                   180                26
## 2 Resample5      2017                   180                26
## 3 Resample4      2018                   180                26
## 4 Resample2      2019                   180                26
## 5 Resample6      2021                   180                26
## 6 Resample3      2022                   181                25
## 7 Resample8      2023                   180                26
## 8 Resample7      2024                   181                25

rf_grid <- expand_grid(
  mtry = c(2, 3, 4, 5, 6),
  min_node_size = c(3, 5, 10)
)

rf_grid

## # A tibble: 15 × 2
##     mtry min_node_size
##    <dbl>         <dbl>
##  1     2             3
##  2     2             5
##  3     2            10
##  4     3             3
##  5     3             5
##  6     3            10
##  7     4             3
##  8     4             5
##  9     4            10
## 10     5             3
## 11     5             5
## 12     5            10
## 13     6             3
## 14     6             5
## 15     6            10

12.2 Televote Random Forest

set.seed(123)

tele_rf_predictions <- vector(
  mode = "list",
  length = nrow(rf_folds)
)

tele_rf_best_params <- vector(
  mode = "list",
  length = nrow(rf_folds)
)

for (i in seq_len(nrow(rf_folds))) {
  
  # -----------------------------
  # Outer split:
  # one whole contest year is held out
  # -----------------------------
  
  outer_split <- rf_folds$splits[[i]]
  
  outer_train <- analysis(outer_split)
  outer_test <- assessment(outer_split)
  
  held_out_year <- unique(outer_test$year)
  
  
  # -----------------------------
  # Inner leave-one-year-out CV
  # using only the outer training data
  # -----------------------------
  
  inner_folds <- group_vfold_cv(
    outer_train,
    group = year,
    v = n_distinct(outer_train$year)
  )
  
  
  # -----------------------------
  # Tune mtry and minimum node size
  # -----------------------------
  
  tuning_results <- rf_grid |>
    mutate(
      rmse = map2_dbl(
        mtry,
        min_node_size,
        function(current_mtry, current_node_size) {
          
          inner_rmse <- map_dbl(
            inner_folds$splits,
            function(inner_split) {
              
              inner_train <- analysis(inner_split)
              inner_test <- assessment(inner_split)
              
              rf_fit <- ranger(
                tele_score_norm ~
                  performance_order +
                  fire +
                  smoke +
                  wind_machine +
                  tempo +
                  energy +
                  danceability +
                  valence,
                data = inner_train,
                num.trees = 1000,
                mtry = current_mtry,
                min.node.size = current_node_size,
                seed = 123
              )
              
              predictions <- predict(
                rf_fit,
                data = inner_test
              )$predictions
              
              sqrt(
                mean(
                  (inner_test$tele_score_norm - predictions)^2
                )
              )
            }
          )
          
          mean(inner_rmse)
        }
      )
    )
  
  
  # -----------------------------
  # Select best hyperparameters
  # -----------------------------
  
  best_params <- tuning_results |>
    slice_min(
      rmse,
      n = 1,
      with_ties = FALSE
    )
  
  
  # -----------------------------
  # Fit final model for this outer fold
  # -----------------------------
  
  final_rf <- ranger(
    tele_score_norm ~
      performance_order +
      fire +
      smoke +
      wind_machine +
      tempo +
      energy +
      danceability +
      valence,
    data = outer_train,
    num.trees = 1000,
    mtry = best_params$mtry,
    min.node.size = best_params$min_node_size,
    seed = 123
  )
  
  
  # -----------------------------
  # Predict the unseen contest year
  # -----------------------------
  
  outer_predictions <- predict(
    final_rf,
    data = outer_test
  )$predictions
  
  
  # Save predictions
  tele_rf_predictions[[i]] <- tibble(
    year = outer_test$year,
    to_country = outer_test$to_country,
    observed = outer_test$tele_score_norm,
    predicted = outer_predictions
  )
  
  
  # Save selected hyperparameters
  tele_rf_best_params[[i]] <- tibble(
    test_year = held_out_year,
    mtry = best_params$mtry,
    min_node_size = best_params$min_node_size,
    inner_rmse = best_params$rmse
  )
}


# Combine all outer-fold predictions
tele_rf_predictions <- bind_rows(
  tele_rf_predictions
)

# Combine selected hyperparameters
tele_rf_best_params <- bind_rows(
  tele_rf_best_params
)

tele_rf_best_params

## # A tibble: 8 × 4
##   test_year  mtry min_node_size inner_rmse
##       <dbl> <dbl>         <dbl>      <dbl>
## 1      2016     2            10      0.211
## 2      2019     2            10      0.210
## 3      2022     2            10      0.210
## 4      2018     2            10      0.219
## 5      2017     2            10      0.207
## 6      2021     2            10      0.209
## 7      2024     2            10      0.200
## 8      2023     2            10      0.209

# Confirm that every RF observation received one out-of-sample prediction
nrow(tele_rf_predictions)

## [1] 206

# Calculate overall out-of-sample predictive performance
tele_rf_performance <- tele_rf_predictions |>
  summarise(
    RMSE = sqrt(
      mean((observed - predicted)^2)
    ),
    MAE = mean(
      abs(observed - predicted)
    ),
    R_squared = 1 -
      sum((observed - predicted)^2) /
      sum((observed - mean(observed))^2)
  )

tele_rf_performance

## # A tibble: 1 × 3
##    RMSE   MAE R_squared
##   <dbl> <dbl>     <dbl>
## 1 0.211 0.163   -0.0647

ggplot(
  tele_rf_predictions,
  aes(
    x = observed,
    y = predicted
  )
) +
  geom_point() +
  geom_abline(
    intercept = 0,
    slope = 1,
    linetype = "dashed"
  ) +
  labs(
  title = "Observed versus Predicted Televote Scores",
  subtitle = "Leave one contest year out predictions",
  x = "Observed normalised televote score",
  y = "Predicted normalised televote score"
) +
  theme_minimal()

12.3 Jury Random Forest

set.seed(123)

jury_rf_predictions <- vector(
  mode = "list",
  length = nrow(rf_folds)
)

jury_rf_best_params <- vector(
  mode = "list",
  length = nrow(rf_folds)
)

for (i in seq_len(nrow(rf_folds))) {
  
  # -----------------------------
  # Outer split:
  # one whole contest year is held out
  # -----------------------------
  
  outer_split <- rf_folds$splits[[i]]
  
  outer_train <- analysis(outer_split)
  outer_test <- assessment(outer_split)
  
  held_out_year <- unique(outer_test$year)
  
  
  # -----------------------------
  # Inner leave-one-year-out CV
  # using only the outer training data
  # -----------------------------
  
  inner_folds <- group_vfold_cv(
    outer_train,
    group = year,
    v = n_distinct(outer_train$year)
  )
  
  
  # -----------------------------
  # Tune mtry and minimum node size
  # -----------------------------
  
  tuning_results <- rf_grid |>
    mutate(
      rmse = map2_dbl(
        mtry,
        min_node_size,
        function(current_mtry, current_node_size) {
          
          inner_rmse <- map_dbl(
            inner_folds$splits,
            function(inner_split) {
              
              inner_train <- analysis(inner_split)
              inner_test <- assessment(inner_split)
              
              rf_fit <- ranger(
                jury_score_norm ~
                  performance_order +
                  fire +
                  smoke +
                  wind_machine +
                  tempo +
                  energy +
                  danceability +
                  valence,
                data = inner_train,
                num.trees = 1000,
                mtry = current_mtry,
                min.node.size = current_node_size,
                seed = 123
              )
              
              predictions <- predict(
                rf_fit,
                data = inner_test
              )$predictions
              
              sqrt(
                mean(
                  (inner_test$jury_score_norm - predictions)^2
                )
              )
            }
          )
          
          mean(inner_rmse)
        }
      )
    )
  
  
  # -----------------------------
  # Select best hyperparameters
  # -----------------------------
  
  best_params <- tuning_results |>
    slice_min(
      rmse,
      n = 1,
      with_ties = FALSE
    )
  
  
  # -----------------------------
  # Fit final model for this outer fold
  # -----------------------------
  
  final_rf <- ranger(
    jury_score_norm ~
      performance_order +
      fire +
      smoke +
      wind_machine +
      tempo +
      energy +
      danceability +
      valence,
    data = outer_train,
    num.trees = 1000,
    mtry = best_params$mtry,
    min.node.size = best_params$min_node_size,
    seed = 123
  )
  
  
  # -----------------------------
  # Predict the unseen contest year
  # -----------------------------
  
  outer_predictions <- predict(
    final_rf,
    data = outer_test
  )$predictions
  
  
  # Save predictions
  jury_rf_predictions[[i]] <- tibble(
    year = outer_test$year,
    to_country = outer_test$to_country,
    observed = outer_test$jury_score_norm,
    predicted = outer_predictions
  )
  
  
  # Save selected hyperparameters
  jury_rf_best_params[[i]] <- tibble(
    test_year = held_out_year,
    mtry = best_params$mtry,
    min_node_size = best_params$min_node_size,
    inner_rmse = best_params$rmse
  )
}


# Combine all outer-fold predictions
jury_rf_predictions <- bind_rows(
  jury_rf_predictions
)

# Combine selected hyperparameters
jury_rf_best_params <- bind_rows(
  jury_rf_best_params
)

jury_rf_best_params

## # A tibble: 8 × 4
##   test_year  mtry min_node_size inner_rmse
##       <dbl> <dbl>         <dbl>      <dbl>
## 1      2016     2            10      0.178
## 2      2019     2            10      0.174
## 3      2022     2            10      0.176
## 4      2018     2            10      0.176
## 5      2017     2            10      0.175
## 6      2021     2            10      0.172
## 7      2024     2            10      0.170
## 8      2023     2            10      0.172

# Confirm that every RF observation received one out-of-sample prediction
nrow(jury_rf_predictions)

## [1] 206

# Calculate overall out-of-sample predictive performance
jury_rf_performance <- jury_rf_predictions |>
  summarise(
    RMSE = sqrt(
      mean((observed - predicted)^2)
    ),
    MAE = mean(
      abs(observed - predicted)
    ),
    R_squared = 1 -
      sum((observed - predicted)^2) /
      sum((observed - mean(observed))^2)
  )

jury_rf_performance

## # A tibble: 1 × 3
##    RMSE   MAE R_squared
##   <dbl> <dbl>     <dbl>
## 1 0.175 0.139   -0.0822

ggplot(
  jury_rf_predictions,
  aes(
    x = observed,
    y = predicted
  )
) +
  geom_point() +
  geom_abline(
    intercept = 0,
    slope = 1,
    linetype = "dashed"
  ) +
  labs(
  title = "Observed versus Predicted Jury Scores",
  subtitle = "Leave one contest year out predictions",
  x = "Observed normalised jury score",
  y = "Predicted normalised jury score"
) +
  theme_minimal()

12.4 Predictive Performance

rf_performance_comparison <- bind_rows(
  
  tele_rf_performance |>
    mutate(
      outcome = "Televote"
    ),
  
  jury_rf_performance |>
    mutate(
      outcome = "Jury"
    )
  
) |>
  select(
    outcome,
    RMSE,
    MAE,
    R_squared
  ) |>
  mutate(
    across(
      c(RMSE, MAE, R_squared),
      ~ round(.x, 3)
    )
  )

rf_performance_comparison

## # A tibble: 2 × 4
##   outcome   RMSE   MAE R_squared
##   <chr>    <dbl> <dbl>     <dbl>
## 1 Televote 0.211 0.163    -0.065
## 2 Jury     0.175 0.139    -0.082

12.5 Variable Importance

set.seed(123)

# Final televote forest for permutation importance
tele_rf_importance_model <- ranger(
  tele_score_norm ~
    performance_order +
    fire +
    smoke +
    wind_machine +
    tempo +
    energy +
    danceability +
    valence,
  data = rf_data,
  num.trees = 1000,
  mtry = 2,
  min.node.size = 10,
  importance = "permutation",
  seed = 123
)


# Final jury forest for permutation importance
jury_rf_importance_model <- ranger(
  jury_score_norm ~
    performance_order +
    fire +
    smoke +
    wind_machine +
    tempo +
    energy +
    danceability +
    valence,
  data = rf_data,
  num.trees = 1000,
  mtry = 2,
  min.node.size = 10,
  importance = "permutation",
  seed = 123
)


# Extract televote variable importance
tele_importance <- tibble(
  variable = names(
    tele_rf_importance_model$variable.importance
  ),
  importance = as.numeric(
    tele_rf_importance_model$variable.importance
  )
) |>
  arrange(desc(importance))

tele_importance

## # A tibble: 8 × 2
##   variable           importance
##   <chr>                   <dbl>
## 1 danceability       0.00233   
## 2 performance_order  0.00154   
## 3 wind_machine       0.000889  
## 4 valence           -0.00000629
## 5 smoke             -0.000338  
## 6 fire              -0.000405  
## 7 energy            -0.000843  
## 8 tempo             -0.00174

# Extract jury variable importance
jury_importance <- tibble(
  variable = names(
    jury_rf_importance_model$variable.importance
  ),
  importance = as.numeric(
    jury_rf_importance_model$variable.importance
  )
) |>
  arrange(desc(importance))

jury_importance

## # A tibble: 8 × 2
##   variable          importance
##   <chr>                  <dbl>
## 1 energy             0.00125  
## 2 danceability       0.000225 
## 3 fire               0.0000455
## 4 performance_order -0.0000481
## 5 wind_machine      -0.000260 
## 6 tempo             -0.000348 
## 7 valence           -0.000387 
## 8 smoke             -0.000538

rf_importance_combined <- bind_rows(
  tele_importance |>
    mutate(
      outcome = "Televote"
    ),
  
  jury_importance |>
    mutate(
      outcome = "Jury"
    )
) |>
  mutate(
    variable = dplyr::recode(
      variable,
      performance_order = "Performance order",
      fire = "Fire",
      smoke = "Smoke",
      wind_machine = "Wind machine",
      tempo = "Tempo",
      energy = "Energy",
      danceability = "Danceability",
      valence = "Valence"
    )
  )


ggplot(
  rf_importance_combined,
  aes(
    x = importance,
    y = reorder(variable, importance)
  )
) +
  geom_col() +
  geom_vline(
    xintercept = 0,
    linetype = "dashed"
  ) +
  facet_wrap(
    ~ outcome,
    scales = "free_y"
  ) +
  labs(
    title = "Random Forest Permutation Variable Importance",
    subtitle = "Televote and jury prediction models",
    x = "Permutation importance",
    y = NULL
  ) +
  theme_minimal()

13. Final Tables and Figures

13.1 Regression Tables

primary_regression_table <- etable(
  "H1: Televote" = model_h1,
  "H2: Jury" = model_h2,
  "H3: Televote" = model_h3,
  "H4: Jury-Televote Gap" = model_h4,
  
  dict = c(
    special_effects_count = "Special-effects count",
    performance_order = "Performance order",
    tempo = "Tempo",
    energy = "Energy",
    danceability = "Danceability",
    valence = "Valence"
  ),
  
  fitstat = ~ n + wr2,
  
  digits = 3,
  signif.code = c(
    "***" = 0.001,
    "**" = 0.01,
    "*" = 0.05,
    "." = 0.10
  )
)

primary_regression_table

##                          H1: Televote         H2: Jury    H3: Televote
## Dependent Var.:       tele_score_norm  jury_score_norm tele_score_norm
##                                                                       
## Special-effects count   0.026 (0.023)                    0.026 (0.023)
## Performance order     0.005** (0.002)                  0.005** (0.002)
## Tempo                                 -0.0002 (0.0005)                
## Energy                                 -0.203* (0.081)                
## Danceability                            -0.079 (0.097)                
## Valence                                 0.0005 (0.071)                
## Fixed-Effects:        --------------- ---------------- ---------------
## year                              Yes              Yes             Yes
## _____________________ _______________ ________________ _______________
## S.E. type             Heteroske.-rob. Heterosked.-rob. Heteroske.-rob.
## Observations                      232              206             232
## Within R2                     0.03614          0.05713         0.03614
## 
##                       H4: Jury-Televo..
## Dependent Var.:       jury_televote_gap
##                                        
## Special-effects count   0.057** (0.021)
## Performance order        0.003. (0.002)
## Tempo                                  
## Energy                                 
## Danceability                           
## Valence                                
## Fixed-Effects:        -----------------
## year                                Yes
## _____________________ _________________
## S.E. type             Heteroskeda.-rob.
## Observations                        232
## Within R2                       0.04594
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

13.2 Hypothesis Figures

# --------------------------------------------------
# Create standardized datasets for visualization only
# --------------------------------------------------

h1_standardized <- sample_h1 |>
  mutate(
    tele_score_z = as.numeric(scale(tele_score_norm)),
    special_effects_z = as.numeric(scale(special_effects_count)),
    performance_order_z = as.numeric(scale(performance_order))
  )

h2_standardized <- sample_h2 |>
  mutate(
    jury_score_z = as.numeric(scale(jury_score_norm)),
    tempo_z = as.numeric(scale(tempo)),
    energy_z = as.numeric(scale(energy)),
    danceability_z = as.numeric(scale(danceability)),
    valence_z = as.numeric(scale(valence))
  )

h3_standardized <- sample_h3 |>
  mutate(
    tele_score_z = as.numeric(scale(tele_score_norm)),
    performance_order_z = as.numeric(scale(performance_order)),
    special_effects_z = as.numeric(scale(special_effects_count))
  )

h4_standardized <- sample_h4 |>
  mutate(
    gap_z = as.numeric(scale(jury_televote_gap)),
    special_effects_z = as.numeric(scale(special_effects_count)),
    performance_order_z = as.numeric(scale(performance_order))
  )


# --------------------------------------------------
# Re-estimate primary models using standardized
# variables for coefficient visualization
# --------------------------------------------------

model_h1_standardized <- feols(
  tele_score_z ~
    special_effects_z +
    performance_order_z |
    year,
  data = h1_standardized,
  vcov = "HC1"
)

model_h2_standardized <- feols(
  jury_score_z ~
    tempo_z +
    energy_z +
    danceability_z +
    valence_z |
    year,
  data = h2_standardized,
  vcov = "HC1"
)

model_h3_standardized <- feols(
  tele_score_z ~
    performance_order_z +
    special_effects_z |
    year,
  data = h3_standardized,
  vcov = "HC1"
)

model_h4_standardized <- feols(
  gap_z ~
    special_effects_z +
    performance_order_z |
    year,
  data = h4_standardized,
  vcov = "HC1"
)


# --------------------------------------------------
# Extract only the focal hypothesis coefficients
# --------------------------------------------------

coefficient_plot_data <- bind_rows(
  
  tidy(model_h1_standardized) |>
    filter(term == "special_effects_z") |>
    mutate(
      hypothesis = "H1",
      predictor = "Special effects"
    ),
  
  tidy(model_h2_standardized) |>
    filter(
      term %in% c(
        "tempo_z",
        "energy_z",
        "danceability_z",
        "valence_z"
      )
    ) |>
    mutate(
      hypothesis = "H2",
      predictor = dplyr::recode(
        term,
        tempo_z = "Tempo",
        energy_z = "Energy",
        danceability_z = "Danceability",
        valence_z = "Valence"
      )
    ),
  
  tidy(model_h3_standardized) |>
    filter(term == "performance_order_z") |>
    mutate(
      hypothesis = "H3",
      predictor = "Performance order"
    ),
  
  tidy(model_h4_standardized) |>
    filter(term == "special_effects_z") |>
    mutate(
      hypothesis = "H4",
      predictor = "Special effects"
    )
  
) |>
  mutate(
    lower_95 = estimate - 1.96 * std.error,
    upper_95 = estimate + 1.96 * std.error,
    label = paste0(
      hypothesis,
      ": ",
      predictor
    )
  )


# --------------------------------------------------
# Plot standardized coefficients and 95% CIs
# --------------------------------------------------

ggplot(
  coefficient_plot_data,
  aes(
    x = estimate,
    y = reorder(label, estimate)
  )
) +
  geom_point(
    size = 2.5
  ) +
  geom_errorbarh(
    aes(
      xmin = lower_95,
      xmax = upper_95
    ),
    height = 0.15
  ) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed"
  ) +
  labs(
    title = "Standardized Coefficients of the Primary Hypothesis Models",
    subtitle = "Points show estimates; horizontal lines show 95% confidence intervals",
    x = "Standardized coefficient",
    y = NULL
  ) +
  theme_minimal()

13.3 Random Forest Tables and Figures

final_rf_performance_table <- rf_performance_comparison |>
  rename(
    Outcome = outcome,
    RMSE = RMSE,
    MAE = MAE,
    `Out-of-sample R²` = R_squared
  )

final_rf_performance_table

## # A tibble: 2 × 4
##   Outcome   RMSE   MAE `Out-of-sample R²`
##   <chr>    <dbl> <dbl>              <dbl>
## 1 Televote 0.211 0.163             -0.065
## 2 Jury     0.175 0.139             -0.082

13.4 Hypothesis Summary

hypothesis_summary <- tibble(
  hypothesis = c(
    "H1",
    "H2",
    "H3",
    "H4"
  ),
  
  wording = c(
    "Performances using more special effects are associated with higher televote scores.",
    "Jury scores are associated with the musical characteristics tempo, energy, danceability, and valence.",
    "Later positions in the Grand Final running order are associated with higher televote scores.",
    "Performances using more special effects are associated with relatively higher televote than jury scores."
  ),
  
  main_result = c(
    "Special effects: beta = 0.026, p = .266",
    "Joint Wald test: p = .049; energy: beta = -0.203, p = .014",
    "Performance order: beta = 0.0045, p = .008",
    "Special effects: beta = 0.057, p = .008"
  ),
  
  conclusion = c(
    "Not supported",
    "Supported with qualification",
    "Supported",
    "Supported"
  )
)

hypothesis_summary

## # A tibble: 4 × 4
##   hypothesis wording                                      main_result conclusion
##   <chr>      <chr>                                        <chr>       <chr>     
## 1 H1         Performances using more special effects are… Special ef… Not suppo…
## 2 H2         Jury scores are associated with the musical… Joint Wald… Supported…
## 3 H3         Later positions in the Grand Final running … Performanc… Supported 
## 4 H4         Performances using more special effects are… Special ef… Supported

14. Final Table Exports

14.1 Main-Text Tables

write_csv(
  descriptive_statistics,
  here("tables", "table_descriptive_statistics.csv")
)

special_effects_export <- tibble(
  effect = c(
    "Fire",
    "Smoke",
    "Wind machine",
    "Confetti",
    "Any special effect"
  ),
  performances = c(
    sum(sample_h1$fire == 1),
    sum(sample_h1$smoke == 1),
    sum(sample_h1$wind_machine == 1),
    sum(sample_h1$confetti == 1),
    sum(sample_h1$special_effects_count > 0)
  )
) |>
  mutate(
    percentage = round(
      performances / nrow(sample_h1) * 100,
      1
    )
  )

write_csv(
  special_effects_export,
  here("tables", "table_special_effects_frequencies.csv")
)

primary_regression_export <- bind_rows(
  
  tidy(model_h1) |>
    mutate(model = "H1"),
  
  tidy(model_h2) |>
    mutate(model = "H2"),
  
  tidy(model_h3) |>
    mutate(model = "H3"),
  
  tidy(model_h4) |>
    mutate(model = "H4")
  
) |>
  select(
    model,
    term,
    estimate,
    std.error,
    statistic,
    p.value
  )

write_csv(
  primary_regression_export,
  here("tables", "table_primary_regressions.csv")
)

write_csv(
  final_rf_performance_table,
  here("tables", "table_rf_performance.csv")
)

write_csv(
  hypothesis_summary,
  here("tables", "table_hypothesis_summary.csv")
)

14.2 Appendix: Correlations and Diagnostics

appendix_h2_correlations <- sample_h2 |>
  select(
    jury_score_norm,
    tempo,
    energy,
    danceability,
    valence
  ) |>
  cor() |>
  round(3) |>
  as.data.frame() |>
  rownames_to_column(
    var = "variable"
  )

write_csv(
  appendix_h2_correlations,
  here("tables", "appendix_h2_correlations.csv")
)

appendix_main_correlations <- sample_h4 |>
  select(
    tele_score_norm,
    jury_score_norm,
    jury_televote_gap,
    special_effects_count,
    performance_order
  ) |>
  cor() |>
  round(3) |>
  as.data.frame() |>
  rownames_to_column(
    var = "variable"
  )

write_csv(
  appendix_main_correlations,
  here("tables", "appendix_main_correlations.csv")
)

h2_vif_export <- car::vif(
  lm(
    jury_score_norm ~
      tempo +
      energy +
      danceability +
      valence +
      factor(year),
    data = sample_h2
  )
) |>
  as.data.frame() |>
  rownames_to_column(
    var = "variable"
  )

write_csv(
  h2_vif_export,
  here("tables", "appendix_h2_vif.csv")
)

14.3 Appendix: Influence Diagnostics

influence_summary_export <- tibble(
  hypothesis = c(
    "H1",
    "H2",
    "H3",
    "H4"
  ),
  observations = c(
    nrow(sample_h1),
    nrow(sample_h2),
    nrow(sample_h3),
    nrow(sample_h4)
  ),
  cooks_threshold = c(
    4 / nrow(sample_h1),
    4 / nrow(sample_h2),
    4 / nrow(sample_h3),
    4 / nrow(sample_h4)
  ),
  influential_observations = c(
    sum(cooks.distance(diag_h1) > 4 / nrow(sample_h1)),
    sum(cooks.distance(diag_h2) > 4 / nrow(sample_h2)),
    sum(cooks.distance(diag_h3) > 4 / nrow(sample_h3)),
    sum(cooks.distance(diag_h4) > 4 / nrow(sample_h4))
  ),
  maximum_cooks_distance = c(
    max(cooks.distance(diag_h1)),
    max(cooks.distance(diag_h2)),
    max(cooks.distance(diag_h3)),
    max(cooks.distance(diag_h4))
  )
)

write_csv(
  influence_summary_export,
  here("tables", "appendix_influence_summary.csv")
)

influential_h1_h3_export <- sample_h3 |>
  mutate(
    cooks_distance = cooks.distance(diag_h3)
  ) |>
  arrange(desc(cooks_distance)) |>
  slice_head(n = 5) |>
  transmute(
    model = "H1/H3",
    year,
    to_country,
    song,
    special_effects_count,
    performance_order,
    outcome = tele_score_norm,
    cooks_distance
  )

influential_h2_export <- sample_h2 |>
  mutate(
    cooks_distance = cooks.distance(diag_h2)
  ) |>
  arrange(desc(cooks_distance)) |>
  slice_head(n = 5) |>
  transmute(
    model = "H2",
    year,
    to_country,
    song,
    special_effects_count = NA_real_,
    performance_order = NA_real_,
    outcome = jury_score_norm,
    cooks_distance
  )

influential_h4_export <- sample_h4 |>
  mutate(
    cooks_distance = cooks.distance(diag_h4)
  ) |>
  arrange(desc(cooks_distance)) |>
  slice_head(n = 5) |>
  transmute(
    model = "H4",
    year,
    to_country,
    song,
    special_effects_count,
    performance_order,
    outcome = jury_televote_gap,
    cooks_distance
  )

influential_observations_export <- bind_rows(
  influential_h1_h3_export,
  influential_h2_export,
  influential_h4_export
)

write_csv(
  influential_observations_export,
  here("tables", "appendix_influential_observations.csv")
)

14.4 Appendix: Robustness Checks

write_csv(
  tidy(model_h1_binary),
  here(
    "tables",
    "appendix_robustness_01_h1_binary_effects.csv"
  )
)

write_csv(
  tidy(model_h4_binary),
  here(
    "tables",
    "appendix_robustness_02_h4_binary_effects.csv"
  )
)

influence_sample_sizes_export <- tibble(
  hypothesis = c(
    "H3",
    "H4"
  ),
  primary_n = c(
    nrow(sample_h3),
    nrow(sample_h4)
  ),
  robustness_n = c(
    nrow(sample_h3_influence),
    nrow(sample_h4_influence)
  ),
  excluded = c(
    nrow(sample_h3) - nrow(sample_h3_influence),
    nrow(sample_h4) - nrow(sample_h4_influence)
  )
)

write_csv(
  influence_sample_sizes_export,
  here(
    "tables",
    "appendix_robustness_03_influence_sample_sizes.csv"
  )
)

write_csv(
  tidy(model_h3_influence),
  here(
    "tables",
    "appendix_robustness_04_h3_influence_model.csv"
  )
)

write_csv(
  tidy(model_h4_influence),
  here(
    "tables",
    "appendix_robustness_05_h4_influence_model.csv"
  )
)

write_csv(
  tidy(model_h2_nonlinear),
  here(
    "tables",
    "appendix_robustness_06_h2_nonlinearity.csv"
  )
)

14.5 Appendix: Random Forest

appendix_robustness_01_h1 <- tidy(
  model_h1_binary
)

write_csv(
  appendix_robustness_01_h1,
  here(
    "tables",
    "appendix_robustness_01_h1_binary_effects.csv"
  )
)

appendix_robustness_02_h4 <- tidy(
  model_h4_binary
)

write_csv(
  appendix_robustness_02_h4,
  here(
    "tables",
    "appendix_robustness_02_h4_binary_effects.csv"
  )
)

appendix_robustness_03_samples <- tibble(
  hypothesis = c(
    "H3",
    "H4"
  ),
  primary_n = c(
    nrow(sample_h3),
    nrow(sample_h4)
  ),
  robustness_n = c(
    nrow(sample_h3_influence),
    nrow(sample_h4_influence)
  ),
  excluded = c(
    nrow(sample_h3) - nrow(sample_h3_influence),
    nrow(sample_h4) - nrow(sample_h4_influence)
  )
)

write_csv(
  appendix_robustness_03_samples,
  here(
    "tables",
    "appendix_robustness_03_influence_sample_sizes.csv"
  )
)

appendix_robustness_04_h3 <- tidy(
  model_h3_influence
)

write_csv(
  appendix_robustness_04_h3,
  here(
    "tables",
    "appendix_robustness_04_h3_influence_model.csv"
  )
)

appendix_robustness_05_h4 <- tidy(
  model_h4_influence
)

write_csv(
  appendix_robustness_05_h4,
  here(
    "tables",
    "appendix_robustness_05_h4_influence_model.csv"
  )
)

appendix_robustness_06_h2 <- tidy(
  model_h2_nonlinear
)

write_csv(
  appendix_robustness_06_h2,
  here(
    "tables",
    "appendix_robustness_06_h2_nonlinearity.csv"
  )
)

14.6 Random Forest Appendix Tables

write_csv(
  rf_fold_summary,
  here(
    "tables",
    "appendix_rf_01_fold_summary.csv"
  )
)

write_csv(
  rf_grid,
  here(
    "tables",
    "appendix_rf_02_tuning_grid.csv"
  )
)

write_csv(
  tele_rf_best_params,
  here(
    "tables",
    "appendix_rf_03_televote_best_parameters.csv"
  )
)

write_csv(
  jury_rf_best_params,
  here(
    "tables",
    "appendix_rf_04_jury_best_parameters.csv"
  )
)

write_csv(
  tele_importance,
  here(
    "tables",
    "appendix_rf_05_televote_variable_importance.csv"
  )
)

write_csv(
  jury_importance,
  here(
    "tables",
    "appendix_rf_06_jury_variable_importance.csv"
  )
)

14.7 Primary Regression Table

primary_regression_export <- etable(
  "H1: Televote" = model_h1,
  "H2: Jury" = model_h2,
  "H3: Televote" = model_h3,
  "H4: Jury-Televote Gap" = model_h4,
  
  dict = c(
    special_effects_count = "Special-effects count",
    performance_order = "Performance order",
    tempo = "Tempo",
    energy = "Energy",
    danceability = "Danceability",
    valence = "Valence"
  ),
  
  fitstat = ~ n + wr2,
  
  digits = 3,
  
  signif.code = c(
    "***" = 0.001,
    "**" = 0.01,
    "*" = 0.05,
    "." = 0.10
  ),
  
  tex = FALSE
)

primary_regression_export_df <- data.frame(
  row = rownames(primary_regression_export),
  primary_regression_export,
  row.names = NULL,
  check.names = FALSE
)

write_csv(
  primary_regression_export_df,
  here(
    "tables",
    "table_primary_regressions.csv"
  )
)

14.8 Final CSV Verification

expected_table_files <- c(
  # Main-text tables
  "table_descriptive_statistics.csv",
  "table_special_effects_frequencies.csv",
  "table_primary_regressions.csv",
  "table_rf_performance.csv",
  "table_hypothesis_summary.csv",
  
  # Appendix correlations and diagnostics
  "appendix_h2_correlations.csv",
  "appendix_main_correlations.csv",
  "appendix_h2_vif.csv",
  "appendix_influence_summary.csv",
  "appendix_influential_observations.csv",
  
  # Appendix robustness
  "appendix_robustness_01_h1_binary_effects.csv",
  "appendix_robustness_02_h4_binary_effects.csv",
  "appendix_robustness_03_influence_sample_sizes.csv",
  "appendix_robustness_04_h3_influence_model.csv",
  "appendix_robustness_05_h4_influence_model.csv",
  "appendix_robustness_06_h2_nonlinearity.csv",
  
  # Appendix random forest
  "appendix_rf_01_fold_summary.csv",
  "appendix_rf_02_tuning_grid.csv",
  "appendix_rf_03_televote_best_parameters.csv",
  "appendix_rf_04_jury_best_parameters.csv",
  "appendix_rf_05_televote_variable_importance.csv",
  "appendix_rf_06_jury_variable_importance.csv"
)

table_export_check <- tibble(
  file = expected_table_files,
  exists = file.exists(
    here(
      "tables",
      expected_table_files
    )
  )
)

table_export_check

## # A tibble: 22 × 2
##    file                                  exists
##    <chr>                                 <lgl> 
##  1 table_descriptive_statistics.csv      TRUE  
##  2 table_special_effects_frequencies.csv TRUE  
##  3 table_primary_regressions.csv         TRUE  
##  4 table_rf_performance.csv              TRUE  
##  5 table_hypothesis_summary.csv          TRUE  
##  6 appendix_h2_correlations.csv          TRUE  
##  7 appendix_main_correlations.csv        TRUE  
##  8 appendix_h2_vif.csv                   TRUE  
##  9 appendix_influence_summary.csv        TRUE  
## 10 appendix_influential_observations.csv TRUE  
## # ℹ 12 more rows

14.9 Final File List

sort(
  list.files(
    here("tables"),
    pattern = "\\.csv$"
  )
)

##  [1] "appendix_h2_correlations.csv"                     
##  [2] "appendix_h2_vif.csv"                              
##  [3] "appendix_influence_summary.csv"                   
##  [4] "appendix_influential_observations.csv"            
##  [5] "appendix_main_correlations.csv"                   
##  [6] "appendix_rf_01_fold_summary.csv"                  
##  [7] "appendix_rf_02_tuning_grid.csv"                   
##  [8] "appendix_rf_03_televote_best_parameters.csv"      
##  [9] "appendix_rf_04_jury_best_parameters.csv"          
## [10] "appendix_rf_05_televote_variable_importance.csv"  
## [11] "appendix_rf_06_jury_variable_importance.csv"      
## [12] "appendix_robustness_01_h1_binary_effects.csv"     
## [13] "appendix_robustness_02_h4_binary_effects.csv"     
## [14] "appendix_robustness_03_influence_sample_sizes.csv"
## [15] "appendix_robustness_04_h3_influence_model.csv"    
## [16] "appendix_robustness_05_h4_influence_model.csv"    
## [17] "appendix_robustness_06_h2_nonlinearity.csv"       
## [18] "table_descriptive_statistics.csv"                 
## [19] "table_hypothesis_summary.csv"                     
## [20] "table_primary_regressions.csv"                    
## [21] "table_rf_performance.csv"                         
## [22] "table_special_effects_frequencies.csv"