GLMをもう少し理解したい① - 統計コンサルの議事メモ

背景
glm
glm.fit
C_Cdqrls
F77_CALL

背景

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

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

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

そんな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に渡す引数を定義しています。これらの引数でよく使われるのはformula、family、dataでしょうか。それぞれ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で指定しているmethodはglmの引数で指定されているものでしたが、デフォルトでは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)

x、yはそれぞれ説明変数と目的変数を指定し、その他の引数は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(...)に渡すcontrolはlist()なので、do.call("glm.control", list())を実行すると以下が返ります；

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

$maxit
[1] 25

$trace
[1] FALSE

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

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

C_Cdqrls

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

> 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;
}

ここまで来ると何がなんだか私にはわかりませんが、returnでansを返しているので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