統計数理（2020）
第 68 巻 第 2 号 287–303©2020 統計数理研究所［原著論文］

グループ正則化に基づく順序ロジットモデルにおける隣接クラスの統合
永沼 瑞穂1 ・吉川 剛平2 ・川野 秀一2（受付 2020 年 3 月 23 日；改訂 6 月 29 日；採択 6 月 29 日）

要

旨

本稿では，多クラス分類を目的とした順序ロジットモデルの枠組みにおいて，クラスを統合する方法を提案する．クラスを統合することによって，モデルの解釈が容易になるとともに，冗長なクラスが存在する場合はその冗長性を排除することができる．クラスの統合は隣接する各クラスの事後確率が等しくなるときに実行し，この目的のために構造的正則化の一つである
グループ正則化を用いる．グループ正則化は微分不可能な点を含むため，モデルに含まれるパラメータの推定値は，交互方向乗数法に基づく推定アルゴリズムにより得る．正則化パラメー
タの値は，グループ正則化により推定されたモデルの自由度をもつベイズ型情報量規準により選択する．モンテカルロ・シミュレーションおよび実データへの適用を通して，提案手法の有効性を検証する．
キーワード：グループ lasso，交互方向乗数法，順序付きカテゴリカルデータ，隣接カ
テゴリーロジットモデル．
1.

はじめに

多クラス分類は，自然科学・社会科学のさまざまな分野で用いられており，それを実行する
ための統計モデルは，クラスを表すカテゴリカルデータを目的変数，それに関係する共変量を説明変数としてモデル化される．多クラス分類においては，クラスの個数が多くなるときには，または冗長なク
注意が必要である．得られた解析結果の解釈が困難になる（Price et al., 2019）
ラスが生じる，といった現象が起きやすくなる. その例として，マーケティング分野における顧客市場のセグメンテーションを挙げることができる．顧客市場のセグメンテーションでは，大規模データ分析を行うことで顧客市場が過剰に細分化されてしまい，クラス毎の解釈が困難になるという現象が生じる．これはユーザーの分析コストを増大させ，マーケティング資源を浪費してしまうという恐れがある．このようなクラスの個数が多い状況においては，クラスを統合することによって，分析コストを削減することができ有用である．近年，そのクラスの統は，多項ロジットモデルと構造的正則化に基づいた統計手合を目的として，Price et al.（2019）（GFMR）を提案している．
法 group fused multinomial logistic regression
一方，カテゴリカルデータの中でも，順序関係を持つデータは順序付きカテゴリカルデータ
1

電気通信大学大学院 情報理工学研究科：〒 182–8585 東京都調布市調布ヶ丘 1–5–1（現 株式会社マクロ

ミル：〒 108–0075 東京都港区港南 2–16–1 品川イーストワンタワー 11F）
2

電気通信大学大学院 情報理工学研究科：〒 182–8585 東京都調布市調布ヶ丘 1–5–1

統計数理

288

第 68 巻 第 2 号 2020と呼ばれている．例えば，企業の格付けや 5 段階で表される学生の成績，アンケートの選好度などがある．このような順序付きカテゴリカルデータと共変量との間をモデリングする統計手法として，順序ロジットモデル（Agresti, 2010）が広く用いられている．前述したクラス過多な状況から生じる問題は，順序付きカテゴリカルデータを目的変数としたときも同様な理由で生じる．したがって，この場合でもクラスを統合することは有用である．Price et al.（2019）では，順序付きカテゴリカルデータを目的変数に用いる場合でも GFMR によってクラス統合を行っている．しかし，GFMR は多項ロジットモデルに基づいており，順序付きカテゴリカル
データを目的変数とするならば，順序ロジットモデルを用いる方が望ましい．
本稿では，順序ロジットモデルの枠組みの下，クラスを統合するための統計手法を提案する．
順序ロジットモデルの中でも，隣接カテゴリーロジットモデル（Goodman, 1984）を用い，クラスの統合は，隣接する各クラスの事後確率が等しくなるときに実行する．この目的のために，
構造的正則化の一つであるグループ正則化，特にグループ lasso（Yuan and Lin, 2006）を用い
る．グループ lasso の正則化項は微分不可能な点を含むため，モデルに含まれるパラメータの推定値は，交互方向乗数法（Boyd et al., 2011）に基づく推定アルゴリズムにより得る．構築したモデルには正則化パラメータが含まれるが，この値をグループ正則化により推定されたモデルの自由度をもつベイズ型情報量規準（Schwarz, 1978）により選択する．
本稿の構成は次の通りである．2 節では順序ロジットモデル，特に隣接カテゴリーロジット
モデルを紹介する．3 節では提案手法とその推定アルゴリズム，正則化パラメータの選択方法，
関連研究の GFMR について述べる．4 節ではモンテカルロ・シミュレーションおよび実データへの適用を通して，提案手法の有効性を検証する．5 節では本稿のまとめを行い，今後の課題について検討する．
2.

順序ロジットモデル

順序ロジットモデルは，各クラスの事後確率の対数オッズを，説明変数の線形結合で表すこと．この対数オッズの設定
によってモデル化される（Ananth and Kleinbaum, 1997; Agresti, 2010）
方法によって，累積ロジットモデル（McCullagh, 1980; Armstrong and Sloan, 1989; Whitehead，連続比ロジットモデル（Fienberg and Mason, 1979; Fienberg, 2007），隣接カテゴ
et al., 2001）（Anderson, 1984; Goodman, 1984;
リーロジットモデル（Adjacent-Categories Logit model, ACL）
Agresti, 1992）などのさまざまなモデルが提案されている．本稿ではこの中でも，隣接カテゴ
リーロジットモデルを用いる．
順序付きカテゴリカル変数 G ∈ {1, . . . , J} と p 次元説明変数 x† = (x1 , . . . , xp ) より，n組のデータ {(gi , x†i ); i = 1, . . . , n} が得られたとする．ここで，順序付きカテゴリカル変数とは，何らかの尺度についてのクラスの順序を示したものである．例えば，アンケート選好度
を変数とした場合，その取り得る選択肢として，1：あてはまる，2：ややあてはまる，3：や
やあてはまらない，4：あてはまらない，の 4 つのクラスは，順序付きカテゴリカル変数で
nある．また，説明変数に関するデータは標準化されているとする．つまり， i=1 x†ij = 0，n(x†ij )2 = n (j = 1, . . . , p) が成り立つものとする．
i=1
説 明 変 数 に 関 す る デ ー タ x†i が 与 え ら れ た 下 で ，j 番 目 の ク ラ ス に 属 す る 事 後 確 率 をPr(gi = j|x†i ) = πj (x†i ) (j = 1, . . . , J) とする．このとき，隣接カテゴリーロジットモデ
ルは
(2.1)

log

πj (x†i )πj+1 (x†i )= β0j + β1j xi1 + · · · + βpj xip = βj xi ,j = 1, . . . , J − 1となる．ここで，xi = (1, x†i ) ，βj = (β0j , β1j , . . . , βpj ) は j 番目のクラスに対する (p + 1)

グループ正則化に基づく順序ロジットモデルにおける隣接クラスの統合

289

次元の回帰係数ベクトルである．つまり，j 番目のクラスに属する確率と (j + 1) 番目のクラスに属する事後確率の対数オッズと，説明変数間をモデル化したものが隣接カテゴリーロジット
モデルである．
いま，i 番目のデータ x†i が j 番目のクラスに属しているとする．このとき，指示変数yij (i = 1, . . . , n; j = 1, . . . , J) を



yij =

(2.2)

1

データ i がクラス j に属しているとき

0

その他と定義すると，データ x†i に対する順序付きカテゴリカル変数 Gi は，J 次元のベクトルJ−1yi = (yi1 , . . . , yiJ ) と表すことができる．ここで，yiJ = 1 − j=1 yij が成り立つことに注意
しておく．
（2.2）
式の構成法から，確率 π1 (x†i ), . . . , πJ (x†i ) を持つ多項分布にしたがうことベクトル yi はがわかる．したがって，隣接カテゴリーロジットモデルの対数尤度関数は，(2.3)

(β1 , . . . , βJ−1 ) =

J−1
n


i=1

j=1

yij

J−1




βk xi − log

k=j

1+

J−1

j=1

exp

J−1




βk xik=jと表すことができる．（2.3）
式を最大にするパラメータは，解析的に陽に表現することが難しい
ので，非線形最適化手法，例えば，ニュートン・ラフソン法によって求める．なお，（2.3）式のを参照されたい．
導出過程については Agresti（2010）
3.

提案手法

本節では，隣接カテゴリーロジットモデルの枠組みの下，クラスを統合するために，モデルに基づに含まれるパラメータを構造的正則化の一つであるグループ lasso（Yuan and Lin, 2006）
いて推定する方法を提案する．その上で，推定アルゴリズムならびに正則化パラメータの選択
方法について述べる．また，関連研究として Price et al.（2019）の GFMR を紹介する．
3.1 クラス統合とグループ正則化対数尤度関数
Price et al.（2019）は，各クラスの事後確率が等しくなるときにクラスを統合することを提案した．本稿でもこの考えにしたがうものとする．しかし，本稿で扱っている隣接カテゴリー
ロジットモデルには，クラスに順序が存在しているため，隣接していないクラスを統合することはせず，隣接しているクラスのみを統合する．すなわち，本稿では，隣接した各クラスに対して，それらの事後確率が等しくなる場合にクラスを統合することを考える．以降，計画行列(x1 , . . . , xn ) の階数は (p + 1)，すなわちフルランクであると仮定する．
隣接カテゴリーロジットモデルにおいて，隣接したクラスの統合と回帰係数の関係について調べる．いま，j 番目のクラスと (j + 1) 番目のクラスを統合することを考える．すなわち，
πj (x†i ) = πj+1 (x†i ) (i = 1, . . . , n) となる状況を考える．このとき，log

πj (x†i )πj+1 (x†i )= 0,

i = 1, . . . , n,となる．この式と（2.1）
式より，βj xi = 0 (i = 1, . . . , n) となり，計画行列がフルランクであるという仮定を用いると βj = 0 となることがわかる．したがって，隣接カテゴリーロジットモ
デルにおいては，πj (x†i ) = πj+1 (x†i ) (i = 1, . . . , n) ならば βj = 0 であることがわかる．一方，上述の手順を逆に考えることにより，βj = 0 ならば πj (x†i ) = πj+1 (x†i ) (i = 1, . . . , n) というこ

統計数理

290

第 68 巻 第 2 号 2020とがわかり，以上より，πj (x†i ) = πj+1 (x†i ) (i = 1, . . . , n) となるための必要十分条件は βj = 0
である．なお，必要条件のみならば，計画行列のフルランク性の仮定は必要ないことに注意されたい．
したがってクラス統合を行うために，βj = 0 と推定することが可能なグループ正則化項が付与された対数尤度関数



(3.1)

(β1 , . . . , βJ−1 ) − λmax

β1 ,...,βJ−1J−1βj 2j=1を最大化し，パラメータを推定する方法を提案する．ここで，(β1 , . . . , βJ−1 ) は（2.3）式の隣
接カテゴリーロジットモデルに対する対数尤度関数，λ は非負の値をとる正則化パラメータ，
 · 2 はベクトルの L2 ノルムである．
（3.1）
式の第 2 項目は，隣接したクラスを統合するためのを用いる．このグルーグループ正則化項であり，本稿ではグループ lasso（Yuan and Lin, 2006）
プ正則化項により，βj (j = 1, . . . , J − 1) が 0 ベクトルか否かが判定され，もし，βj = 0 ならば
j 番目のクラスと (j + 1) 番目のクラスを統合する．なお，pj を βj の要素数と定義した場合，
J−1 √（3.1）式においては p1 = · · · = pJ−1
グループ lasso は通常 j=1 pj βj 2 と定式化されるが，√のため， pj は含めずに定式化していることに注意されたい．
3.2 推定アルゴリズム（3.1）式に含まれるグループ正則化項は，微分不可能な点をもつため，ニュートン・ラフソン
法等の非線形最適化手法を用いて推定値を求めることは困難である．本稿では，交互方向乗数（3.1）式に対するパラメータ推定値を法（Boyd et al., 2011）に基づいた推定アルゴリズムにより得る．
まず，（3.1）式を最小化問題に書き換えると



min
β

−(β) + λJ−1βj 2

j=1

となる．ここで，β = (β1 , . . . , βJ−1
) である．次に，微分不可能な正則化項を分離するために (p + 1) 次元の変数 zj (j = 1, . . . , J − 1) を用いて，以下の等式制約付きの最小化問題に変形
する．


min
β,z

−(β) + λJ−1zj 2

subject to

β = z.

j=1
ここで，z = (z1 , . . . , zJ−1
) である．したがって，この最小化問題に対する拡張ラグランジア
ンは
J−1
J−1
J−1


ρ
Lρ (β, z, u) = −(β) + λzj 2 +
u
(β
−
z
)
+
βj − zj 22
j
j
j
2
j=1

j=1

j=1



となる．ここで，uj (j = 1, . . . , J −1) は (p+1) 次元のラグランジュ乗数で，u = (u
1 , . . . , uJ−1 )である．また，ρ は正の値を取る調整パラメータである．
交互方向乗数法では，まず各パラメータ (β, z) について，拡張ラグランジアンを最小化するパラメータを推定する．次に，その最小となるパラメータを用いてラグランジュ乗数を勾配法（最急上昇法）により更新する．つまり，t (t = 1, 2, . . .) 回目の更新値をそれぞれ β t , zjt , ut とすると，

グループ正則化に基づく順序ロジットモデルにおける隣接クラスの統合(3.2)

β t+1 = argmin Lρ (β, z t , ut ),(3.3)

zjt+1 = argmin Lρ (β t+1 , zj , ut ),
zj

291

β

j = 1, . . . , J − 1,

ut+1 = ut + ρ(β t+1 − z t+1 )のように更新していく．
式の具体的な更新式を述べる．まず β の更新式
（3.2）について，右辺の最小（3.2）式と（3.3）
化は，

J−1
J−1

ρ
t
min −(β) +
ut
(β
−
z
)
+
βj − zjt 22
j
j
j
β
2
j=1

j=1

と同値である．この最小化問題はパラメータ β に対して微分可能であるため，更新値 β t+1 はニュートン・ラフソン法により得ることができる．ニュートン・ラフソン法を実行するための1 次導関数および 2 次導関数は順に，
∂Lρ
∂
=−
+ ρβ − ρz + u,
∂β
∂β
∂2
∂ 2 Lρ
=−+ ρI(J−1)(p+1)∂β∂β 
∂β∂β 

となる．ここで，I(J−1)(p+1) は (J − 1)(p + 1) × (J − 1)(p + 1) 型の単位行列であり，対数尤度関数の 1 次導関数と 2 次導関数の具体的な式は付録 A に載せている．
（3.3）
式について，右辺の最小化は，
次に z の更新式について述べる．



(3.4)

min λ
z

J−1


zj 2 +

j=1

J−1


ρ  t+1
βj − zj 22
2
J−1

t+1
ut
− zj ) +
j (βj

j=1

j=1

と同値である．ここで，等式
t+1ut− zj ) +j (βjρ
βjt+1 − zj
2

2

=
2

ρ
1
βjt+1 − zj + utj
2
ρ

に注意すると，問題（3.4）の z に関する最小化は，
J−1

j=1



1
ρ
zj −
2

2

λ
+ zj 2ρ2

1
βjt+1 + utj
ρ

2

−
2

ρ
2

1 t
uj
ρ

2

2



の z に関する最小化と同値となる．またこれは，各 j (j = 1, . . . , J − 1) 毎に問題が分離されているため，更新式は(3.5)

zjt+1 = S λρ1
βjt+1 + utjρ

,

j = 1, . . . , J − 1となる．ここで，Sv (a) はベクトルに対する軟閾値作用素，
Sv (a) =

1−

v
a2



a
+

であり，(·)+ は任意の実数 x に対して，(x)+ = max(0, x) と定義される. なお，更新式（3.5）の導出については川野 他（2018）
を参照されたい．

統計数理

292

第 68 巻 第 2 号 2020

3.3 正則化パラメータの選択
パラメータの推定値は，調整パラメータ ρ と正則化パラメータ λ を含んでいる．調整パラメータ ρ は Boyd et al.（2011）にしたがい ρ = 1 と設定する．一方，正則化パラメータ λ は情に基づいて選択する．
報量規準 BIC（Schwarz, 1978）
BIC はその第 2 項目にモデルの自由度を必要とする．Yuan and Lin（2006）
は，線形回帰モデ
ルにおいて，グループ lasso により推定したときの推定値 ξ̂ を持つモデルの自由度をJI(ξ̂j 2 > 0) +

j=1

J

ξ̂j 2j=1

ξ̂jLS 2(pj − 1)と定義した．ここで ξ̂jLS は最小二乗推定値である．本稿ではこの考えに基づいて，提案手法のモデルの自由度を
(3.6)

df =

J−1


I(β̂j 2 > 0) + p

j=1

J−1

β̂j 2j=1

β̂jridge 2とする．ここで，β̂j は最大化問題（3.1）の解，β̂jridge はリッジ推定値を表している．リッジ推定値に含まれる正則化パラメータの値は 10−5 に固定した．最尤推定値ではなくリッジ推定値を用いた理由は，4 節の数値実験において最尤推定値が求まらなかったことが多々あったためである．また，p1 = · · · = pJ−1 = p + 1 となる関係性に注意されたい．
（3.6）式の自由度を用いて以上より，本稿で用いる情報量規準 BIC は
BIC = −2(β̂) + df log(n)となる．BIC の値が最小となる正則化パラメータ λ の値を最適な値として採用する．
3.4 関連研究
Price et al.（2019）
は，順序関係のないカテゴリカルデータを目的変数とした多項ロジットモ
デルに対し，クラス統合を実行する方法として GFMR を提案した．GFMR は以下の最小化問
題として定式化される．

minβ1 ,...,βJ−MLM (β1 , . . . , βJ−1 ) + λβj − βm 2subject toβJ = 0.

(j,m)∈L

ここで，MLM (β1 , . . . , βJ−1 ) は多項ロジットモデルに対する対数尤度関数，λ は正の値をとる
正則化パラメータである．第 2 項目はクラス統合を行うための正則化項であり，候補集合 L は
S = {(a, b) ∈ J × J : a < b} からクラスの組み合わせを任意に設定する．ここで J = {1, . . . , J}である．例えば，1 番目のクラスと 3 番目のクラスが統合する可能性がある場合は，L = {(1, 3)}とする．クラスの関係性に何も事前の情報がなければ，L = S とする．βJ = 0 は多項ロジット
モデルに関する制約条件である．なお，クラスを統合する方法として，提案手法ではグループ
lasso，GFMR ではグループ連結 lasso（group fused lasso）（Bleakley and Vert, 2011; Alaíz et al.,2013）を用いていることに注意されたい．
推定アルゴリズムは，交互方向乗数法により得ている．更新式の詳細については Price et al.
（2019）
を参照されたい．正則化パラメータの値は，K 分割交差検証法により選択している．具
体的には，K 分割したデータセットのうち k 番目のデータセットを Vk としたとき，評価関数
(3.7)

K
J


k=1 i∈Vk l=1(−Vk )

yil log(π̂l(xi , λ))グループ正則化に基づく順序ロジットモデルにおける隣接クラスの統合

293を最大にする正則化パラメータの値を採用している．ここで，yil は（2.2）と同様にデータ i がクラス l に属しているときに 1，その他のときに 0 をとる指示変数である．π̂l(−Vk ) (xi , λ) は，k
番目のデータセットを除いたデータセットにより推定したモデルの事後確率を表す．
4.

数値実験

本節では，モンテカルロ・シミュレーションと実データへの適用を通して，提案手法の有効性を検証する．
4.1 モンテカルロ・シミュレーション
真のモデルを(4.1)

log

πj∗ (x†i )∗
πj+1(x†i )= βj∗ xi ,(i = 1, . . . , n, j = 1, . . . , J − 1)∗∗∗と設定した．ここで，回帰係数パラメータ βj∗ = (βj0
, βj1
, . . . , βj10
) の設定は次の 8 通りとした．

• 設定 1：(β1∗ , β2∗ , β3∗ ) = (0, δ, 0)• 設定 2：(β1∗ , β2∗ , β3∗ ) = (0, δ, δ)• 設定 3：(β1∗ , β2∗ , β3∗ , β4∗ ) = (0, δ, δ, δ)• 設定 4：(β1∗ , β2∗ , β3∗ , β4∗ ) = (0, δ, δ, 0)⎛

0• 設定 5：(β1∗ , β2∗ , β3∗ ) = ⎝ 1
0

0
0.5
0

0
• 設定 6：(β1∗ , β2∗ , β3∗ ) = ⎝ 1
1

0
1
1

⎛

0
1
0

0
0.5
0

0
1
0

0
0.5
0

0
1
0

0
0.5
0

0
1
0

⎞

0
0.5
0

0
1 ⎠
0

⎞

0
1
1

0
1
0.5

0
1
0.5

0
1
0.5

0
1
0.5

0
1
0.5

0
0.5
0.5

0
0.5
0.5

0
0.5 ⎠
0.5

0
⎜
1
• 設定 7：(β1∗ , β2∗ , β3∗ , β4∗ ) = ⎜
⎝ 1
1

0
1
1
0.5

0
1
1
1

0
1
0.5
0.5

0
1
0.5
1

0
1
0.5
0.5

0
1
0.5
1

0
1
0.5
0.5

0
0.5
0.5
1

0
0.5
0.5
0.5

0
⎜ 1
∗
∗
∗
∗

• 設定 8：(β1 , β2 , β3 , β4 ) = ⎜
⎝ 1
0

0
1
1
0

⎛

⎛

0
1
1
0

0
1
0.5
0

0
1
0.5
0

0
1
0.5
0

0
1
0.5
0

0
1
0.5
0

0
0.5
0.5
0

0
0.5
0.5
0

⎞

0
0.5 ⎟
⎟
0.5 ⎠
1

⎞

0
0.5 ⎟
⎟
0.5 ⎠
0

ただし，δ = (δ, . . . , δ ) であり，δ = 0.5, 1 の 2 通りを考えた．

  
11

ここで，設定したパラメータの意味について，設定 1 を例に取って説明する．設定 1 では
β1∗ = β3∗ = 0 であることから，クラス 1 とクラス 2，クラス 3 とクラス 4 はそれぞれ事後確率が等しいことがわかる．すると，3.1 節の議論より，クラス 1 とクラス 2，クラス 3 とクラス 4はそれぞれ同じクラスであり，実際には 2 クラスしかないことがわかる．つまり，設定 1 は，
見かけ上は 4 クラスであるが，実際は 2 クラスしかないシミュレーション設定である．これより，すべての設定における実際のクラスの個数と見かけ上のクラスの個数は，表 1 のようにまとめることができる．

294

統計数理

第 68 巻 第 2 号 2020

表 1．実際のクラスの個数と見かけ上のクラスの個数．

説明変数に関するデータ x†i は N (0, Σ) から発生させた．ここで，Σ は 10 × 10 型の分散共分散行列であり，その (i, j) 成分を ζ |i−j| とし，ζ = 0, 0.5 の 2 通りを考えた．（4.1）式より計算される真の事後確率 πj∗ (x†i ) (j = 1, . . . , J − 1) から，乱数を用いて各データ x†i (i = 1, . . . , n) が所属するクラスを決定した．サンプルサイズは n = 200, 300, 500 とし，シミュレーション回数は100 回とした．
による GFMR と 2 節で述べた ACL を用いた．ACL に含比較手法として，Price et al.（2019）
まれるパラメータはリッジ推定により求め，正則化パラメータの値は 10−5 に固定した．提案
手法と GFMR に含まれる正則化パラメータの値は，3.3 節の BIC により選択した．正則化パラを用い，提案手法では 0.1 から 7，GFMR
メータの値の候補はベイズ最適化（Snoek et al., 2012）
では 0.1 から 100 の範囲で探索した．ベイズ最適化の計算には統計ソフトウェア R のパッケー
ジ rBayesianOptimization を用いた．GFMR の候補集合 L については，順序付きカテゴリカ
で推奨されている L = {(1, 2), (2, 3), . . . , (J − 1, J)}
ルデータを解析する際に Price et al.（2019）
を用いた．
クラス統合後のクラスの個数，統合正答回数，真の事後確率と推定した事後確率間のカルの 3 つの観点からそれぞれの手法を評価した．100 回のシミュ
バック・ライブラー情報量（KL）
レーションにおいて，クラスが正しく統合された回数を，統合正答回数とした．ここで，各シ
ミュレーションにおいてクラスが正しく統合されるとは，見かけ上のクラスが実際のクラスに統合されることをいう．例えば，設定 1 においては，クラス 1 とクラス 2，クラス 3 とクラス 4がそれぞれ統合されるときは正しく統合されているといえるが，クラス 1 とクラス 2，クラス
2 とクラス 3 がそれぞれ統合されるときは正しく統合されているとはいえない．また，KL は∗KL(π̂, π ) =

J

j=1



log

π̂j (x†i , λ)πj∗ (x†i )π̂j (x†i , λ)と定義される．ただし，π̂j (x†i , λ) は推定したモデルから求めた事後確率である．100 回のシ
ミュレーション毎にテストデータを 100 個発生させて KL を計算し，得られた KL の平均と標準偏差を計算した．
表 2，表 3，表 4 はそれぞれ ζ = 0 とした場合の設定 1 と設定 2，設定 3 と設定 4，設定 5 から設定 8 に対する計算結果である．また，表 5，表 6，表 7 はそれぞれ ζ = 0.5 とした場合の設定 1 と設定 2，設定 3 と設定 4，設定 5 から設定 8 に対する計算結果である．なお，設定 3 の
n = 200, δ = 1, ζ = 0.5 においては，リッジ推定値が発散して計算することができなかった．
まず，ζ = 0 のときの結果を考察する．表 2，表 3，表 4 から，ほぼすべての場合において，提案手法は ACL よりも小さい KL の値を与えており，よりあてはまりの良いモデルを構築できていることがわかる．また，GFMR と比較しても，同様に，提案手法がよりあてはまりの良いモデルを構築できていることがわかる．統合正答回数の観点から見ると，サンプルサイズが200 の場合には GFMR が他の手法よりも多いが，サンプルサイズが大きくなるにつれ提案手法が多くなる．
次に，ζ = 0.5 のときの結果を考察する．表 5，表 6，表 7 から，ζ = 0 のときに比べ GFMR

グループ正則化に基づく順序ロジットモデルにおける隣接クラスの統合

295

表 2．設定 1 から設定 2 に対する結果（ζ = 0）．実験の各設定において，KL が一番小さい値，
統合正答回数が一番多い値を太字で表している．

はより良い結果を与えているが，総合的に見て提案手法の方が KL の値は小さい．また，統合
正統回数もサンプルサイズが大きくなるにつれ提案手法が GFMR より大きい値をとっている
ことがわかる．以上より，提案手法は従来手法よりも KL を小さくし，さらにサンプルサイズが大きい場合は正確にクラス統合を実行することがわかる．
提案手法において BIC により選択されるモデルと K 分割交差検証法（CV）により選択される式を用い，分割数は K = 5 とした．KL と統合正答モデルとを比較した．交差検証法には（3.7）
回数の 2 つの観点から評価し，KL については BIC による結果を分子，CV による結果を分母とし，統合正答回数については CV による結果を分子，BIC による結果を分母としてそれらの割合を計算した．これらより，KL と統合正答回数のどちらの指標においても，1 より小さい値の場合は本稿で用いた BIC の方がより良いといえる．
設定 1 から設定 4 の結果を表 8，設定 5 から設定 8 の結果を表 9 に載せている．KL については，ときおり 1 より大きい値をとっているが，基本的には 1 より小さい．統合正答回数については，ほとんどすべての場合において，1 より小さい値をとっている．したがって，BIC の方が KL を小さくし，統合正答回数を多くするモデルを選択することがわかる．
4.2 実データへの適用
赤ワインの品質ランクデータ（Cortez et al., 2009）に対し提案手法を適用した．本データは，

296

統計数理

第 68 巻 第 2 号 2020

表 3．設定 3 から設定 4 に対する結果（ζ = 0）．実験の各設定において，KL が一番小さい値，
統合正答回数が一番多い値を太字で表している．

UCI Machine Learning Repository（https://archive.ics.uci.edu/ml/index.php）から取得できる．赤ワインの品質は 6 段階で評価されており，段階が上がるにつれて赤ワインの品質の評
価は高くなる．つまり，1 段階目の赤ワインの品質が最も評価が低く，6 段階目に属する赤ワインの品質が最も評価が高い．この評価された段階を目的変数とし，以降段階をクラスと呼ぶ．
説明変数には，酒石酸濃度や残留糖分濃度などの赤ワインの成分に対する 11 種類の変数を用いた．各クラスのサンプルサイズは表 10 の通りである．
（3.6）
比較手法として GFMR を考え，BIC により正則化パラメータの値を選択した．ここでは
J−1
式の自由度に含まれるリッジ推定値が発散したため，自由度 df = (p + 1)
I(
β̂

> 0)
j
2
j=1

をもつ BIC を用いた．正則化パラメータの値の候補は，提案手法では 0.1 から 7，GFMR では
0.1 から 200 の範囲からベイズ最適化を用いて探索した．GFMR の候補集合 L については，4.1
節と同様に L = {(1, 2), (2, 3), . . . , (J − 1, J)} を用いた．
ブートストラップ法（Efron and Tibshirani, 1993）を用いてクラス統合の安定性を検証した．
まずクラス 1 をクラス 1a とクラス 1b に，サンプルサイズが同じになるようランダムに 2 分割
し，このデータセットからブートストラップ標本を発生させて提案手法と GFMR を当てはめた．この操作を 100 回繰り返し，クラス 1a とクラス 1b が統合された回数（統合正答回数）を数え上げた．なお，統合正答回数の数え上げでは，各ブートストラップ標本においてクラス 1a とクラス 1b だけを統合するときに限って数え上げることに注意しておく．クラス 2 からクラス

グループ正則化に基づく順序ロジットモデルにおける隣接クラスの統合

297

表 4．設定 5 から設定 8 の結果（ζ = 0）．実験の各設定において，KL が一番小さい値，統合
正答回数が一番多い値を太字で表している．

6 に対しても同様の手順でブートストラップ標本を発生させ，統合正答回数を数え上げた．
結果を表 11 にまとめている．クラス 2 とクラス 5 以外の結果については，提案手法は統合
正答回数が 80 回以上であり，GFMR より多くなっている．特に，クラス 1 とクラス 6 の結果
においては，提案手法の統合正答回数は GFMR の 2 倍以上である．一方，クラス 2 やクラス 5
では，そもそも統合正答回数が少ない．したがって，このときの結果についてはあまり議論すべきではないかも知れないが，クラス 2 においては提案手法の統合正答回数は GFMR より多くなっている．以上より，統合正答回数の観点から，提案手法は GFMR よりも多くの場合において安定してクラス統合が行われていることがわかる．
5.

まとめ

隣接カテゴリーロジットモデルとグループ lasso により，順序付きカテゴリカルデータを対象とするモデルにおけるクラス統合の方法を提案した．推定アルゴリズムを交互方向乗数法により構成し，グループ lasso により推定されたモデルの自由度をもつ情報量規準 BIC により正則化パラメータの値を選択した．数値実験の結果をいくつかの観点から考察することにより，従来手法に比べて提案手法は多くの点で有用であることがわかった．
4.1 節のモンテカルロ・シミュレーションにおいて，サンプルサイズが小さいときに，提案

298

統計数理

第 68 巻 第 2 号 2020

表 5．設定 1 から設定 2 の結果（ζ = 0.5）．実験の各設定において，KL が一番小さい値，統合正答回数が一番多い値を太字で表している．

手法が GFMR に統合正答回数で劣ることがあった．この理由としては，次の 2 点が考えられ
る．1 点目は，サンプルサイズとクラス間の順序との関係性である．一般に，サンプルサイズが大きい場合には各クラスの性質が明確になるため，クラス間の順序が強調されることが知られている．一方，サンプルサイズが小さい場合には各クラスの性質が明確にならないため，クラス間の順序が強調されず，提案手法の統合正答回数は少なくなったと考えられる．2 点目は，正則化パラメータの選択方法である．サンプルサイズが小さい場合には，統合後のクラスの個数が真のクラスの個数よりも大きくなり，その結果統合正答回数が少なくなる状況が見受けられた．これは，正則化パラメータの値が小さめに選択されていることを表しており，BIC を用
いた選択方法がうまく機能していないことを示唆している．
4.2 節の実データへの適用において，クラス 2 とクラス 5 の統合正答回数が非常に少なかった．これは次の 2 つの理由が重なったため生じたと考えられる．1 つ目は，クラス 2 とクラス 5はそれぞれ 2 つのクラスに挟まれており，統合する隣接クラスの候補が多いということ．2 つ目は，クラス 3 とクラス 4 に比べて，サンプルサイズが小さいということである．このような状況において，統合正答回数を改善する統計モデルの構築に取り組むことは大変意義がある．
本稿で用いた BIC におけるモデルの自由度は Yuan and Lin（2006）の考えに基づいて導出しで導出されている自由度は線形回帰モデルかつ計画行列の各列が直たが，Yuan and Lin（2006）
交するという仮定が課されている．この点において提案モデルと乖離しており，提案モデルの

グループ正則化に基づく順序ロジットモデルにおける隣接クラスの統合299

表 6．設定 3 から設定 4 の結果（ζ = 0.5）．実験の各設定において，KL が一番小さい値，統合正答回数が一番多い値を太字で表している．

自由度についてはまだ多くの問題が残っている．以上は今後の研究課題としたい．
謝辞

本稿の改訂に当たって，編集委員および査読者の方から大変貴重なご指摘と適切なご意見をいただきました．ここに記して御礼申し上げます．本研究は日本学術振興会科学研究費補助金（19K11854）の助成を受けたものです．
付
A.

録

隣接カテゴリーロジットモデルの対数尤度関数の 1 次導関数および 2 次導関数

（2.3）式の隣接カテゴリーロジットモデルの対数尤度関数は，(A.1)

(β1 , . . . , βJ−1 ) =

n

i=1



∗∗ ∗
1
J−1 Yi Xi β − log



1+

J−1


exp

 

Ij Xi∗ β





j=1

と書き換えることができる．ここで，1J−1 は要素がすべて 1 の (J − 1) 次元ベクトル，β = (β1 , . . . , βJ−1) は (J − 1)(p + 1) 次元ベクトル，Ij = (I(1 > j − 1), . . . , I(J − 1 > j − 1))

統計数理

300

第 68 巻 第 2 号 2020

表 7．設定 5 から設定 8 の結果（ζ = 0.5）．実験の各設定において，KL が一番小さい値，統合正答回数が一番多い値を太字で表している．

は (J − 1) 次元ベクトル，Yi∗∗ は

⎛

yi1
⎜0
⎜
Yi∗∗ = ⎜ .
⎝ ..
0

yi1
yi2
..
.
0

···
···
..
.
···

⎞

yi1
yi2 ⎟
⎟
.. ⎟
. ⎠
yiJ−1

となる (J − 1) × (J − 1) 型の正方行列，Xi∗ は⎛x
i⎜0 · · · 0⎜
Xi∗ = ⎜ .
⎝ ..
0···0

0···0
x
i
..
.
0···0

···
···
..
.
···

⎞

0
0 ⎟
⎟
.. ⎟
. ⎠
x
i

となる (J − 1) × (J − 1)(p + 1) 型の行列である．このとき，（A.1）式の対数尤度関数の 1 次導関数および 2 次導関数は，順に


∂
=
∂β
n

i=1


 

1n Yi∗∗ Xi∗



−

1

1+

J−1

J−1


η
j=1 ij j=1ηij γij,

グループ正則化に基づく順序ロジットモデルにおける隣接クラスの統合

301

表 8．設定 1 から設定 4 における正則化パラメータ選択法の比較．実験の各設定において，1より小さい値を太字で表している．

表 9．設定 5 から設定 8 における正則化パラメータ選択法の比較．実験の各設定において，1より小さい値を太字で表している．

表 10．各クラスのサンプルサイズ．

表 11．各クラスのブートストラップ選択回数．各クラスにおいて統合正答回数が多い方を太字にしている．

⎡
⎢
∂ 
=
⎣
∂β∂β 
2

n

i=1

1

1+

J−1
j=1

ηij

 

2

J−1

j=1


∗

ηij γij J−1

j=1

⎤



ηij γij−

1

1+となる．ここで，ηij = exp Ij Xi β , γij = Ij Xi∗ と置いている．

J−1


J−1

η
j=1 ij j=1⎥ηij γij γij
⎦

302

統計数理

参

第 68 巻 第 2 号 2020

考

文

献

Agresti, A. (1992). Analysis of ordinal paired comparison data, Journal of the Royal Statistical Society:
Series C (Applied Statistics), 41(2), 287–297.
Agresti, A. (2010). Analysis of Ordinal Categorical Data (Second Edition), John Wiley & Sons, NewYork.
Alaíz, C. M., Barbero, A. and Dorronsoro, J. R. (2013). Group fused lasso, International Conferenceon Artificial Neural Networks, 66–73, Springer, Berlin, Heidelberg.
Ananth, C. V. and Kleinbaum, D. G. (1997). Regression models for ordinal responses: A review ofmethods and applications, International Journal of Epidemiology, 26(6), 1323–1333.
Anderson, J. A. (1984). Regression and ordered categorical variables, Journal of the Royal StatisticalSociety: Series B, 46(1), 1–22.
Armstrong, B. G. and Sloan, M. (1989). Ordinal regression models for epidemiologic data, American
Journal of Epidemiology, 129(1), 191–204.
Bleakley, K. and Vert, J.-P. (2011). The group fused lasso for multiple change-point detection, arXiv
preprint, arXiv:1106.4199.
Boyd, S., Parikh, N., Chu, E., Peleato, B. and Eckstein, J. (2011). Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends ® in
Machine Learning, 3(1), 1–122.
Cortez, P., Cerdeira, A., Almeida, F., Matos, T. and Reis, J. (2009). Modeling wine preferences bydata mining from physicochemical properties, Decision Support Systems, 47(4), 547–553.
Efron, B. and Tibshirani, R. J. (1993). An Introduction to the Bootstrap, Chapman & Hall/CRC, NewYork.
Fienberg, S. E. (2007). The Analysis of Cross-classified Categorical Data, Springer, New York.
Fienberg, S. E. and Mason, W. M. (1979). Identification and estimation of age-period-cohort models
in the analysis of discrete archival data, Sociological Methodology, 10, 1–67.
Goodman, L. A. (1984). The Analysis of Cross-classified Data Having Ordered Categories, Harvard
University Press, Cambridge, Massachusetts.
川野秀一, 松井秀俊, 廣瀬慧 (2018).『スパース推定法による統計モデリング』, 共立出版, 東京.
McCullagh, P. (1980). Regression models for ordinal data, Journal of the Royal Statistical Society:
Series B, 42(2), 109–127.
Price, B. S., Geyer, C. J. and Rothman, A. J. (2019). Automatic response category combination in
multinomial logistic regression, Journal of Computational and Graphical Statistics, 28(3), 758–
766.
Schwarz, G. (1978). Estimating the dimension of a model, The Annals of Statistics, 6(2), 461–464.
Snoek, J., Larochelle, H. and Adams, R. P. (2012). Practical Bayesian optimization of machine learningalgorithms, Advances in Neural Information Processing Systems 25, 2951–2959.
Whitehead, A., Omar, R. Z., Higgins, J. P., Savaluny, E., Turner, R. M. and Thompson, S. G. (2001).
Meta-analysis of ordinal outcomes using individual patient data, Statistics in Medicine, 20(15),2243–2260.
Yuan, M. and Lin, Y. (2006). Model selection and estimation in regression with grouped variables,
Journal of the Royal Statistical Society: Series B, 68(1), 49–67.

Proceedings of the Institute of Statistical Mathematics Vol. 68, No. 2, 287–303 (2020)

303

Fusing Adjacent Classes in an Ordinal Logistic Model
via Group Regularization
Mizuho Naganuma1 , Kohei Yoshikawa2 and Shuichi Kawano2
1 Graduate School of Informatics and Engineering, The University of Electro-Communications;

Now at Macromill, Inc.
2 Graduate School of Informatics and Engineering, The University of Electro-Communications

This paper aims to fuse adjacent classes in an ordinal logistic model in light of themulti-class classification problem. Fusing the classes enables us to easily interpret theconstructed model and remove irrelevant classes. Fusion of classes is performed when two
adjacent classes have the same posterior probability. To this end, we developed an ordinal logistic model with group regularization for fusing adjacent classes. We establishedan estimation algorithm based on the alternating direction method of multipliers, andused Monte Carlo simulations and real data analysis to investigate the usefulness of our
proposed method.

Key words: Adjacent-categories logit model, alternating direction method of multipliers, group lasso, ordinal
categorical data.