Tutorial Statystyki Bayesowskiej

Część IV: Wprowadzenie do modeli wielopoziomowych

Post będzie aktualizowany.

Dzisiaj zajmiemy się bayesowskimi wielopopoziomowymi (mieszanymi/hierarchicznymi) modelami liniowymi. Na ten temat poświęca się całe książki, (choć z drugiej strony czasami mam wrażenie, że traktuje się ten temat jak czarną magię). Zobaczmy o co całe zmieszanie.

Motywacja

Zamiast zaczynać od wprowadzenia definicji mieszanego modelu liniowego, przyjrzyjmy się poniższym sytuacjom.

Niewykrycie prawdziego efektu

Załóżmy, że wykonujemy przez pewien czas pomiary zadowolenia z obsługi i wysokości napiwków u stałych klientów pewnej restauracji.

Widzimy tu brak zależności liniowej. Teraz zróbmy tak by na wykresie kolor kropki odpowiadał konkretnemu klientowi.

Widzimy, że u każdego z klientów występuje związek liniowy pomiędzy napiwkiem i zadowoleniem. Jednak ponieważ, nasze pomiary zagnieżdone są w klientach, co objawia się tym, że każdy klient średnio daje inną wartość napiwku, przeoczylibyśmy ten związek, nie uwzględniając tego w modelu.

Wykrycie nieprawdziwego efektu

Czy może nam się zdarzyć sytuacja odwrotna? Używając prostej analizy liniowej wykryjemy efekt, którego nie ma? Załóżmy, że chcemy ocenić wpływ picia napojów kofeinowych na wyniki ucznia w nauce. Zebraliśmy dane z trzech klas liceum pytając o średnie spożycie napojów kofeinowych dziennie i średnią ocen.

Zaznaczmy na wykresie kolorami uczniów przynależących do danej klasy.

Jak widzimy, w poszególnych klasach spożycie kofeiny nie ma wpływu na wyniki ucznia w nauce. Kiedy możemy zaobserwować taką zależność? Załóżmy na przykład, że każda z klas ma innego wychowawcę i ich wpływ przedkłada się na średnią ocen. Załóżmy też, że zupełnie przypadkowo w klasach z lepszym wychowawcą, średnie spożycie napojów kofeinowych jest większe. Otrzymamy wtedy sytuację jak na powyższych wykresach.

Oczywiście, może istnieć jakaś zależność pomiędzy średnim spożyciem kofeiny na klasę a jej średnią ocen lecz to nie odpowiada na pytanie, które zadaliśmy tzn. o wpływ kofeiny na osiągnięcia ucznia. Tutaj warto wprowadzić sobie rozróżnienie między efektem węwnątrzobiektowym i międzyobiektowym. Zrobiłem to w jednym z poprzednich wpisów, można przeczytać tutaj.

Z czym się mierzymy?

Problem możemy zdefiniować na kika sposobów. W pierwszej sytuacji chcemy zbadać związek pomiędzy napwkami a zadowoleniem z obsługi. Na przeszkodzie jednak staje nam efekt kienta - każdy klient ma inny bazowy średni napiwek. Klient jest tu zmienną zakłócającą. By poprawnie wyestymować wpływ zadowolenia na napiwki, moglibyśmy kontrolować wpływ poszczególnego klienta. W istocie, na tym z grubsza polega mechanizm działania mieszanych modeli liniowych.

Przyjrzyjmy się teraz problemowi z innej perspektywy. Gdy chcemy zbadać wpływ zmiennej \(X\) na Y, na przykład licząc współczynnik korelacji Pearsona, musimy wykonać następujące kroki (książkowo): Losujemy obserwacje do naszej próby. Następnie, ponieważ każda obserwacja miała takie samo prawdopobieństwo dostania się do próby, traktujemy te obserwacje jako niezależne zmienne losowe o identycznym rozkładzie (iid, independent and identically distributed random variables). Na chłopski rozum, oznacza to, że jeśli stworzymy sobie histogram zmiennej X, otrzymamy z grubsza estymatę rozkładu zmiennej \(X\) w populacji. W przypadku powyższych danych, nie jest to prawda. Przyjrzyjmy sie dokładniej tej definicji.

Definicja składa się z dwóch zasadniczych elementów. Niezależność oznacza, że obserwacje są niepowiązane (wylosowanie \(x_1\) nie mówi nam nic o prawdopobieństwie wylosowania \(x_2\)), co często nie jest spełnione w danych pogrupowanych w klastry, ponieważ obserwacje wewnątrz klastra są do siebie bardziej podobne. Innymi słowy obserwacje w klastrze są ze sobą skorelowane.

Drugą częścią definicji jest fakt, że obserwacje zostały wylosowane z tego samego rozkładu. Przykładowo, jeśli wykonujemy pojedynczy pomiar zadowolenia u różnych osób odwiedzających restaurację, to możemy powiedzieć, że wszystkie te obserwacje zostały wylosowane z tego samego rozkładu zadowolenia (by być dokładnym, z tego samego rozkładu w danym klastrze, czyli w tej samej restauracji). Jednakże, w praktyce wielokrotnie wykonujemy pomiary zadowolenia u tych samych osób, co oznacza, że zmienna “zadowolenie” ma inny rozkład u klienta A, a inny u klienta B.

By więc poprawnie przeanalizować takie dane, musimy zastosować model mieszany.

Czym jest więc model mieszany?

Wikipedia opisuje go tak: Mieszany model liniowy to model statystyczny, który uwzględnia zarówno efekty stałe (fixed effects), jak i losowe (random effects) w celu analizy zgrupowanych (clustered) danych.

Jest to definicja częstościowa, niemniej jednak pewnie każdy, kto miał styczność z modelami mieszanymi, słyszał o efektach stałych i losowych. Bez wchodzenia zbytnio w tę terminologię, efekty stałe są traktowane jako efekty, które mają określoną wartość, którą staramy się oszacować, na przykład wpływ zadowolenia na wysokość napiwku. Efekt losowy to efekt, który jest opisany jakimś rozkładem statystycznym, a my napotykamy na jakąś jego realizację, na przykład hojność danego klienta. Jednak gdy przeprowadzimy kolejne badanie na innych klientach, wciąż estymujemy ten sam związek między zadowoleniem a wysokością napiwku, ale nasi klienci będą mieli inną bazową wartość hojności.

My nie będziemy korzystać z tej terminologii. Według mnie czasami wporowadza ona zamieszanie w rozumieniu tej klasy modeli. Na szczęście mamy na to bardzo dobre uzasadnienie. Przypomnijcie sobie, że w statystyce Bayesowskiej wszystkie parametry traktujemy jako rozkłady statytyczne. Więc rozróżnienie na efekty stałe i losowe nie ma większego sensu, ponieważ w modelach Bayesowskich wszystkie efekty są losowe.

Modele wielopoziomowe

Model z efektem klastra

Jak więc będzie wyglądał bayesowski model, którym poprawnie możemy przeanalizować dane z powyższych przykładów?

W poprzednich częściach estymowaliśmy regresję liniową. Miała ona miała postać:

\[y_{i} \sim N(\mu = \alpha + \beta * x_i, \sigma)\] A wszystkie paramtery w tym równaniu miały rozkłady a priori ze zdefiniowanymymi przez nas hiperparametrami (czyli takimi parametrami, które są zdefiniowane arbitralnie, a nie estymowane z danych).

Teraz musimy uwzględnić w modelu różnice pomiędzy obserwacjami wynikajace z przynależności do klastra:

\[y_{i} \sim N(\mu = \beta * x_i + k_j , \sigma)\] \[k_{j} \sim N(\alpha, \tau)\] \[\alpha \sim U(-10^{3}, 10^{3})\] \[ \tau \sim N(10^{-6}, 10^{3}) \] \[ \beta \sim N(0, 10^{2}) \] \[ \sigma \sim N(10^{-6}, 10^{3}) \] Żeby dokładnie zrozumieć co się dzieje, spójrzmy na diagram naszego nowego modelu.

Żeby zrozumieć mechanizm tego modelu, warto zwizualizować sobie jak jest estymowany (mi to osobiscie pomagało zrozumieć jak działają bayesowskie modele). Przy użyciu MCMC, najpierw losujemy \(\alpha\) i \(\tau\) z ich rozkładów a priori, a następnie te wartości są używane do wylosowania wartości \(k_j\) dla każdego klastra, korzystając z rozkładu normalnego. Na koniec, mając wylosowane także wartości \(\beta\) i \(\sigma\), sprawdzamy, jak prawdopodobne jest, że model o takich parametrach wyprodukował zaobserwowane dane \(y\).

Parametr \(k_j\), który występuje w funckcji wiarygodności ma rozkład a priori opisany rozkładem normalnym o średniej \(\alpha\) i wariancji \(\tau\), który także jest estymowany z danych, a jego parametry mają własne rozkłady a priori. Reszta modelu pozostaje taka sama jak w przypadku standardowej regresji. Zamiast pojedynczego interceptu, model ma teraz \(j\) interceptów, gdzie \(j\) to liczba klastrów.

Rozkład a priori \(P(k_j|\alpha,\tau)\) określa prawdopodobieństwo wystąpienia intercepta o danej wartości, gdzie \(\alpha\) oznacza średni intercept dla całego zbioru danych, a \(\tau\) to wariancja interceptów.

Implementacja

Posłużymy się danymi (do pobrania tutaj) z drugiego przykładu i zaimplementujemy model użyciu JAGS.

library(rjags)
library(coda)

colnames(d) = c("Zadowolenie", "Klient", "Napiwek")

mod = "model {
  # Priors
  a ~ dunif(-1000, 1000)
  b ~ dnorm(0, 100^-2)
  sigma ~ dunif(0.000001,1000)
  tau ~ dunif(0.000001,1000)
  
  for (i in 1:5){
  k[i] ~ dnorm(a, tau^-2)
  }
  
  # Likelihood
  for (i in 1:length(Napiwek)) {
    mu[i] <- b * Zadowolenie[i] + k[Klient[i]]
    Napiwek[i] ~ dnorm(mu[i], sigma^-2)
  }}"

params = c("a","b","k","tau","sigma")
n.adapt = 100
ni = 3000
nb = 6000
nt = 1
nc = 6
jmod = jags.model(file = textConnection(mod), data = d, n.chains = nc, inits = NULL, n.adapt = n.adapt)
update(jmod, n.iter=ni, by=1)
post = coda.samples(jmod, params, n.iter = nb, thin = nt)

summary(post)

Iterations = 3101:9100
Thinning interval = 1 
Number of chains = 6 
Sample size per chain = 6000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

         Mean       SD  Naive SE Time-series SE
a     16.5100  7.87697 0.0415153      0.0765400
b      0.2777  0.02614 0.0001378      0.0003779
k[1]   6.2376  0.25082 0.0013220      0.0028611
k[2]  11.4953  0.23383 0.0012324      0.0024677
k[3]  16.2927  0.24670 0.0013002      0.0026946
k[4]  21.6225  0.22811 0.0012023      0.0023383
k[5]  26.5587  0.23112 0.0012181      0.0024175
sigma  0.8541  0.06148 0.0003240      0.0003498
tau   13.1714 11.67397 0.0615272      0.2492754

2. Quantiles for each variable:

         2.5%     25%     50%     75%   97.5%
a      2.3227 13.2364 16.4631 19.6847 30.4752
b      0.2266  0.2602  0.2776  0.2951  0.3298
k[1]   5.7422  6.0706  6.2380  6.4056  6.7278
k[2]  11.0347 11.3398 11.4964 11.6500 11.9617
k[3]  15.8088 16.1289 16.2932 16.4570 16.7785
k[4]  21.1711 21.4704 21.6255 21.7746 22.0685
k[5]  26.1064 26.4034 26.5576 26.7139 27.0142
sigma  0.7443  0.8110  0.8502  0.8936  0.9839
tau    5.2418  8.0037 10.5506 14.7879 35.9838

Jeśli spojrzymy na kwantyle parametru \(\beta\) zauważymy, że 0 znajduje się poniżej 2.5% percentyla. Innymi słowy, wykryliśmy efekt pomiędzy napiwkami a zadowoleniem z obsługi.

Zauważmy też, że model moglibysmy zaimplementować w następującej formie:

\[y_{i} \sim N(\mu =\alpha + \beta * x_i + k_j , \sigma)\] \[k_{j} \sim N(0, \tau)\]

Otrzymalibyśmy dokłokładnie to samo z tą różnicą, że teraz nie \(k_j\) nie byłoby interceptem dla klastra, a wartością o jaką różni się intercept dla klastra od \(\alpha\) (czyli, \(\alpha\) + \(k_j\) byłby intercepem dla klastra \(j\).)

Wspomniałem wyżej, że modele mieszane zasadniczo można interpretować w kategorii kontrolowania zmiennych zakłócających. Policzmy zwykłą regresję (częstościową, ale dotyczy to też bayesowskiego odpowiednika), w której uwzględnimy predyktory kategorialne wskazujące przynależność do klastra.

library(broom)
tidy(lm(Napiwek ~ Zadowolenie + Klient, d))

# A tibble: 6 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    6.23     0.246       25.3 4.42e-45
2 Zadowolenie    0.278    0.0256      10.8 1.67e-18
3 Klient2        5.26     0.261       20.1 9.36e-37
4 Klient3       10.1      0.260       38.6 1.61e-61
5 Klient4       15.4      0.263       58.6 1.23e-78
6 Klient5       20.3      0.262       77.6 1.84e-90

Otrzymaliśmy prawie taki sam współczynnik regresji! Różny sposób ujęcia problemu prowadzi do tego samego rozwiązania.

Ot, cała tajemnica modeli mieszanych. W identyczny sposób możemy rozszerzyć model, by uwzględniał różnice we współczynniku regresji pomiędzy klastrami. Bądź zmodyfikować go tak by uwzględniał więcej poziomów i.e. klastry wenątrz klastrów. I taki przykład sobie zrobimy gdy będziemy omawiać case studies, ale najpierw opowiedzmy sobie o właściwościach modeli wielopoziomowych.

Niezależność warunkowa

Łączny rozkład prawdopodobieństwa post priori naszych parametrów to:

\[P(\beta,k_j,\sigma,\alpha,\tau|y) \propto P(y|\beta,k_j,\sigma) P(k_j|\alpha,\tau) P(\beta)P(\alpha)P(\tau)P(\sigma)\]

To na co warto zwrócić uwagę, to relacje niezależności.

Warunkowa niezależność

Niech \(A\), \(B\) i \(C\) będą zdarzeniami. \(A\) i \(B\) są warunkowo niezależne jeśli:

\[P(A|B,C)= P(A|C) \] Ekwiwalentnie \[P(A,B|C) = P(A|C)P(B|C)\]

Czyli, jak widzimy na diagramie, bezpośrednio zależne są ze sobą tylko elementy połączone strzałkami. Na przykład \(P(y|\beta,k_j,\sigma)\) nie zależy bezpośrednio od \(P(\alpha)\) ponieważ:

\[P(y|\beta,k_j,\sigma,\alpha) = P(y|\beta,k_j,\sigma)\] Innymi słowy, cała informacja jaką \(\alpha\) ma o prawdopodobieństwie \(y\), jest już zawarta w \(k_j\).

Możemy zadać sobie pytanie, skoro ta klasa modeli nazywana jest (też) modelami wielopoziomowymi, to gdzie te poziomy. Jeśli mamy funkcję wiargodności i rozkłady a priori dla parametrów, które zarządzane sa przez wybrane przez badacza wartości hiperparametrów, to modelujemy dane na 1 poziomie. Jeśli jednak parametry rozkładów a priori same mają rozkłady a priori wtedy modelujemy dane na 2 (lub więcej) poziomach.

Bayesowskie dzielenie informacji

Rozważmy teraz następujący problem: zmierzyliśmy zmienną \(X\) w różnych podgrupach, a jedyne, co nas interesuje, to średnia wartość zmiennej \(X\) w danej podgrupie. Jednak podgrupy różnią się liczebnością. Mamy trzy możliwości postępowania:

• Policzyć średnią dla każdej grupy, ignorując pozostałe grupy w tym procesie oraz fakt, że dokładność średnich będzie mniejsza dla mniej licznych grup (no pooling).

• Policzyć średnią zmiennej \(X\) , uznając, że lepiej jest skorzystać z najpewniejszej estymaty i zignorować podział na grupy (full pooling).

• Użyć modelu mieszanego, który pozwoli nam wykorzystać wszystkie dostępne dane, poprzez wykorzystanie informacji z całej naszej próby (partial pooling, borrowing information).

W jaki sposób model mieszany pozwoli nam na to? Dla uproszczenia zdefiniujmy model, w którym interesuje nas prawdopodobieństwo otrzymania średniej \(y_{ij}\) w próbie dla danej podgrupy \(j\). Jako funkcję wiarygodności wybierzmy rozkład normalny:

\[P(y|\theta_j, \sigma_j) \sim N(\theta_j , \sigma_j)\] Parametry tego rozkładu to prawdziwa średnia w grupie \(\theta_j\) oraz odchylenie standardowe w grupie \(\sigma_j\).

Jednocześnie, rozkład prawdziwych średnich w grupach zamodelujmy następująco:

\[\theta_{j} \sim N(\mu, \tau)\] gdzie \(\mu\) to ogólna średnia zmiennej \(X\) , a \(\tau\) to odchylenie standardowe rozkładu ogólnej średniej.

Dobrze, to co nas naprawdę interesuje, to policzenie najlepszej estymaty prawdziwej średniej w grupie \(\theta_j\). Do tego potrzebujemy rozkładu \(P(\theta_j|\bar{x_{i}})\). Aby obliczyć szukany przez nas rozkład, możemy zastosować twierdzenie Bayesa:

\[P(\theta_j|y) \propto P(y|\theta_j)P(\theta_j)\]

Rozkład Normalny

Rozkład normalny dany jest wzorem: \[x \sim N(\mu,\sigma); \; \; P(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left( -\frac{(x-\mu)^2}{2\sigma^2}\,\right)\]
Gdzie \(\mu\) i \(\sigma\) to kolejno średnia i wariancja \(x\). Ponieważ w statystyce bayesowskiej często wystarczą nam relacje proporcjonalności możemy usnąć wszystko, co tylko skaluje nasz rozkład:
\[P(x) \propto \exp\left( -\frac{(x-\mu)^2}{2\sigma^2}\,\right)\]

Funkcja wiarygodności naszych danych to iloczyn prawdopdobieństwa otrzymania każdej obserwacji:

\[P(y|\theta_j) = \prod_{i=1}^{n} \exp\left( -\frac{(y_{ij}-\theta_j)^2}{\sigma^2}\,\right) = \exp\left( -\frac{\sum_{i=1}^{n}(y_{ij}-\theta_j)^2}{2\sigma^2}\,\right)\]

A więc nasz rozkład post priori będzie proporcjonalny do:

\[P(\theta_i|y) \propto \exp\left(-\frac{\sum_{i=1}^{n}(y_{ij} - \theta_j)^2}{2\sigma_j^2}\right)\exp\left(-\frac{(\theta_j-\mu)^2}{2\tau^2}\right)\]

Gdy po wykonaniu kilku przekształceń otrzymamy:

\[P(\theta_i|y) \propto \exp\left(-\frac{(\theta_j - \hat{\theta_j})^2}{2S_i^2}\right)\] Rozkład normalny jest rozkładem zgodym (conjugate prior) dla samego siebie. A więc rozkład a posteriori jest rozkładem normalnym o średniej \(\hat{\theta_i}\) i wariancji \(S^2\).

Gdzie

\[\hat{\theta_j} = \frac{\tau^2 \bar{x_j} + \frac{\sigma_j^2}{n}\mu}{\tau^2 + \frac{\sigma_j^2}{n_j}} = \left(\frac{\tau^2}{\tau^2 + \frac{\sigma_j^2}{n_j}}\right)\bar{x_j} + \left(\frac{\frac{\sigma_j^2}{n_j}}{\tau^2 + \frac{\sigma_j^2}{n}}\right)\mu\]
\[\bar{x_j} = \sum_{i=1}^{n_j}\frac{y_{ij}}{n_j}\]
\[\frac{1}{S_j^2} = \frac{1}{\tau^2} + \frac{1}{\frac{\sigma_j^2}{n}}\]

Czyli nasz rozkład posteriori \(\theta_j\) ma rozkład o średniej \(\hat{\theta}_{j}\) i wariancji \(S^2\). Zwykle, jeśli chcemy otrzymać punktową estymatę parametru z rozkładu posteriori, bierzemy jego średnią (w przypadku rozkładu normalnego jego średnia ma najwyższą wartość gęstości prawdopodobieństwa). W tym przypadku będzie to \(\hat{\theta_i}\). Możemy zauważyć, że wzór na \(\hat{\theta_i}\), który otrzymaliśmy, to średnia ważona średniej empirycznej w próbie w danej grupie \(\bar{x_j}\) i ogólnej średniej w próbie \(\mu\). Możemy zapisać to następująco:

\[\hat{\theta}_j = \lambda_j \bar{x_j} + (1-\lambda_j)\mu\]

gdzie:

\[\lambda_j = \frac{\tau^2}{(\tau^2 + \frac{\sigma_j^2}{n})}\]

Widzimy, że najlepsza estymata średniej w grupie będzie się różniła od średniej w próbie dla grupy \(j\) w zależności od tego, jak bardzo \(\lambda_i\) będzie mniejsze od 1. Będzie się tak działo, gdy \(\frac{\sigma_j^2}{n_j}\) będzie duży. A kiedy będzie duży?

Im mniejsza wariancja w danej grupie i im większa liczebność grupy, tym bardziej estymata średniej będzie się zbliżała do średniej empirycznej w grupie. Natomiast, jeśli wariancja w grupie jest duża i/lub grupa ma małą liczebność, estymata \(\hat{\theta}_i\) będzie bardziej przesunięta w stronę \(\mu\), ponieważ prawdopodobieństwo, że średnia w grupie \(j\) znacznie odbiega od średniej ogólnej, jest mniejsze (przypomnijmy sobie przykład z księgarnią z pierwszej części tutorialu).

\(\lambda_i\) nazywamy współcznynnikiem ściągania (shrinkage coefficient). Choć jego formuła nie będzie zawsze wyglądała tak elegancko, jak w tym przypdaku (gdzie rozkład a priori i funkcja wiarygodności dane są rozkładami normalnymi), widzimy, że model mieszany pozwala nam na estymację średnich w grupach z większą pewnością, ponieważ wykorzystujemy do tego obserwacje ze wszystkich podgrup!

Ta własność pozwala nam nie tylko na uzyskanie większej pewności naszych estymat, ale także na nie martwienie się problemem wielokrotnych porównań. Model w naturalny sposób przybliża średnie w grupach do średniej ogólnej gdy wariancja grup \(\tau^2\) zbliża sie do 0. Nie będziemy omawiać tego dokładniej, ale Andrew Gelman opisał to w tym artykule (Gelman et al., 2012), który wam polecam.

Modelujmy dalej

Zobaczmy teraz nieco bardziej skomplikowane przykłady modelu wielopoziomowego.

Radon jest radioaktywnym gazem szlachetnym, drugą po paleniu tytoniu, przyczyną raka płuc. Zanieczyszcza on wnętrza budynków przedostając się do środka przez małe szczeliny i otwory, stwarzając zagrożenie dla ludzi.

Mamy dane z 919 domów z 85 hrabstw, pomiar radonu w mieszkaniach oraz infromacje czy dom ma piwnicę (0) czy nie (1). Posiadamy także pomiar średniego stężenia uranu w glebie hrabstwa (Dane możecie je pobrać tutaj). Zobaczmy je.

data = read.csv("Radon_Data.csv")

head(data)

  county floor     radon    uranium
1      1     1 0.7884574 -0.6890476
2      1     0 0.7884574 -0.6890476
3      1     0 1.0647107 -0.6890476
4      1     0 0.0000000 -0.6890476
5      2     0 1.1314021 -0.8473129
6      2     0 0.9162907 -0.8473129

Modelowanie interceptu i wpspółczynnika regresji dla klastra

Radon uwalnia się z gleby, więc chcemy dowiedzieć się, jak posiadanie piwnicy wpływa na stężenie radonu. Nasz model będzie miał postać:

\[y_i = N(a_j + b_jx_i,\sigma)\] Indeks \(i\) oznacza numer obserwacji, a indeks \(j\) - numer hrabstwa. W naszym modelu uwzględniamy zarówno wyraz wolny (intercept), jak i współczynnik regresji dla danego hrabstwa. W ten sposób modelujemy korelacje pomiędzy obserwacjami w hrabstwach, które wynikają nie tylko z bazowego poziomu stężenia radonu w hrabstwie, ale także z różnic w poziomie stężenia radonu wynikających z posiadania piwnicy.

Warto zwrócić uwagę na fakt, że teraz każde hrabstwo ma swój własny intercept \(a_j\) i współczynnik regresji \(b_j\). Prawdopodobnie istnieje pewna zależność między tymi parametrami. Na przykład w hrabstwach o wysokim bazowym poziomie radonu (\(a_j\)) posiadanie piwnicy może mieć niewielki wpływ na poziom stężenia (\(b_j\)), a w hrabstwach o niskim bazowym poziomie radonu posiadanie piwnicy może znacząco wpłynąć na poziom stężenia. Taka zależność również wpływa na korelacje obserwacji w grupach (hrabstwach), dlatego należy ją uwzględnić w modelu.

Jak to zrobić? Aby to zrozumieć, przyjrzyjmy się wielowymiarowemu rozkładowi normalnemu. Podobnie jak opisujemy rozkład zmiennej \(x\) za pomocą rozkładu normalnego, tak możemy opisać łączny rozkład zmiennych \(x\) i \(y\) za pomocą wielowymiarowego rozkładu normalnego. Wielowymiarowy rozkład normalny jest opisany dwoma parametrami: wektorem średnich zmiennych \(\mu\) oraz macierzą kowariancji \(\boldsymbol{\Sigma}\), która opisuje związki między zmiennymi. Jeśli dwie zmienne są niezależne, to odpowiednie elementy macierzy kowariancji są równe zero.”

Wielowymiarowy Rozkład Normalny

Wielowymiarowy rozkład normalny opisuje sposób, w jaki zachowuje się \(m\) zmiennych losowych, podobnie jak jednowymiarowy rozkład normalny opisuje sposób, w jaki zachowuje się jedna zmienna losowa. Jest on opisany przez wektor średnich \(\boldsymbol{\mu}\) oraz macierz kowariancji \(\boldsymbol{\Sigma}\). Wielowymiarowy rozkład normalny oznaczamy symbolem \(\mathcal{N}_m(\boldsymbol{\mu}, \boldsymbol{\Sigma})\), gdzie \(m\) określa wymiarowość przestrzeni, w której zmienne losowe mają swoje wartości.

Funkcja gęstości prawdopodobieństwa dla \(m\)-wymiarowego rozkładu normalnego ma postać: \[f(\boldsymbol{x}) = \frac{1}{(2\pi)^{\frac{m}{2}}\sqrt{|\boldsymbol{\Sigma}|}} \exp\left(-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right)\]
gdzie \(\boldsymbol{x}\) jest \(m\)-elementowym wektorem zmiennych losowych. Symbol \(|\boldsymbol{\Sigma}|\) oznacza wartość wyznacznika macierzy \(\boldsymbol{\Sigma}\).

A więc rokład a priori dla \(a_j\), \(b_j\) będzie wyglądał tak:

\[ \begin{aligned} \begin{pmatrix}a_j \\ b_j \end{pmatrix} &\sim \mathcal{N_m}(\boldsymbol{\mu}, \boldsymbol{\Sigma}) \\ \end{aligned}\]

Gdzie

\[\boldsymbol{\mu} = \begin{pmatrix} \boldsymbol{\alpha} \\ \boldsymbol{\beta} \end{pmatrix}, \qquad \boldsymbol{\Sigma} = \begin{pmatrix} \mathrm{Var}(a_j) & \mathrm{Cov}(a_j,b_j) \\ \mathrm{Cov}(a_j,b_j) & \mathrm{Var}(b_j) \end{pmatrix}\] Oczywiście wszystkie parametry z powyższego wzoru mają swoje rozkłady a priori (możecie zobaczyć je w kodzie poniżej).

Nasz model w JAGS będzie wyglądał następująco:

model =  "model {
#Likehood
for (i in 1:length(radon)){
radon[i] ~ dnorm(radon.hat[i], tau.y)
radon.hat[i] <- a[county[i]]+ b[county[i]]*floor[i]
}

#Prior dla sigma
sigma.y ~ dunif (0, 1000)
tau.y <- pow(sigma.y, -2)

#Prior dla alpha
alpha ~ dnorm (0, .0001)

#Prior dla beta
beta ~ dnorm (0, .0001)

#Priory dla a_j i b_j
for (j in 1:max(county)){
mu[j,1] <- alpha 
mu[j,2] <- beta 
g[j,1:2] ~ dmnorm(mu[j,1:2],tau.g[1:2,1:2])
a[j] <- g[j,1]
b[j] <- g[j,2]
}

#Prior dla wariancji a_j
sigma.a ~ dunif (0, 100)
sigma.g[1,1] <- pow(sigma.a, 2)

#Prior dla wariancji b_j
sigma.g[2,2] <- pow(sigma.b, 2)
sigma.b ~ dunif (0, 100)

#Prior dla wspóczynnika korelacji a_j i b_j
rho ~ dunif (-1, 1)

#Kowariancja a_j i b_j
sigma.g[1,2] <- rho*sigma.a*sigma.b
sigma.g[2,1] <- sigma.g[1,2]

#Zamiana na macierz precyzji
tau.g[1:2,1:2] <- inverse(sigma.g[,])
}"

O ile nie interesuje nas relacja między piętrem a radonem w konkretnym hrabstwie, możemy śledzić jedynie wektor \(\boldsymbol{\mu}\) zawierający ogólny intercept (\(\alpha\)) i współczynnik regresji \(\beta\). Wyestymujmy także \(\rho\), czyli nasz współczynnik korelacji pomiędzy \(a_j\) i \(b_j\).

params = c("alpha","beta", "rho")
n.adapt = 100
ni = 3000
nb = 6000
nt = 1
nc = 4
jmod = jags.model(file = textConnection(model), data = data, n.chains = nc, inits = NULL, n.adapt = n.adapt)
update(jmod, n.iter=ni, by=1)
post = coda.samples(jmod, params, n.iter = nb, thin = nt)

summary(post)

Iterations = 3101:9100
Thinning interval = 1 
Number of chains = 4 
Sample size per chain = 6000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

         Mean      SD  Naive SE Time-series SE
alpha  1.4630 0.05470 0.0003531      0.0008216
beta  -0.6866 0.08453 0.0005457      0.0024429
rho   -0.2220 0.40279 0.0026000      0.0213931

2. Quantiles for each variable:

         2.5%     25%     50%       75%   97.5%
alpha  1.3563  1.4260  1.4625  1.499288  1.5720
beta  -0.8456 -0.7444 -0.6896 -0.630995 -0.5148
rho   -0.8533 -0.5100 -0.2832 -0.009378  0.8514

Tak jak się spodziewaliśmy, im większy indywidualny intercept, tym mniejszy współczynnik regresji dla posiadania piwnicy. Taki model, który właśnie stworzyliśmy, jest konceptualnym odpowiednikiem modelu z random interceptem i random współczynnikiem w statystyce częstościowej.

Predyktory na wyższych poziomach

Zauważmy, że w danych mamy jeszcze informacje o średnim stężeniu uranu w glebie w danym hrabstwie. Radon powstaje w wyniku rozpadu uranu, więc jest to istotny predyktor stężenia radonu w mieszkaniu. Jednakże, posiadamy tylko dane dotyczące hrabstwa, a nie mieszkania, co wymaga od nas zamodelowania związku na poziomie klastra. Jak będzie wyglądał nasz model? Funkcja wiarygodności pozostaje taka sama:

\[y_i = N(a_j + b_jx_i,\sigma)\]

ale nasz rozkład a priori \(a_j\) i \(b_j\) będzie wyglądał tak:

\[ \begin{aligned} \begin{pmatrix}a_j \\ b_j \end{pmatrix} &\sim \mathcal{N_m}(\begin{pmatrix} \boldsymbol{\alpha + \lambda u_j} \\ \boldsymbol{\beta} \end{pmatrix} \begin{pmatrix} \mathrm{Var}(a_j) & \mathrm{Cov}(a_j,b_j) \\ \mathrm{Cov}(a_j,b_j) & \mathrm{Var}(b_j) \end{pmatrix}) \\ \end{aligned}\]

Jak możemy zauważyć, intercept dla danego hrabstwa jest teraz modelowany jako zmienna losowa z rozkładu o średniej \(\alpha + \lambda u_j\). Widzimy też dosyć wyraźnie, że modelujemy relację pomiędzy stężeniem radonu a średnim stężeniem uranu na innym poziomie niż relację zmienną oznaczającą posiadanie piwnicy. Jeśli nazwiemy predyktorem dla pierwszego poziomu zmienną floor, to znaczy zmienną na poziomie każdej obserwacji (mieszkania), to średnie stężenie uranu wpływa tak samo na każde mieszkanie w hrabstwie, a więc jest predyktorem drugiego poziomu, to znaczy predyktorem interceptu dla klastra.

Możemy dopełnić skomplikowanie modelu w następujący sposób.

\[ \begin{aligned} \begin{pmatrix}a_j \\ b_j \end{pmatrix} &\sim \mathcal{N_m}(\begin{pmatrix} \boldsymbol{\alpha + \lambda u_j} \\ \boldsymbol{\beta + \gamma u_j} \end{pmatrix} \begin{pmatrix} \mathrm{Var}(a_j) & \mathrm{Cov}(a_j,b_j) \\ \mathrm{Cov}(a_j,b_j) & \mathrm{Var}(b_j) \end{pmatrix}) \\ \end{aligned}\]

Czym w takim wypadku będzie \(\gamma\)? Będzie nam mówił jak zmienia się współczynnik dla piętra w zależności od średniego poziomu uranu w hrabstwie. Innymi słowy będzie współczynnikiem dla interakcji dla predyktora pierwszego poziomu - piętra i predyktora drugiego poziomu - uranu.

Zaimplementujmy taki model.

model =  "model {
#Likehood
for (i in 1:length(radon)){
radon[i] ~ dnorm(radon.hat[i], tau.y)
radon.hat[i] <- a[county[i]]+ b[county[i]]*floor[i]
}

#Prior dla sigma
sigma.y ~ dunif (0, 1000)
tau.y <- pow(sigma.y, -2)

#Prior dla alpha
alpha ~ dnorm (0, .0001)

#Prior dla beta
beta ~ dnorm (0, .0001)

#Prior dla lambda
lambda ~ dnorm (0, .0001)

#Prior dla gamma
gamma ~ dnorm (0, .0001)

#Priory dla a_j i b_j
for (j in 1:max(county)){
mu[j,1] <- alpha + lambda*uranium[j]
mu[j,2] <- beta + gamma*uranium[j]
g[j,1:2] ~ dmnorm(mu[j,1:2],tau.g[1:2,1:2])
a[j] <- g[j,1]
b[j] <- g[j,2]
}

#Prior dla wariancji a_j
sigma.a ~ dunif (0, 100)
sigma.g[1,1] <- pow(sigma.a, 2)

#Prior dla wariancji b_j
sigma.g[2,2] <- pow(sigma.b, 2)
sigma.b ~ dunif (0, 100)

#Prior dla wspóczynnika korelacji a_j i b_j
rho ~ dunif (-1, 1)

#Kowariancja a_j i b_j
sigma.g[1,2] <- rho*sigma.a*sigma.b
sigma.g[2,1] <- sigma.g[1,2]

#Zamiana na macierz precyzji
tau.g[1:2,1:2] <- inverse(sigma.g[,])
}"

By użyć tego modelu musimy wyciągnąć uranium z naszych danych i stworzyć z niego wektor tylu wartości, ile jest hrabstw.

uranium = unique(data %>% select(county,uranium)) %>% arrange(county) %>% select(uranium)
data = as.list(data)
data[[4]] = uranium$uranium

Zaimplementujmy model

params = c("alpha","beta", "rho", "lambda","gamma")
n.adapt = 100
ni = 3000
nb = 6000
nt = 1
nc = 4
jmod = jags.model(file = textConnection(model), data = data, n.chains = nc, inits = NULL, n.adapt = n.adapt)
update(jmod, n.iter=ni, by=1)
post = coda.samples(jmod, params, n.iter = nb, thin = nt)

summary(post)

Iterations = 3101:9100
Thinning interval = 1 
Number of chains = 4 
Sample size per chain = 6000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean      SD  Naive SE Time-series SE
alpha   1.4662 0.03733 0.0002409       0.001365
beta   -0.6720 0.08772 0.0005662       0.003427
gamma  -0.4356 0.22992 0.0014841       0.008282
lambda  0.8169 0.09538 0.0006157       0.003355
rho     0.3093 0.43181 0.0027873       0.030863

2. Quantiles for each variable:

          2.5%      25%     50%     75%    97.5%
alpha   1.3962  1.44042  1.4654  1.4915  1.54106
beta   -0.8353 -0.73408 -0.6755 -0.6171 -0.48892
gamma  -0.8884 -0.58505 -0.4392 -0.2794  0.02015
lambda  0.6301  0.75404  0.8161  0.8799  1.00167
rho    -0.5815  0.00579  0.3332  0.6672  0.96624

Podsumujmy, co dowiedzieliśmy się z tego modelu. Średnia beta jest ujemna, co oznacza, że jeśli dom nie posiada piwnicy, to średnie stężenie radonu jest mniejsze. Średnia lambda jest dodatnia, co sugeruje, że im wyższy średni poziom radonu w danym hrabstwie, tym wyższe stężenie radonu w mieszkaniu w tym hrabstwie.

Jednak średnia gamma jest ujemna, co oznacza, że im wyższe stężenie uranu w danym hrabstwie, tym mniejsze stężenie radonu w domach bez piwnicy. To wydaje się być nieco dziwne, ale warto zwrócić uwagę na przedział wiarygodności. Dla gamma 95% przedział wiarygodności zawiera 0, co sugeruje, że ta estymata może nie być zbyt wiarygodna.

Podsumowanie

Zobaczyliśmy, jak tworzy się modele hierarchiczne w podejściu bayesowskim. Niemniej jednak, ponieważ ten post skupiał się tylko na wybranych aspektach, możliwe, że pozostajecie z niedosytem. W takim przypadku polecam książkę Gelmana i Hill o modelach hierarchicznych, w której w sposób systematyczny pokazują, jak modelować dane wielopoziomowe (Gelman & Hill, 2006).