統計コンサルの議事メモ

統計や機械学習の話題を中心に、思うがままに

統計数理研究所公開講座「統計の哲学を理解するために」参加記録

1/31に開催された統計数理研究所公開講座「統計の哲学を理解するために」に参加してきましたのでそのメモを共有しておきます。全体的にはエリオット・ソーバーという尤度主義者から見た頻度主義・ベイズ主義に対する批判的観点の紹介という構成で、それぞれの立場が答えようとしている問いが浮き彫りになるような内容でした。

なおこの講義は「科学と証拠」という書籍に基づいています。とても面白い本です。

科学と証拠―統計の哲学 入門―

科学と証拠―統計の哲学 入門―

また最後に、あまり時間がなくて駆け足での紹介となっていましたが、Deborah Mayoによる「誤り統計学」について簡単に紹介がありました(松王先生曰く世界初の資料、とのこと)。

Statistical Inference as Severe Testing: How to Get Beyond the Statistics Wars

Statistical Inference as Severe Testing: How to Get Beyond the Statistics Wars


日時

2020/1/31 10:00@立川

講師

北海道大学 松王政浩

エリオット・ソーバー

  • アメリカの代表的な科学哲学者
  • 尤度主義 + AICを支持
  • 著書
    • 科学と証拠
    • OCKHAM'S RAZORSの翻訳がでる予定

統計学論争は終わっていない

  • Statistical Inference as Severe Testing
    • Deborah Mayo
      • 頻度主義
  • 有意検定論争、レシピ的な統計学の食い止め

統計の哲学

  • 統計の基礎をめぐる議論の総体
  • 哲学は常に論争で成り立ち、結論には到達しない
    • ベイズ主義
      • 信念度合い
        • 主観
        • 客観
    • 頻度主義
      • ネイマン-ピアソン
        • 無限回の施行を前提とした考えた
      • フィッシャー
        • 相対頻度、推測確率
  • 統一的見解を示したものはないので各主義を学ぶしかない
    • ベイズ
      • 実態が捉えにくい
      • 意思決定中心、科学的仮説の確証中心(科学哲学)、ベイズ統計の流儀
      • バーガー&ウォルパート(1988)、ゲルマン(一連)
    • 頻度主義
      • Mayoのものは包括的だがわかりにくい
    • 尤度主義
      • コンパクトでわかりやすい
      • ソーバーやロイヤルは話が短く論点が明快だが排他主義

3つの主義

  • ロイヤルの3つの問い
    • 証拠をもとに何がわかるか → 尤度主義
      • 静的現在主義、外在主義
      • データが直接示している「仮説に関する情報(証拠)」をキャッチせよ
      • いったんデータが得られたら他のデータの可能性は捨象せよ
    • 証拠をもとに何を信じるか → ベイズ主義
      • 動的現在主義、内在主義
      • 他のデータの可能性は一切捨象せよ
      • データが「(可能な)仮説に関する情報」をどう変化させるかキャッチせよ
    • 証拠をもとに何をするべきか → 頻度主義
      • 反事実主義、規約主義
      • この場で生じていることだけでなく、生じた可能性がある事柄もすべて併せて仮説について判断せよ
      • 低い確率を棄却のサインとみなそう、同じルールを何度も適用すれば誤りの可能性が低いのだから今回もそのルールを適用する

尤度「原理」について

  • ソーバーの議論の核となる2つの原理
    • 穏当な原理
      • もしEが真であると知ることによって命題Pを棄却することが正当化され、かつ、この情報を得てはじめてPの棄却が正当化されたのであれば、EはPに反する証拠とされねばならない
    • 全証拠の原理
      • 実験によって得られたデータはすべて証拠の判断に用いられなければならない
  • 尤度の法則とは異なる
    • 尤度主義、ベイズ主義の両方に共通する源泉
  • 尤度原理(LP)
    • xが観測されたあと、θについて推論(決定)する際に、実験について関係のあるすべての情報は、観測されたxに対する尤度関数に含まれている。さらに2つの尤度関数がθの関数として互いに比例の関係にあるなら2つの関数はθについて同じ情報を含んでいる
    • フィッシャーに由来
    • これを満たすか満たさないかで、陣営が分かれる
      • 満たす
      • 満たさない
        • 頻度主義
    • いまだ論争の火種
    • もしLPがパラメータ推定の根本原理であるなら、
      • P(x' | θi) = kP(x | θi) が成り立つときEv(E, x) = Ev(E', x')となるはず。Evは証拠、Eは実験
      • これが成り立たなければそのような結果をもたらす推論は不適切
      • バーンバウムが証明
        • 証明の是非は未決着
    • さらに遡ると
      • 十分性の原理
      • 条件付けの原理
      • これら2つの原理は統計学者なら誰でも受け入れられるはずの原理
      • これら2つの原理と尤度原理の等価性を「証明」(バーンバウム)、2つの原理を受け入れるなら尤度原理も受け入れないといけない
    • バーンバウムの問題意識
      • 統計学的に導ける証拠」が「実験における証拠」になっている
        • (E, x) と Ev(E, x) は区別される
          • (E, x)
            • パラメータ空間Ωについての記述
            • Eの可能な結果xのサンプル空間についての記述
          • Ev(E, x):実験的証拠
            • どう評価するか?
      • 解明のポイント
        • 2つの統計的証拠(E, x) と (E', y) が関係するあらゆる点で等しいと言える条件?
          • 統計的証拠(E, x) と (E', y) が関係するあらゆる重要な点で等しい時、 Ev(E, x) = Ev(E', y)
  • 尤度原理からの重要な帰結
    • サンプルスペース(可能だが実際には得られなかった確率変数の値)の無関係性
      • 重要な争点
      • 反事実の重要性を否定

ベイズ主義

  • デフィネッティ
    • 確率、賭け、信念の関係
      • 「合理的な賭けが成立するための条件」「確率の3つの規則」同じ条件
    • D.ギリース「確率の哲学理論」
  • サヴェッジ
    • Inductive inference と Inductive behavior
    • 後者がより重要
      • Inference
        • 意見を変えること
      • Behavior
        • 分布と最終的な行為の経済的事実を用いて最も期待効用の高いものを選ぶ
  • 主観的 vs 客観的
    • 信念変化の合理性
    • 基礎付け主義
    • 客観的事前確率
  • 哲学的ベイズ主義
  • 実用への転換
    • Lindley & Smithによる階層ベイズ
    • Gelfand & SmithによるMCMCの導入

尤度主義とは

  • グレムリン仮説
    • 哲学ではしばしば「説明可能性」と「確率」が結びつく
  • 尤度主義の限界
    • 区間を持った仮説(複合仮説)は考察の対象とならない
    • 尤度の平均を求めることは可能だが、事前分布が必要となる
    • ソーバーの立場
      • 頻度主義の統計テストは用いたくない
      • ベイズ主義はできるだけ控えたい
      • 結果
  • 尤度もモデルに依存するため恣意性が混じるが、なぜ尤度主義者はベイズ主義を批判できるのか?
    • ミスリーディング確率というものを用いることで客観性を保証できる、というのがソーバーの立場
  • ソーバーのAICの法則、どれぐらい使われているのか?
    • あまり

対頻度主義

  • フィッシャーの有意性検定
    • 単独の仮説に対する検定
    • 帰無仮説は正しいと証明されることはない
    • 確率論的モーダストレンスへの批判
  • ネイマン - ピアソン
    • 行為選択の規則

GAMをもう少し理解したい④

前回の続きです。過去記事はこちらから。

ushi-goroshi.hatenablog.com

ushi-goroshi.hatenablog.com

ushi-goroshi.hatenablog.com

GAMの実装

bakfit

さて bakfit ですが、このサブルーチンはバックフィッティングというアルゴリズムを定義しています。このアルゴリズムは、簡単に言えば「関心のある平滑化対象変数を順番に残差に対してフィッティングしていく」方法です。Wikipediaバックフィッティングのページからアルゴリズムの部分について抜粋しましょう:

f:id:ushi-goroshi:20191125142457p:plain

ポイントは4行目の (a) のところで、


{\displaystyle \mathrm{Smooth} \left[\left\{ y_{i} - \hat{\alpha} - \sum_{k \neq j}\hat{f_{k}}(x_{ik})\right\}_{1}^{N}\right]}

とあるように、目的変数から定数部(線形部分)および関心のある平滑化対象変数"以外"(非線形部分)を引いています。すなわち「平滑化対象変数で説明される部分」と「確率的な誤差」で構成される部分のみを残し、平滑化を行なっていることがわかります。これを非線形としたい変数全てについて収束するまで順次繰り返す、というのがバックフィッティングの流れです。

ちなみにこのバックフィッティング、Wikipediaページでは以下のような説明があります:

the backfitting algorithm is equivalent to the Gauss–Seidel method algorithm for solving a certain linear system of equations.

さらにGauss–Seidel法について調べると(こちら)、これはn元の連立一次方程式を反復的に解くための手法であることが解説されているのですが、その中で以下の記述があります:

The Gauss–Seidel method now solves the left hand side of this expression for x, using previous value for x on the right hand side.

で、これは以下のように解くことができます:


 {\displaystyle x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum _{j=i+1}^{n}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\dots ,n.}

ここで括弧の中に着目すると、$b$から$ax$が順次引かれている(二項目と三項目)のですが、総和記号のインデックスには$i$が含まれていません。つまり更新の対象である$i$番目の$x$以外を順番に引いており、バックフィッティングと同じような計算になっていることがわかります。なお括弧の中の三項目は添え字が (k) ですが、二項目では (k+1) になっており、更新後の$x$を直接使っていることがわかります。さらにちなみに、ここで逐次更新するのではなく一度全ての解を求めてから一斉に更新する方法をヤコビ法といいます。

では上記のアルゴリズムbakfit でどのように実装されているのかを確認してみましょう。と言いつつ実は bakfit 自体は基本的に以下のように backf1 を呼び出しているだけです。

subroutine bakfit(x,npetc,y,w,which,spar,dof,match,nef,
                  etal,s,eta,beta,var,tol,
                  qr,qraux,qpivot,effect,work)
implicit double precision(a-h,o-z)
logical ifvar
integer npetc(7),iter
integer n,p,q,which(*),match(*),nef(*),nit,maxit,qrank,qpivot(*)
double precision x(*),y(*),w(*),spar(*),dof(*),
etal(*),s(*),eta(*),beta(*),var(*),tol,
qr(*),qraux(*),effect(*),work(*)
n=npetc(1)
p=npetc(2)
q=npetc(3)
ifvar=.false.
if(npetc(4)==1)ifvar=.true.
maxit=npetc(6)
qrank=npetc(7)
do i=1,q{work(i)=dof(i)}
### ここで backf1 を呼び出している
call backf1(x,n,p,y,w,q,which,spar,dof,match,nef,
            etal,s,eta,beta,var,ifvar,tol,nit,maxit,
            qr,qraux,qrank,qpivot,effect,work(q+1),work(q+n+1),
            work(q+2*n+1),work(q+3*n+1),work(q+4*n+1))

npetc(7)=qrank
return
end

なので backf1 を確認することにしましょう。

backf1

前回同様まずは backf1 の処理全体をさっと眺めて見通しを良くします。 backf1 ではざっくりと

  1. 前処理
  2. バックフィッティング
  3. eta を再定義

という処理が行われます。

subroutine backf1(x,n,p,y,w,q,which,spar,dof,match,nef,
                  etal,s,eta,beta,var,ifvar,tol,nit,maxit,
                  qr,qraux,qrank,qpivot,effect,z,old,sqwt,sqwti,work)
implicit double precision(a-h,o-z)
logical ifvar
integer n,p,q,which(q),match(n,q),nef(q),nit,maxit,qrank,qpivot(p)
double precision x(n,p),y(n),w(n),spar(q),dof(q),
etal(n),s(n,q),eta(n),beta(p),var(n,q),tol,
qr(n,p),qraux(p),effect(n),work(*)
double precision z(*),old(*),dwrss,ratio
double precision sqwt(n),sqwti(n)
logical anyzwt
double precision deltaf, normf,onedm7
integer job,info

onedm7=1d-7
job=1101;info=1
if(q==0)maxit=1
ratio=1d0

## 1. 前処理
### データに重みをつける
anyzwt=.false.
do i=1,n{
  if(w(i)>0d0){
    sqwt(i)=dsqrt(w(i))
    sqwti(i)=1d0/sqwt(i)
  }
  else{
    sqwt(i)=0d0
    sqwti(i)=0d0
    anyzwt=.true.
  }
}

### QR分解する
if(qrank==0){
  do i=1,n{
    do j=1,p{
      qr(i,j)=x(i,j)*sqwt(i)
    }
  }
  do j=1,p{qpivot(j)=j}
  call dqrdca(qr,n,n,p,qraux,qpivot,work,qrank,onedm7)
}

### eta に 非線形項 s(i,j) を加算する
do i=1,n{
  eta(i)=0d0
  for(j=1;j<=q;j=j+1){
    eta(i)=eta(i)+s(i,j)
  }
}


## 2. バックフィッティング
nit=0
while ((ratio > tol )&(nit < maxit)){
  deltaf=0d0
  nit=nit+1
  
  ### まずは Linear part + error 部分に最小二乗法を当てはめる
  do i=1,n{
    z(i)=(y(i)-eta(i))*sqwt(i) 
    old(i)=etal(i)
  }
  call dqrsl(qr,n,n,qrank,qraux,z,work(1),effect(1),beta,
             work(1),etal,job,info)
  
  ### 重みの逆数で戻し、 etal を更新
  do i=1,n{
    etal(i)=etal(i)*sqwti(i)
  }
  
  ### 次に平滑化対象の変数についてループ
  for(k=1;k<=q;k=k+1){
    j=which(k)
    do i=1,n{
      old(i)=s(i,k)
      z(i)=y(i)-etal(i)-eta(i)+old(i)
    }
    if(nit>1){dof(k)=0d0}
    
    ### splsmで平滑化
    call splsm(x(1,j),z,w,n,match(1,k),nef(k),spar(k),
               dof(k),s(1,k),s0,var(1,k),ifvar,work)
    do i=1,n{ 
      eta(i)=eta(i)+s(i,k)-old(i)
      etal(i)=etal(i)+s0
    }
    
    deltaf=deltaf+dwrss(n,old,s(1,k),w)
  }
  
  ### whileループの打ち切り判定
  normf=0d0
  do i=1,n{
    normf=normf+w(i)*eta(i)*eta(i)
  }
  if(normf>0d0){
    ratio=dsqrt(deltaf/normf)
  }
  else {ratio = 0d0}
}

### 3. eta を再定義
do j=1,p {work(j)=beta(j)}
do j=1,p {beta(qpivot(j))=work(j)}
if(anyzwt){
  do i=1,n {
    if(w(i) <= 0d0){
      etal(i)=0d0
      do j=1,p{
        etal(i)=etal(i)+beta(j)*x(i,j)
      }
    }
  }
}

do i=1,n
eta(i)=eta(i)+etal(i)

do j=1,q {
  call unpck(n,nef(j),match(1,j),var(1,j),old)
  do i=1,n {var(i,j)=old(i)}
}

return
end

3つ目の工程で eta の再定義とあるのですが、これまでの処理(〜 s.wam)において eta は「 y をリンク関数で変換したもの」として扱われてきました。一方 backf1 の中では 「平滑化対象変数によって構成される部分」(非線形/平滑化部分)として扱われており、混乱しやすいので注意が必要です。つまり backf1 の中では y を「線形部分」と「非線形(平滑化)部分」に分割し、前者を etal 、後者を eta として計算が行われます。

では一つずつ見ていきましょう。

1. 前処理

まずは前処理です。ここでは

  • データに重みをつける
  • QR分解する
  • eta非線形部分)を定義する

という処理が行われます。

### weight の平方根をとって逆数にする
anyzwt=.false.
do i=1,n{
  if(w(i)>0d0){
    sqwt(i)=dsqrt(w(i)) ## Fortranの組み込み関数。sqrt()
    sqwti(i)=1d0/sqwt(i) ## weight の逆数
  }
  else{
    sqwt(i)=0d0
    sqwti(i)=0d0
    anyzwt=.true.
  }
}

### qrank が 0 のときは、x に重みを加えた上で QR 分解する
if(qrank==0){
  do i=1,n{
    do j=1,p{
      ### 重みを乗じた x を生成
      qr(i,j)=x(i,j)*sqwt(i) # qr は n*p の行列
    }
  }
  do j=1,p{qpivot(j)=j} ### qvivot は長さ p の integer1~p の数字を格納
  
  ### QR分解。dqrdca のソースは以下
  ### https://github.com/cran/gam/blob/master/inst/ratfor/linear.r
  call dqrdca(qr,n,n,p,qraux,qpivot,work,qrank,onedm7)
}

### eta に 非線形項 s(i,j) を加算する
### ただしここで eta は fitted value ではなく 0 スタート
### s は線形部分とは独立な平滑化部分で、バックフィッティングで更新する
do i=1,n{
  eta(i)=0d0 # 0.0 にリセット 
  for(j=1;j<=q;j=j+1){ # 平滑化対象の変数について
    eta(i)=eta(i)+s(i,j) # 平滑化部分を順次加算していく
  }
}

上のブロックでQR分解を行なっていますが、これは次の工程で線形部分に対して通常の回帰を実行するためです。以前の記事で glm がどのように実装されているかを調べた際、 glm では、

  • dqrdc2 というサブルーチンを用いてQR分解
  • dqrsl というサブルーチンで最小二乗法を当てはめ

という処理で解が得られることを示しましたが、 glm 同様に gam においても線形部分のフィッティングではQR分解によって最小二乗解を求める、という流れになっているようです。

ushi-goroshi.hatenablog.com

2. バックフィッティング

では続いてバックフィッティングの実装を見ていきましょう。ここが、ある意味で gam の本体とも言えるルーチンになります。

### while ループ開始
nit=0
while ((ratio > tol )&(nit < maxit)){ # ratio のデフォルトは 1.0 。 tol は多分 0.0005
  # first the linear fit
  deltaf=0d0
  nit=nit+1 # イテレータに +1
  
  ### まずは Linear part + error 部分に最小二乗法を当てはめる
  ### y から非線形項を減じたもの(Linear part + error)に重みを乗じて z を更新(これを y として dqrsl に渡す)
  # etal を old に格納。ループ1回目の時点では etal は Nulldo i=1,n{
    z(i)=(y(i)-eta(i))*sqwt(i) 
    old(i)=etal(i)
  }
  
  ### https://github.com/wch/r-source/blob/trunk/src/appl/dqrsl.f
  ### dqrsl の6個目の引数として z を渡しているが、これは y になる
  ### etal は11番目の引数として渡しているが、これは xb 。
  call dqrsl(qr,n,n,qrank,qraux,z,work(1),effect(1),beta,
             work(1),etal,job,info)
  
  ### 重みの逆数で戻し、 etal を更新
  ### etal は y の予測値の線形部分
  do i=1,n{
    etal(i)=etal(i)*sqwti(i)
  }

まずは上のブロックで線形部分のフィッティングが行われます。 dqrsl の直前の処理に注目すると、 y から eta を減じたものを z として定義し、これを dqrsl の引数として渡しています。すなわち目的変数から eta非線形部分)を減じた、線形部分ですね。これを、先ほど解説したように dqrsl というサブルーチンに渡し、最小二乗法による解を得ます。引数には eta ではなく etal が入っていることに注意しましょう。この etal は直後の処理で重みを乗じた上で後続の処理に引き渡されます。

続いて平滑化対象変数の処理に入ります。

  ### ここで平滑化対象の変数についてループ
  for(k=1;k<=q;k=k+1){ # 平滑化対象の変数ごとに
    j=which(k)
    do i=1,n{
      ### old に k 番目の平滑化対象変数の値を入れる
      old(i)=s(i,k)
      ### y から etal (Linear part)と eta (Nonlinear part) を引き、old (k 番目の Nonlinear part) を足している
      ### 残差 + k 番目の平滑化変数の値になっている
      z(i)=y(i)-etal(i)-eta(i)+old(i)
    }
    ### df は 0 にリセットされてしまう
    if(nit>1){dof(k)=0d0}
    
    ### splsm を呼びだす
    ### 切片と k 番目の平滑化変数を x として渡す
    ### 9番目の引数 s が splsm 内での平滑化変数 smo なので s が更新される
    call splsm(x(1,j),z,w,n,match(1,k),nef(k),spar(k),
               dof(k),s(1,k),s0,var(1,k),ifvar,work)
    
    do i=1,n{
      ### eta を更新
      ### 古い eta にk番目の非線形部分を加算して、古い非線形部分を減じる)
      eta(i)=eta(i)+s(i,k)-old(i)
      ### etal を更新(古い etal に s0 を乗じる。)
      ### s0 は weighted mean of y
      etal(i)=etal(i)+s0
    }
    ### deltaf を更新(while ループの打ち切り判定)
    deltaf=deltaf+dwrss(n,old,s(1,k),w)
  } # ここまで for ループ

細かな処理はコメントとして書いてありますが、このブロックでは平滑化の対象となる変数について平滑化を実行して s を更新します。ポイントは

  z(i)=y(i)-etal(i)-eta(i)+old(i)

のところで、この記事の冒頭でバックフィッティングのアルゴリズムを解説した際に「関心のある平滑化対象変数を順番に残差に対してフィッティングしていく」と述べましたが、そのために残差を求めているのが上の式になります。 etal は線形部分、 eta非線形部分なので両者を引くと確率的な誤差しか残りませんが、そこに old(i) つまり s(i, k) を加算することで、k番目の平滑化対象変数 + 確率的な誤差のみが残ります。

続いてこの誤差に対して平滑化を行なっているのが

splsm(x(1,j),z,w,n,match(1,k),nef(k),spar(k),
      dof(k),s(1,k),s0,var(1,k),ifvar,work)

の部分です。この splsm は内部でさらに sbart というCのサブルーチンを呼び出しているのですが、この解説を始めるとキリがないのでまた別の機会に取っておきましょう。

平滑化により s が更新されたら、 eta にもk番目の平滑化対象変数の更新を反映させ、最後に while ループの打ち切り判定のための指標を計算します。

### normf の更新(while ループの打ち切り判定)
normf=0d0
do i=1,n{
  normf=normf+w(i)*eta(i)*eta(i)
}
### ratio(while ループの判定に使われる)を計算
if(normf>0d0){
  ratio=dsqrt(deltaf/normf)
}
else {ratio = 0d0}
#         call DBLEPR("ratio",-1,ratio,1)
}
3. eta の再定義

以上の工程で十分に収束したと判断されたらバックフィッティングアルゴリズムは終了し、 s.wam に値を返す処理に入ります。冒頭で解説したように backf1 内では eta非線形部分として扱われていたので、これを線形部分と加算して元の eta に戻します。なお beta は線形部分の各変数の回帰係数です。

### etal を再定義
do j=1,p {work(j)=beta(j)}
do j=1,p {beta(qpivot(j))=work(j)}
if(anyzwt){
  do i=1,n {
    if(w(i) <= 0d0){
      etal(i)=0d0
      do j=1,p{
        etal(i)=etal(i)+beta(j)*x(i,j)
      }
    }
  }
}

### 非線形項 + 線形項として eta を再定義
do i=1,n
eta(i)=eta(i)+etal(i)

### 結果をまとめる    
do j=1,q {
  call unpck(n,nef(j),match(1,j),var(1,j),old)
  do i=1,n {var(i,j)=old(i)}
}

return
end

さて、これで backf1 および bakfit の処理は終了したので s.wam に戻ることになります。 s.wam の残りの工程では、前回記事で確認したように結果を格納したオブジェクトである rl を生成するのでした。それと元々のデータを順位で置き換える形で生成された smooth.frame をまとめて返します。

さらに s.wam を抜けると gam.fit 内でのループに戻ってくるわけですが、これは前々回記事で解説した通り、残差からデビアンスを計算した上でループの打ち切りを判定します。

gam.fit

収束が十分であれば、 gam.fit のループを抜けます。残りの工程は基本的に結果をまとめあげて gam に返すだけのようです。

### gam.fit に戻ってきた
fitqr <- fit$qr
xxnames <- xnames[fitqr$pivot]
nr <- min(sum(good), nvars)
if (nr < nvars) {
  Rmat <- diag(nvars)
  Rmat[1:nr, 1:nvars] <- fitqr$qr[1:nr, 1:nvars]
}
else Rmat <- fitqr$qr[1:nvars, 1:nvars]
Rmat <- as.matrix(Rmat)
Rmat[row(Rmat) > col(Rmat)] <- 0
dimnames(Rmat) <- list(xxnames, xxnames)
names(fit$residuals) <- ynames
names(mu) <- ynames
names(eta) <- ynames

### eta and mu
fit$additive.predictors <- eta
fit$fitted.values <- mu
names(fit$weights) <- ynames
names(fit$effects) <- c(xxnames[seq(len = fitqr$rank)], rep.int("", 
                                                                sum(good) - fitqr$rank))
if (length(fit$smooth) > 0) 
  fit$smooth.frame <- smooth.frame[smooth.labels]
wtdmu <- if (a$intercept) 
  sum(weights * y)/sum(weights)
else linkinv(offset)
nulldev <- sum(dev.resids(y, wtdmu, weights))
n.ok <- nobs - sum(weights == 0)
nulldf <- n.ok - as.integer(a$intercept)
rank <- n.ok - fit$df.residual
aic.model <- aic(y, nobs, mu, weights, new.dev) + 2 * rank
if (!is.null(fit$smooth)) {
  nonzeroWt <- (wz > 0)
  nl.chisq <- gam.nlchisq(fit$qr, fit$residuals, wz, fit$smooth)
}
else nl.chisq <- NULL
fit <- c(fit, list(R = Rmat, rank = fitqr$rank, family = family, 
                   deviance = new.dev, aic = aic.model, null.deviance = nulldev, 
                   iter = iter, prior.weights = weights, y = y, df.null = nulldf, 
                   nl.chisq = nl.chisq))
fit
}

gam

ようやく gam に帰ってきました。 gam.fit を呼び出した以降の処理は以下のブロックだけです。特定の条件下においてNull Deviance(切片のみ考慮されたモデルによる逸脱度)を計算し、結果をまとめあげて返します。

### offset が指定されており intercept 項がある場合
if (length(offset) && attr(mt, "intercept") > 0) {
  fit$null.dev <- glm.fit(x = X[, "(Intercept)", drop = FALSE], 
                          y = Y, weights = weights, offset = offset, family = family, 
                          control = control[c("epsilon", "maxit", "trace")], 
                          intercept = TRUE)$deviance
}
if (model) 
  fit$model <- mf
fit$na.action <- attr(mf, "na.action")
if (x) 
  fit$x <- X
if (!y) 
  fit$y <- NULL
fit <- c(fit, list(call = call, formula = formula, terms = mt, 
                   data = data, offset = offset, control = control, method = method, 
                   contrasts = attr(X, "contrasts"), xlevels = .getXlevels(mt, 
                                                                           mf)))
class(fit) <- c("Gam", "glm", "lm")
if (!is.null(fit$df.residual) && !(fit$df.residual > 0)) 
  warning("Residual degrees of freedom are negative or zero.  This occurs when the sum of the parametric and nonparametric degrees of freedom exceeds the number of observations.  The model is probably too complex for the amount of data available.")
fit
}

終わりに

以上で終わりです!長くなりましたが、ようやく gam の内部で何が行われているかを確認することができました。振り返ると、 gam という関数のコアとなる部分では:

  1. バックフィッティングによる非線形項のフィッティング( bakfit および backf1
  2. バックフィッティングに渡すためにデータを順位に置き換える( s.wam
  3. {general, s, lo}.wam のいずれのエンジンに渡すかを特定し、 bf.call を作成( gam.fit
  4. gam.fit に渡すための諸々の設定( gam

といった処理が行われているようでした。またバックフィッティングというアルゴリズムは、端的には「関心のある変数による影響のみを残してフィッティングすることを繰り返す」というものであることがわかりました。

一方で、今回の一連の記事では平滑化の仕組みを追いかけることがまだ出来ていません。また平滑化の関数にしても s.wam しか見ておらず、 lo.wamgeneral.wam の挙動を確認していません。その辺りもいずれ調べたいですね。

樋口先生退任記念シンポジウム参加記録

11/5に開催された元統計数理研究所の所長である樋口先生の退任記念シンポジウムに参加してきましたので、そのメモを共有しておきます。ご参考まで。やっぱり長らく研究をされてきた方のお話というのは面白いですね。

ご挨拶

統数研 椿広計 氏

バイオインフォマティクスベイズモデリング

東京大学 井本清哉 氏

  • 樋口先生との共同研究
    • ベイジアンネットワークによる遺伝子ネットワークの推定
      • 遺伝子間の因果関係(制御関係)
      • グラフ構造の推定
      • ヒトゲノム、21000個の遺伝子
        • それらが複雑に関係し合う
        • たかだか有限のデータでは推定できない
      • ネットワーク構造の事後確率
      • 生物学的知識、他のデータで事前確率を構成
        • 結合配列
        • タンパク質の相互作用
        • 進化情報
        • データベース
        • 遺伝子ノックダウン
      • バイオインフォマティクスベイズモデリングの格好の応用対象であり、問題の宝庫
      • マイクロバイオームを含んで「私」
      • 健康や病気をデータ化できる時代が来た
    • データ同化による細胞内分子シミュレーションモデルの構築
    • インフルエンザ感染拡大エージェントベースシミュレーション
      • インフルエンザの拡大を抑えるためには会社員のワクチン接種が有効
  • 予測だけを目的としたら、AI(機械学習技術)は人を超えたかも
    • Kaggle、AutoMLがコンペで2位になった
  • ゲノム医療、理由が説明できないと医療に入っていくことは難しい現状
    • 医療分野で AI は Augumented Intelligence
    • 人の能力を向上させるもの
  • 医科学研究所、2900万件の医学論文を学習
  • 思考のリミッターを外しなさい by 樋口先生

消費者行動理解のためのベイジアンモデリング

筑波大学 佐藤忠彦 氏

  • マーケティングに関わる領域も最近は樋口先生の表芸になってきている
  • 能書き
    • マーケティング、現実の市場で発生する現象
    • 問題を解決するために以下の問いに答える必要がある
      • マーケティングを取り囲む環境は?
      • 市場はどこにある?
      • 敵は誰?
      • 顧客はどこに?
      • 顧客の嗜好は?
      • 消費者の価格反応は?
      • 広告、販促への反応は?
      • ものやサービスを選ぶ基準は?
    • 上記は経験や感覚で評価できるか?
      • 科学的アプローチが必要
    • 消費者を理解する = 仮説を発見する
    • 消費者を理解する ≠ 真の姿を解明する
    • 消費者理解のプロセス
      • 複数のモデルを立て
      • 全てのモデルを推定・比較し
      • この時点でもっとも妥当と考えられる構造の推定
      • 情報の変換
      • 結果の活用
    • 現在持ちうる知識で記述能力・汎化能力が高いモデルを選ぶ
    • 消費者理解のキーワード
      • 消費者異質性
      • 時間的異質性
      • 潜在変数
    • パラメータがいっぱいで普通のモデルでは対応できない
    • 事後分布が関心の対象
    • 樋口流ベイジアンモデリング
      • 理論のみならずドメイン知識に基づく
      • こんな構造はないか、あんな構造はないか、という構成的なモデリング
  • マーケティングの事例①:佐藤、樋口(2009)
    • スーパーで牛乳を購買する確率が、潜在変数である家庭内在庫量/消費量の推移に基づき変化する構造
    • 一般状態空間
    • 店に来て牛乳を買う状況
    • n時点での在庫量 = 前日開始時点での在庫量 + 購買量 - 消費量
    • 消費量 = 残っている量 * 消費率
    • 店に来る効用、牛乳を買う効用
    • 動的に変化する
    • 在庫量が変化すると購入確率がどのように変化するか
  • マーケティングの事例②:宮津、佐藤(2015)
    • スーパーにおける購買個数が、心理的財布の状況により変化する構造
    • 給料日直前、給料日
    • 25、10、17日を起点とした3つの累積購買金額の線形結合で柔軟に構成する
    • 閾値ポアソン回帰
      • 心理的財布が閾値に到達しているか否かで行動が変わる
  • まとめ
    • Supply Side Thinking ではなく Demand Side Thinking
    • ぐにゃぐにゃモデルのすすめ
    • ベイズモデリングが有効

固体地球科学とデータ同化

東京大学 長尾大道 氏

  • 北川源四郎 時系列解析プログラミング
  • データ同化
    • 数値シミュレーションと観測データの融合
    • 大規模な状態空間モデルと大規模データ
  • 四次元変分法
  • データ同化入門

樋口先生とブリジストンと私

ブリジストン 花塚泰史 氏

  • 樋口先生と10年以上の関係がある
  • タイヤ事業が8割、多角化事業が2割
  • 売り切り型からソリューションプロバイダーとして
  • 状態空間モデルとカーネルマシンを併用した・・・
  • ブリジストンでの事例
    • タイヤ、ほぼ唯一の電子化されていない部品
    • CAISによる安全・安心なモビリティ社会への貢献
      • インテリジェント化したタイヤ
      • 空気圧、温度を拾うのに加え、タイヤ内のセンサーから
        • 磨耗状態推定
        • 路面状態推定 → 自動車の制御へフィードバック
          • 加速度
          • 路面状態に応じて波形が特徴的に変化
            • ウェット、凍結
            • 周波数領域において特徴的な波形
          • Windowをシフトしながら
          • カーネルSVM
            • (速度によって)系列長が異なるものを比較
            • グローバルアライメントカーネル
              • 計算量がでかい
              • リアルタイムでやりたい
            • 窓同士の類似度
            • ドライとウェットを分類する判別問題
          • GAカーネルの高速化
            • 類似度を計算する窓を制限する
            • タイヤの波形は特徴的な2つのピークが出るので・・・
          • 元のGA、93.9%で判別できるがタイヤ100回転分の時間が必要
          • 工夫により96.9%、タイヤ一回転分ぐらいで判断できるようになった
            • 7種類に判別
          • 冬季道路管理の最適化への適用
            • 路面のメンテナンス
              • 巡回車に乗っている人の目で判断
            • 北海道で実用化されている
              • ISCOS
    • 次世代のデータサイエンティスト育成
      • 社内研修
  • 将来的な発展
    • センシング技術
    • 解析、予測
    • デジタルサービス

ベイズモデリングと歩んだ30年

中央大学 樋口知之 氏

  • 平成元年に入所、平成31年に退任
    • 200MBのデータ解析に大型計算機が必要だった時代
  • 人間の未来予測は当たらない
  • モデルとモデリング
    • 何が違うのか
      • 常にモデルを改良する姿勢
      • ものの見方、捉え方を柔軟に変化
      • 思考のプロセス自体を科学する
    • モデルを比較するための羅針盤としての情報量基準
  • 赤池先生との出会い
    • 多変量自己回帰を応用したセメント工場、火力制御
    • 船のオートパイロット(北川先生)
  • 当時
    • 支配方程式がない、計算コストが大
  • 生成モデルと識別モデル
    • 生成モデルはデータが生成される過程
    • 識別モデルはデータやパラメータを与えたときにどのように識別されるか
    • 後者が簡単(機械学習コミュニティで盛ん)
    • 統計の人たちはデータがどんな機序で生まれたのかに強い関心がある
  • 頻度主義 vs ベイズ主義
  • ベイズの定理がなぜ役立つか?
    • コンピュータの性能向上
      • 周辺尤度の計算
    • 高精度センサーのコモディティ化
      • 尤度関数
    • ストレージの廉価化
      • 細かい情報を保存
  • 1980年代の統数研
  • 統数研
    • 統数研の助手は世界でもっとも恵まれている
    • 人材採用方針
      • 数理統計、確率論
      • 物理、計算科学、計算機科学
      • データの質、量の観点で大きく変化のある分野から
    • 赤池先生
  • アカデミアから産業界への人材流出は欧米でも同じ
    • クロスアポイント
  • メゾスコピックモデリング
    • 数理表現された法則
    • 第一原理法則
  • 演繹と帰納
  • 逐次データ同化技術
    • KFとPFの間にある
    • 状態ベクトルの次元数、アンサンブルメンバー数
  • 統計的推論法のパラダイムシフトを促す外的要因
    • ビッグデータ
    • データ駆動サービス社会
      • 因果から相関へ
    • 目的特化型計算機
  • 分析者の技術から計算パワーへ
  • 力学的ダウンスケーリング
  • エミュレーション、仮想計測
  • やっていること
    • GANを使ってエミュレータを作る
    • 深層学習と状態空間モデルの融合
    • 高等教育
  • 機械学習と統計
    • Data driven vs Science driven
    • Predictive model vs Descriptive model
    • Correlation vs Causality
    • 認識科学から設計科学
    • 対象理解から帰納の最適化
    • 「真のモデルなど存在しない」という見方は、当面の問題のイメージ能力の欠如

AI社会での数理モデリングと現場主義の意義と価値

東京大学 合原一幸 氏 東北大学 早瀬敏幸 氏 大阪大学 鷲尾 隆 氏 ブリジストン 花塚泰史 氏 中央大学 樋口知之 氏

  • 分野横断的なアプローチで豊富な共同研究実績の持ち主
  • 合原氏
    • 数理工学
      • 統計とは兄弟みたいな
      • 暮らしを変える驚きの数理工学
      • ニューラルコンピュータ
      • アルファ碁解体新書
    • AIと人間の共同作業、協働
      • Pair Go
      • ファッション
        • GANでデザイナーっぽいデザイン
  • 早瀬氏
    • 流体力学、制御工学、生体工学
    • データ同化
  • 鷲尾氏
    • 演繹的なモデリングの自動化
      • 要素還元主義・普遍主義
      • 外挿に利用可能
      • 深い理解が可能
    • 帰納的なモデリング
      • 高精度な近似はできるが外挿できない
      • XAIの観点から、浅い理解しかできない
    • 演繹的モデルと帰納的モデルのより広範な融合
  • 三者は生成モデルにこだわりがある
  • 3つのお題(6つの事前質問から)
    • 現場主義
      • 1.現場主義とKKDの違いは? 現場主義に重要性について
      • 2.学生や研究者に現場主義をどう教育していくべきか?
        • 樋口
          • 現場に出向く研究者魂
          • 現場の人と厳しく議論できるコミュニケーション力、分野跳躍力
          • 本質的な課題を数理の問題にしたてるセンス
          • しつこさ
        • 鷲尾
          • 向こうの分野の研究者に統計や機械学習を教えた方が良いという経験
          • 現場を知っている人の方がよく使える
        • 合原
          • 中小企業において、社長が現場を見に行く
            • そうすると階層構造にショートカットが生まれる
        • 早瀬
          • 研究時間が減っていて、研究に割ける時間がない
        • 花塚
          • ブリジストンでの取り組み
            • 工場のオペレーター、タイヤ館の店員
            • 何に困っているか、ソリューションが役にたつかを話を聞きに行く
            • 障害
              • 相手が自分のナレッジを理解していない
              • KKDが正しいかの見極め
            • ドメインの知識を作り手側が持っていないといけない
      • コミュニケーションが取りづらい相手は?
        • 無茶振りする(合原)
          • 専念しても一年かかるようなもの
          • データ化されていないものを大量に持ってくる
          • 割ける時間の中でうまく成果が出るテーマ選び
          • 他の人に無茶振りしてみる
            • その人だったら無茶ではないかもしれない
            • 人脈づくり、人的ネットワーク
              • 有名な人を知っているだけではダメ
              • あの人はいま時間がある
        • 分析・テーマの振り返りが繰り返されるので、じっくり取り組む必要がある(鷲尾)
    • 人とAI
      • 3.複雑な現象を数理モデルで表現する意義は?
      • 4.今後の研究スタイルはどう変わっていくか?
        • 合原
          • Data drivenは実用上は問題ない
          • そこから法則をスッキリ抜き出すことができれば面白い
          • シンプルな非線形性で抜き出したい
          • ニューラルネットは完全にTransparent
            • 言葉にできないだけ
        • 鷲尾
          • 機械学習の研究者が精度を競っているので精一杯
          • 最近XAIが成長してきた
          • 中身を理解しやすいモデルはどうやって作るか
          • 現在の機械学習は、現実との引き合わせができていない
        • 花塚
          • 製造業においてはケースバイケース
        • 早瀬
          • 順モデルなのか逆モデルなのか
          • どのような過程で生まれたのか説明される必要がある
          • 現実のシステムがもつ因果
    • 教育システム
      • 5.分野横断的な学術体系をどのようにして教育システムに組み込むか
      • 6.10年後に国立大学はどのような姿になっているか
        • 東北大学では全員にデータサイエンス系の講義を受けさせることになる(早瀬)
        • オーストラリアの某大学は数学科を三等分している(合原)
        • 若い人を早めに独立させる(合原)
          • 教授、准教授、助教みたいな階層構造をやめる
        • 数学が苦手な人たちにいかに脱落させないか(鷲尾)
          • OJTとして授業より前に使ってもらう
        • 古典的製造業は縦割りがピシッとしている(花塚)
          • 物売りからこと売り、ソリューションやるなら横串が必要
          • サプライチェーンの初めから終わりまでを担う
            • それがデータサイエンス組織
          • コミュニケーションを密にして話をしにいく
          • 現場ごとにデータサイエンスを理解している人が必要
          • CoEとしての機能
          • 育成
            • CoE
              • 高度な手法を理解できる
            • 現場
              • 話を理解できる人を多く育てる
              • 裾野を広げる
          • 知っている人をいかにつなげていくか
      • 統数研に期待すること、文科省に頑張ってほしいこと
        • 予算(合原)
        • 日本のトップの研究所として世界をリードする研究(早瀬)
        • 地道に学問を究めるところが魅力(鷲尾)
          • 外にアピールして外部資金を獲得
          • 基礎研究を地道にコツコツと続けられる
        • リカレント教育(花塚)
          • DSを外から取ってくるのは至難
          • 社内だけでの教育も困難
            • アカデミアとの連携

GAMをもう少し理解したい③

前回の続きです。過去記事はこちらから。

ushi-goroshi.hatenablog.com ushi-goroshi.hatenablog.com

GAMの実装

s.wam

s.wamgam.fit の中で bf.call として定義されたフィッティングのための関数で、 s.wamwamweighted additive modelの意味のようです。名前からしgam の本体のような気がしてきますが、先に言ってしまうとこの関数内で bakfit というFortranの関数を呼び出しており、少なくともRの「 gam というライブラリの本体」という意味では合っているのではないでしょうか。

s.wam の説明に入る前に、見通しを良くするために関数全体をさっと眺めてみましょう。 s.wam はざっくりと以下の4つの処理に分けることができそうです:

  1. 平滑化の対象変数について元のデータを順位に置き換え、smooth.frameを再定義する
  2. 後でFortranに渡すために必要な指定を行う
  3. Fortran で書かれた bakfit を呼び出す(本体)
  4. 後処理
### gam::s.wam
function (x, y, w, s, which, smooth.frame, maxit = 30, tol = 1e-07, 
          trace = FALSE, se = TRUE, ...) 
{
  ### 1. 平滑化の対象変数について元のデータを順位に置き換え、smooth.frameを再定義する
  if (is.data.frame(smooth.frame)) {
    first <- TRUE
    data <- smooth.frame[, names(which), drop = FALSE]
    
    smooth.frame <- gam.match(data)
    dx <- as.integer(dim(x))
    smooth.frame$n <- dx[1] 
    smooth.frame$p <- dx[2] 
    oldClass(data) <- NULL
    smooth.frame$spar <- unlist(lapply(data, attr, "spar"))
    smooth.frame$df <- unlist(lapply(data, attr, "df"))
  }
  else first <- FALSE
  
  ### 2. 後でFortranに渡すために必要な指定を行う
  storage.mode(tol) <- "double"
  storage.mode(maxit) <- "integer"
  which <- unlist(which)
  storage.mode(which) <- "integer"
  storage.mode(y) <- "double"
  storage.mode(w) <- "double"
  p <- smooth.frame$p
  n <- smooth.frame$n
  
  for (ich in which) x[, ich] = signif(x[, ich], 6)
  
  ### 3. Fortran で書かれた bakfit を呼び出す(本体)
  fit <- .Fortran("bakfit", x, npetc = as.integer(c(n, p, length(which), 
                                                    se, 0, maxit, 0)), y = y, w = w, which, spar = as.double(smooth.frame$spar), 
                  df = as.double(smooth.frame$df), as.integer(smooth.frame$o), 
                  as.integer(smooth.frame$nef), etal = double(n), s = s, 
                  eta = double(n), beta = double(p), var = s, tol, qr = x, 
                  qraux = double(p), qpivot = as.integer(1:p), effects = double(n), 
                  double((10 + 2 * 4 + 5) * (max(smooth.frame$nef) + 2) + 
                           15 * n + 15 + length(which)), PACKAGE = "gam")
  
  ### 4. 後処理
  #### 収束に関する警告
  nit <- fit$npetc[5]
  qrank <- fit$npetc[7]
  if ((nit == maxit) & maxit > 1) 
    warning(paste("s.wam convergence not obtained in ", maxit, 
                  " iterations"))
  
  #### ループの2回目以降は序盤の処理をスキップさせる
  if (first) {
    smooth.frame$spar <- fit$spar
    first <- FALSE
  }
  
  #### 最終的な返り値となる rl に情報を追加していくための準備
  names(fit$df) <- dimnames(s)[[2]]
  names(fit$beta) <- labels(x)[[2]]
  qrx <- structure(list(qr = fit$qr, qraux = fit$qraux, rank = qrank, 
                        pivot = fit$qpivot, tol = 1e-07), class = "qr")
  effects <- fit$effects
  r1 <- seq(len = qrx$rank)
  dn <- colnames(x)
  if (is.null(dn)) 
    dn <- paste("x", 1:p, sep = "")
  names(effects) <- c(dn[qrx$pivot[r1]], rep.int("", n - qrx$rank))
  
  #### rl を生成
  rl <- list(coefficients = fit$beta, residuals = fit$y - fit$eta, 
             fitted.values = fit$eta, effects = effects, weights = w, 
             rank = qrank, assign = attr(x, "assign"), qr = qrx, smooth = fit$s, 
             nl.df = fit$df - 1)
  rl$df.residual <- n - qrank - sum(rl$nl.df) - sum(fit$w == 
                                                      0)
  if (se) 
    rl <- c(rl, list(var = fit$var))
  c(list(smooth.frame = smooth.frame), rl)
}

以下、一つずつ見ていきましょう。

1. 平滑化の対象変数について元のデータを順位に置き換え、smooth.frameを再定義する

s.wam では始めに smooth.frame に対する処理を行います。具体的には

  • 平滑化の対象変数を抽出
  • 指定のdigitsで丸める
  • ソートした上でユニークなデータ数を得る
  • 元データをラベル(数値の大小で並べかえたときの順位)に変換
  • 行数・列数を付与
  • 平滑化パラメータおよび自由度を付与

という加工を、下記の処理にて実施します。

### smooth.frame はループ初回は data.frame なのでここが評価される
if (is.data.frame(smooth.frame)) {
  first <- TRUE
  data <- smooth.frame[, names(which), drop = FALSE] # 平滑化対象変数を抽出
  
  smooth.frame <- gam.match(data) # 下で解説
  dx <- as.integer(dim(x)) # 行列のサイズ
  smooth.frame$n <- dx[1] # number of row
  smooth.frame$p <- dx[2] # number of column
  oldClass(data) <- NULL
  ### 各列に対する平滑化パラメータおよび自由度
  smooth.frame$spar <- unlist(lapply(data, attr, "spar"))
  smooth.frame$df <- unlist(lapply(data, attr, "df"))
}
else first <- FALSE

上のブロックのポイントは gam.match という関数で、 gam.match(data) とすると以下の結果が返ります:

> gam.match(data)
$o
s(year, 4) s(age, 5)
[1,]          4         1
[2,]          2         7
[3,]          1        28
[4,]          1        26
[5,]          3        33
[6,]          6        37

### 中略

$nef
s(year, 4)  s(age, 5) 
7         61 

ちなみに data は元データです:

> head(data)
s(year, 4) s(age, 5)
231655       2006        18
86582        2004        24
161300       2003        45
155159       2003        43
11443        2005        50
376662       2008        54

これは gam.match という関数において下記の部分が評価された結果です。

### gam.match から
xr <- signif(as.vector(x), 6) # signif 関数は指定の有効桁で丸める。デフォルトは6
sx <- unique(sort(xr)) # それをソートしてユニークにする
nef <- as.integer(length(sx)) # ユニークなデータ数
if (nef <= 3) 
  stop("A smoothing variable encountered with 3 or less unique values; at least 4 needed")
o <- match(xr, sx, nef + 1) # 元のデータを順位に変更する。 +1 は欠損用。
o[is.na(o)] <- nef + 1

match の部分が少しわかりにくいかもしれませんが、以下のような結果が得られます:

> head(xr) # 元の値
[1] 2006 2004 2003 2003 2005 2008

> match(head(xr), sx, nef + 1) # 順位に変換
[1] 4 2 1 1 3 6

これは xr の値が sx の何番目に位置するかを表わしていますが、参照先の sx は元の値をソートしてユニークにしたものなので、順位に変換していることになります。 また gam.match の結果を smooth.frame に代入されているので、この時点で smooth.frame は数値の大小関係のみを持つことになり、また一連の処理で smooth.framedata.frame ではなくなるため、ループの2回目以降ではこの処理はスキップされます。

2. 後でFortranに渡すために必要な指定を行う

次のブロックは後ろの工程でFortranの関数に渡すための準備となります。

### 後でFortranに渡すために必要な指定
storage.mode(tol) <- "double"
storage.mode(maxit) <- "integer"
which <- unlist(which)
storage.mode(which) <- "integer"
storage.mode(y) <- "double"
storage.mode(w) <- "double"
p <- smooth.frame$p
n <- smooth.frame$n

また以下のコードでは、先程の smooth.frame 同様に有効桁を6桁に丸めます。

### 平滑化対象の変数を6桁に丸める
for (ich in which) x[, ich] = signif(x[, ich], 6)
3. Fortran で書かれた bakfit を呼び出す(本体)

以上で準備が整いました。以下が s.wam の本体になります。 .FortranFortranで書かれた bakfit というモジュールを呼び出します。

### ここが本体。Fortran で書かれた bakfit を呼び出す
fit <- .Fortran("bakfit", x, npetc = as.integer(c(n, p, length(which), 
                                                  se, 0, maxit, 0)), y = y, w = w, which, spar = as.double(smooth.frame$spar), 
                df = as.double(smooth.frame$df), as.integer(smooth.frame$o), 
                as.integer(smooth.frame$nef), etal = double(n), s = s, 
                eta = double(n), beta = double(p), var = s, tol, qr = x, 
                qraux = double(p), qpivot = as.integer(1:p), effects = double(n), 
                double((10 + 2 * 4 + 5) * (max(smooth.frame$nef) + 2) + 
                         15 * n + 15 + length(which)), PACKAGE = "gam")

さてこの bakfitソースコードとしてはおそらくこちら が正しいのかなぁと思うのですが、 gam をインストールする際に同時にダウンロードされる以下のファイル:

/Users/hogehoge/Library/R/3.6/library/gam/ratfor/backfit.r

の方が説明もあってわかりやすいので、こちらベースに進めたいと思います。なお前回記事でも述べたことですが、 s.wam は第二引数として y を受け取ります:

### s.wam
function (x, y, w, s, which, smooth.frame, maxit = 30, tol = 1e-07, 
          trace = FALSE, se = TRUE, ...) 

ただしこの s.wam にはYそのものではなく、以下のように z として渡されるのでした:

### 以下は gam.fit から再掲
z <- eta - offset
z[good] <- z[good] + (y - mu)[good]/mu.eta.val[good]

### 中略 

# eval(bf.call) ← これは下の表現となる
s.wam(x, z, wz, fit$smooth, which, fit$smooth.frame, bf.maxit, 
      bf.epsilon, trace)

要するにリンク関数で変換したあとのデータを渡している、ということです。

では bakfit に移りたいのですが、その前にまずは残りの工程を見ておきましょう。

4. 後処理

s.wam の最後の部分です。

nit <- fit$npetc[5] # number of iteration
qrank <- fit$npetc[7] # qrank

### nit が maxit と同じであった場合、収束していないことを警告
if ((nit == maxit) & maxit > 1) 
  warning(paste("s.wam convergence not obtained in ", maxit, 
                " iterations"))

### first == T の場合、spar を加えて FALSE にして返す
### ループ2回目以降は評価されない
if (first) {
  smooth.frame$spar <- fit$spar
  first <- FALSE
}

names(fit$df) <- dimnames(s)[[2]] # 平滑化対象変数の列名
names(fit$beta) <- labels(x)[[2]] # 説明変数
qrx <- structure(list(qr = fit$qr, qraux = fit$qraux, rank = qrank, 
                      pivot = fit$qpivot, tol = 1e-07), class = "qr") # 属性を指定しながらオブジェクトを作成
effects <- fit$effects
r1 <- seq(len = qrx$rank)
dn <- colnames(x)
if (is.null(dn)) 
  dn <- paste("x", 1:p, sep = "")
names(effects) <- c(dn[qrx$pivot[r1]], rep.int("", n - qrx$rank))

#### 最終的な返り値である rl を生成
rl <- list(coefficients = fit$beta, residuals = fit$y - fit$eta, 
           fitted.values = fit$eta, effects = effects, weights = w, 
           rank = qrank, assign = attr(x, "assign"), qr = qrx, smooth = fit$s, 
           nl.df = fit$df - 1)
rl$df.residual <- n - qrank - sum(rl$nl.df) - sum(fit$w == 
                                                    0)
if (se) 
  rl <- c(rl, list(var = fit$var))

#### return
c(list(smooth.frame = smooth.frame), rl)

基本的には s.wam という関数の返り値となる rl を生成するための工程と言えそうです。 bakfit はまた次回。

応用統計学フロンティアセミナー「データサイエンスと応用統計学」参加記録

10/19に開催された応用統計学フロンティアセミナーに参加してきましたので、そのメモを共有しておきます。話を聴きながらのメモなので単語しか書けず意味がわかりにくいところもありますが、ご参考まで。なおセミナーの様子は以下のtogetterでまとめられていますので、そちらも合わせてご覧ください。

togetter.com

データサイエンスにおける応用統計人材の育成

横浜市立大学 岩崎学 氏

  • これまでとこれからの統計的データ解析

    1. 研究目的の設定
    2. データ収集法の立案(実験、観察、調査)
    3. データの収集(モニタリング)
    4. データの電子化
    5. データのチェック・マージ
    6. データの集計とグラフ化(記述統計)
    7. 統計的推測ないしは予測(推測統計)
    8. プレゼン・文章化、意思決定(終了もしくは最初に戻る) → 役に立たないと意味がない
  • いまのデータサイエンス、データありきになっているかもしれない

    • 推測統計の部分にフォーカスが当たっている
    • 全体の流れを見ることが大事
  • データサイエンスとは?

    • 統計学 + 情報科学)* 社会展開
    • 理系的な要素と文系的な要素の両面
      • 文理融合は当然
  • 横浜市立大学での状況

    • 線形代数の学部生の試験成績、文系と理系でほとんど差がない(120点満点で76 vs 72とか)
    • 2020年、大学院
      • ヘルスデータサイエンス専攻(医学部があるので)
      • データサイエンス専攻(M20名、D3名)

楽天技術研究所におけるデータサイエンスおよび統計の様々な応用事例

楽天 森正弥 氏

  • 1.3Bユーザー、48Kの店舗(ビジネスパートナー)

  • 250Mのアイテムについて需要予測

    • 価格、在庫の調整
    • 非線形回帰。基盤を作ったロシア人が手作業でフィッティング。DLに勝つ。
  • 金融市場予測

    • Now-Cast
      • Forecastに対してNowcast
      • 今の状態を推測する
      • Hal Varian
      • 独自の景気動向指数
        • Google
        • 東大
        • 48Kのビジネスパートナーのデータを使って
        • 8000カテゴリの売上データ
        • LASSO回帰
          • 2521のカテゴリ
          • MAE0.4%
          • 宝石、エアコン、ワークステーション、コメディのブルーレイが売れるとCIが良くなる
          • 交互作用は見てない
        • DL派が増えた
          • データぶっこもうぜ
          • 丁寧な分析はしなくなった
    • 会社の株価
      • 売上と株価、2ヶ月のラグ
      • 相関が0.96
      • でも予測精度は芳しくない
    • 化粧品会社の株価
    • 多段のDoc2Vecでメーカー名のリンキング
    • 楽天トラベルのデータ、インバウンド

自動車とデータサイエンス

日産自動車 上田哲郎 氏

  • Connected Carのエキスパート、ITのエキスパート
  • 自動車、AI/ITに関しては新卒の人たちの方が専門性高い
  • 今はリアルタイムで車の車種まで特定できる
    • 二週間でインターンが片付ける
    • 正答率93.6%
    • アテンション層4000次元で車を分類できる知識構造ができているのでは
      • メーカーらしさを表すベクトルができている?
        • 今まではアンケートしかなかった
        • アテンション層4000次元を可視化
          • tSNE
          • 三次元:
            • メーカーらしさ?
            • 時間変化?
          • 二次元:
            • 車の向きが揃っているクラスターができる
            • ドアが開いている
          • 海外メーカーはクラスターが分かれる
          • 日本車は特徴ない?
      • このような手法は客観的に思える
  • BigData、 AI、IoT
    • Corporate Value
    • Customer Value:Autonomous Drive、Connected Car, Electric Car, MaaS
  • 箱根登って下るとバッテリー回復するので沼津まで行ける
  • シミュレーションではなく過去のリーフのデータから、どこまで行けるか示せる
  • ProPilot 2.0
    • マジでやばい
    • ハンドルから手を離して良いと言っている(高速に限る)
    • ドライバーを監視している、前方を見ている限り
    • Adaptive Cruise Control
      • 単眼可視光カメラで前方の車との距離を測る
      • ナンバープレートのサイズがわかっているのでカメラの画素数から車間距離がわかる
  • AI
    • データ型
      • 言語:✖︎
      • 時間:▲
      • 空間:●
    • 能力
      • 生成:次元拡張
      • 認識次元圧縮
    • 上の軸をマトリックスにする
  • オイルフィルターの真贋判定
  • 空気抵抗を計算で求める
    • スパコンでも数日かかる
    • CADとCdの関係性
    • Voxel dataとCdの関係性を学習
    • 精度はまだまだ
      • R値で0.7(R2のこと?)
  • GANで車のデザインを生成
  • 画像の生成を三次元に拡張する
    • 内装とかもできる
    • フランスから来たインターンの学生
    • VRに展開
      • VRの中で白板に書くと目の前に現れる
    • Voxel Data との関係性がわかってるので、手書きで書いた絵の空気抵抗がわかる
  • データアナリスト
    • 問題のモデル化
  • 人材育成
    • 自動車業界、サプライヤーが持ってくる
    • うまくいってない
    • データを使っていて車に興味を持っている人に来てもらいたい
  • 自動車がなくなっていく?
    • MaaS、CASE
    • 痛し痒しと思っている
    • 体力のあるうちに全部やっておこう

企業ビッグデータから捉える企業活動と未来の活用可能性

帝国データバンク 中川みゆき 氏

  • 帝国データバンク保有するデータ
    • 調査員が現地に足を運んで入手した、ネットでは得られないデータベース
    • 100近い項目
    • 定量、定性
    • 中小企業が99%
    • 年間60万件超の調査
  • 活用事例

早稲田大学におけるデータサイエンス人材育成への取り組み

早稲田大学 松嶋敏泰 氏

  • 大隈重信統計学
    • 統計院を設立
      • 統計センターの前身
  • DSとは
    • Interdisciplinary Field
    • concept to unify
    • 日本だとSexiest Jobで広まったので、ビジネスの面が強い
    • Fourth paradigm of science
      • メタサイエンス
      • 知覚と思考の拡大
  • 早稲田のフォーカス領域
    • 専門性(応用先)
    • データサイエンス
  • 教育
    • 学内が中心
      • 理念
        • コンピュータにぶち込んで答えを出す、というのは×
        • 理論を押さえる
        • 座学だけでもダメ
        • 理論、専門、スキル
      • e-learningで実施
        • モジュール化されたコンテンツを組み合わせて
        • 様々なバックグラウンドに対応
    • 学外
      • 高校生
      • 社会人
      • 他大学
      • 企業
    • データサイエンスコンペティション
      • 参議院選
      • 精度がよかったチームは3人しか当落を外してなかった(最優秀ではない)
      • 政経のチームがPythonでコードを書いている

学校における「統計」教育の課題 我々は木に縁りて魚を求めてはいないか

代々木ゼミナール 西岡康夫 氏

  • 初等・中等教育の現状
  • サッカー、世界で活躍できるようになった理由は多くの子供がサッカーをやるようになったから
    • 棟梁レベルのDSを多く育てるなら裾野を広げないといけない
    • 国によってはエリート制を採用しているが、日本のカルチャーに合わないのでは
  • 教育改革
    • センター試験の廃止
      • 混沌の極み
      • 大学入学選抜は暗記、再生
        • ドーリットル、甚兵衛、小ピピン
      • 理系尚もて数学す、況や文系をや
    • 統計が入ってきた
      • 産業界の要請
        • 教員研修や教職課程のカリキュラム改革
          • まったく進んでない
          • 統計のイロハのイもダメ
    • 若者の強い自己否定感
      • 大学入試結果による序列付けが原因であるとされた
        • できません、を前提に話はじめる
        • ダニングクルーガー効果
          • できない人に限って自己肯定感が強い
      • 生産性
        • AIからIA
          • Intelligence amplifier
          • AIをうまく使う
      • 必要とされる理由
        • 第四次産業革命
          • アロ・ポイエティック(他者制作的)
          • オート・ポイエティック(自己制作的)
        • 論理国語、文学国語(非論理国語?)
        • STEM、STEAM教育
          • AIの苦手項目
            • 発問、小情報からの創意、独自指標の創出、規範の再構築、定義不明瞭
          • 数学の教師が統計を担当というのは物理の教師が生物を担当するようなもの
            • 演繹的推論、帰納的推論
          • やりたくない理由
            • 自信がないのでやりたくない
              • 統計的仮説検定を背理法として解説
            • 統計の諸概念をビルドアップ型で授業したいが、カリキュラムが対応してない
              • ネイピア数の存在が、積分の前に出てくる
              • 統計がベクトルを追い出した、と言われている
              • 理論的な扱いに深入りせず、という要領が出てる
            • そもそも生徒は統計は選択されないのでは
              • センターの2Bで選択問題、統計を選択するのは3.5%で、かつ平均点も低い
              • ベクトルが苦手な学生が選択?
              • 問題が長い
              • 捨象力は試されるので良い
            • 9/74校(2014年)しか推測統計を出さない

岩崎先生あいさつ

  • Technologyは変わっていくがPrincipleは変わらない

GAMをもう少し理解したい②

前回の続きです。 gam.fit() から。 前回記事はこちら。

ushi-goroshi.hatenablog.com

GAMの実装

gam.fit()

それでは gam.fit の中身を覗いてみましょう。

### gam.fit は x, y に加えて smooth.frame を受け取る。これは gam で作った mf で、中身は平滑化に関する情報を持った data.frame 
function (x, y, smooth.frame, weights = rep(1, nobs), start = NULL, 
    etastart = NULL, mustart = NULL, offset = rep(0, nobs), family = gaussian(), 
    control = gam.control()) 

始めに、受け取ったオプションをもとにデータのサイズ、 controlfamily に関するパラメータの取り出しを行います。

### 列名
ynames <- if (is.matrix(y)) 
    dimnames(y)[[1]]
else names(y)
xnames <- dimnames(x)[[2]]

### データの行数、列数
nobs <- NROW(y)
nvars <- ncol(x)

## その他の gam.control オプション
maxit <- control$maxit
bf.maxit <- control$bf.maxit
epsilon <- control$epsilon
bf.epsilon <- control$bf.epsilon
trace <- control$trace

### digits, weigths, offset
digits <- -log10(epsilon) + 1
if (is.null(weights)) 
    weights <- rep.int(1, nobs)
if (is.null(offset)) 
    offset <- rep.int(0, nobs)

### family に関するパラメータ
variance <- family$variance
dev.resids <- family$dev.resids
aic <- family$aic
linkinv <- family$linkinv
mu.eta <- family$mu.eta
if (!is.function(variance) || !is.function(linkinv)) 
    stop("illegal `family' argument")
valideta <- family$valideta
if (is.null(valideta)) 
    valideta <- function(eta) TRUE
validmu <- family$validmu
if (is.null(validmu)) 
    validmu <- function(mu) TRUE
eval(family$initialize)
if (is.null(mustart)) {
    eval(family$initialize)
}
else {
    mukeep <- mustart
    eval(family$initialize)
    mustart <- mukeep
}
### eta の初期値
eta <- if (!is.null(etastart)) 
    etastart

### エラーチェック
else if (!is.null(start)) 
    if (length(start) != nvars) 
        stop("Length of start should equal ", nvars, " and correspond to initial coefs for ", 
            deparse(xnames))
    else {
        coefold <- start
        offset + as.vector(if (NCOL(x) == 1) 
            x * start
        else x %*% start)
    }
else family$linkfun(mustart)

### mu の初期値
mu <- linkinv(eta)
if (!(validmu(mu) && valideta(eta))) 
    stop("Can't find valid starting values: please specify some")

### デビアンス(残差平方和)
new.dev <- sum(dev.resids(y, mu, weights))

ここまでは関数に渡されたオプションを元に処理を進めていますが、以降のプロセスで gam 特有の処理が行われています。まずは smoothers の取り出しです。

### smoothers が指定されている場合に smoother を抽出する
### 今回のケースでは s が入る
a <- attributes(attr(smooth.frame, "terms"))
smoothers <- a$specials

オブジェクトの属性を取り出すにあたり、 attr は属性名を指定する必要がありますが、 attributes では全ての属性をリスト形式で取り出します。 a というオブジェクトには、具体的には以下のようなリストが格納されます。

> a
$variables
list(wage, s(year, 4), s(age, 5), education)

$factors
           s(year, 4) s(age, 5) education
wage                0         0         0
s(year, 4)          1         0         0
s(age, 5)           0         1         0
education           0         0         1

$term.labels
[1] "s(year, 4)" "s(age, 5)"  "education" 

$specials
$specials$s
[1] 2 3

$specials$lo
NULL

$specials$random
NULL


$order
[1] 1 1 1

$intercept
[1] 1

$response
[1] 1

$class
[1] "terms"   "formula"

$.Environment
<environment: R_GlobalEnv>

$predvars
list(wage, s(year, 4), s(age, 5), education)

$dataClasses
      wage s(year, 4)  s(age, 5)  education 
 "numeric"  "numeric"  "numeric"   "factor" 

a$specials$s には 2 3 という数字が入っていますが、これは smooth.frame の2・3列目が平滑化の対象であるということだと思います。

今回は平滑化が含まれるので以下が処理されます:

if (length(smoothers) > 0) {
    ### NA ではない要素を取り出す
    smoothers <- smoothers[sapply(smoothers, length) > 0]
    
    ### 以下の処理はちょっとよくわからない
    for (i in seq(along = smoothers)) {
        tt <- smoothers[[i]]
        ff <- apply(a$factors[tt, , drop = FALSE], 2, any)
        smoothers[[i]] <- if (any(ff)) 
            seq(along = ff)[a$order == 1 & ff]
        else NULL
    }
}

smoothers に格納されている数値が 1 2 になりました。

> smoothers
$s
[1] 1 2

次に、以下の処理ではパラメータの推定に使われるエンジンを決めます。平滑化が含まれるケースでは general.wam または {s, lo}.wam が、そうでなければ lm.wfit が選ばれます。二種類以上の平滑化関数(s , lo, random)が指定されている、または s lo 以外の平滑化関数が指定されている場合には general.wam が選ばれるようです。

条件にしたがってエンジンを決めたあとは、後ろの工程で評価するための bf.call を作成しています。

### smoother が指定されている場合、ここが処理される
if (length(smoothers) > 0) {
    gam.wlist = gam.smoothers()$wlist
    smooth.labels <- a$term.labels[unlist(smoothers)]
    assignx <- attr(x, "assign")
    assignx <- assign.list(assignx, a$term.labels)
    which <- assignx[smooth.labels]

    ### 2つ以上の smoothers が指定されている場合は general.wam が指定される
    ### wam は weighted additive model
    if (length(smoothers) > 1) 
        bf <- "general.wam"

    ### 今回のケース(平滑化の関数として s を指定)はこちらで、 s.wam が指定される。
    else {
        sbf <- match(names(smoothers), gam.wlist, FALSE)
        bf <- if (sbf) 
            paste(gam.wlist[sbf], "wam", sep = ".")
        else "general.wam"
    }

    ### bf にオプション部分を文字列で結合
    bf.call <- parse(text = paste(bf, "(x, z, wz, fit$smooth, which, fit$smooth.frame,bf.maxit,bf.epsilon, trace)", 
        sep = ""))[[1]]
    s <- matrix(0, length(y), length(which))
    dimnames(s) <- list(names(y), names(which))
    fit <- list(smooth = s, smooth.frame = smooth.frame)
}
### 平滑化しない場合、通常の lm.wfit に渡す。なお general.wam でも lm.wfit が使われる。
else {
    bf.call <- expression(lm.wfit(x, z, wz, method = "qr", 
        singular.ok = TRUE))
    bf <- "lm.wfit"
}

なお上記コードブロックの最後に lm.wfit が使われていますが、LMやGLMの中身がどうなっているのかについては過去記事を参照してください。

ushi-goroshi.hatenablog.com

bf.call が作られたので、反復しながら解を求めにいきます。ここからが本体となるブロックです。

### デビアンス(ループの打ち切り基準)
old.dev <- 10 * new.dev + 10

### ここから反復に入る
for (iter in 1:maxit) {

    ### weight が 0 のデータは除外する
    good <- weights > 0
    varmu <- variance(mu)
    if (any(is.na(varmu[good]))) 
        stop("NAs in V(mu)")
    if (any(varmu[good] == 0)) 
        stop("0s in V(mu)")
    mu.eta.val <- mu.eta(eta)
    if (any(is.na(mu.eta.val[good]))) 
        stop("NAs in d(mu)/d(eta)")
    good <- (weights > 0) & (mu.eta.val != 0)
    
    ### z を生成。ただし bf.call で二番目の引数は y なので、この z は目的変数の意味
    z <- eta - offset
    z[good] <- z[good] + (y - mu)[good]/mu.eta.val[good]

    ### wz を生成。重み。今回のケースでは無視して良い
    wz <- weights
    wz[!good] <- 0
    wz[good] <- wz[good] * mu.eta.val[good]^2/varmu[good]

上のブロックでは weights を元に分析対象を定義し、 z および wz を生成しています。

ところで、今回のケースでは bf.call は以下のようになるのでした:

> bf.call
s.wam(x, z, wz, fit$smooth, which, fit$smooth.frame, bf.maxit, 
    bf.epsilon, trace)

z は第二引数となっていますが、後ほど紹介するようにこの s.wam という関数の中では二番目の引数は y なので、この z は目的変数となります(そのものではありませんが)。

続いて以下の処理が行われますが、これが gam.fit の本体となります:

    ### ここで bf.call が評価される。 s.wam の場合、bakfit が呼ばれる。
    ### bf.call で指定されている smooth.frame はこの時点では単なる data.frame 
    fit <- eval(bf.call)

eval によって bf.call が評価されるので、 s.wam が実行されます。 s.wam を見に行くまえに、このループ内での残りの処理を簡単に見ておきましょう。

    ### 予測値にオフセットを加算する
    eta <- fit$fitted.values + offset

    ### eta から mu に変換する
    mu <- linkinv(eta)

    ### デビアンスを更新
    old.dev <- new.dev
    new.dev <- sum(dev.resids(y, mu, weights))

    ### ここの trace は対角和ではなく、 gam のオプションでループごとのデビアンスをモニターできる
    if (trace) 
        cat("GAM ", bf, " loop ", iter, ": deviance = ", 
            format(round(new.dev, digits)), " \n", sep = "")

    ### ループの打ち切り判定
    ### デビアンスが NA となった場合、警告を出して打ち切る
    if (is.na(new.dev)) {
        one.more <- FALSE
        warning("iterations terminated prematurely because of singularities")
    }
    ### 差が十分に小さければ終了
    else one.more <- abs(old.dev - new.dev)/(old.dev + 0.1) > 
        epsilon
    if (!one.more) 
        break
}

評価結果を受けとり、デビアンスをもとにループを継続するかを判定しています。 gam.fit の残りの工程は一旦無視して、 s.wam を見に行きましょう。また次回。

GAMをもう少し理解したい

とても久しぶりの更新です。

背景

業務でモデリングを行うとき、私は大抵の場合GLMから始めます。目的変数に合わせて柔軟に分布を選択することが可能で、回帰係数という極めて解釈性の高い結果を得ることができるというのが理由です。

一方でGLMを使っていて不満に感じることの一つが、( \eta の世界で)非線形な効果を表現できないという点です。もちろん2次・3次の項や交互作用項を追加することである程度それらの不満は解消できるのですが、もう少しデータからそれらの特徴を学習したいと思うことがあります。

今回取り上げる一般化加法モデル(Generalized Additive Model, GAM)は、そのような複雑な関連性を表現できるよう説明変数に非線形な変換を行うもので、GLMを拡張したものとなります。ちょっと古いですが、2015年にMicrosoft RearchがKDDに出した論文(PDF)では、GAMを指して「the gold standard for intelligibility when low-dimensional terms are considered」と言っており、解釈性を保ちつつ高い予測精度を得ることができるモデルとしています。なおこの論文ではGAMに2次までの交互作用を追加するGA2Mという手法を提案しています。

このGAMの実装について調べた内容を書き留めておきます。GAMがどういうものかとか、平滑化に関する説明は、他に良いページがありますのでそちらを参照してください。例えば:

GAMの実行結果

始めに、GAMを使うとどのような結果を得ることができるのか確認しましょう。ちなみに検証した環境は以下の通りです。

> sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] compiler_3.6.0 tools_3.6.0    knitr_1.23     xfun_0.7 

GAMの実装としてRでは {mgcv}が使われることが多いようですが、今回は「Rによる統計的学習入門」を参考に{gam}を使用しました。ちなみに、この本は統計・機械学習の主要な手法を網羅的に押さえつつ章末にRでの実行方法が紹介されており、大変勉強になる良書です。

Rによる 統計的学習入門

Rによる 統計的学習入門

  • 作者: Gareth James,Daniela Witten,Trevor Hastie,Robert Tibshirani,落海浩,首藤信通
  • 出版社/メーカー: 朝倉書店
  • 発売日: 2018/08/03
  • メディア: 単行本(ソフトカバー)
  • この商品を含むブログ (1件) を見る

同書で用いているサンプルデータ Wage を使用するため、{ISLR}を同時に読み込みます。この{ISLR}パッケージは上記の本の原著であるIntroduction of Statistical Learning with Applications in Rから来ているようです。またこのデータは、アメリカの大西洋岸中央部における男性3000人の賃金、および年齢や婚姻状況、人種や学歴などの属性が記録されています。

library(gam)
library(ISLR)

データの中身を見てみましょう。

> head(Wage)
       year age           maritl     race       education             region       jobclass
231655 2006  18 1. Never Married 1. White    1. < HS Grad 2. Middle Atlantic  1. Industrial
86582  2004  24 1. Never Married 1. White 4. College Grad 2. Middle Atlantic 2. Information
161300 2003  45       2. Married 1. White 3. Some College 2. Middle Atlantic  1. Industrial
155159 2003  43       2. Married 3. Asian 4. College Grad 2. Middle Atlantic 2. Information
11443  2005  50      4. Divorced 1. White      2. HS Grad 2. Middle Atlantic 2. Information
376662 2008  54       2. Married 1. White 4. College Grad 2. Middle Atlantic 2. Information
               health health_ins  logwage      wage
231655      1. <=Good      2. No 4.318063  75.04315
86582  2. >=Very Good      2. No 4.255273  70.47602
161300      1. <=Good     1. Yes 4.875061 130.98218
155159 2. >=Very Good     1. Yes 5.041393 154.68529
11443       1. <=Good     1. Yes 4.318063  75.04315
376662 2. >=Very Good     1. Yes 4.845098 127.11574

このデータを用いて早速フィッティングしてみましょう。 glm と同様に関数 gam でモデルを当てはめることができます。

res_gam <- gam(wage ~ s(year, 4) + s(age, 5) + education, data = Wage)

ここで s(age, 5) は、説明変数 year を平滑化した上でモデルに取り込むことを意味し、5は平滑化の自由度です。なお s() の時点ではまだ平滑化は行われておらず、平滑化に必要な情報を属性として付与しているだけのようです。

> head(s(Wage$age, 5))
[1] 18 24 45 43 50 54
attr(,"spar")
[1] 1
attr(,"df")
[1] 5
attr(,"call")
gam.s(data[["s(Wage$age, 5)"]], z, w, spar = 1, df = 5)
attr(,"class")
[1] "smooth"

ではフィッティングした結果を見てみましょう。

> summary(res_gam)

Call: gam(formula = wage ~ s(year, 4) + s(age, 5) + education, data = Wage)
Deviance Residuals:
    Min      1Q  Median      3Q     Max 
-119.43  -19.70   -3.33   14.17  213.48 

(Dispersion Parameter for gaussian family taken to be 1235.69)

    Null Deviance: 5222086 on 2999 degrees of freedom
Residual Deviance: 3689770 on 2986 degrees of freedom
AIC: 29887.75 

Number of Local Scoring Iterations: 2 

Anova for Parametric Effects
             Df  Sum Sq Mean Sq F value    Pr(>F)    
s(year, 4)    1   27162   27162  21.981 2.877e-06 ***
s(age, 5)     1  195338  195338 158.081 < 2.2e-16 ***
education     4 1069726  267432 216.423 < 2.2e-16 ***
Residuals  2986 3689770    1236                      
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Anova for Nonparametric Effects
            Npar Df Npar F  Pr(F)    
(Intercept)                          
s(year, 4)        3  1.086 0.3537    
s(age, 5)         4 32.380 <2e-16 ***
education                            
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

パッと見ると、 glm を当てはめたときと同様の結果が得られるようです。実際、 gam の本体(と思われる) gam.fit の中では stats::lm.wfit が呼ばれます。2019/10/17修正 今回のケースでは lm.wfit は使われないので正しくありませんでした。

しかし、 glm の結果と大きく異なる点として、各説明変数の回帰係数が出ていません。これはどういうことでしょう。以下のように説明変数の効果をプロットしてみます。

par(mfrow = c(1, 3))
plot(res_gam, se = TRUE, col = "blue")

f:id:ushi-goroshi:20191016112558p:plain

真ん中の age が分かりやすいですが、この変数は34まではy軸の値が0未満となっており、 age は( yeareducation を固定した下で)35歳程度まで wage の平均を下回っていることがわかります。その後45歳をピークに減少傾向に入り、65歳を過ぎると再び平均を下回るようになり、以降は急激に減少していきます。このような現象は、年齢とともに役職が上がることで賃金が増加し、退職によって減少することを考えれば非常に納得感のあるものだと思います。

なお上のプロットに必要なデータは以下のように取れるので、説明変数の各点における目的変数に対する影響を数値でも確認できます(必要な関数がエクスポートされていないので gam::: で直接呼び出しています)。

tmp <- gam:::preplot.Gam(res_gam, terms = gam:::labels.Gam(res_gam))
age_sm <- cbind(tmp$`s(age, 5)`$x, tmp$`s(age, 5)`$y) 
age_sm_uniq <- unique(age_sm[order(age_sm[, 1]), ])

プロットしてみましょう。同じ線が描けます。

plot(age_sm_uniq, type = "l", xlab = "age", ylab = "s(age, 5)")

f:id:ushi-goroshi:20191016113801p:plain

さて、上記の結果は例えば age の二次の項を含めることで lm でも再現できるかもしれません。例えば以下のようになります:

### lmでフィッティング
res_lm <- lm(wage ~ year + poly(age, 2) + education, Wage) ## poly(age, 2)で二次の多項式とする

### 予測用のデータ作成。ageだけを変化させ、yearとeducationは固定する。
x_lm <- seq(min(Wage$age), max(Wage$age), length.out = 100)
nd <- data.frame(year = rep(2003, 100),
               age = x_lm,
               education = rep("2. HS Grad"))
prd_lm <- predict(res_lm, nd, se.fit = T)

### 予測値から平均を引いてからプロット
prd_m <- 90
plot(x_lm, prd_lm$fit - prd_m, type = "l", col = "blue", ylim = c(-40, 10))
lines(x_lm, prd_lm$fit - prd_m + 2*prd_lm$se.fit, lty = "dashed")
lines(x_lm, prd_lm$fit - prd_m - 2*prd_lm$se.fit, lty = "dashed")

f:id:ushi-goroshi:20191016114041p:plain

大体同じようなプロットを作成する事ができました。しかしピークを過ぎてからの緩やかな減少は表現できていませんし、一つ一つの変数について何次までの多項式を含めるかを検討していくのは少し手間がかかります。GAMを使えばデータに存在する細やかな変化を自動的に捉えることができます(もちろん代償もあります)。

GAMの実装

gam()

それでは gam がどのようにフィッティングを行っているのかを見ていきましょう。本体は最後の方に出てくる gam.fit なのですが、途中も少し細かく追ってみます。コンソールで gam を実行すると、以下のように関数の中身を見ることができます。

まず以下のブロックでは、 gam にオプションとして指定した内容に沿った処理を実行しています。

function (formula, family = gaussian, data, weights, subset, 
        na.action, start = NULL, etastart, mustart, control = gam.control(...), 
        model = TRUE, method = "glm.fit", x = FALSE, y = TRUE, ...) 
{
### 関数の引数を名前付きで確定。gam(wage ~ s(year, 4) + education, Wage)として与えた場合、
### formula = と data = がそれぞれ保持される。
### match.call returns a call in which all of the specified arguments are specified by their full names.
call <- match.call()

### familyの判定
if (is.character(family)) 
  family <- get(family, mode = "function", envir = parent.frame())
if (is.function(family)) 
  family <- family()
if (is.null(family$family)) {
  print(family)
  stop("`family' not recognized")
}

### データが指定されていない場合
if (missing(data)) 
  data <- environment(formula)

### 指定されている引数の取り出し
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "weights", "etastart", 
             "mustart", "offset"), names(mf), 0L)
mf <- mf[c(1L, m)]

### 指定されていない引数を指定し、stats::model.frame()の形式に仕立てる
mf$na.action = quote(na.pass)
mf$drop.unused.levels <- TRUE
mf[[1L]] <- quote(stats::model.frame)

次に、ここで一つ gam ならではの処理として平滑化に使う関数を取り出しています。

### 平滑化の関数を取ってくる(s, lo, random)
gam.slist <- gam.smoothers()$slist

gam.smoothers として新しい平滑化関数を指定することも出来るようですが、デフォルトは slo および random で、それぞれ 平滑化スプライン局所回帰ランダム効果としての指定を意味しているようです。 3つめの random がわからなかったので調べてみたところ、これはカテゴリ変数に対する指定で、パラメータの推定においていわゆる縮小推定を行なうもののようでした。

https://www.rdocumentation.org/packages/gam/versions/1.16.1/topics/random

これらに接頭語として gam を加えた(e.g. gam.s)関数が平滑化のための関数として実行されます。 gam.s については後述します。

次のブロックですが、 call クラスであった mf を評価することで data.frame に変換しています。

### term を mf$formula に渡す
mt <- if (missing(data)) 
  terms(formula, gam.slist)
else terms(formula, gam.slist, data = data)
mf$formula <- mt

### ここで mf 、つまり model.frame が実行されて data.frame になる
### ただし平滑化は実行されず、平滑化のパラメータは attribute として持っている
mf <- eval(mf, parent.frame())
if (missing(na.action)) {
  naa = getOption("na.action", "na.fail")
  na.action = get(naa)
}
mf = na.action(mf)
mt = attributes(mf)[["terms"]]

ここまでは mfcall クラス、すなわち未評価の関数およびその引数を要素に持つオブジェクトでした。それが eval で評価されたため stats::model.frame が実行され、 formula にしたがい data.frame が生成された、という流れのようです(間違っていたらすみません)。

### method の指定によって処理を分ける。 glm.fit または glm.fit.null 以外の場合はエラー
switch(method, model.frame = return(mf), glm.fit = 1, glm.fit.null = 1, 
       stop("invalid `method': ", method))

ここでは method によって処理を分けていますが、現状は glm.fit または model.frame のみ受け付けているようです。 なおhelpを参照すると、 model.frame を指定した場合フィッティングは行われないようですね。以下のようになります。

> gam(wage ~ s(year, 4) + s(age, 5) + education, data = Wage, method = "model.frame")
wage s(year, 4) s(age, 5)          education
231655  75.04315       2006        18       1. < HS Grad
86582   70.47602       2004        24    4. College Grad
161300 130.98218       2003        45    3. Some College
155159 154.68529       2003        43    4. College Grad
11443   75.04315       2005        50         2. HS Grad

この通りデータが返ってきます。

以下では、YおよびXをそれぞれ抽出し、さらにフィッティングに必要なオプションを指定しています。

### Y を取り出す
Y <- model.response(mf, "any")

### X を matrix で取り出す。
### gam を実行したときのエラーメッセージ( `non-list contrasts argument ignored` )はここで出ている。 contrasts の指定が良くない様子。
X <- if (!is.empty.model(mt)) 
  model.matrix(mt, mf, contrasts) 
else matrix(, NROW(Y), 0)

### その他パラメータ(weights, offset, mustart, etastart)
weights <- model.weights(mf)
offset <- model.offset(mf)
if (!is.null(weights) && any(weights < 0)) 
  stop("Negative wts not allowed")
if (!is.null(offset) && length(offset) != NROW(Y)) 
  stop("Number of offsets is ", length(offset), ", should equal ", 
       NROW(Y), " (number of observations)")
mustart <- model.extract(mf, "mustart")
etastart <- model.extract(mf, "etastart")

以上で準備が完了し、 gam.fit で当てはめを行います。

### ここが本体。 gam.fit を呼び出している
fit <- gam.fit(x = X, y = Y, smooth.frame = mf, weights = weights, 
               start = start, etastart = etastart, mustart = mustart, 
               offset = offset, family = family, control = control)

いったん後続の部分は無視して gam.fit に移りたいと思いますが、長くなったので今回はここまでにします。