統計コンサルの議事メモ

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

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 を見に行きましょう。また次回。