GAMをもう少し理解したい②
前回の続きです。 gam.fit()
から。
前回記事はこちら。
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())
始めに、受け取ったオプションをもとにデータのサイズ、 control
や family
に関するパラメータの取り出しを行います。
### 列名 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の中身がどうなっているのかについては過去記事を参照してください。
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
を見に行きましょう。また次回。