GAM
とその周辺
田中祐輔
1∗伊庭克拓
2竹澤邦夫
3辻谷将明
4 1大阪電気通信大学大学院 工学研究科
2東京ＣＲＯ (株) DM 統計解析部
3独立行政法人 農業・食品産業技術総合研究機構 中央農業総合研究センター
4大阪電気通信大学 情報通信工学部
1
はじめに
近年のコンピュータ技術の飛躍的進歩や多様なデー タの収集と蓄積によって，非線形性をもつ統計的多変 量データ解析が脚光を浴びつつある．特に，平滑化ス プライン，カーネルマシン，サポートベクターマシン (SVM:Support Vector Machine)などのマシン・ラーニ ング [6] が有力な武器になりつつある．本稿では，GAM の一種である平滑化スプライン による判別モデルを取 上げ [2, 4, 23]，シミュレーション研究を試みる [8]． 平滑化スプラインによるアプローチの利点は，i) 共 変量の非線形性を摘出し，ii) グラフ化により非線形構 造を視覚的に把握できる．更に，iii) モデルの未知パラ メータの最尤解が一意で，計算に要する時間が短い，な どが挙げられる．しかし，平滑化スプラインは未知パラ メータを推定するだけでは不十分であり，予測精度の高 い平滑化スプラインを得るためには，最適な平滑化パ ラメータを決定しなければならない．局所評点化法を用 いて，最適な平滑化パラメータをグリッド検索すると， ill-conditionに落込み，解が収束しないという欠点があ る．その難点を回避するため，R では Wood[20, 22, 23] のアルゴリズムを採用している．本稿では，(１例消去) クロス・バリデーション (CV) 法による平滑化パラメー タの決定について考察する． 2群判別問題の適用例として，共変量間に交互作用 効果があると考えられるマサバ卵データ (表 1) を解析 する [24]．これは，北大西洋のビスケー湾におけるマ サバ資源に関する 634 例の調査データである． 応答変数には egg(卵の有無) を取上げ，共変量として， lon(経度)，lat(緯度)，b.depth(√海底までの深さ;m)， c.dist(200m等深線からの距離;degree)，temp.surf (海 面温度) を考える．経緯度は調査地点を表す．海の状態 (深さ，水温など) は経緯度に対して逐一変化し，加法的 な変化はせず，凹凸のあるねじれた構造をもつ．すなわ ち，経緯度に 2 因子交互作用があると考えられる．また， 200m等深線からの距離 (c.dist) は調査地点から大陸棚 ∗連絡先： 大阪電気通信大学大学院 工学研究科 〒 572-8530 大阪府寝屋川市初町 18-8 E-mail: so3bs [email protected]端までの距離の代用であり，(緯度)1 度=60 分，(緯度)1 分=1 海里=1852m から，距離は c.dist× 60 × 1852(m) と算出できる．
表 1: マサバ卵データ
No. lon lat b.depth c.dist temp.surf 群
1 -4.65 44.57 65.89 0.84 15.00 0 2 -4.48 44.57 65.83 0.86 15.40 0 . . . . . . . . . . . . . . . . . . . . . 334 -7.27 48.25 13.15 0.37 13.40 1 . . . . . . . . . . . . . . . . . . . . . 633 -3.72 47.25 10.00 1.00 16.40 1 634 -3.25 47.25 8.00 1.26 16.70 0
2
ロジスティック判別
2.1
平滑化スプライン
(1)
局所評点化法
2値応答 yiがベルヌーイ分布 B(1, πi)，すなわち， E[yi] = πi, V ar[yi] = πi(1− πi) に従うとき，n 組のデータを (x(i), yi) , i = 1,· · · , n とする．x(i) = (x1(i),· · · , xI(i)) T は I 個の要素から成 る共変量ベクトルである．この x(i) を用い共変量の データ行列 x を x = (x(1),· · · , x(n))T とした, 平滑化 スプライン ηi ∆ = ln ( πi 1− πi ) =
α + s1(x1(i)) + s2(x2(i)) +· · · + sI(xI(i)),
i = 1, 2, . . . , n (1)
に関する局所評点化法は，次の通りである ([4] の 2.3 節，[5] の 6.3 節)．
人工知能学会研究会資料 SIG-DMSM-A803-07 (3/3)
[ステップ 1(初期化)] ˆ sj(xj(i))≡ 0, j = 1, 2, . . . , I, i = 1, 2, . . . , n ˆ α = ln ( ¯ y 1− ¯y ) , y =¯ n ∑ i=1 yi / n [ステップ 2] ˆ ηi= ˆα, πˆi= 1 1 + exp(−ˆηi) [ステップ 3(反復)] zi= ˆηi+ yi− ˆπi ˆ πi(1− ˆπi) :調整従属変数 wi = ˆπi(1− ˆπi) :重み として，(x(1), z1), (x(2), z2), . . . , (x(n), zn)に後 退当てはめ法を適用し，i = 1, 2, . . . , I について ˆ siを算出する．そして，ˆπiを推定する． [ステップ 4(収束判定)] 逸脱度 Dev(y, ˆπ) =−2 n ∑ i=1 {yiln ˆπi+ (1−yi) ln(1− ˆπi)} (2) が収束するまで [ステップ 3] を繰返す．ここに， y = (y1,· · · , yn) ′ とする． [ステップ 3] で用いられている後退当てはめ法は，I = 2の場合，以下の通りである： [ステップ 1(初期化)] zi= α + s1(x1(i)) + s2(x2(i))に ついて，ˆs[0]2 (x2(i))≡ 0 とする． [ステップ 2] x1(i)に対して，zi−ˆs [0] 2 (x2(i))− ˆα を平滑 化する．すなわち，yiを yi ∆ = zi− ˆs [0] 2 (x2(i))− ˆα とし， “ x1(1), z1− ˆs[0]2 (x2(1))− ˆα ” ,“x1(2), z2− ˆs[0]2 (x2(2))− ˆα ” , . . . ,“x1(n), zn− ˆs[0]2 (x2(n))− ˆα ” について，ˆs[1]1 を算出する． [ステップ 3] x2(i)に対して，zi− ˆs [1] 1 (x1(i))− ˆα を平 滑化し，ˆs[1]2 を算出する． [ステップ 4] x1(i)に対して，zi− ˆs [1] 2 (x2(i))− ˆα を平 滑化し，ˆs[2]1 を算出する． [ステップ 5] ˆs[k]1 (x1(i)), ˆs[k]2 (x2(i))が収束するまで繰 返す． なお，3 次の平滑化スプラインを用いると局所評点 化法は，ペナルティ付き対数尤度 n ∑ i=1 {yiln πi+ (1− yi) ln(1− πi)} −n 2 I ∑ j=1 λj { s′′j(t)}2dt (3) の最大化と等価である．SAS や S-Plus では，局所評点 化法が採用されている．
(2) Wood
法
Wood[20, 22, 23]は，最適な平滑化パラメータ λ を 求めるアルゴリズムを提案している．B-スプライン関 数 f1(x1i), f2(x2i), . . .に基づくロジスティック・モデル ηi ∆ = ln ( πi 1− πi ) = X β n×m m×1 (4) を考える．例えば，I=2 なら s1(x) = q1 ∑ j=1 βjb1j(x), s2(x) = q2 ∑ j=1 βq1+jb2j(x) β = (β1, . . . , βq1, βq1+1, βq1+2, . . . , βq1+q2) X = b11(x11)· · · b1,q1(x1n) b2,q1+1(x21)· · · b2,q1+q2(x2n) b11(x12)· · · b1,q1(x1n) b2,q1+1(x22)· · · b2,q1+q2(x2n) • · · · • • · · · • b11(x1n)· · · b1,q1(x1n) b2,q1+1(x22)· · · b2,q1+q2(x2n) となる．ここに，q1, q2 は節点の個数である。対数尤 度は l(y, X, β) = n ∑ i=1 {yiln πi+ (1− yi) ln(1− πi)} (5) で与えられる．ここに， πi= 1 1 + exp(−Xiβ) とする．ペナルティ付き対数尤度は lp= l(y, X, β)− 1 2 m ∑ i=1 λiβTSiβ (6) となる ([23] の p.168)．ただし，Siは (非負定) ペナル ティ行列である．(6) 式のペナルティ付き対数尤度を最 大にする β をスコア法で求めることができる [20, 21]． (6)式は 最小化 1 n n ∑ i=1 { (yi− πi)2 V ar[yi] } + m ∑ i=1 λiβTSiβ (7) と同等である ([23] の 4.3 節)1．スコア法の第 k 反復に おいて，平滑化スプラインを用いた GAM の場合， { 重み : W[k]= Cov[y]の逆行列 調整従属変数 : z = Xβ + W[k]−1(y− π) 1(6)式より−2l p=−2Pni=1li(ti,Xi,˛) +Pmi=1λi˛TSi˛ となる．よって，(6) 式の最大化は，(2) 式より，−2lp= Dev(y, ˆı)+ Pm i=1λi˛TSi˛ の最小化と同等になる．とし，付録 A の (A.8) 式にペナルティ項を付けた 最小化 °°°√W[k](z− Xβ)°°°2 +∑Ii=1λiβTSiβ w.r.t. β 制約式 Cβ = 0 (8) を解けばよい ([23] の p.185)．このとき，制約式の係 数 C は，∑mj=1βj= 0, ∑m j=1βjxj= 0を意味する．こ の制約式に QR 分解 CT = U ( P 0 ) を施し，U を U ≡ (D : B) に分解する．β = BβBとおくと Cβ = ( PT0) ( DT BT ) BβB= 0を満たす ([23] の 1.8.1 節)． (8)式は ˜z =√W[k]z, ˜X =√W[k]XB, ˜S i= BTSiB とおくと 最小化 °°°˜z − ˜XβB°°° 2 +∑λiβBTS˜iβB w.r.t. βB となり，ハット行列は ˜ Hλ= ˜X ( ˜ XTX +˜ ∑λiS˜i )−1 ˜ XT (9) と書ける [20]． ˜y =√W[k]yとおくと，GCV 規準は Vg= n∥˜y − ˜Hλy˜∥2 { n− tr( ˜Hλ) }2 (10) で与えられる．y の予測値 ˆyは ˆ y =√W[k]−1H˜λy˜ (11) から推定すればよい．なお，R では (10) 式の残差平方和 を逸脱度で置換えた GCV 逸脱度 (Generalized Cross-Validated deviance) Vg= nDev(y, ˆπ) { n− tr( ˜Hλ) }2 (12) が採用されている ([5] の p.159，[23] の p.178，[13] の p.376)． ペナルティ付き対数尤度を最大化する λ = (λ1, λ2, . . . , λm)を求める際，(9) 式の ˜Hλに未知パラメータ λ が 含まれており，(9) 式が最小になる λ をニュートン・ラフ ソン法で求める ([23] の 4.6 節)．ς = (ς1, ς2, . . . , ςm)≡ (ln λ1, ln λ2, . . . , ln λm)とおくと，更新式 ς[k]= ς[k−1]− ( ∂2Vg ∂ςi∂ςj )−1 ς=ς[k−1] ( ∂Vg ∂ςi ) ς=ς[k−1] (13) を得る．∂Vg ∂ςi および ∂ 2V g ∂ςi∂ςj の導出については，[15] に 詳しい. 次に，予測変数が 2 個の場合につて，B-スプライン 関数の有意性検定について述べる． ηi ∆ = ln ( πi 1− πi ) = s1(x1i) + s2(x2i) = X β n×q q×1 (14) において，β = (β1; β2) = (β1, β2,· · · , βq1; βq1+1,· · · , βq1+q2)の例えば，最初の予測変数に関する 帰無仮説 H0: β1= 0(i.e., β1= β2=· · · = βq1 = 0) を検定しよう．ベルヌーイ分布の場合，overdispersion は発生しないから，逸脱度に基づく尤度比検定を推奨 する (overdispersion については，付録 B を参照せよ)． (14)式のモデルに関する逸脱度を Dev(1)，モデルの残 差自由度を ν1，帰無仮説のもとでのそれらを Dev(0) および ν0 とすると Dev(0)− Dev(1) ∼ χ2ν0−ν1 (15) となる ([2] の p.30). なお，R では，(14)，(15) 式に基 づく F 検定 (Wald 検定) が採用されている2．
2.2
薄板平滑化スプライン
n組のデータ (x(i), yi) , i = 1,· · · , n について，共変 量間に交互作用がある場合 ln ( πi 1− πi ) =∑ α<β sαβ(xα(i), xβ(i)), i = 1, 2,· · · , n (16) を考える [3, 18, 19]．ここで，(16) 式の右辺を s(x)≡∑ α<β sαβ(xα(i), xβ(i)), i = 1, 2,· · · , n (17) とおく．sαβ(xα(i), xβ(i))は薄い弾性体の板を曲げた 形に似ていることに由来し，薄板スプライン(thin plate spline)と呼ばれている． 2値応答ベクトル y = (y1, y2,· · · , yn) T について，ペ ナルティー付き対数尤度 n ∑ i=1 { yis(x(i))− ln(1 + es(x(i))) } −n 2λJmI(η) (18) が最大になる s を推定する．ここに，∑ni=1{yis(x(i))− ln(1 + es(x(i)))} は対数尤度である3．J mI(η)は揺動 (wiggliness)を測るペナルティー関数，λ は平滑化パラ 2Rのマニュアルはウェヴサイト http://gd.tuwien.ac.at/languages/R/doc/packages/mgcv.pdf から得られる 3(3)式の対数尤度はPn i=1{yiln πi+ (1− yi) ln(1− πi)} = Pn i=1{yiln “ π i 1−πi ” + ln(1− πi)} =Pni=1{yis(x(i))− ln(1 + πi 1−πi)} と書ける．メータである．揺動のペナルティーは JmI = ∫ · · ·∫ℜI ∑ ν1+···+νI=m m! ν1!···νI! ( ∂ms ∂xν11 ···∂xνII )2 dx1· · · dxI (19) と定義する．揺動のペナルティに関する偏微分の階数 mは 2m > I を満たす整数値である．このとき，(18) 式が最大になる関数は s (x) = n ∑ i=1 ciηmI(∥x − x(i)∥) + M ∑ j=1 djϕj(x) (20) と表せる．ここに，c と d は T′c = 0, Tij = ϕj(x(i)) を満たす未知パラメータで，M = ( m + I− 1 I ) 個 の関数 ϕi(x)は線形独立な多項式である．更に ηmI(r) = (−1)m+1+I/2 22m−1πI/2(m− 1)!(m − I/2)!r 2m−Ilog(r) r : 偶数 Γ(I/2− m) 22mπI/2(m− 1)!r 2m−I r : 奇数 (21) となる ([3] の 4.4 節，[23] の 4.1.5 節)． 例えば、m = I = 2 の場合，(19) 式の JmIは， J22= ∫+∞ −∞ ∫+∞ −∞ {( ∂2s(x1,x2) ∂x2 1 )2 + 2 ( ∂2s(x1,x2) ∂x1∂x2 )2 + ( ∂2s(x1,x2) ∂x2 2 )2} dx1dx2 (22) となり，(20) 式において，T′c = 0より，ϕ1(x) = 1, ϕ2(x) = x1, ϕ3(x) = x2を得，(18) 式を最大にす る解は η (x1, x2) = d0+ d1x1+ d2x2 + n ∑ i=1 ci· Gi(x1− x1(i), x2− x2(i)) (23) と表せる． (21)式より Gi(x1− x1(i), x2− x2(i)) = {t2 i 8π· log(ti) , t > 0 0 , t = 0 (24) である. ただし， ti= √ (x1− x1(i))2+ (x2− x2(i))2 (25) である．データ数が n のとき，(23) 式には，n + 3 個 の未知パラメータが含まれている． 3個の不足するデータは 3 つの条件 n ∑ i=1 ci= 0, n ∑ i=1 cix1(i) = 0, n ∑ i=1 cix2(i) = 0 (26) に対応する ([13] の 4.3 節)． ちなみに，m = 2, I = 1 の場合，(19) 式の JmIは， J21= ∫ +∞ −∞ ( ∂2s (x) ∂x2 )2 dx (27) を得る．(20) 式において，TTc = 0より ϕ1(x) = 1, ϕ2(x) = xを得，(18) 式を最大にする解は s (x) = d0+ d1x + n ∑ i=1 ci|x − xi|3 (28) と表せる．データ数が n のとき，(28) 式には，n + 2 個の未知パラメータが含まれている．2 個の不足する データは 2 つの条件 n ∑ i=1 ci= 0, n ∑ i=1 cixi= 0 (29) に対応する [12].
2.3
リサンプリング法
グループ化されていない 2 値データ (すなわち，ベル ヌーイ分布) の場合，(2) 式の逸脱度に対する漸近カイ二 乗性は成立たない ([1] の 3.8.3 節，[10]，[14],[15],[16])． そこで，ブートストラップ法 [9] を援用し，(2) 式の棄 却点を算出する: [ステップ 1] 初期標本 X = (X<1>, X<2>, . . . , X<n>) からリサンプリング（ペア・サンプリング）によ り，B 個のブートストラップ標本 X∗= (X<1>∗, X<2>∗, . . . , X<n>∗) を生成する． [ステップ 2] b(= 1, . . . , B) 番目のブートストラップ標 本を Xb∗とし，逸脱度 Dev(b) =−2 ln L ( Xb∗;β(Xˆb∗) ) (30) を計算する．ここに， ˆβ(Xb∗)は b 番目のブート ストラップ標本 Xb∗から推定される回帰係数で， その推定値から算出される Xb∗の対数尤度が ln L ( Xb∗;β(Xˆb∗) ) である．[ステップ 3] Dev≥ Dev∗なら，有意水準 α でモデル は妥当でないとみなす．ここに，Dev は初期標本 に対する (2) 式の逸脱度で， Dev∗= Dev(b)を小さい順に並べたときの 第 j 番目の値，α = 1− j/(B + 1) とする．
2.4
適用例
2群判別問題の実際例として，5 個の共変量から成る 表 1 のマサバ卵データを取上げる．このデータについ て，交互作用スプラインモデル ln ( p 1− p )= s(lon, lat) + s(b.depth)
+s(c.dist) + s(temp.surf ) を当てはめる．すなわち，lon と lat とに薄板平滑化ス プラインを用いる． このモデルについて，(１例消去)CV 法による厳密解 法を採用し、自由度 (平滑化パラメータ) の最適選択を 行う．Wood 法 (GCV 規準) との比較結果を表 2 に示 す．なお，厳密解法で最適化を行った結果より GCV を，Wood 法で最適化を行った結果より CV をそれぞ れ再計算している．同表から，Wood 法による結果で は，厳密解法に比べて GCV が小さくなり，逸脱度も かなり小さくなる．よって，データへの当てはまりが 良いといえる．しかし，厳密解法に比べて CV が大き い．図 1 に s(c.dist) の推定値を示す．図 1(b) は Wood 法による結果であり，図 1(a) の厳密解法と比べて，推 定値曲線の凹凸が多く，GCV 規準による最適 df 値の 選択ではデータへの過剰当てはめの傾向が見られる． 次に，CV 法による厳密解法の結果について考察を 加える．図 2 は，b.depth に対する平滑化スプラインの 推定値 (実線) である．実線の周りの領域は，推定値か らその標準偏差の 2 倍離れた位置を示している．同図 では，b.depth に対する非線形性が顕著に現れており， マサバが生息する最適な海域が，海底 900m 程度の海 域であることが分かる． 図 3(a) に経緯度に対する薄板平滑化スプラインの推 定値 (曲面) を示す．また，図 3(b) にその等高線図を示 す．図 3(a) より，薄板平滑化スプラインは共変量が 2 (a) CV法による厳密解法 (b) Wood法 図 1: s(c.dist) の推定値 !" #$ % &' ( )* ' 図 2: s(b.depth) の推定値 個の場合，滑らかな曲面をもってデータ構造を表現す る．次に，図 3(b) は，白い部分に近づくほどマサバ卵 が捕れやすく，黒い部分に近づくほど捕れにくいこと を表す．ここから，フランス西岸のビスケー湾ではマ サバがあまり生息しておらず，そこから北西の沖合に 多く生息している傾向が見られる． 平滑化スプラインを適用した効果を調べるため，各 共変量に関する尤度比検定を行なう． • s(lon, lat) Dev(0)− Dev(1) = 521.23 − 438.68 = 82.54 df = 625.59− 614.51 = 11.08 82.54 > 24.85 = χ2 11.08(0.01) 表 2: 最適 df 値の比較
s(lon, lat) s(b.depth) s(c.dist) s(temp.surf ) CV GCV 逸脱度 CV法による厳密解法 12.22 2.69 1.89 1.69 0.76 0.74 438.68 Wood法 (GCV 規準) 20.74 2.40 8.04 3.86 0.83 0.72 406.54
(a) 3次元プロット ! " # $% (b)等高線図 図 3: s(lon, lat) の推定値 • s(b.depth) Dev(0)− Dev(1) = 447.08 − 438.68 = 8.40 df = 616.72− 614.51 = 2.21 8.40 > 6.38 = χ2 2.21(0.05) • s(c.dist) Dev(0)− Dev(1) = 446.83 − 438.68 = 8.15 df = 615.95− 614.51 = 1.44 8.15 > 7.85 = χ21.44(0.01) • s(temp.surf) Dev(0)− Dev(1) = 442.66 − 438.68 = 3.98 df = 615.91− 614.51 = 1.40 3.98 > 3.51 = χ2 1.40(0.10) より，各共変量は有意である． これらの結果より，経度 (-15∼-10 度)，緯度 (48∼55 度) の座標において，海底までの深さ 900m 程度の海域 で産卵される可能性が高い，つまり，マサバの生息範 囲がこの辺りに集中していることが分かる． このように平滑化スプラインを適用することによっ て，母集団構造の非線形性を摘出し，その非線形構造 を視覚的に把握するグラフ表現が可能になる． ペア・サンプリングに基づく逸脱度のブートストラッ ピングに対するヒストグラムを図 4 に与えておく．同 図より，Dev∗ = 454.70 > 438.68 = Devを得，モデ ルは妥当とみなす. 次に一つの観測値を除去することによるモデル適合度 への影響を調べるため逸脱度に関する DIFDEV(Differe nce of Deviance)
∆Dev[d]= Dev− Dev[d] (31)
! " # $ 図 4: ブートストラップされた逸脱度のヒストグラム を用いる ([7] の 5.3 節)．ここに，Dev と Dev[d]は，そ れぞれ初期標本および d 番目の個体を除去したときの 逸脱度である．Dev の自由度 df = n− p, Dev[d] の 自由度 df = n− p − 1 より，∆Dev[d]の自由度 df = n− p − (n − p − 1) = 1 となる．よって，∆Dev[d]は 自由度 1 のカイ二乗分布に従う．図 5 に DIFDEV をプ ロットしておく．同図から，有意水準 1% で No.601 が 影響観測値であることが分かる．この観測値を除去後 の逸脱度は 429.92 と減少する．
図 5: DIFDEV
3
シミュレーション
3.1
シミュレーション設計
本稿のシミュレーションでは，各観測値の共変量を x = (x1, x2)Tの 2 個とし，それぞれ [−1, +1] 区間の一 様乱数として生成する．そして，生成した x = (x1, x2)T を p = 1 1 + exp{−f(x1, x2)} (32) に代入して，各観測値を 2 つの群に分ける．f (x1, x2) は任意の非線形関数を用いることができる．f (x1, x2) として，Vach モデル [17] を変型した f (x1, x2) = sin(2πx1) + αx1x2+ sin(2πx2) (33) を採用する．(33) 式の α は，第 2 項を共変量間の積に よる交互作用項と見なし，交互作用効果の大きさを調 整する．ここでは，α = 1, 10 とする．p は 0≤ p ≤ 1 と なり，p の 3 次元プロット (x1, x2, p)は図 6,7 のように なる．同図より，交互作用効果が大きくなるほど，デー タ構造にはねじれた形状が顕著になる． 各観測値は，ベルヌーイ乱数 y∼ Bin(1, p), y = { 1(第 1 群) 0(第 2 群) (34) に基づいて 2 つの群に分割される．よって，シミュレー ション・データは共変量 x = (x1, x2)Tと 2 値応答 y と をもつ観測値から成る． 本稿ではモデルへのサンプル・サイズの影響を調べ るため，サイズを 100, 200, 400, 600, 800, 1000 例と変化 させた訓練標本を生成する．次に，構築したモデルの 未知の標本に対する判別精度を評価するため，同サイ ズの検証標本を生成する． 図 6: 3 次元プロット (α = 1) 図 7: 3 次元プロット (α = 10) 図 8 に，α = 10 における 1000 例のシミュレーショ ン・データを示す．図 8(a) は，図 7 の 3 次元プロット を p に関する等高線図に置き直したものである．同図 の白の部分は確率 p が高い部分を，黒の部分は低い部 分をそれぞれ表す．散布図は図 8(b) のように表わされ， 図 8(a) の等高線図の白の部分，すなわち，確率 p が高 い右上と左下部分に第 1 群に属する個体が集まる．ま た，灰色の確率 p が 0.5 付近では，群がばらつく傾向 が見られる．3.2
性能評価
生成したシミュレーション・データを用いて，以下 の 2 つのモデルの非線形性をもつデータへの有効性を 比較検討する． • 主効果スプラインモデル ln ( p 1−p ) = s(x1) + s(x2) • 交互作用 (薄板) スプラインモデル ln ( p 1−p ) = s(x1, x2) 他のマシン・ラーニングの手法として，ニューラル ネットワーク，および従来手法として，線形判別モデ(a)等高線図 (b)散布図 図 8: シミュレーション・データ (α = 10) ルとロジスティック判別モデルを取上げる． シミュレーション・データは，1000 例の訓練標本と検 証標本とを 1 組とし，各標本に対して，直前のサイズを 含む形で 100 例から 1000 例までの 6 サイズを設定する． そして，交互作用効果を調整した α = 1, 10 について， 個体を変えた各 100 組のデータを生成する．α = 1, 10 の 2 パターンについて，それぞれ 100 試行のシミュレー ションを行った結果について考察する．なお，各試行 における最適な平滑化パラメータの決定には Wood 法 を用いる．また，ニューラルネットワークの最適な隠 れユニット数は，すべての試行において，100,200 例で は 5 個，400,600,800,1000 例では 6 個とする．
(1)
α = 1
α = 1において，図 9 にサンプル・サイズを変えた ときの，種々のモデルの下での誤判別率を示す．各サ イズについて，図の左から主効果スプライン，交互作 用スプライン，ニューラルネットワーク，線形判別，ロ ジスティック判別の結果である． 図 9(a) はモデルの構築に用いたデータ (訓練標本) に 対する判別結果であり，図 9(b) は未知のデータ (検証 標本) の結果である．また，箱ひげ図のひげの長さは， 箱の長さの 1.5 倍とし，その上端・下端はひげの長さ 以内の最大値・最小値である． 図 9(b) より，検証標本において，主効果スプライン， 交互作用スプラインの誤判別率がもっとも小さく，将 来のデータに対する判別精度がもっとも良い．また， ニューラルネットワークにおいても誤判別率は小さく なり，これらのモデルはデータの非線形構造を近似で きている．それに対して，従来手法の線形判別，ロジ スティック判別は，全体的に他のモデルと比べて判別 精度に劣る．よって，非線形性をもつデータに対して， 従来手法よりも有効であることが確認できる． 次に，サンプル・サイズについて考察を加える．主 効果スプライン，交互作用スプライン，ニューラルネッ トワークでは，100,200 例において，図 9(a) の訓練標 本の誤判別率が他のサイズと比べて小さい．それに対 して，図 9(b) の検証標本では誤判別率が大きくなる． よって，サイズが小さい場合には，モデル構築に利用す るデータへの過剰当てはめの傾向が見られる．しかし， このときスプラインモデルにおいて，逆に非常に強い (a)訓練標本 (b)検証標本 図 9: 誤判別率に関する箱ひげ図 (α = 1)平滑化，すなわち，平滑化パラメータに大きな値が選 ばれるケースも見られる．その際の判別精度は，線形 判別，ロジスティック判別の結果と類似する．これらの 傾向は，100,200 例における訓練標本の誤判別率のばら つきが大きいことに表れており，サイズが小さい場合 には，モデル構築がデータによる影響を強く受ける． 400例では，検証標本の誤判別率が減少する．その 後，サイズが大きくなるに従い，緩やかに減少し，判別 精度の向上が確認できる．また，誤判別率のばらつき も小さくなり，安定したモデルが得られている．更に， 訓練標本の誤判別率は，サイズが大きくなるに従い，検 証標本の結果へと近づく．よって，サンプル・サイズ が十分に大きい場合には，モデル構築に用いるデータ から，構築したモデルの将来のデータに対する判別精 度を評価できる．
(2)
α = 10
α = 10において，図 10 にサンプル・サイズを変えた ときの，種々のモデルの下での誤判別率を示す．同図 より，2 因子交互作用効果が大きくなると，主効果スプ ラインの誤判別率は大きくなる．その判別精度は，線 形判別，ロジスティック判別と類似している．一方，交 互作用スプライン，ニューラルネットワークの判別精 度は良い．この結果より，この 2 つのモデルは，交互 作用効果を考慮したモデルを構築できることが分かる． また，サンプル・サイズによる傾向は，α = 1 の場 合と同様である．3.3
モデル評価
構築したモデルの真のモデルへの近似精度を比較検 討する．ここでは，確率 p と推定値 ˆpとの平均二乗誤 差に基づく比較を行う．また，各モデルについて，推 定値 ˆpに関する 3 次元プロットを取上げ，視覚的に考 察する．評価対象として，α = 1, 10 から，それぞれ 1 組のシミュレーション・データを用いて構築したモデ ルを取上げる．(1)
α = 1
α = 1について，図 11 はサンプル・サイズに対す る種々のモデルの平均二乗誤差を示す．また，図 13 に 100,400,1000例における推定値 ˆpに関する 3 次元プロッ トを示す．図 11 より，主効果スプライン，交互作用ス プライン，ニューラルネットワークは，サイズが大き くなると共に真のモデルとの誤差が小さくなる．ここ では主効果スプラインがもっとも近似精度が良いとい える．それに対して，線形判別，ロジスティック判別 は，誤差がサイズが大きくなると共に緩やかに大きく なる． また，図 13 より，図 6 の真のモデルと比べて，主効 果スプライン，交互作用スプライン，ニューラルネッ トワークは，100 例ではデータのもつ非線形な構造を 表現できていないが，400,1000 例では非線形構造を的 確に表現している．一方，線形判別，ロジスティック判 別は，真のモデルを近似できておらず，サイズに対し てモデルが大きく変化する様態が見られない． (a)訓練標本 (b)検証標本 図 10: 誤判別率に関する箱ひげ図 (α = 10)図 11: サンプル・サイズに対する平均二乗誤差 (α = 1)
(2)
α = 10
α = 10において，図 12 にサンプル・サイズに対す る種々のモデルの平均二乗誤差を示す．また，図 14 に 100,400,1000例における推定値 ˆpに関する 3 次元プロッ トを示す． 図 12 より，2 因子交互作用効果が大きくなるため， 主効果スプラインでは，真のモデルとの誤差が大きく なり，その結果は線形判別，ロジスティック判別と類似 する．一方，交互作用スプライン，ニューラルネット ワークは，真のモデルとの誤差が小さく，交互作用を 含むデータ構造を近似している． また，図 14 より，図 7 の真のモデルと比べて，主 効果スプラインは交互作用を含むデータ構造を表現で きず，平面に近いモデルとなる．交互作用スプライン， ニューラルネットワークは，このようなデータ構造を 近似している．2 変量間の積による交互作用効果は，ね じれた構造をもって表現される．また，100 例におい て推定値 ˆpが 0 または 1 に近い値となり，過学習の傾 向が見られる． 図 12: サンプル・サイズに対する平均二乗誤差 (α = 10)4
結び
本稿では，2 群判別問題を例とし，一般化加法モデル (GAM)の一種である平滑化スプラインについて，自由 度 (平滑化パラメータ) の最適選択とシミュレーション を行った． 平滑化スプラインは最適な自由度の決定が必要とな るが，データのもつ非線形な構造を抽出し，グラフ化に より視覚的に把握することができる．この非線形性の 可視化は，他のマシン・ラーニングの手法であるニュー ラルネットや SVM などでは不可能であり，平滑化ス プラインを用いる利点である． シミュレーションにおいては，シミュレーション・デー タを用いて，非線形性をもつデータに対して，従来の 解析手法よりも有効であることを示した．また，サン プル・サイズについて，サイズが大きくなると共に未 知のデータに対する判別精度は向上する．しかし，サ イズが小さい場合には安定したモデルが得られないこ とに留意する必要がある．更に，薄板平滑化スプライ ンおよびニューラルネットワークを適用することによ る 2 因子交互作用のモデル化が有効であることを確認 した．(a)主効果スプライン (100 例) (b)主効果スプライン (400 例) (c)主効果スプライン (1000 例) (d)交互作用スプライン (100 例) (e)交互作用スプライン (400 例) (f)交互作用スプライン (1000 例) (g)ニューラルネット (100 例) (h)ニューラルネット (400 例) (i)ニューラルネット (1000 例) (j)線形判別 (100 例) (k)線形判別 (400 例) (l)線形判別 (1000 例) (m)ロジスティック判別 (100 例) (n)ロジスティック判別 (400 例) (o)ロジスティック判別 (1000 例) 図 13: 推定値 ˆpに関する 3D プロット (α = 1)
(a)主効果スプライン (100 例) (b)主効果スプライン (400 例) (c)主効果スプライン (1000 例) (d)交互作用スプライン (100 例) (e)交互作用スプライン (400 例) (f)交互作用スプライン (1000 例) (g)ニューラルネット (100 例) (h)ニューラルネット (400 例) (i)ニューラルネット (1000 例) (j)線形判別 (100 例) (k)線形判別 (400 例) (l)線形判別 (1000 例) (m)ロジスティック判別 (100 例) (n)ロジスティック判別 (400 例) (o)ロジスティック判別 (1000 例) 図 14: 推定値 ˆpに関する 3D プロット (α = 10)
謝辞
本稿をまとめるに当たり，日本科学技術連盟 DE・O 部会におきましては，大変貴重なご意見を頂きました． ここに記して厚く御礼申し上げます．
参考文献
[1] Collett,D.(2003): Modelling Binary Data, 2nd ed., Chapman and Hall, London.
[2] Faraway,J.J.(2006) Extending the Linear Model
with R Generalized linear,Mixed Effects and Nonparametric Regression Models, New York,
Chapman & Hall.
[3] Gu, C.(2002): Smoothing Spline ANOVA
Mod-els, Springer, New York.
[4] Hastie, T.J. and Tibshirani, R.J. (1987): Gen-eralized additive models: Some applications, J.
Amer. Stat. Assoc. 82, 371-386
[5] Hastie, T.J. and Tibshirani, R.J.(1990):
Gener-alized Additive Models, Chapman and Hall,
Lon-don.
[6] Hastie, T., Tibshirani, R. and Freidman, J. (2001): The Elements of Statistical Learn-ing: Data Mining, Inference, and Prediction, Springer, New York.
[7] Hosmer,D.W. and Lemshow,S. (2000):Applied Logistic Regression, John Wiley, New York. [8] 伊庭克拓，辻谷将明 (2006): ニューロ判別モデル
に関するシミュレーション研究，大阪電気通信大 学研究論集 (自然科学編), 41, 13-39
[9] Ishiguro, M., Sakamoto, Y. and Kitagawa, G. (1997): Bootstrapping log likelihood and EIC, An extension of AIC. Ann. Inst. Stat. Math. 49, 411-434.
[10] Landwehr, J.M, Pregibon,D. and Shoe-maker,A.C.(1984): Graphical methods for assessing logistic regression models. J. Amer.
Stat. Assoc. 79, 61-71.
[11] Ruppert,D., Wand, M.P., and Caroll,R.J.(2003): Semiparametric Regression, Cambridge Univ. Press,Cambridge.
[12] 桜井明 (1993):C によるスプライン関数，東京電機 大学出版局．
[13] Takezawa, K.(2006): Introduction to Nonpara-metric Regression, John Wiley.
[14] Tsujitani, M. and Koshimizu,T.(2000)：Neural Discriminant Analysis. IEEE Trans. Neural
Net-works 11, 1394-1401.
[15] 辻谷将明，外山信夫 (2007):R による GAM 入門, 行動計量学, 34, 111-131.
[16] M.Tsujitani and M.Sakon(2009): Analysis of survival data having time-dependent covari-ates.IEEE transactions on Neural Networks, in Press. ”
[17] Vach, W., Roβner, R., and Schumacher, M.(1996): Neural betworks and logistic regres-sion:Part II, Comput.Statist. & Data Analy., 21, 683-701.
[18] Wang, Y., Wahba, G., Gu, C., and Klein, R., and Klein, B.(1997): Using smoothing spline ANOVA to examine the relation of risk factors to the in-cidence and progression of diabetic retinopathy,
Statist. Medicine, 16, 1357-1376.
[19] Wang, Y.(1997): GRKPAK: Fitting smoothing spline ANOVA models for exponential families,
Commun. Statist.-Simulat. Comput., 26,
765-782.
[20] Wood,S.N. (2000): Modelling and smoothing parameter estimation with multiple quadratic penalties,J. Roy. Statist. Soc. B, 62, 413-428. [21] Wood,S.N. and Nicole H.Augustin(2002): GAMs
with integrated model slection using pe-nalized regression splines and applications to environmental modelling, Ecological
Mod-elling, 157, 157-177.
[22] Wood,S.N.(2004):Stable and efficient multiple smoothing parameter estimation for generalized additive models, J.Amer.Stat.Asoc.,99,673-686. [23] Wood,S.N.(2006) Generalized Additive Models
An Introduction with R, New York, Chapman &
Hall.
[24] Wood,S.N.(2008):Fast stable direct fitting and smoothness selection for generalized additive models, J.Roy. Statsit.Soc. B,70,495-518.
