• Graf: zbiór wierzchołków + zbiór krawędzi łączących pary wierzchołków.

• W modelach grafowych każdy wierzchołek odpowiada zmiennej losowej, a graf stanowi graficzną prezentację łącznego rozkładu całego zbioru zmiennych losowych.

• Graf nieskierowany: krawędzie nie mają kierunku.

• W nieskierowanych modelach graficznych (sieć Markova, pole losowe Markova) brak krawędzi ma szczególne znaczenie: niepołączone zmienne losowe są warunkowo niezależne, pod warunkiem pozostałych zmiennych.

• Rzadkie grafy (z niewielką liczbą krawędzi) są wygodne w interpretacji i bardzo użyteczne, np. w biologii molekularnej.

• Graf \(G\) składa się z pary (V,E), gdzie V - zbiór wierzchołków, E - zbiór krawędzi definiowany przez parę wierzchołków.

• Wierzchołki X i Y nazywamy sąsiadującymi, jeśli istnieje łącząca je krawędź (ozn. \(X \sim Y\)).

• Ścieżka \(X_1, X_2, \ldots, X_n\) to zbiór połączonych wierzchołków (ozn. \(X_{i-1}\sim X_i, \ i=2,\ldots,n.\))

• Graf pełny to graf, w którym każda para wierzchołków jest połączona krawędzią.

• Podgrafem \(U \in V\) nazywamy podzbiór wierzchołków wraz odpowiadającymi im krawędziami.

• Rozważmy graf \(G\) o zbiorze wierzchołków \(V\) odpowiadającym zbiorowi zmiennych losowych o łącznym rozkładzie \(P\).

• Zmienne X i Y są warunkowo niezależne pod warunkiem Z, jeżeli:

\[P(X|Y,Z)=P(X|Z)\]

• Brak krawędzi w grafie \(G\) oznacza warunkową niezależność dwóch niepołączonych zmiennych pod warunkiem pozostałych zmiennych. Oznaczamy:

\[X \perp Y|V',\ gdzie\ V'=V \setminus \{X,Y\}.\]

• Niech A, B oraz C będą podgrafami grafu \(G\). Wówczas C rozdziela A i B, jeśli każda ścieżka pomiędzy A i B przechodzi przez wierzchołek należący do grafu C.

• Grafy rozdzialające dzielą graf \(G\) na warunkowo niezależne części. Dla grafu \(G\) o podgrafach A, B, C oznaczamy:

\[C\ rozdziela\ A\ i\ B \implies A\perp B|C.\]

• Własność ta pozwala na dekompozycję grafu na na mniejsze, łatwiejsze w interpretacji i obliczeniach części - kliki (kompletne podgrafy).

• Ze względu na mnogość analitycznych własności, zakładamy że obserwacje mają wielowymiarowy rozkład normalny \(N(\mu, \Sigma).\)

• Wszystkie rozkłady warunkowe są również normalne.

• Odwrotność macierzy kowariancji, \(\Sigma^{-1},\) zawiera informację o kowariancji pomiędzy parami zmiennych \(i\) oraz \(j\), pod warunkiem pozostałych zmiennych.

• Jeśli \(\Theta = \Sigma^{-1},\) oraz \(\Theta_{ij}=0\), to \(i\) oraz \(j\) są warunkowo niezależne pod warunkiem pozostałych zmiennych.

Definicja:

Niech \(X\) będzie zmienną z rozkładu wielowymiarowego normalnego o średniej \(\mu\) i kowariancji \(\Sigma\) i niech \(x_i\), \(i = 1, \ldots, N\) będzie realizacją zmiennej \(X\). Empiryczną macierzą kowariancji nazwiemy macierz \(S\) taką, że \[S = \dfrac{1}{N}\sum_{i = 1} ^ N (x_i - \overline{x})(x_i - \overline{x})^T,\] gdzie \(\overline{x}\) jest próbkowym wektorem średnich. Łatwo pokazać, że \(S\) jest estymatorem największej wiarygodności macierzy \(\Sigma\).

\[l(\Theta) = \log \det\Theta - \text{trace}(S\Theta)\]

• Gdy znamy strukturę grafu zakładamy, że w przypadku braku krawędzi odpowiadający jej element macierzy \(\Theta = \Sigma^{-1}\) jest zerowy

• Maksymalizujemy logarytm wiarygodności przy założeniu, że niektóre z parametrów są równe zero

Logarytm wiarygodności ograniczony za pomocą stałych Lagrange’a:

\[l_C(\Theta) = \log\det\Theta - \text{trace}(S\Theta) - \sum_{(j, k)\notin E} \gamma_{jk}\theta_{jk}\] Otrzymujemy równanie gradientowe:

\[\Theta^{-1} - S - \Gamma = 0,\] gdzie \(\Gamma\) jest macierzą parametrów Lagrange’a z niezerowymi wartościami dla każdej pary wierzchołków niepołączonych krawędzią.

Z równania

\[\Theta^{-1} - S - \Gamma = 0\] Niech \(W = \Theta^{-1}\). Dla uproszczenia wyodrębniamy ostatni wiersz i ostatnią kolumnę. Powyższe równanie zapisujemy jako \[w_{12} - s_{12} - \gamma_{12} = 0.\] Z założenia, że \(W = \Theta^{-1}\) zapisujemy \[ \left( \begin{array}{cc} W_{11} & w_{12} \\ w_{12}^T & w_{22} \\ \end{array} \right) \left( \begin{array}{cc} \Theta_{11} & \theta_{12} \\ \theta_{12}^T & \theta_{22} \\ \end{array} \right) = \left( \begin{array}{cc} I & 0 \\ 0^T & 1 \\ \end{array} \right) \]

Stąd dostajemy \(w_{12} = - \frac{\theta_{12}}{\theta_{22}}W_{11}\).

Niech \(\beta = - \frac{\theta_{12}}{\theta_{22}}\), wtedy \(w_{12} = \beta W_{11}\).

Stąd dostajemy \[\beta W_{11} - s_{12} - \gamma_{12} = 0.\]

• Zakładamy, że mamy \(q\) krawędzi w grafie. Stąd, niech \(p - q\) będzie liczbą niezerowych elementów w \(\gamma_{12}\).

• Usuwamy \(p - q\) wierszy, które nie niosą żadnej informacji

• Redukujemy \(\beta\) do \(\beta^*\) usuwając zerowe elementy

Otrzymujemy układ równań wymiaru \(q\times q\):

\[W^*_{11}\beta^* - s^*_{12} = 0\]

Ze wzoru \[W^*_{11}\beta^* - s^*_{12} = 0\] możemy łatwo wyznaczyć \(\beta^* = W_{11}^{*-1}s^*_{12}\), co odpowiada \(q\) niezerowym elementom. Po uzupełnieniu do pełnego wymiaru \(p-q\) zerami otrzymujemy \(\hat{\beta}\).

Z równań \(\beta = - \frac{\theta_{12}}{\theta_{22}}\) i \(w_{12} = \beta W_{11}\) dostajemy

\[\dfrac{1}{\theta_{22}} = w_{22} - w_{12}^T \beta\]

1. Inicjujemy \(W = S\)
2. Powtarzamy dla \(j = 1, 2, \ldots, p, 1, \ldots\) aż do zbieżności:

• wyodrębniamy j-ty wiersz i j-tą kolumnę
• Rozwiązujemy \(W^*_{11}\beta^* - s^*_{12} = 0\) dla \(\beta^*\) używając zredukowanego układu równań. Otrzymujemy \(\hat{\beta}\) uzupełniając zerami w odpowiednich miejscach.
• zamieniamy wiersz i kolumnę \(w_{12} = W_{11}\hat{\beta}\)

3 . Rozwiązujemy (dla każdego j) \(\hat{\theta_{12}} = -\hat{\beta}\hat{\theta}_{22}\) z \(\dfrac{1}{\hat{\theta}_{22}} = s_{22} - w_{12}^T \hat{\beta}\)

\[ \hat{\Sigma} = \left( \begin{array}{cccc} 10.00 & 1.00 &1.31 &4.00 \\ 1.00& 10.00 &2.00 &0.87 \\ 1.31 &2.00& 10.00& 3.00 \\ 4.00& 0.87& 3.00& 10.00 \\ \end{array} \right)\]

\[\hat{\Sigma}^{-1} = \left( \begin{array}{cccc} 0.12 &−0.01& 0.00 &−0.05\\ −0.01 &0.11& −0.02& 0.00\\ 0.00 &−0.02& 0.11 &−0.03\\ −0.05 &0.00 &−0.03& 0.13\\ \end{array} \right) \]

• W większości praktycznych zastosowań struktura grafu nie jest znana, tj. nie wiemy które krawędzie powinniśmy pominąć.

• Estymacja w oparciu o regresje: lasso i SLOPE.

• Zamiast estymacji całej macierzy \(\bf{\Sigma}\) lub \(\bf{\Theta}=\bf{\Sigma}^{-1}\), szukamy jedynie jej niezerowych elementów.

• Regresja lasso: każda zmienna jest zmienną odpowiedzi, a pozostałe predyktorami (p modeli).

• Element \(\theta_{ij}\) jest niezerowy, jeżeli współczynniki przy zmiennej \(i\) (gdy \(j\) jest zmienną odpowiedzi) LUB przy zmiennej \(j\) (gdy \(i\) jest zmienną odpowiedzi) są niezerowe.

• Alternatywna wersja: koniunkcja zamiast alternatywy.

• Procedura asymptotycznie estymuje zbiór niezerowych elementów \(\bf{\Theta}\).

• Rozważamy problem maksymalizacji regularyzowanego logarytmu funkcji wiarygodności (penalized log-likelihood): \[l_{\lambda}(\bf{\Theta})=log\ det\ \bf{\Theta} - trace(\bf{S \Theta})-\lambda||\bf{\Theta}||_1,\] gdzie:
  • \(||\bf{\Theta}||_1\) jest normą \(L_1\) - sumą wartości bezwzględnych elementów macierzy \(\bf{\Sigma}^{-1},\)
  • S - estymator macierzy kowariancji \(\bf{\Sigma}\),
  • \(\bf{\Theta}\) - macierz precyzji,
  • \(\lambda\) - parametr regularyzacji.

• \(-l_{\lambda}(\bf{\Theta})\) jest wypukłą funkcją macierzy precyzji \(\bf{\Theta}.\)

• Maksimum możemy wyliczyć metodą lasso.

• Gradient: \[\bf{\Theta}^{-1}-\bf{S}-\lambda \cdot Sign(\bf{\Theta}) = 0,\] gdzie: \[Sign(\theta_{jk}) = \left\{\begin{matrix} sign(\theta_{jk}) & gdy & \theta_{jk} \neq 0 \\ \in [-1,1] & gdy & \theta_{jk} = 0 \end{matrix}\right.\]

• Dzielimy macierze \(\bf{\Theta}\) i \(\bf{W}=\bf{\Theta}^{-1}\) na dwie części: pierwsze p-1 wierszy i kolumn oraz p-te wiersz i kolumnę:

\[\begin{pmatrix} \bf{W}_{11} & w_{12} \\ w_{12}^{T} & w_{22} \end{pmatrix} \begin{pmatrix} \bf{\Theta}_{11} & \theta_{12} \\ \theta_{12}^{T} & \theta_{22} \end{pmatrix}= \begin{pmatrix} \bf{I} & 0 \\ 0 & 1 \end{pmatrix}\]

• Stąd dostajemy:

\[w_{12}=-\bf{W}_{11} \dfrac{\theta_{12}}{\theta_{22}} = \bf{W}_{11} \beta,\]

gdzie \(\beta = -\dfrac{\theta_{12}}{\theta_{22}}.\)

• Wstawiamy do równania gradientu i dostajemy (górny blok):

\[\bf{W}_{11} \beta-s_{12}+\lambda \cdot Sign(\beta)=0.\]

• Dla zmiennej odpowiedzi \(\bf{y}\) oraz macierzy predyktorów \(\bf{Z}\), lasso minimalizuje wyrażenie:

\[\dfrac{1}{2}(\bf{y}-\bf{Z}\beta)^{T} (\bf{y}-\bf{Z}\beta) + \lambda \cdot ||\beta||_{1}.\]

• Gradient tego wyrażenia jest postaci:

\[\bf{Z}^{T} \bf{Z}\beta-\bf{Z}^{T}\bf{y} + \lambda \cdot Sign(\beta) = 0.\]

• Zamiana \(\bf{Z}^{T}\bf{y}\) na \(s_{12}\) oraz \(\bf{Z}^{T} \bf{Z}\) na \(\bf{W}_{11}\) pozwala zastosować metodę lasso. Tak powstały problem optymalizacyjny nazywamy grafowym lasso.

Źródło: T.Hastie, R.Tibshirani, J.Friedman, ‘The Elements of Statistical Learning. Data Mining, Inference, and Prediction’.

• Powyższy algorytm jest bardzo efektywny. Przy umiarkowanej rzadkości danych pozwala rozwiązać problem o 1000 wierzchołków poniżej minuty.

• Parametr \(\lambda\) można dobrać tak, aby usunąć z grafu konkretną krawędź i zostawić te o silniejszych zależnościach. Wynika to z faktu, że dla \(\lambda_{jk}=\infty\) mamy \(\hat{\theta}_{jk}=0.\)

Źródło: T.Hastie, R.Tibshirani, J.Friedman, ‘The Elements of Statistical Learning. Data Mining, Inference, and Prediction’.

• Dobierając wartość \(\lambda\) chcemy kontrolować FWER.
• Metoda Banerjee (w modelach gaussowskich metoda Banerjee pozwala kontrolować FWER na poziomie \(\alpha\)):

\[\lambda^{Banejree}(\alpha)=\underset{i<j}{max}(s_{ii}, s_{jj})\dfrac{qt_{n-2}(1-\dfrac{\alpha}{2p^2})}{\sqrt{n-2+qt_{n-2}(1-\dfrac{\alpha}{2p^2})}},\] gdzie \(s_{ii}\) - wariancja i-tej zmiennej, \(qt_{n-2}(\phi)\) - kwantyl rzędu \(\phi\) rozkładu Studenta z \(n-2\) stopniami swobody.

• Wartości przy niektórych wierzchołkach mogą być niewidoczne (ukryte lub brakujące).
• Jeśli brakuje tylko kilku wartości, można je uzupełnić algorytmem EM.
• Jeśli przy wierzchołku brakuje wszystkich wartości, można je uzupełnić modelem opartym o obserwowane wierzchołki.

• Dla modelu liniowego \(y=X\beta+e,\) SLOPE ma postać problemu minimalizacji wyrażenia:

\[\underset{b}{min}\dfrac{1}{2}||y-Xb||^2_2+\lambda_1|b_{(1)}|+\lambda_2|b_{(2)}|+\ldots+\lambda_p|b_{(p)}|,\] gdzie:

• \(\lambda_1 \geq \lambda_2 \geq \ldots \geq \lambda_p,\)
• \(|b|_{(1)} \geq |b|_{(2)} \geq \ldots \geq |b|_{(p)},\)
• \(b_{(i)}\) oznacza i-tą statystykę pozycyjną.

• Dla \(\beta \in \mathbb{R}^p\) oraz \(\lambda \in \mathbb{R}^p,\) \(\lambda_1 \geq \ldots \geq \lambda_p,\) OL1 definiujemy jako: \[J_{\lambda}(\beta)=\sum_{i=1}^{p}\lambda_i|\beta|_{(i)}.\]

• W celu zastosowania \(J_{\lambda}\) w modelach grafowych musimy go zdefiniować w przestrzeni \(p\times p.\)

• Dla symetrycznej, dodatnio półokreślonej (ozn. \(S^p_{+}\)) macierzy precyzji \(\bf{\Theta} \in \mathbb{R}^{p\times p}\) oraz \(\lambda \in \mathbb{R}^{p^2},\) \(\lambda_1 \geq \ldots \geq \lambda_{p^2},\) definiujemy: \[J_{\lambda}(\bf{\Theta})=\sum_{i}\lambda_i|\theta|_{(i)}.\]

• Rozważamy logarytm funkcji wiarygodności:

\[l_{\lambda}(\bf{\Theta},\bf{X})=log\ det\ \bf{\Theta} - trace(\bf{S\Theta})-J_{\lambda}(\bf{\Theta}).\]

• Estymatorem macierzy precyzji jest zatem rozwiązanie problemu optymalizacyjnego:

\[\bf{\hat{\Theta}} \in \underset{\bf{\Theta \in S^p_{+}}}{arg\ max} \{log\ det\ \bf{\Theta} - trace(\bf{S\Theta})-J_{\lambda}(\bf{\Theta}) \}\] - Problem można rozwiązać przy pomocy algorytmu ADMM.

• Wybór ciągu \(\lambda_i\): poprawka Bonferonniego, metoda Holma, procedura Benjaminiego-Hochberga.