Test permutacyjny i “Dzień świstaka”

Czy świstak Phil z filmu “Dzień świstaka” rzeczywiście przewiduje pogodę? 2 lutego wychodzi ze swej nory i jeśli zobaczy swój cień (czyli poranek jest słoneczny), zima ma trwać jeszcze 6 tygodni, a jeśli nie zobaczy (pochmurno), szybko nadejdzie wiosna.

Podejdźmy do tego na poważnie.

Dane

Dane o tym, czy świstak widział cień, czy nie, można pobrać stąd (są do 2016 roku). Mamy 12 przypadków, w których informacja, co widział świstak, nie została zarejestrowana, oraz w jednym roku świstak widział swój cień “częściowo” — wszystkie te przypadki pomijam.

Niestety nie jest łatwo określić, co oznacza “szybkie nadejście wiosny”. Będę analizował średnią temperaturę w marcu w Pensylwanii, bo właśnie w tym stanie urzęduje świstak. Porównam, czy w latach, w których świstak widział swój cień, te średnie temperatury są niższe, niż gdy tego cienia nie widział.

Główny cel tego artykułu to pokazanie, jak wykonać prosty test permutacyjny, także udostępniam również kod R. Poniżej wyczytuję dane, pomijam problematyczne przypadki i zamieniam skalę Fahrenheita na Celsjusza. Oczywiście ta modyfikacja nie będzie mieć wpływu na wyniki, kwestia wygody. W pierwszych latach (1886-1894) nie jest znana temperatura, więc te obserwacje też pomijam.

library("tidyverse")
theme_set(ggpubr::theme_pubr(base_size = 13))

phil <- read_csv("phil.csv") %>% 
  rename(
    Shadow = `Punxsutawney Phil`,
    Temperature = `March Average Temperature (Pennsylvania)`
  ) %>% 
  filter(Shadow %in% c("Full Shadow", "No Shadow")) %>% 
  mutate(Shadow = as.factor(Shadow)) %>%
  mutate(Temperature = (Temperature - 32) * 5/9) %>% 
  drop_na(Temperature)

count(phil, Shadow)

# A tibble: 2 × 2
  Shadow          n
  <fct>       <int>
1 Full Shadow   100
2 No Shadow      15

Ostatecznie mamy 115 obserwacji, w 15 przypadkach świstak nie widział cienia, w 100 widział.

Rozkład temperatur

Sprawdźmy, jak rozkładają się średnie temperatury w marcu w zależności od tego, co “twierdzi” świstak, oraz policzmy podstawowe statystyki (średnie ze średnich).

ggplot(phil, aes(Temperature, fill = Shadow)) + 
  geom_density(alpha = 0.6)

summarise(phil, n = n(), M = mean(Temperature), SD = sd(Temperature),
  .by = Shadow)

# A tibble: 2 × 4
  Shadow          n     M    SD
  <fct>       <int> <dbl> <dbl>
1 Full Shadow   100  2.19  2.41
2 No Shadow      15  2.4   1.97

Różnic praktycznie brak (0,2 stopnia), a zatem nawet jeśli świstak coś “wie”, to można by jego wskazania zignorować.

Ale jak pisałem, chcemy do tego podejść poważnie — zbyt poważnie. Załóżmy, że ta różnica 0,2 stopnia to jednak coś. Ponieważ jest “w dobrą stronę” (niższa temperatura, gdy świstak widział cień, więc zima miała być dłuższa), to sprawdźmy, czy jest istotna statystycznie.

Permutacja

Zacznijmy od zapisania zaobserwowanej różnicy do zmiennej real_diff (bo chcę mieć dokładny wynik, a nie przybliżony).

real_diff <- phil %>% 
  summarise(M = mean(Temperature), .by = Shadow) %>% 
  pull(M) %>%
  diff()

Pytanie brzmi, czy jeśli świstak jest losowym prognostą, to różnice na poziomie 0,2 powinny nas dziwić, czy nie? Są zgodne z założeniem, że świstak nic nie wie, czy nie dają się w ten sposób obronić?

Żeby to sprawdzić, możemy np. zapomnieć, w jakich sytuacjach była mierzona temperatura: czy wtedy, gdy świstak widział cień, czy nie. Technicznie da się to rozwiązać tak, że losowo zamieniamy miejscami informacje o cieniu — czyli dokonujemy permutacji. Zauważmy, że zachowany zostanie rozkład zarówno zmiennej Shadow (wciąż w 15 przypadkach świstak nie widzi cienia), jak i Temperature. Jedynie związek (rzekomy) między nimi zostanie utracony: nie wiemy już, której sytuacji odpowiadała która temperatura.

Wykonajmy jedną taką permutację i policzymy różnicę między temperaturami.

set.seed(11)
phil %>% 
  mutate(Shadow = sample(Shadow)) %>% 
  summarise(M = mean(Temperature), .by = Shadow) %>% 
  pull(M) %>%
  diff()

[1] 0.4968519

Uzyskaliśmy różnicę 0,5 stopnia, MIMO ŻE świstak nic nie wie o pogodzie (bo przecież zniszczyliśmy związek między cieniem a temperaturą). Czyli to 0,2 wygląda przy tym naprawdę marnie.

Test permutacyjny

Ale może tak się nieszczęśliwie to losowanie (permutacja) potoczyło? Jest zdecydowanie za wcześnie, by wyciągać wnioski. Wykonajmy 10 tysięcy takich permutacji i za każdym razem policzmy różnicę średnich temperatur. I sprawdźmy, w którym miejscu to nasze skromne 0,2 się znajduje. Dokładnie: policzmy, w jakiej proporcji przypadków otrzymamy różnice równe lub większe od 0,2 — czyli p-wartość.

N <- 1e4
perm_diffs <- numeric(N)
for (i in 1:N) {
  perm_diffs[i] <- phil %>% 
    mutate(Shadow = sample(Shadow)) %>% 
    summarise(M = mean(Temperature), .by = Shadow) %>% 
    pull(M) %>%
    diff()
}
mean(perm_diffs >= real_diff)

[1] 0.3398

Otrzymujemy p = 0,34, także cóż, może i świstak coś wie, ale na podstawie danych nie da się tego udowodnić.