EU:n lainsäädäntö jakaumana — Osa 4: Bayesilainen muutospiste-analyysi

library(tidyverse)
#library(eurlex)
library(here)
library(lubridate)
library(scales)
library(ggtext)
library(patchwork)
library(bcp)        # Bayesian Change Point analysis

col_red    <- "#e63946"
col_green  <- "#2a9d8f"
col_orange <- "#f4a261"
col_navy   <- "#1d3557"
col_blue   <- "#457b9d"

theme_set(
  theme_minimal(base_size = 14) +
    theme(
      plot.title       = element_markdown(face = "bold", size = 15),
      plot.subtitle    = element_markdown(color = "grey40"),
      plot.caption     = element_text(color = "grey55", size = 9),
      axis.title       = element_text(color = "grey30"),
      panel.grid.minor = element_blank()
    )
)

Miksi muutospisteanalyysi — eikä pelkkä trendiviiva?

Osassa 1 piirrettiin EU:n vuosittainen säädöstuotanto ja lisättiin LOESS-käyrä. Se näytti suuntaa. Mutta LOESS-käyrä ei vastaa oikeaan kysymykseen:

Milloin täsmälleen EU:n lainsäädäntötahti muuttui rakenteellisesti — ja kuinka varma voimme siitä olla?

Trendiviiva on visuaalinen apuväline, ei tilastollinen testi. Se ei kerro epävarmuudesta muutoksen ajankohdasta. Se ei erota rakenteellista murrosta satunnaisesta vuosivaihtelusta.

Bayesilainen muutospiste-analyysi tekee tämän oikein.

bcp-paketti (Barry & Hartigan 1993) laskee MCMC-simulaatiolla posterioritodennäköisyyden sille, että kussakin ajankohdassa tapahtuu rakenteellinen muutos tason tai varianssin suhteen. Tuloksena ei ole yksi pistemäinen “muutospiste” — tuloksena on jakauma mahdollisista muutospisteistä, joka kvantifioi epävarmuuden eksplisiittisesti.

Tämä on juuri se mitä “maailma on jakauma” tarkoittaa sovellettuna aikasarja-analyysiin.

---

Data: osa 1:n säädöstuotanto

data_path <- here("data/eu/eu_saadanto_raw.rds")

if (!file.exists(data_path)) {
  stop(
    "Tiedostoa ", data_path, " ei löydy.\n",
    "Aja ensin osa 1 lokaalisti tallentaaksesi datan."
  )
}

raw <- readRDS(data_path)
message("Data ladattu: ", nrow(raw), " riviä")

# Siivotaan ja lasketaan vuosittaiset tuotantoluvut
df <- raw |>
  filter(!is.na(date)) |>
  mutate(
    date         = as.Date(date),
    vuosi        = year(date),
    saadostyyppi = case_when(
      resource_type == "regulation"     ~ "Asetus",
      resource_type == "directive"      ~ "Direktiivi",
      resource_type == "decision"       ~ "Päätös",
      resource_type == "recommendation" ~ "Suositus"
    )
  ) |>
  filter(vuosi >= 1960, vuosi <= 2023)

# Vuosittainen tuotanto tyypeittäin
vuosi_df <- df |>
  count(vuosi, saadostyyppi) |>
  # Varmistetaan täydellinen vuosisarja — puuttuvat vuodet saavat nollan
  complete(vuosi = 1960:2023, saadostyyppi, fill = list(n = 0))

# Kokonaistuotanto (kaikki tyypit yhteensä)
vuosi_yht <- vuosi_df |>
  group_by(vuosi) |>
  summarise(n_yht = sum(n), .groups = "drop")

cat("Vuosia analyysissa:", nrow(vuosi_yht), "\n")

Vuosia analyysissa: 64 

cat("Säädöksiä yhteensä:", sum(vuosi_yht$n_yht), "\n")

Säädöksiä yhteensä: 208164 

---

Vaihe 1: Visuaalinen esitutkimus

Ennen muutospisteanalyysiä katsotaan aikasarjaa; ei LOESS:n kanssa vaan ilman, jotta silmä ei totu odottamaan tiettyä muotoa.

p1 <- ggplot(vuosi_yht, aes(vuosi, n_yht)) +
  geom_col(fill = col_blue, alpha = 0.75, width = 0.8) +
  scale_x_continuous(breaks = seq(1960, 2023, 5)) +
  scale_y_continuous(labels = comma_format(big.mark = " ")) +
  labs(
    title    = "**EU:n vuosittainen säädöstuotanto** 1960–2023",
    subtitle = "Kaikki neljä tyyppiä yhteensä — ilman trendikäyrää",
    x = NULL, y = "Säädösten lukumäärä",
    caption  = "Lähde: EUR-Lex SPARQL via eurlex (R). Kristian Vepsäläinen / kristianvepsalainen.com"
  )

p2 <- vuosi_df |>
  filter(saadostyyppi %in% c("Asetus", "Direktiivi", "Päätös")) |>
  ggplot(aes(vuosi, n, color = saadostyyppi)) +
  geom_line(linewidth = 0.9, alpha = 0.85) +
  scale_color_manual(values = c(col_red, col_green, col_blue)) +
  scale_x_continuous(breaks = seq(1960, 2023, 5)) +
  scale_y_continuous(labels = comma_format(big.mark = " ")) +
  labs(
    title    = "Tyypeittäin: **muutokset eivät ole synkronisia**",
    subtitle = "Asetusten, direktiivien ja päätösten tuotanto on eri prosesseja",
    x = NULL, y = "Säädösten lukumäärä", color = NULL,
    caption  = ""
  ) +
  theme(legend.position = "top")

p1 / p2

EU:n vuosittainen säädöstuotanto 1960–2023 ilman trendikäyrää. Mitä rakennetta näet?

---

Vaihe 2: bcp-malli kokonaistuotannolle

Analyysiin käytetään R:n funktiota bcp(), joka saa syötteenä numerovektorin (vuosittaiset lukumäärät järjestyksessä) ja palauttaa jokaiselle ajankohdalle posterioritodennäköisyyden muutokselle. Käytämme oletusprioria p0 = 0.2 (20 % todennäköisyys muutokselle ennen dataa) — tämä on kohtalaisen heikko priori joka antaa datan puhua.

# Valmistetaan syötedata: järjestetty aikasarja
n_vec <- vuosi_yht |>
  arrange(vuosi) |>
  pull(n_yht)

vuodet <- vuosi_yht |>
  arrange(vuosi) |>
  pull(vuosi)

# bcp-malli
set.seed(42)
bcp_fit <- bcp(
  y    = as.numeric(n_vec),
  p0   = 0.2,      # priori muutospisteen todennäköisyydelle
  burnin = 500,
  mcmc   = 5000
)

# Tarkistus
cat("bcp-malli ajettu.\n")

bcp-malli ajettu.

cat("Korkein posterioritodennäköisyys muutokselle:\n")

Korkein posterioritodennäköisyys muutokselle:

max_prob <- max(bcp_fit$posterior.prob, na.rm = TRUE)
max_vuosi <- vuodet[which.max(bcp_fit$posterior.prob)]
cat("  Vuosi:", max_vuosi, "— posteriori p:", round(max_prob, 3), "\n\n")

  Vuosi: 1972 — posteriori p: 0.998 

# Kootaan tulokset data frameksi
bcp_df <- tibble(
  vuosi         = vuodet,
  n_yht         = n_vec,
  post_prob     = bcp_fit$posterior.prob,
  post_mean     = bcp_fit$posterior.mean[, 1],
  post_var      = bcp_fit$posterior.var[, 1]
) |>
  mutate(
    # Muutospisteen kynnysarvo: p > 0.5 = vahva evidenssi
    muutospiste_vahva    = post_prob > 0.5,
    muutospiste_kohtalainen = post_prob > 0.2
  )

cat("Vahvat muutospisteet (p > 0.5):\n")

Vahvat muutospisteet (p > 0.5):

bcp_df |> filter(muutospiste_vahva) |> select(vuosi, post_prob) |> print()

# A tibble: 12 × 2
   vuosi post_prob
   <dbl>     <dbl>
 1  1966     0.972
 2  1967     0.989
 3  1968     0.857
 4  1972     0.998
 5  1979     0.998
 6  1980     0.817
 7  1981     0.961
 8  1985     0.861
 9  1986     0.958
10  1994     0.958
11  2001     0.522
12  2006     0.968

cat("\nKohtalaisen vahvat muutospisteet (p > 0.2):\n")

Kohtalaisen vahvat muutospisteet (p > 0.2):

bcp_df |> filter(muutospiste_kohtalainen) |> select(vuosi, post_prob) |> print()

# A tibble: 13 × 2
   vuosi post_prob
   <dbl>     <dbl>
 1  1961     0.320
 2  1966     0.972
 3  1967     0.989
 4  1968     0.857
 5  1972     0.998
 6  1979     0.998
 7  1980     0.817
 8  1981     0.961
 9  1985     0.861
10  1986     0.958
11  1994     0.958
12  2001     0.522
13  2006     0.968

# Institutionaaliset tapahtumat kontekstia varten
tapahtumat <- tribble(
  ~vuosi, ~teksti,
  1973,   "UK/DK/IE\nliittyminen",
  1986,   "Yhtenäis-\neuroppa-asiak.",
  1993,   "Sisämarkkinat\n(Maastricht)",
  1995,   "AT/FI/SE\nliittyminen",
  2004,   "Itä-Eurooppa\nliittyminen",
  2009,   "Lissabon"
)

p_ts <- ggplot(bcp_df, aes(vuosi, n_yht)) +
  # Taustamerkinnät institutionaalisille tapahtumille
  geom_vline(data = tapahtumat, aes(xintercept = vuosi),
             linetype = "dotted", color = "grey70", linewidth = 0.4) +
  geom_col(fill = col_blue, alpha = 0.5, width = 0.8) +
  # Posteriori-keskiarvo = mallin estimoima "taso" kullakin hetkellä
  geom_line(aes(y = post_mean), color = col_red,
            linewidth = 1.2) +
  geom_text(data = tapahtumat,
            aes(x = vuosi, y = Inf, label = teksti),
            vjust = 1.3, size = 2.4, color = "grey50", linewidth = 0.3) +
  scale_x_continuous(breaks = seq(1960, 2023, 5)) +
  scale_y_continuous(labels = comma_format(big.mark = " ")) +
  labs(
    title    = "EU:n säädöstuotanto ja **posteriori-keskiarvo**",
    subtitle = "Punainen viiva = bcp-mallin estimoima taso. Hypyt kertovat rakenteellisista muutoksista.",
    x = NULL, y = "Säädösten lukumäärä"
  )

p_prob <- ggplot(bcp_df, aes(vuosi, post_prob)) +
  geom_hline(yintercept = c(0.2, 0.5), linetype = "dashed",
             color = c(col_orange, col_red), linewidth = 0.6) +
  geom_col(aes(fill = post_prob > 0.5), width = 0.8, show.legend = FALSE) +
  scale_fill_manual(values = c(col_blue, col_red)) +
  annotate("text", x = 1962, y = 0.52, label = "Vahva evidenssi (p > 0.5)",
           color = col_red, size = 3, hjust = 0) +
  annotate("text", x = 1962, y = 0.22, label = "Kohtalainen evidenssi (p > 0.2)",
           color = col_orange, size = 3, hjust = 0) +
  scale_x_continuous(breaks = seq(1960, 2023, 5)) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  labs(
    title    = "**Posterioritodennäköisyys muutospisteelle** vuosittain",
    subtitle = "Tämä on se, mitä trendiviiva ei näytä: epävarmuus muutoksen ajankohdasta",
    x = NULL, y = "P(muutospiste)",
    caption  = "Malli: bcp (Barry & Hartigan 1993), p0 = 0.2, MCMC = 5000. Lähde: EUR-Lex via eurlex (R).\nKristian Vepsäläinen / kristianvepsalainen.com"
  )

p_ts / p_prob + plot_layout(heights = c(2, 1))

Bayesilainen muutospiste-analyysi EU:n kokonaistuotannolle. Yläkuva: aikasarja posteriori-keskiarvolla. Alakuva: muutospisteen posterioritodennäköisyys per vuosi.

---

Vaihe 3: Tyyppikohtaiset muutospisteen analyysit

Kokonaistuotanto peittää alleen eri prosessit. Asetuksilla, direktiiveillä ja päätöksillä on todennäköisesti omat rakenteelliset murroskohtansa — eikä niiden tarvitse osua samaan vuoteen.

# Ajetaan bcp erikseen kullekin säädöstyypille
tyypit_analyysissa <- c("Asetus", "Direktiivi", "Päätös")

aja_bcp <- function(tyyppi) {
  sarja <- vuosi_df |>
    filter(saadostyyppi == tyyppi) |>
    arrange(vuosi) |>
    pull(n)

  set.seed(42)
  fit <- bcp(
    y      = as.numeric(sarja),
    p0     = 0.2,
    burnin = 500,
    mcmc   = 5000
  )

  tibble(
    vuosi        = 1960:2023,
    n            = sarja,
    post_prob    = fit$posterior.prob,
    post_mean    = fit$posterior.mean[, 1],
    saadostyyppi = tyyppi
  )
}

bcp_tyypit <- map_dfr(tyypit_analyysissa, aja_bcp)

# Yhteenveto tyypeittäin
cat("Vahvimmat muutospisteet tyypeittäin:\n")

Vahvimmat muutospisteet tyypeittäin:

bcp_tyypit |>
  group_by(saadostyyppi) |>
  slice_max(post_prob, n = 3) |>
  select(saadostyyppi, vuosi, post_prob) |>
  mutate(post_prob = round(post_prob, 3)) |>
  print()

# A tibble: 9 × 3
# Groups:   saadostyyppi [3]
  saadostyyppi vuosi post_prob
  <chr>        <int>     <dbl>
1 Asetus        1967     0.996
2 Asetus        1979     0.993
3 Asetus        1981     0.992
4 Direktiivi    2009     0.603
5 Direktiivi    2014     0.516
6 Direktiivi    2011     0.368
7 Päätös        1982     0.747
8 Päätös        1974     0.67 
9 Päätös        1962     0.576

ggplot(bcp_tyypit, aes(vuosi, post_prob, fill = post_prob > 0.5)) +
  geom_col(width = 0.8, show.legend = FALSE) +
  geom_hline(yintercept = 0.5, linetype = "dashed",
             color = col_red, linewidth = 0.5) +
  scale_fill_manual(values = c(col_blue, col_red)) +
  scale_x_continuous(breaks = seq(1960, 2023, 10)) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  facet_wrap(~saadostyyppi, ncol = 1) +
  labs(
    title    = "**Muutospisteen posterioritodennäköisyys tyypeittäin**",
    subtitle = "Osuvatko murrokset samaan aikaan? Vai onko kullakin tyypillä oma dynamiikkansa?",
    x = NULL, y = "P(muutospiste)",
    caption  = "Malli: bcp (Barry & Hartigan 1993), p0 = 0.2, MCMC = 5000. Lähde: EUR-Lex via eurlex (R).\nKristian Vepsäläinen / kristianvepsalainen.com"
  )

Muutospisteen posterioritodennäköisyys tyypeittäin. Eri tyypit reagoivat eri institutionaalisiin tapahtumiin.

# Poimitaan vahvat muutospisteet pystyviivoja varten
vahvat_mp <- bcp_tyypit |>
  filter(post_prob > 0.5) |>
  select(saadostyyppi, vuosi, post_prob)

ggplot(bcp_tyypit, aes(vuosi, n)) +
  geom_col(fill = col_blue, alpha = 0.5, width = 0.8) +
  geom_line(aes(y = post_mean), color = col_red, linewidth = 1.1) +
  geom_vline(data = vahvat_mp, aes(xintercept = vuosi),
             color = col_red, linetype = "dashed",
             linewidth = 0.6, alpha = 0.7) +
  scale_x_continuous(breaks = seq(1960, 2023, 10)) +
  scale_y_continuous(labels = comma_format(big.mark = " ")) +
  facet_wrap(~saadostyyppi, ncol = 1, scales = "free_y") +
  labs(
    title    = "**Aikasarjat ja posteriori-keskiarvot** tyypeittäin",
    subtitle = "Punainen viiva = mallin estimoima taso. Pisteviivat = vahvat muutospisteet.",
    x = NULL, y = "Säädösten lukumäärä",
    caption  = "Malli: bcp (Barry & Hartigan 1993), p0 = 0.2, MCMC = 5000. Lähde: EUR-Lex via eurlex (R).\nKristian Vepsäläinen / kristianvepsalainen.com"
  )

Aikasarja ja posteriori-keskiarvo tyypeittäin. Punaiset pystyviivat = vahvat muutospisteet (p > 0.5).

---

Vaihe 4: Prioriherkkyyden tarkastelu

Hyvä bayesilainen analyysi raportoi, kuinka herkästi tulokset muuttuvat priorin suhteen. Ajetaan sama malli kolmella eri p0-arvolla kokonaistuotannolle:

• p0 = 0.1 — konservatiivinen priori: muutospisteitä on harvoin
• p0 = 0.2 — perusmalli
• p0 = 0.4 — liberaali priori: muutospisteitä voi olla usein

priorit <- c(0.1, 0.2, 0.4)
priorit_nimet <- c("Konservatiivinen (p0=0.1)",
                   "Perusmalli (p0=0.2)",
                   "Liberaali (p0=0.4)")

aja_bcp_priori <- function(p0_arvo, nimi) {
  set.seed(42)
  fit <- bcp(
    y      = as.numeric(n_vec),
    p0     = p0_arvo,
    burnin = 500,
    mcmc   = 5000
  )
  tibble(
    vuosi     = vuodet,
    post_prob = fit$posterior.prob,
    post_mean = fit$posterior.mean[, 1],
    priori    = nimi
  )
}

herkkyys_df <- map2_dfr(priorit, priorit_nimet, aja_bcp_priori) |>
  mutate(priori = factor(priori, levels = priorit_nimet))

ggplot(herkkyys_df, aes(vuosi, post_prob, color = priori)) +
  geom_line(linewidth = 0.9, alpha = 0.85) +
  geom_hline(yintercept = 0.5, linetype = "dashed",
             color = "grey50", linewidth = 0.5) +
  scale_color_manual(values = c(col_green, col_navy, col_red)) +
  scale_x_continuous(breaks = seq(1960, 2023, 5)) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  labs(
    title    = "**Prioriherkkyyden tarkastelu** — robustit muutospisteet pysyvät",
    subtitle = "Muutospisteet, jotka näkyvät kaikilla kolmella priorilla, ovat datavahvistettuja",
    x = NULL, y = "P(muutospiste)", color = NULL,
    caption  = "Malli: bcp (Barry & Hartigan 1993), MCMC = 5000. Lähde: EUR-Lex via eurlex (R).\nKristian Vepsäläinen / kristianvepsalainen.com"
  ) +
  theme(legend.position = "top")