#はじめに
川野秀一先生・他著「スパース推定法による統計モデリング」共立出版の第5章のGraphical lassoの内容をまとめました。
#Graphical lassoとは
複数の確率変数間の統計的な独立性に着目し、ガウシアングラフィカルモデル$N(\mu,\Omega)$のネットワーク構造を推定することを考えます。
この時に、変数間の関係をスパースモデリングの考えを用いて推定する手法がGraphical lassoです。
##ガウシアングラフィカルモデルについて
多変量正規分布における$\Omega$は共分散逆行列であり、第$(i,j)$成分を$\omega_{ij}$とする
この時、着目する変数$X_j$を目的変数に、それ以外の変数$X_{-j}$を説明変数とした回帰モデルが成り立ちます。
X_j = \sum_{k≠j} \beta_{jk} X_{k} + \epsilon_j \tag{1}
ここで、 $\beta_{ij}=-\omega_{jk}/\omega_{jj}$、$\epsilon_j \sim
N(0,1/\omega_{jj})$
$\omega_{jk}=\omega_{kj}=0$の時、$X_j$と$X_k$は相互に影響を与えないため、$X_j$と$X_k$は他の変数を与えたもとでの条件付き独立となります。
#パラメータ推定
##共分散選択による推定
共分散選択と呼ばれる逐次的に変数を独立=偏相関を０に設定していく方法。
変数の数が少ないと良いが、変数の数が増えると計算時間が増えてしまうことになる。
##lasso回帰による推定
着目する変数$X_j$を目的変数に、それ以外の変数$X_{-j}$を説明変数とした(1)式に対して、lasso回帰による推定を行う方法が提案されている。
しかし、この推定方法では共分散逆行列の対称性が失われる場合がある。
##正則化最尤法による推定
多変量正規分布$N(0,\Omega^{-1})$を仮定した対数尤度関数に基づく損失関数に、$\Omega$に関数にL1正則化項を加えた関数を最小化することにより共分散逆行列を推定する。
l_{\lambda}(\Omega)=-log|\Omega|+tr(\Omega S)+ \lambda \sum_{j=1}^p \sum_{k=1}^p |\omega_{ik}| \tag{2}
ここで、$S=(s_{ij})$は、標本共分散行列である。
このように推定することの利点は3つあげられる。
- 共分散逆行列の対称性を保てる
- 正則化パラメータ$\lambda$が十分に大きい場合に計算が高速に行える
- 標本分散共分散行列を得られれば、直接グラフを推定できるため、順位相関係数を使った推定も可能で、データが正規分布に従っていなかったり、外れ値がある場合でも推定できる
##Graphical lassoによるパラメータ推定
実際は、ブロック座標効果法に用いて、$\Omega$のみを推定するのではなく、$\Sigma$と$\Omega$の両方を推定することになる。
はじめに、2式を$\Omega$に関して微分する
2(-\Omega^{-1}+S+\lambda \Gamma)-diag(-\Omega^{-1}+S+\lambda \Gamma) \tag{3}
ここで$\Gamma=(\gamma_{jk})$であり、 $\gamma \in sign^*(\omega_{jk})$となる。
$sign^*(・)$は、実数から実数の部分集合族への写像である。
sign^*(x)=\left\{
\begin{array}{ll}
{sign(x)} & (x \neq 0) \\
[-1,1] & (x = 0)
\end{array}
\right.
この時、次のような推定方程式が得られる。
\Omega^{-1}-S-\lambda \Gamma = O \tag{4}
$\Sigma$の対角成分は、$\Gamma$の対角成分は常に正であるので、$\sigma_{jj}=s_{jj}+\gamma$で得ることができる。
実際の推定においては、次のように$\Sigma$を初期化する形となる。
\Sigma=S+\gamma I_p \tag{5}
そして、$\Sigma$の非対角要素と$\Omega$を次に推定する。
まず、$\Sigma$、$\Omega$、$S$を次のように区分けする。
\Sigma=\begin{pmatrix}
\hat \Sigma & \sigma_1 \\
\sigma_1^T & \sigma
\end{pmatrix}
,
\Omega = \begin{pmatrix}
\hat \Omega & \omega_1 \\
\omega_1^T & \omega
\end{pmatrix}
,
S = \begin{pmatrix}
\hat S & s_1 \\
s_1^T & s
\end{pmatrix}
\tag{6}
ここで、$\hat \Sigma,\hat \Omega,\hat S$は$(p-1)\times(p-1)$行列であり、$ \sigma_1,\omega_1, s_1$は$(p-1)$次元ベクトル、$ \sigma_1,\omega_1, s_1$はスカラーである。
また、次の関係が成り立つ。
\begin{pmatrix}
\hat \Sigma & \sigma_1 \\
\sigma_1^T & \sigma
\end{pmatrix}
\begin{pmatrix}
\hat \Omega & \omega_1 \\
\omega_1^T & \omega
\end{pmatrix}
=
\begin{pmatrix}
I_{p-1} & 0 \\
0^T & 1
\end{pmatrix} \tag{7}
(4)式の右上のブロックから次の関係が得られる。
\sigma_1 - s_1 - \lambda\gamma = 0 \tag{8}
また、(7)式からは次の関係が得られる。
\omega_1=\omega \hat \Sigma^{-1} \sigma_1 = -\omega \beta_0 \tag{9}
以上より、次の関係が得られる。
\hat \Sigma \beta_0 - s_1 + \lambda \gamma = 0 \tag{10}
ここで、$\gamma_{j} \in sign^*(\beta_{0j})$。
この関係を満たす$\beta_0$は以下で得られ、lasso回帰により推定することになる。
\beta_0 = argmin_\beta \left\{\frac{1}{2}\|y - \hat \Sigma^{\frac{1}{2}} \beta \|^2_2 + \lambda \|\beta\|_1 \right\}, y=\hat \Sigma^{-\frac{1}{2}}s_1 \tag{11}
実際の推定では、$\beta_0$の推定を行った後に、次のように各値を更新していく。
\sigma_1 = \hat \Sigma \beta_0 \\
\omega = \frac{1}{\sigma - \sigma^T_1 \beta_0}\\
\omega_1 = - \omega\beta_0
\tag{12}
これを$j=1,...,p$に対して値が収束するまで操作を繰り返すことで、$\Sigma$の非対角要素と$\Omega$を推定することができる。
##交互方向乗数法によるパラメータ推定
Graphical lasso推定では、経時データにおける分散逆行列の推定は行うことができない。
そのため、Graphical lasso推定と結合lassoを組み合わせた Joint Graphical lassoを用いて推定することになる。
これについては、次の記事として考えているGraphical lassoによる異常検知についてでまとめようと思います。
#RによるGraphical lasso
##サンプルデータに対してのGraphical lasso
RではglassoパッケージでGraphical lassoを提供しています。
生成したサンプルデータに対して、真のネットワーク構造を推定することができるのかを検討します。
基本的なコードはこちらの記事を参考にさせてもらいました。
真のネットワーク構造を次のように設定しました。
library(qgraph)
n <- 200
p <- 10
trueGamma <- diag(p)
for(j in 5:p){
trueGamma[j,j-2] <- 0.5
trueGamma[j-2,j] <- 0.5
}
trueS <- solve(trueGamma)
true_graph <- ifelse(trueGamma!=0 & row(trueGamma)!=col(trueGamma),1,0)
qgraph(true_graph, layout = "circle", theme = "Borkulo")
設定した共分散行列からサンプルデータを生成します。
生成したサンプルデータの概形をpair plotで確認します。
library(dplyr)
# サンプルデータの生成
x <- matrix(rnorm(n*p),n,p) %*% chol(trueS)
#pairplotでデータの確認
library(GGally)
ggpairs(x)
結構傾向がわかりやすいデータになっていますね。
次に、BICを用いて正則化の強さを設定する$\lambda$の最適値を求めます。
library(glasso)
Omega <- cov(x)
lambda <- seq(0.1,1,length=n)
bic <- lambda
for(i in 1:n){
tmp_model <- glasso(Omega,lambda[i])
p_off_d <- sum(tmp_model$wi != 0 & col(Omega) < row(Omega))
bic[i] <- -2*(tmp_model$loglik) + p_off_d*log(n)
}
best <- which.min(bic)
lambda_and_bic <- data.frame(lambda, bic)
print(lambda_and_bic[best,])
gp <- ggplot(lambda_and_bic, aes(x = lambda, y = bic))
gp <- gp + geom_point()
plot(gp)
最もデータに当てはまりが良い構造が得られるのは、 $\lambda$を0.1678に設定した時のようです。
この値を採用した時のモデルを求め、変数間の関係性を可視化します。
gmodel <- glasso(Omega, lambda[best])
P <- gmodel$wi %>% round(4)
diag(P) <- 0
qgraph(P,layout = "circle", theme = "Borkulo", edge.labels=T)
設定した以外の変数間にも弱い関係性が出ていますが、概ね真の構造が推定できたと考えられます。
##attitudeデータに対するGraphical lasso
データ attitude は無作為に抽出した30部門に所属している35人の事務員のアンケート回答から得られたデータです。
pair plotでデータの概形を確認します。
x <- attitude %>% scale
n <- nrow(x)
ggpairs(x)
総合評価(rating)と従業員の苦情の取り扱い(complaints)は、高い正の相関を示していますが、その他の変数同士でもやや相関があるようです。
以下の処理は基本的にサンプルデータに対する処理と同じなので、省略します。
最も当てはまりの良いのは、$\lambda=0.6897$、$BIC=78.06$の時のようです。
推定されたグラフィカルモデルは次のようなものでした。
総合評価(rating)と従業員の苦情の取り扱い(complaints)のみが関係しているようです。
今回の例としてあまり良いデータではなかったかもしれません。
グラフィカルlassoは、今後、色々な場面で利用していこうかなと思います。
#参考
http://sun.econ.seikei.ac.jp/~shinmura/archive.html
https://rpubs.com/kur0cky/313430
https://www4.stat.ncsu.edu/~reich/BigData/code/glasso.html
https://bigdata.nii.ac.jp/eratokansyasai4/wp-content/uploads/2017/09/8f33ab8d22ce8cbc29f4a50ebc661982.pdf