298. Symulacja Monte Carlo jako sposób rozumienia rozkładów

W poprzednich lekcjach analizowaliśmy konkretne rozkłady prawdopodobieństwa: normalny, t-Studenta, chi-kwadrat, F, Poissona czy wykładniczy. Teraz wykonamy krok w stronę myślenia symulacyjnego. Zamiast tylko korzystać z gotowych wzorów, zobaczymy, jak rozkład może powstać w wyniku wielokrotnego losowania i powtarzania tego samego mechanizmu.

Symulacja Monte Carlo jest sposobem badania zjawisk za pomocą dużej liczby losowych, a dokładniej pseudolosowych, realizacji modelu. Jej główna idea jest prosta: jeśli umiemy opisać mechanizm losowy, możemy wielokrotnie go odtworzyć na komputerze, a następnie obserwować częstości, średnie, kwantyle i inne charakterystyki wyników. Dzięki temu rozkłady przestają być wyłącznie abstrakcyjnymi krzywymi, a stają się obrazem wielu możliwych powtórzeń tego samego procesu.

W tej lekcji zobaczymy, czym różni się pojedyncze losowanie od pełnego badania symulacyjnego, jak z symulacji uzyskuje się przybliżenie prawdopodobieństwa oraz dlaczego wynik symulacji również ma własną niepewność. Przykłady będą osadzone w obszarze medycyny, zdrowia publicznego i metodologii badań. Temat można przećwiczyć także w ramach szkolenia na statystyk.edu.pl.

Na czym polega metoda Monte Carlo?

Załóżmy, że interesuje nas zmienna losowa Y. Może ona oznaczać na przykład:

• liczbę pacjentów, u których wystąpi działanie niepożądane;
• czas oczekiwania na przyjęcie do szpitala;
• różnicę średnich wyników między dwiema grupami;
• wartość kosztu leczenia w modelu ekonomicznym;
• odsetek osób, które uzyskają wynik dodatni w teście diagnostycznym.

Jeżeli znamy rozkład zmiennej Y albo potrafimy wygenerować jej wartości zgodnie z określonym modelem, możemy powtórzyć losowanie wiele razy. Otrzymujemy wtedy ciąg:

Y1, Y2, Y3, ..., YB

Symbol B oznacza liczbę powtórzeń symulacji. Dla każdego powtórzenia zapisujemy wynik, a następnie analizujemy cały zbiór wartości. Histogram wyników przybliża rozkład zmiennej Y, średnia z wyników przybliża wartość oczekiwaną, a odsetek wyników spełniających określony warunek przybliża prawdopodobieństwo.

W klasycznym ujęciu metoda Monte Carlo została opisana jako statystyczne podejście do rozwiązywania problemów, w których bezpośrednie obliczenia analityczne są trudne lub niepraktyczne (Metropolis i Ulam, 1949). Współcześnie metoda ta jest używana zarówno do obliczeń numerycznych, jak i do oceny właściwości metod statystycznych.

Od prawdopodobieństwa do częstości

Najważniejsza intuicja wynika z prawa wielkich liczb. Jeżeli wielokrotnie powtarzamy doświadczenie losowe w tych samych warunkach, względna częstość badanego zdarzenia zwykle stabilizuje się w pobliżu jego prawdopodobieństwa.

Jeżeli w każdej symulacji występuje zdarzenie A, możemy zdefiniować zmienną wskaźnikową:

Ii = 1, gdy w i-tej symulacji wystąpiło A; Ii = 0, gdy A nie wystąpiło.

Przybliżenie prawdopodobieństwa zdarzenia ma wtedy postać:

p̂ = (1/B) × Σ Ii

gdzie:

• p̂ — symulacyjne przybliżenie prawdopodobieństwa zdarzenia;
• B — liczba powtórzeń symulacji;
• I_i — wskaźnik wystąpienia zdarzenia w i-tym powtórzeniu;
• Σ — suma po wszystkich powtórzeniach.

Jeżeli w 10 000 symulacji zdarzenie wystąpiło 2370 razy, otrzymujemy p̂ = 2370/10000 = 0,237. Nie znaczy to, że prawdziwe prawdopodobieństwo wynosi dokładnie 0,237. Jest to wynik uzyskany w skończonej liczbie losowań i dlatego podlega błędowi Monte Carlo.

Błąd Monte Carlo i liczba powtórzeń

Wynik symulacji nie jest całkowicie stabilny. Gdybyśmy uruchomili ten sam program z innym ziarnem generatora liczb pseudolosowych, otrzymalibyśmy nieco inne wartości. Ta zmienność nie musi świadczyć o błędzie programu. Jest naturalną konsekwencją skończonej liczby powtórzeń.

Dla przybliżenia częstości zdarzenia o prawdopodobieństwie p przybliżony błąd standardowy można zapisać jako:

SEMC(p̂) ≈ √[p̂(1 − p̂)/B]

Wzór oznacza, że:

• im większe B, tym mniejszy błąd Monte Carlo;
• największa niepewność występuje zwykle w pobliżu p = 0,5;
• zwiększenie liczby powtórzeń czterokrotnie zmniejsza błąd około dwukrotnie, a nie czterokrotnie.

To ostatnie spostrzeżenie jest praktycznie ważne. Przejście z 1000 do 4000 powtórzeń może przynieść podobną redukcję błędu jak przejście z 4000 do 16 000, ale koszt obliczeń rośnie. Liczbę powtórzeń należy więc dobierać do celu analizy i wymaganej dokładności, a nie kierować się wyłącznie dużą, efektowną liczbą.

W badaniach symulacyjnych należy odróżniać niepewność Monte Carlo od niepewności wynikającej z samego modelu. Można wykonać milion powtórzeń, ale jeśli założono nierealistyczny rozkład wejściowy, bardzo dokładnie oszacujemy niewłaściwy model.

Jak powstaje rozkład w symulacji?

Wyobraźmy sobie, że chcemy zrozumieć rozkład liczby powikłań w grupie 20 pacjentów. Zakładamy, że u każdego pacjenta prawdopodobieństwo powikłania wynosi 0,10, a wyniki pacjentów są niezależne. W jednym powtórzeniu symulujemy 20 obserwacji typu „powikłanie: tak/nie”, a następnie liczymy ich sumę.

Jedno powtórzenie może dać wynik 1, inne 3, jeszcze inne 0 albo 5. Gdy wykonamy 10 000 takich powtórzeń, otrzymamy 10 000 możliwych wartości liczby powikłań. Ich histogram powinien przypominać rozkład dwumianowy z parametrami n = 20 i p = 0,10.

W ten sposób można zrozumieć, że rozkład nie jest „dodatkową substancją” istniejącą poza doświadczeniem. Jest matematycznym opisem tego, jak często pojawiają się różne wyniki przy wielokrotnym powtarzaniu określonego mechanizmu losowego.

Przykład 1: symulacja wyniku testu diagnostycznego

Założenia

Rozważmy hipotetyczne badanie przesiewowe. W populacji częstość choroby wynosi 5%. Czułość testu wynosi 90%, a swoistość 95%. Chcemy oszacować, jaki odsetek wyników dodatnich będzie prawdziwie dodatni, czyli jaka będzie dodatnia wartość predykcyjna.

Dla każdego symulowanego pacjenta wykonujemy trzy kroki:

1. Losujemy, czy pacjent rzeczywiście ma chorobę, z prawdopodobieństwem 0,05.
2. Jeżeli choroba występuje, losujemy wynik testu z czułością 0,90.
3. Jeżeli choroby nie ma, losujemy wynik testu ujemny ze swoistością 0,95.

Następnie powtarzamy procedurę dla 100 000 osób. Obliczamy:

PPV = liczba prawdziwie dodatnich / liczba wszystkich dodatnich

Wartość predykcyjna dodatnia będzie w przybliżeniu wynosić:

PPV = (0,05 × 0,90) / [(0,05 × 0,90) + (0,95 × 0,05)] ≈ 0,486

W symulacji wynik może być na przykład równy 0,485 albo 0,488, zależnie od losowania. Najważniejsza jest interpretacja: nawet przy czułości 90% i swoistości 95% dodatni wynik nie musi oznaczać choroby z prawdopodobieństwem 90%. Zależy także od częstości choroby w badanej populacji.

Symulacja pomaga zobaczyć tę zależność na poziomie konkretnych osób. Możemy zmienić częstość choroby na 1%, 10% albo 20% i sprawdzić, jak zmienia się struktura wyników dodatnich. Jest to dydaktycznie użyteczne, ponieważ pozwala odróżnić czułość i swoistość od wartości predykcyjnych.

Przykład 2: niepewność oszacowania odsetka szczepień

Załóżmy, że w badaniu ankietowym 240 spośród 400 osób deklaruje wykonanie określonego szczepienia. Estymowany odsetek wynosi:

p̂ = 240/400 = 0,60

Możemy potraktować 0,60 jako punktową ocenę prawdopodobieństwa, a następnie symulować wielokrotne próby po 400 osób z prawdopodobieństwem sukcesu równym 0,60. W każdym powtórzeniu zapisujemy odsetek osób zaszczepionych. Rozkład 10 000 otrzymanych odsetków pokazuje, jak bardzo wynik może się różnić wyłącznie z powodu losowego doboru próby.

Jeżeli 95% symulowanych odsetków mieści się na przykład między 0,552 a 0,648, zakres ten obrazuje niepewność związaną z losowym próbkowaniem przy przyjętym modelu. Nie jest to automatycznie pełnoprawny przedział ufności dla wszystkich możliwych źródeł błędu. Symulacja nie uwzględnia sama z siebie błędu pomiaru, odmów udziału, błędu pokrycia populacji ani tendencyjnego raportowania.

Ten przykład pokazuje ważną zasadę: symulacja jest tak dobra, jak jej założenia. Jeżeli model obejmuje wyłącznie losowanie proste, wynik dotyczy przede wszystkim losowania prostego. Nie wolno bez dodatkowych argumentów przenosić go na złożony sondaż z warstwowaniem, klastrami i wagami analitycznymi.

Symulacja rozkładu statystyki

Monte Carlo można wykorzystać nie tylko do generowania pojedynczych obserwacji, lecz także do badania zachowania statystyki. Przykładowo chcemy sprawdzić, jak wygląda rozkład różnicy średnich w małych próbach.

W jednym powtórzeniu losujemy 15 wyników z grupy kontrolnej i 15 wyników z grupy interwencji, zgodnie z przyjętym modelem. Obliczamy różnicę średnich. Następnie powtarzamy cały proces 5000 razy. Otrzymujemy rozkład 5000 różnic.

Możemy wtedy ocenić:

• średni wynik estymatora, czyli jego obciążenie względem wartości założonej;
• odchylenie standardowe wyników, czyli zmienność estymatora;
• częstość, z jaką przedział ufności obejmuje prawdziwą wartość;
• częstość odrzucenia hipotezy zerowej w określonych warunkach;
• wpływ odstających obserwacji, skośności lub braków danych.

W ten sposób symulacja staje się eksperymentem metodologicznym. Nie odpowiada bezpośrednio na pytanie, jaki jest efekt w populacji, lecz pozwala sprawdzić, jak zachowuje się wybrana metoda analizy w sytuacjach, które sami zdefiniowaliśmy.

Planowanie dobrej symulacji

W literaturze metodologicznej zaleca się, aby symulację planować podobnie jak eksperyment. Burton, Altman, Royston i Holder podkreślają znaczenie jasno określonego celu, mechanizmu generowania danych, liczby powtórzeń i sposobu oceny wyników (Burton et al., 2006). Morris, White i Crowther proponują uporządkowanie opisu za pomocą schematu ADEMP: Aims, Data-generating mechanisms, Estimands, Methods, Performance measures (Morris et al., 2019).

1. Cel: określ, co ma zostać sprawdzone.
2. Mechanizm generowania danych: opisz rozkłady, zależności, liczebności i braki danych.
3. Estymand: wskaż dokładną wielkość, którą chcesz ocenić.
4. Metody: ustal, jakie procedury będą porównywane.
5. Miary wyników: wybierz np. obciążenie, błąd średniokwadratowy, pokrycie lub moc.

Warto także ustalić ziarno generatora liczb pseudolosowych, zachować kod i opisać liczbę powtórzeń. Ziarno nie czyni wyników „prawdziwszymi”, ale umożliwia ich odtworzenie. Jeśli symulacja jest wykonywana równolegle na kilku procesorach, trzeba dodatkowo zadbać o niezależne strumienie liczb losowych.

Najczęstsze błędy

• Traktowanie symulacji jako dowodu empirycznego. Symulacja bada konsekwencje przyjętego modelu; nie zastępuje obserwacji rzeczywistej populacji.
• Zbyt mała liczba powtórzeń. Przy kilku lub kilkudziesięciu powtórzeniach histogram może być bardzo niestabilny, a wnioski przypadkowe.
• Brak informacji o założeniach. Sam wykres wyników nie wystarcza, jeśli nie wiadomo, z jakich rozkładów i zależności wygenerowano dane.
• Mylenie błędu Monte Carlo z błędem próbkowania. Pierwszy wynika ze skończonej liczby uruchomień symulacji, drugi z losowania jednostek do badania.
• Zmiana wielu elementów bez planu. Jeżeli jednocześnie zmieniamy liczebność, wariancję, poziom braków i korelacje, trudno ustalić, co odpowiada za wynik.
• Raportowanie wyłącznie wartości p. Symulacja powinna opisywać wielkość błędu, trafność przedziałów, obciążenie i znaczenie praktyczne, a nie tylko odsetek wyników z wartością p poniżej ustalonego progu.
• Założenie, że bardzo duża liczba powtórzeń naprawi zły model. Nie naprawi. Zwiększy tylko precyzję obliczeń wykonanych w ramach błędnych założeń.

Przykład praktyczny: symulacja mocy badania

Zespół pielęgniarski planuje badanie porównujące dwie formy edukacji pacjentów z cukrzycą. Głównym wynikiem ma być odsetek osób, które po trzech miesiącach osiągną założony poziom samokontroli. Badacze przewidują wynik 55% w grupie standardowej i 70% w grupie rozszerzonej.

Zamiast polegać wyłącznie na przybliżeniu analitycznym, można przeprowadzić symulację:

1. Ustalić liczebność obu grup, na przykład po 80 osób.
2. Wygenerować wynik binarny dla każdego pacjenta zgodnie z odpowiednim prawdopodobieństwem.
3. Przeprowadzić zaplanowany test lub model regresji.
4. Zarejestrować, czy wynik jest statystycznie istotny przy wcześniej określonym poziomie alfa.
5. Powtórzyć procedurę np. 5000 razy.

Odsetek symulacji, w których wykryto różnicę, jest przybliżeniem mocy w warunkach przyjętych przez badaczy. Można osobno sprawdzić scenariusz z mniejszą różnicą, większą utratą obserwacji albo nierównymi grupami.

Interpretacja musi być ostrożna. Symulowana moc nie jest gwarancją, że rzeczywiste badanie osiągnie określony wynik. Jest oceną działania planowanej procedury, jeśli założenia dotyczące odsetków, liczebności i braków danych okażą się rozsądne.

Co warto zapamiętać?

• Symulacja Monte Carlo polega na wielokrotnym odtwarzaniu losowego mechanizmu i analizowaniu uzyskanych wyników.
• Histogram wyników symulacji może przybliżać rozkład zmiennej losowej lub rozkład statystyki.
• Częstość zdarzenia w wielu powtórzeniach przybliża jego prawdopodobieństwo, ale wynik ma błąd Monte Carlo.
• Zwiększanie liczby powtórzeń zmniejsza niepewność obliczeniową, lecz nie usuwa błędów modelu.
• W medycynie symulacja pomaga rozumieć testy diagnostyczne, planować badania, oceniać moc i badać zachowanie metod statystycznych.
• Nie należy utożsamiać wyniku symulacji z dowodem przyczynowym ani z bezpośrednim wynikiem badania empirycznego.
• Dobra symulacja wymaga jawnych założeń, określonego celu, adekwatnych miar oceny i przejrzystego raportowania.

Najważniejsza intuicja jest następująca: rozkład można oglądać jako zapis wielu możliwych powtórzeń procesu. Monte Carlo pozwala te powtórzenia wygenerować, policzyć i zobaczyć. Dzięki temu abstrakcyjne pojęcia prawdopodobieństwa, zmienności i niepewności stają się możliwe do prześledzenia krok po kroku.

Źródła i dalsza lektura

1. Metropolis, N.; Ulam, S. (1949). The Monte Carlo Method. Journal of the American Statistical Association, 44(247), 335–341. DOI: 10.1080/01621459.1949.10483310.
2. Concato, J.; Feinstein, A. R. (1997). Monte Carlo methods in clinical research: applications in multivariable analysis. Journal of Investigative Medicine, 45(6), 394–400. Pełny opis źródła.
3. Burton, A.; Altman, D. G.; Royston, P.; Holder, R. L. (2006). The design of simulation studies in medical statistics. Statistics in Medicine, 25(24), 4279–4292. DOI: 10.1002/sim.2673.
4. Morris, T. P.; White, I. R.; Crowther, M. J. (2019). Using simulation studies to evaluate statistical methods. Statistics in Medicine, 38, 2074–2102. DOI: 10.1002/sim.8086.
5. Spiegelhalter, D. J.; Abrams, K. R.; Myles, J. P. (2004). Bayesian Approaches to Clinical Trials and Health-Care Evaluation. 1st edition. John Wiley & Sons. Informacja wydawcy / źródło.
6. Spiegelhalter, D. J.; Myles, J. P.; Jones, D. R.; Abrams, K. R. (1999). An introduction to Bayesian methods in health technology assessment. BMJ, 319(7208), 508–512. DOI: 10.1136/bmj.319.7208.508.