Group Lasso に基づく因子分析モデリング九州大学大学院数理学府 廣瀬
九州大学大学院数理学研究院 小西1

慧貞則

はじめに

因子分析モデルとは，観測される多変量データの相関構造から，観測されない共通因子を見出すモデルであり，心理学，社会科学，生命科学をはじめとする諸科学の様々な分野で応用されている．因子分析モデルは通常最尤法によって推定されるが，観測変数の数や因子数が多くなると，モデルに含まれるパラメータ数が膨大となるため，推定が不安定になりやすい．さらに，因子分析を行う上で不要な観測変数が存在するとき，その変数が原因で推定が不安定になり，時には不適解（独自分散の推定値が負になる現象）が発生し，不適切なモデルを推定することがある（たとえば，狩野，1998）．また，因子分析モデルの目的は，推定されたモデルから観測されない共通因子を解釈することにあるが，因子数が変わると，推定されたモデルの解釈も変わる．
本報告では，これらの問題点に着目し，次の 3 つの項目を同時に実行できる新たな因子分析モデリング手法を提案する．

• 安定した推定量を求める．
• 因子分析を行う上で不要な観測変数が存在した時にその変数を削除する．
• 適切な因子数を選択する．
まず，対数尤度関数に罰則項を加えた罰則付き対数尤度関数の最大化によってモデルを推定することにより，推定の安定化を図る．ここで，どのような罰則項を用いるかが重要となるが，因子分析を行う上で不要な観測変数が存在した時にその変数を削除するため，Group Lasso (Yuanand Lin, 2006) の考え方に基づいて，因子負荷行列の各行ベクトルに L1 ノルムの罰則項を導入する．また，このモデルを EM アルゴリズム (Rubin and Thayer, 1982) と 2 次近似 (Fanand Li, 2001) を組み合わせたアルゴリズムを用いて推定する．さらに，このモデルに含まれる平滑化パラメータの値および因子数を客観的に選択するために，モデル評価基準 (Konishi andKitagawa, 2008) を導出する．提案した一連のモデリング手法を実データに適用し，その特徴と有効性を検証する．

2

因子分析モデルp 次元観測変数ベクトルを X = (X1 , · · · , Xp ) ，k 次元潜在変数ベクトルを F = (F1 , · · · , Fk )とする．因子分析モデルでは，観測変数と潜在変数の間に次のような線形関係を仮定する．
X − μ = LF + ε.

(1)

ただし，μ = (μ1 , · · · , μp ) は平均ベクトル，ε = (ε1 , · · · , εp ) は独自因子ベクトル，L = (ij )は p × k 因子負荷行列とする．以下，一般性を失うことなく μ = 0 の場合を考える．ここで，1

潜在変数ベクトル F と独自因子ベクトル ε は互いに独立とし，それぞれ次の多変量正規分布に従うとする．
F ∼ Nk (0, I), ε ∼ Np (0, Ψ).
(2)
ただし，Ψ = diag(ψ12 , · · · , ψp2 ) とする．この仮定のもとで，観測変数ベクトル X は平均ベク
トル 0，分散共分散行列 LL + Ψ の p 次元正規分布

X ∼ Np (0, LL + Ψ)(3)に従う．これより，第 i 番目の観測変数 Xi の分散はkVar(Xi ) =2ij + ψi2(4)j=1と表される．ここで， kj=1 2ij は共通性と呼ばれ，第 i 番目の観測変数 Xi の分散の中で，共通因子で表される分散の大きさを表す．(4) 式は，第 i 番目の観測変数 Xi の分散が共通性と独自分散の和で書けることを意味する．
N 個の p 次元データ x1 , · · · , xN が観測されたとき，モデルのパラメータ L, Ψ を最尤法によって推定することを考える．(3) 式より，対数尤度関数 l(L, Ψ) はl(L, Ψ) = −Np log(2π) + log |Σ| + tr(Σ−1 S)2(5)で与えられる．ただし，Σ = LL + Ψ とし，S は標本分散共分散行列

1 
S=xn xnNN(6)n=1とする．対数尤度関数 l(L, Ψ) の最大化によって最尤推定値 L̂, Ψ̂ を得る．ただし，(5) 式を最大にするパラメータ L, Ψ を解析的に求めることは困難であるので，反復法によって推定する必要がある．本報告では EM アルゴリズム (Rubin and Thayer, 1982) を適用することにより推定する．

3Group Lasso に基づく因子分析モデリング(1) 式で与えられる因子分析モデルの第 i 番目の観測変数を成分表示すると，
Xi =k

ir Fj + εij=1となる．ここで，i1 , · · · , ik がすべて 0 であるという状況を考える．このとき，

Xi = εiとなり，第 i 変数が独自因子のみの変数となる．すなわち，第 i 変数は共通因子で表されないと考えることができ，因子分析を行う上で不要な変数とみなすことができる．
本報告では，因子分析を行う上で不要な観測変数 Xi が含まれるとき，i1 , · · · , ik のすべてを 0 に縮小推定できるモデリング手法を提案する．

2

3.1

罰則項付き最尤法

Group Lasso (Yuan and Lin, 2006) は，線形回帰モデルにおいて，回帰係数ベクトルをグループに分けたとき，そのグループに含まれるパラメータの推定値をすべて 0 に縮小推定することのできる手法である．本報告では，Group Lasso の基本的な考え方に基づいて，対数尤度関数に因子負荷行列 L の各列に L1 ノルムの罰則項を加えた罰則付き対数尤度関数lλ (L, Ψ) = l(L, Ψ) − N λpi (7)i=1の最大化によりモデルを推定することを考える．ただし，i は因子負荷行列 L の第 i 行ベクトル，λ (> 0) は平滑化パラメータ，i  は k
i  = 
2ijj=1とする．この罰則項は，i1 , · · · , ik のすべてを 0 に縮小推定することができるという特徴を持つ．
しかし，(7) 式を用いてパラメータを推定すると，平滑化パラメータの値が小さい時は 0 に推定されず，ある程度大きな値になると同時にほとんど 0 に推定するという現象が起こった．この問題に対処するため，次の罰則項に重みをつけた罰則付き対数尤度関数の最大化を考える．

lλ (L, Ψ) = l(L, Ψ) − N λp

wi i .

(8)

i=1

ただし，wi は正数とする．これは，Zou (2006) や，Tateishi, Matsui and Konishi (2009) の提案した重み付きの Lasso 推定に基づいて構成されている．Tateishi, Matsui and Konishi (2009)
は，ガウス型基底関数に基づく非線形回帰モデルにおいて，重みをつけた Lasso を適用することにより，より安定した推定量を求めることを可能にした．
ここで，どのような重み wi を用いるかが重要となる．本報告では，変数間の相関に着目し，次のような重みを用いる．
p11  hh· ii .
wi =s(9)psh=1

ただし，sii は，分散共分散行列の逆行列の第 (i, i) 成分とする．sii は分散拡大係数とよばれ，変数相互間の相関の高さを表し，この値が大きいほど他の変数との相関が高いと考えられる．
ゆえに，この重みを採用することにより他の変数と相関の低い変数の共通性が 0 に推定されるという特徴を持つ．

3.2

推定

因子分析モデルを最尤法で推定するアルゴリズムは数多く提案されており，よく用いられるアルゴリズムとして，行列 Ψ−1/2 SΨ−1/2 または S −1/2 ΨS −1/2 の固有値問題を解くことによりL と Ψ を交互に推定するという方法がある (たとえば，Jennrich and Robinson, 1982)．
しかし，この方法を罰則付き最尤法の推定に拡張することは困難である．そこで，EM アルゴリズム (Rubin and Thayer, 1982) を用いて罰則付き対数尤度関数の最大化を行うことを考える．

3まず，潜在変数 F を不観測変数とみなし，観測データ X = x が与えられたもとでの F の事
後分布 F |x を計算する．次に，完全対数尤度関数に罰則項を加えた罰則付き完全対数尤度関数を計算し，その最大化を行う．しかし，罰則項をそのまま用いると M ステップでの因子負荷行
列 L の更新式を陽に表すことが困難となる．そこで，Fan and Li (2001) に基づいて，罰則項を次のように 2 次近似することを考える．
1(old)(old) (old)(i  ≈ i  +−i )for i = 0.
(10)iii(old)2i (old)ただし，iはパラメータを更新する前の値とする．この 2 次近似を適用することにより，L, Ψの更新式を求めることができ，次で与えられる．
−1 1 −1  −1ˆi = (M −1 + M −1 L Ψ−1 SΨ−1 LM −1 )/ψi2 + (λwi /i )IkM L Ψ Si. ,ψi2Ψ̂ = diag S − 2SΨ−1 LM −1 L̂ + L̂M −1 L̂ + L̂M −1 L Ψ−1 SΨ−1 LM −1 L̂ .
ただし，Si. は行列 S の第 i 列，M = L Ψ−1 L + Ik とする．なお，EM アルゴリズムにおいてステップを更新している途中に i  の値が十分小さな値に達したとき，i = 0 と推定する．

3.3

モデル評価基準Group Lasso に基づく因子分析モデルにおいて，どの観測変数を不要な変数と見なすかということが本質的となる．この不要な変数の選択は，平滑化パラメータ λ の選択と対応する．また，因子分析モデルにおけるもう一つの重要な問題として，適切な因子数をどのように選択するかということがある．本報告では，モデル評価基準 GBIC (Konishi et al., 2004) を Group
Lasso に基づく因子分析モデルの枠組みで導出し，因子数と平滑化パラメータの値を同時に選択する方法を提案する．
まず，(8) 式の最大化は，パラメータ L に次の事前分布を仮定して，事後分布のモードで推定することと同等である．
pπλ (L) = C
exp (−N λwi i ) .
(11)
i=1

ただし，C は基準化定数

C = (Ck )p−kk−1

Cj ,Cj =j=1(N λwj )j
j−12j π 2 Γ j+12(12)とする．ここで，Γ(·) はガンマ関数とする．
一般化ベイズ型モデル評価基準 GBIC は，ベイズ型モデル評価基準 BIC (Schwarz, 1978) を拡張することによって導出された評価基準である．BIC は最尤法で推定されたモデルを評価する基準であるが，GBIC は罰則付き最尤法で推定されたモデルを評価することのできる基準である．GBIC の基本的な考え方は，モデルが生起する確からしさを表す周辺尤度を最大とするようなモデルを選ぶというものである．一般に，周辺尤度の積分計算は困難であるため，ラプラス近似を適用することにより次の GBIC の一般式を得る．
GBIC = −2 log
f (x1 , · · · , xN |θ)π(θ|λ)dθ(13)≈ −p∗ log(2π) + p∗ log N + log |Jλ (θ̂p∗ )| − 2 log f (x1 , · · · , xN |θ̂p∗ ) + log π(θ̂p∗ |λ) .

(14)4ただし，p∗ はパラメータ数p∗ = p(k + 1) −k(k − 1)
2(15)とし，θ̂p∗ は事後分布のモード，Jλ (θ̂p∗ ) は
∂21log f (x1 , · · · , xN |θ) + log π(θ|λ)
Jλ (θ̂p∗ ) = −N ∂θ∂θ 

(16)θ̂p∗とする．
(11) 式の事前分布を (14) 式に代入することにより，提案した Group Lasso に基づく因子分析モデルの枠組みで GBIC を導出できる．

GBIC ≈ −p∗ log(2π) + p∗ log N + log |Jλ (θ̂p∗ )| + N p log(2π) + log |Σ̂| + tr(Σ̂−1 SN )−2 log C + 2N λpwi i .

(17)i=1

ただし，Σ̂, ψ̂i2 は Σ, ψi2 の事後分布のモードとする．なお，Jλ (θ̂p∗ ) に含まれる罰則付き対数尤
度関数の 2 階微分は次で与えられる．

∂ 2 lλ (θ)=N∂ab ∂cd(Σ−1 )(a,c) (L Σ−1 L)(b,d) + (Σ−1 L)(a,d) (Σ−1 L)(c,b)−(Σ−1 SΣ−1 )(a,c) (L Σ−1 L)(b,d) − (Σ−1 L)(a,d) (Σ−1 SΣ−1 L)(c,b)−(Σ−1 )(a,c) (L Σ−1 SΣ−1 L)(b,d) − (Σ−1 SΣ−1 L)(a,d) (Σ−1 L)(c,b)−(Σ−1 )(a,c) (Ik )(b,d) + (Σ−1 SΣ−1 )(a,c) (Ik )(b,d)(L)(a,b) (L)(c,d) (Ip )(a,c)(Ip )(a,c) (Ik )(b,d)+λwc− λwc 
3(LL )(c,c)(LL )(c,c)∂ 2 lλ (θ)=N∂ψi2 ∂cd(18)(Σ−1 )(c,i) (Σ−1 L)(i,d) − (Σ−1 SΣ−1 )(c,i) (Σ−1 L)(i,d)−(Σ−1 )(c,i) (Σ−1 SΣ−1 L)(i,d)∂ 2 lλ (θ)
N=222∂ψi ∂ψj,,(Σ−1 )2(i,j) − 2(Σ−1 )(i,j) (Σ−1 SΣ−1 )(i,j)

(19)
.

(20)

ただし，A(α,β) は，行列 A の (α, β) 成分とする．
ここで，共通性が 0 に推定されると，GBIC をそのまま計算することが困難となるという問題が生じる．実際，共通性が 0 に推定されると，(18) 式は発散する．この問題に対処するため，共通性が 0 となった観測変数に対して，推定する前に共通性が 0 であるという制約のもとで推定したと解釈する．すなわち，第 i 番目の変数の共通性が 0 と縮小推定されたら，i1 , · · · , ikはあらかじめ 0 という条件のもとで推定されたと見なす．また，それに伴い，パラメータ数も変わる．一般に，u 個の観測変数の共通性が 0 に推定されたら，パラメータ数はp∗ = p(k + 1) −k(k − 1)− uk2(21)となる．
なお，因子負荷行列の回転の不変性があるため，GBIC を導出する際に，一意性の条件が必要となる．本報告では，因子負荷行列の一意性の条件である ij = 0（i < j ）という制約 (たとえば，Anderson and Rubin, 1956 参照) のもとで GBIC を導出する．

5表 1: 因子数が 1 ∼ 6 のときの AIC，BIC，GBIC と平滑化パラメータ λ の値．AIC，BIC は最尤法によって推定されたモデルを評価しており，GBIC は Group Lasso に基づく因子分析モデルを評価している．
因子数1234564

AIC
1660.541605.39
1562.461532.431525.34
1521.25

BIC
1716.68
1687.72
1669.11
1661.541675.04
1689.66

不適解の発生提案手法（GBIC）

発生しなかった
発生しなかった
発生しなかった
発生した
発生した
発生した

1718.50
1705.18
1699.33
1708.89
1725.95
1748.31

λ0.0363
0.05130.08290.07960.12870.1286

実データへの適用提案したモデリング手法の有効性を検証するため，Kendall のデータ（Kendall, 1980）に適用した．このデータは，ある求人に対する 48 人の応募者についての 15 変数スコア，すなわち，N = 48, p = 15 である．その具体的なスコアは以下のとおりである．

1．履歴書
6．明晰さ
11．抱負

2．容姿
7．誠実さ
12．理解力

3．学力
8．技量
13．将来性

4．好感度
9．経験
14．志望程度

5．自信
10．志望動機
15．会社向き

このデータに対し，因子数が 1 ∼ 6 のとき，最尤法および提案手法を適用した．まず，最尤法で推定し，推定されたモデルを AIC，BIC で評価した．また，提案した Group Lasso に基づく因子分析モデルを推定し，推定されたモデルを GBIC で評価した．なお，不適解が発生した
ときは，独自分散が 0.005 以上という制約のもとで推定した．これらの結果を表 1 にまとめた．
表 1 より，AIC，BIC はともに不適解の発生したモデルを選択している．一方，提案手法は因子数 3 を選択し，かつ「学力」という変数を不要な変数と見なしている（表 3 参照）．
図 1 は，因子数を 3，4 としたときの横軸を平滑化パラメータ λ，縦軸を GBIC としたときのグラフである．これらのグラフから，GBIC は凸になっておらず，所々急激に小さな値を取るところがある．この境目が，不要な変数を取り除いた部分であり，たとえば，因子数 3 のときは，最初の変化点が学力という変数を取り除いたときとなっている．
因子数が 4 の時の最尤推定値と，提案したモデリング手法が選択したモデルを比較してみる．
因子数が 4 の時の最尤推定値を表 2 に，提案手法が選択したモデルを表 3 に示す． 表 2 の 4 因子モデルでは不適解が発生したが，第 4 因子をみると，因子負荷量は「学力」と「志望程度」のみ大きな値を取っており，そのほかの値は 0 に近い．よって，第 4 因子を「学力」と解釈することはできる．
一方，表 3 の結果を見ると，「学力」という変数を不要と見なしている．4 因子モデルのときには，「学力」という変数が重要な役割を果たしていたが，3 因子モデルでは不要であると見なされるのは興味深い．Group Lasso に基づいて推定された共通因子は，第一因子が「やる気と能力」，第 2 因子が「経歴と適性」．第三因子が「人柄」と適切な解釈ができる．
また，3 因子モデルで，提案手法が学力を「不要」とみなしたのは，λ = 0.0478 の時であった．一方，罰則項に重みをつけない (7) 式の最大化によりモデルを推定すると，学力を「不要」

6

1800
1780
1760
1700

1720

1740

GBIC

0.050.100.15

1780
1720

1740

1760

GBIC

1800

1820

1840

lambda0.05

0.100.150.20

lambda

図 1: 因子数 3（上図），4（下図）のときの，平滑化パラメータを動かしたときのグラフ．横軸が平滑化パラメータ λ，縦軸が GBIC の値を表す．

7

表 2: 4 因子モデルの最尤推定値

履歴書
容姿
学力
好感度
自信
明晰さ
誠実さ
技量
経験
志望動機
抱負
理解力
将来性
志望程度
会社向き

Factor1
0.13
0.45
0.07
0.23
0.92
0.84
0.25
0.89
0.09
0.77
0.90
0.78
0.73
0.42
0.36

Factor2
0.72
0.14
0.12
0.24
-0.10
0.12
-0.22
0.24
0.78
0.39
0.18
0.28
0.35
0.39
0.77

Factor3
0.11
0.24
0.00
0.83
0.15
0.30
0.74
0.08
-0.05
0.18
0.11
0.36
0.44
0.56
0.05

Factor4
-0.12
0.16
0.68
-0.05
-0.09
0.06
-0.02
-0.07
0.17
-0.06
-0.06
0.16
0.25
-0.59
0.14

uniqueness
0.44
0.69
0.52
0.20
0.11
0.19
0.34
0.13
0.36
0.22
0.14
0.15
0.09
0.00
0.25

表 3: 提案手法の選択したモデル（λ = 0.0478）

履歴書
容姿
学力
好感度
自信
明晰さ
誠実さ
技量
経験
志望動機
抱負
理解力
将来性
志望程度
会社向き

Factor1
0.08
0.27
0.00
0.10
0.79
0.69
0.18
0.74
0.05
0.63
0.75
0.65
0.58
0.32
0.26

Factor2
0.45
0.12
0.00
0.17
-0.12
0.08
-0.17
0.18
0.60
0.30
0.10
0.21
0.28
0.19
0.70

8

Factor3
0.10
0.17
0.00
0.86
0.11
0.25
0.56
0.08
-0.05
0.12
0.10
0.24
0.34
0.480.05

uniqueness
0.57
0.72
1.00
0.01
0.11
0.19
0.45
0.14
0.41
0.23
0.16
0.21
0.19
0.42
0.16

とみなしたのは，λ = 0.2 の時であった．すなわち，提案手法は，平滑化パラメータの値が小さくても不要な因子を削除することができるという特徴を持つことが分かった．

5

まとめと今後の課題

本報告では，安定した推定量を求め， 因子分析を行う上で不要な観測変数が存在した時にその変数を削除し，さらに適切な因子数を選択することのできる新たな因子分析モデリング手法を提案した．まず，推定の不安定性という問題に対処するために，対数尤度関数に罰則項を加えた罰則付き対数尤度関数の最大化によりモデルを推定した．ここで，どのような罰則項を用いるかということが重要となるが，Group Lasso (Yuan and Lin, 2006) の基本的な考え方に基づいて，因子負荷行列の各行ベクトルに L1 ノルムの重み付き罰則項を導入した．このモデルを推定するため，EM アルゴリズムと 2 次近似を組み合わせるアルゴリズムを用いた．さらに，このモデルに含まれる平滑化パラメータの値および因子数を同時に選択するために，モデル評価基準 (Konishi et al., 2004) を導出した．
また，提案したモデリング手法を Kendall のデータに適用し，その有効性を検証した．まず，最尤法で推定し，推定されたモデルを BIC で評価すると，因子数 4 が選択された．しかし，因子数 4 では不適解が発生しており，推定されたモデル解釈が困難となる．一方，提案手法は，「学力」という変数を取り除いた因子数 3 のモデルが選択された．このようなモデルは，最尤法では決して推定できないモデルであり，提案手法の有用性を実証できた．

参考文献[1] Anderson, T. W. and Rubin, H. (1956). Statistical inference in factor analysis, In J.
Neyman (Ed.), Proceedings of the Third Berkeley Symposium on Mathematical Statisticsand Probability, 5, (pp. 111–150). Berkeley: University of California Press.
[2] Efron, B., Hastie, T., Johnstone, I. and Tibshirani, R. (2004). Least angle regression(with discussion). Ann. Statist., 32, 407–499.
[3] Fan, J., and Li, R. (2001). Variable selection via nonconcave penalized likelihood and its
oracle properties. J. Amer. Statist. Assoc., 96, 1348–1360.
[4] Jennrich, R. I. and Robinson, S. M. (1969). A Newton-Raphson algorithm for maximum
likelihood factor analysis. Psychometrika, 34, 111–123.
[5] 狩野裕 (1998)．不適解の原因と処理．探索的因子分析．大阪大学人間科学部紀要.
[6] Kendall, M. G. (1980). Multivariate Analysis (2nd.ed.). Charles Griﬃn.
[7] Konishi, S., Ando, T. and Imoto, S. (2004). Bayesian information criteria and smoothing
parameter selection in radial basis function networks. Biometrika, 91, 27–43.
[8] Konishi, S. and Kitagawa, G. (2008). Information Criteria and Statistical Modeling.
Springer.

9[9] Rubin, D. B. and Thayer, D. T. (1982). EM algorithms for ML factor analysis. Psychometrika, 47, 69–76.
[10] Schwarz, G. (1978). Estimating the dimension of a model. Ann. Statist., 6, 461-464.
[11] Shimamura, T., Imoto, S., Yamaguchi, R. and Miyano, S. (2007). Weighted lasso ingraphical Gaussian modeling for large gene network estimation based on microarraydata. Genome Inform. 19, 142–153.
[12] Tateishi, S., Matsui, H. and Konishi, S. (2009). Nonlinear regression modeling via thelasso-type regularization. Preprint, MI 2009-8, Kyushu University.
[13] Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. Roy. Statist.
Soc. Ser. B, 58, 267–288.
[14] Yuan, M and Lin, Y. (2006). Model selection and estimation in regression with groupedvariables. J. Roy. Statist. Soc. Ser. B, 68, 49–67.
[15] Zou, H. (2006). The adaptive Lasso and its oracle properties. J. Amer. Statist. Assoc.,
101, 1418–1429.

10