統計コンサルの議事メモ

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

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 はまた次回。