統計コンサルの議事メモ

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

GLMをもう少し理解したい①

背景

一般化線形モデル(GLM)は、一般に線形回帰モデルを正規分布を含む指数分布族に拡張したものだと捉えられています。アイディアとしてはシンプルである割に非常に有用で、GLMによって

  • 整数値(ポアソン回帰)
  • 二値(ロジスティック回帰)
  • 0〜1の実数(ベータ回帰) ※2020/4/8追記 ベータ回帰は一般的でなかったので消しておきます
  • 0以上の実数(ガンマ回帰) ※2020/4/8追記 代わりにガンマ回帰を追加

などを扱うことができ、しかも回帰係数という非常に解釈性の高い結果を得ることができます1

そんなGLMですが、よく使う割には内容を今ひとつ理解できていないなと思うことがあったので、もう少しだけGLMを理解したいと思いRのglmの中身を見てみました。その内容をメモしておきます。

ちなみにこの検証を行っている環境は以下の通りです:

> sessionInfo()
R version 3.3.3 (2017-03-06)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS  10.13.3

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] tools_3.3.3  yaml_2.1.13  knitr_1.15.1

glm

まずはRのglmがどのように定義されているかを見てみましょう。コンソールでglmと入力することで、以下のようにglmという関数の定義を見ることができます。

> glm
function (formula, family = gaussian, data, weights, subset, 
          na.action, start = NULL, etastart, mustart, offset, control = list(...), 
          model = TRUE, method = "glm.fit", x = FALSE, y = TRUE, contrasts = NULL, 
          ...) 

まずはここでglmに渡す引数を定義しています。これらの引数でよく使われるのはformulafamilydataでしょうか。それぞれglmに渡す線形予測子(数式)、Yの従う分布、モデリングに用いるデータを指定しています。その他、データポイントの一つ一つの重みを変えたい場合にはweights、データの一部を使用する場合にはsubsetを指定したりします。

関数の定義は以下より始まりますが、細かい話は飛ばしてglmの本体に向かいましょう。

{
   
   call <- match.call()
   if (is.character(family)) 
      family <- get(family, mode = "function", envir = parent.frame())   
   
   ##
   ## 中略
   ##
   
   ## 本体はココのようです
   fit <- eval(call(if(is.function(method)) "method" else method,
                    x = X, y = Y, weights = weights, start = start, etastart = etastart, 
                    mustart = mustart, offset = offset, family = family, 
                    control = control, intercept = attr(mt, "intercept") > 0L))

glmの本体と呼べそうな部分はどうやらこのfitを定義している部分です。最初のevalは与えられた文字列をスクリプトとして解釈するための関数なので、call以降を実行するようです。またcallここによると「与えられた名前の関数の、与えられた引数への適応からなる未評価の表現式である」とのことなので、callに続くmethodおよび残りの引数がmethod(...)の形でevalに与えられ、関数として評価されます。

つまり

eval(call(if(is.function(method)) "method" else method, ...

method(...)

と同じとなるはずで、以下の例では同じように動いていることが確認できました。

# 関数を定義
return_cube <- function(x) x^3

# 普通に呼び出す
> return_cube(3)
[1] 27

# eval(call(...))で呼び出す
> eval(call("return_cube", x = 3))
[1] 27

# match.call
> eval(match.call(return_cube, call("return_cube", x = 3)))
[1] 27

> eval(match.call(return_cube, call("return_cube", 3)))
[1] 27

さて、callで指定しているmethodglmの引数で指定されているものでしたが、デフォルトではglm.fitが入力されています。したがってmethod(...)glm.fit(...)となるはずです。そこで今度はglm.fitの定義を確認してみましょう。

glm.fit

glmと同じく、glm.fitについてもコンソールに直接打ち込むことで関数の定義を表示することができます。まずは引数から見てみましょう。

> glm.fit
function (x, y, weights = rep(1, nobs), start = NULL, etastart = NULL, 
    mustart = NULL, offset = rep(0, nobs), family = gaussian(), 
    control = list(), intercept = TRUE) 

xyはそれぞれ説明変数と目的変数を指定し、その他の引数はglmから引き継がれるようですね。

またメインとなるのは以下のループ部分のようです。

   control <- do.call("glm.control", control)
   x <- as.matrix(x)
   xnames <- dimnames(x)[[2L]]

   ##
   ## 中略
   ##

      for (iter in 1L:control$maxit) {

         ##
         ## 中略
         ##

         z <- (eta - offset)[good] + (y - mu)[good]/mu.eta.val[good]
         w <- sqrt((weights[good] * mu.eta.val[good]^2)/variance(mu)[good])
         fit <- .Call(C_Cdqrls, x[good, , drop = FALSE] * 
             w, z * w, min(1e-07, control$epsilon/1000), check = FALSE)
          
         ##
         ## 中略
         ##
      }

   ##
   ## 中略
   ##
      

do.call(...)に渡すcontrollist()なので、do.call("glm.control", list())を実行すると以下が返ります;

> do.call("glm.control", list())
$epsilon
[1] 1e-08

$maxit
[1] 25

$trace
[1] FALSE

maxitが25なので、このループは最大で25回実行されます。ではこのループ内で何が行われているかというと、zwを新たに定義し、それをxおよび互いに乗じた形でC_Cdqrlsに渡しています。また.CallはCで書かれたルーチンを呼び出すための関数なので、ここではC_Cdqrlsという関数にx * wz * wといった引数を渡しているようです。

ではこのC_Cdqrlsはどこにあるのでしょうか?今度はC_Cdqrlsを探してみましょう。

C_Cdqrls

実はこのC_Cdqrlsstatsパッケージの関数として定義されています。しかしエクスポートされていないため、そのままコンソールに打ち込んでも表示されません。そのような場合には:::を使います。

> stats:::C_Cdqrls
$name
[1] "Cdqrls"

$address
<pointer: 0x101a2cdd0>
attr(,"class")
[1] "RegisteredNativeSymbol"

$dll
DLL name: stats
Filename:
         /Library/Frameworks/R.framework/Versions/3.3/Resources/library/stats/libs/stats.so
Dynamic lookup: FALSE

$numParameters
[1] 4

attr(,"class")
[1] "CallRoutine"      "NativeSymbolInfo"

しかしstats:::C_Cdqrlsと打ち込んでも、これまでと異なり関数の定義が表示されません。これは先ほど書いた通り、C_CdqrlsがCで書かれた関数であり、.callで呼び出されるためです。

ではどこから呼び出されるのかと言うと、私の環境では上記のFilenameで指定されている場所のようなのですが、これ自体は実行ファイル(stats.so)となっていてソースが見当たりません。それじゃどこにあるのかということで色々とググってみたところ、どうやらここで見れそうです。ファイル名を見てわかる通り、これはlmを定義しているコードです。glmは深く潜っていくとlmにたどり着くようです

C_Cdqrlsを定義している部分を見てみると:

SEXP Cdqrls(SEXP x, SEXP y, SEXP tol, SEXP chk)
{
   SEXP ans;

   ###
   ### 中略
   ###

   work = (double *) R_alloc(2 * p, sizeof(double));
   F77_CALL(dqrls)(REAL(qr), &n, &p, REAL(y), &ny, &rtol,
           REAL(coefficients), REAL(residuals), REAL(effects),
           &rank, INTEGER(pivot), REAL(qraux), work);
   SET_VECTOR_ELT(ans, 4, ScalarInteger(rank));
   for(int i = 0; i < p; i++)
    if(ip[i] != i+1) { pivoted = 1; break; }
   SET_VECTOR_ELT(ans, 8, ScalarLogical(pivoted));
   UNPROTECT(nprotect);

   return ans;
}

ここまで来ると何がなんだか私にはわかりませんが、returnansを返しているのでansを定義している箇所に着目すると、どうもF77_CALLが怪しい感じです。F77_CALLこのページによるとCからFortranを呼び出すための関数のようです。

F77_CALL

ではFortranで書かれたdqrlsソースコードはどこで見れるのかと言うと、ここのようです。重要そうなところを抜き出すと:

subroutine dqrls(x,n,p,y,ny,tol,b,rsd,qty,k,jpvt,qraux,work)
      integer n,p,ny,k,jpvt(p)
      double precision x(n,p),y(n,ny),tol,b(p,ny),rsd(n,ny),
     .                 qty(n,ny),qraux(p),work(p)
      integer info,j,jj,kk
      
      ### Householder transformation
      call dqrdc2(x,n,n,p,tol,k,qraux,jpvt,work)
      
      ### 
      if(k .gt. 0) then
         do 20 jj=1,ny
            call dqrsl(x,n,n,k,qraux,y(1,jj),rsd(1,jj),qty(1,jj),
     1           b(1,jj),rsd(1,jj),rsd(1,jj),1110,info)
   20       continue
      else
         do 35 i=1,n
            do 30 jj=1,ny
                rsd(i,jj) = y(i,jj)
   30       continue
   35   continue
      endif

となっており、dqrdc2dqrsl(紛らわしいけどdqrlsではない)を呼んでいます。これらはそれぞれ、

  • Householder変換を行う関数
  • そのアウトプットに対して加工および最小二乗解を与える関数

となっています。

随分かかりましたが、ここに来てようやく解を得ることができました。ここまでを振り返ると、glmという関数のコアの部分の役割はそれぞれ:

  1. Householder法による最小二乗解の推定(C_Cdqrls
  2. 上記の反復による収束判定(glm.fit
  3. もろもろの条件設定など(glm

となっているようでした。つまりglmは最尤法によって解を推定していると思われている(そしてそれは正しい)のだけれど、実際には最小二乗解をHouseholder法によって得ているのだということがわかりました。

長くなってしまったので、一旦切ります。次回はこの意味についてもう少し追いかけてみたいと思います。


  1. そのため個人的にはGLMをモデリングのベースラインとすることが多く、ここで十分な精度が得られるかでその後の対応を決めたりしています