🔍

R の線形回帰 lm 関数の実装を辿る

2024/04/28に公開

この記事は？

この記事では、R にデフォルトで入っている stats パッケージの線形回帰関数 lm を題材に、R のパッケージ・関数がどのように実装されているかを辿っていく。

はじめに

R のパッケージ・関数の使い方は、? コマンドや Web 上でドキュメントを読めば一通り知ることができる。

しかし、時としてどうしてもドキュメントには書かれていないような詳細を知りたくなるケースに遭遇することがある。
このような時には、実装コードを読みにいくことで疑問が解決することが多い。

また、定評のあるパッケージの実装を知ることで、それまで知らなかった書き方や言語仕様に触れることができ、コードを書く際の参考にもなる。

そこで、この記事では線形回帰の lm 関数に着目して、実装の辿り方を紹介するとともに、R パッケージの実装がどのように行われているかを見ていく。

実装を辿るうえでの着目ポイント

実装を辿るとはいっても、細かく辿っていくととても一本の記事に収まりきらないので、今回は以下の3つのポイントについて見ていくことにする：

実装どのようにして辿ればよいか？
モデル推定はどのように実装されているか？
print 関数や summary 関数で推定結果を表示できるのはなぜか？

`lm` 関数の基本的な使い方

まず具体的な実装を読み始める前に、lm 関数をどのように使うかを紹介する。

サンプルデータとして、以下を用いる：

set.seed(11)

x <- c(30, 20, 15, 45, 25)
y <- 2 + 0.3 * x1 +  1.0 * rnorm(5)
df <- data.frame(x, y)

モデル作成

上記のデータをもとに、以下のようにして lm 関数を使ってモデルを作成する：

model <- lm(y ~ x, data=df)

入力データ点と予測結果の直線を表したのが、次の図である：

モデル情報の表示

上記で作成したモデル model を引数に取ってprint 関数ないし summary 関数を実行することで、モデルの情報を表示することができる。

print 関数では、モデル内容と係数・切片の値を確認できる：

print(model)

Call:
lm(formula = y ~ x, data = df)

Coefficients:
(Intercept)            x
     2.0893       0.2799

一方、summary 関数では上記に加えて残差や係数の標準誤差と検定結果など、より詳細な情報を表示することができる。

summary(model)

Call:
lm(formula = y ~ x, data = df)

Residuals:
       1        2        3        4        5
-0.07774  0.33903 -1.30455 -0.54808  1.59135

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.08928    1.56363   1.336    0.274
x            0.27991    0.05411   5.173    0.014 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.246 on 3 degrees of freedom
Multiple R-squared:  0.8992,	Adjusted R-squared:  0.8656
F-statistic: 26.76 on 1 and 3 DF,  p-value: 0.01402

lm 関数の実装をどうやって辿るか？

lm 関数の基本的な使い方がわかったところで、いよいよ実装を見ていこうと思う。

だが、そもそも lm 関数のソースコードはどこで読めるのだろうか？
そして、どうやって読めば良いだろうか？

まず、どこにあるかについてだが、多くのオープンソースのソフトウェアと同様に GitHub 上で公開されている。

しかし、GitHub 上でコードを頭から順に読んでいくのは苦行なので、RStudio 上で debug 実行しながら実装を理解する方法も紹介する。

GitHub のコードを読む

lm 関数は、stats という R にデフォルトで入っているパッケージに含まれており、その実装は R 自身のソースコードと同じ場所にある。

以下の GitHub の repository で R のソースコードを読むことができる：

https://github.com/wch/r-source

ただし、実際には R のソースコードは GitHub ではなく Subversion の以下の repository で管理・開発が行われており、上記の GitHub の repository はその read-only のミラーとのことである：^[1]

https://svn.r-project.org/R/

GitHub の方でコードを探すと、stats パッケージは以下の directory に見つかる：

そして今回見たい lm 関数は、以下のファイルで定義されている。
[https://github.com/wch/r-source/blob/master/src/library/stats/R/lm.R]

根気があればこれを愚直に読んでいけば実装が理解できる。

RStudio から debug 実行しながら処理を辿る

GitHub 等のソースコードでは実装全体が見えるものの、実行時に関数内部でどのような処理が行われているか分かりづらいことも多い。

例えば、以下は lm 関数の冒頭である：

これだけ読んでもよほど慣れている人でないとどのような処理を行っているかわからないであろう。

そこで、このような時は RStudio で Debug モードで実行してみることをオススメする。

Debug モードで実行すると、実行の最中に変数にどのような値が入っているか1行1行立ち止まって確認することができる。

やり方は簡単。
まず、以下のようなコードを R Script として作成する（RMarkdown では不可）：

x <- c(30, 20, 15, 45, 25)
y <- 2 + 0.3 * x1 + 1.0 * rnorm(5)

df <- data.frame(x, y)

plot(x, y)

model <- lm(y ~ x, data = df)
summary(model)

エディタ上で実行を止めたい行の行番号あたりをクリックすると赤い印がつくので、その状態でコードを実行すればよい。

すると、赤い印を付けた行にカーソルが現れハイライトされるとともに、console に以下のような表示が出る：

この状態ではハイライトされた行を実行する直前で止まっており、console からスコープ内で定義されている変数を入力するとその中身を表示することができるほか、特殊なコマンドを実行することで debug 箇所を移動できる。
コマンドを含めた debug の詳細は以下を参照いただきたい：

今回は lm 関数の中で行われている処理を知りたいので、lm 関数を呼び出している行にカーソルが置かれた状態で console から関数の中に入るコマンド s を実行する。
すると、表示されているコードが lm 関数のものに切り替わる。

この状態で console で n を入力すると、1行1行停止しながら実行していくことができる。

例えば8行目（GitHub のコードでは28行目）の mf <- match.call(expand.dots = FALSE) という処理で変数 mf に何が入っているかを知りたい場合は、コマンド n で9行目まで進んだうえで console に mf と入力すればよい。
すると、以下のように表示される：

Browse[3]> mf
lm(formula = y ~ x, data = df)

なんとも不思議な結果が表示されたが、そもそも match.call 関数とはなんだろうか？
このような時は、console に ?match.call と入力し Help を表示させるとよい。

どうやら match.call 関数は現在実行中の関数とその引数を引数名の省略なしで表した call クラスオブジェクトというものを返すらしい。

なお、call クラスオブジェクトは、eval 関数の引数にとることによって実行できるとのこと。

同様の手順で1行1行実行しながら確かめていくと、最終的に lm 関数実行時に引数として指定した内容のうち formula, data, subset, weights, na.action, offset の内容を取り出して stats の model.frame にまとめていることがわかった。
15行目（GitHub 上では36行目）までカーソルを進めてから console で mf の中身を表示した結果が以下である：

Browse[3]> mf
          y  x
1 10.408969 30
2  8.026594 20
3  4.983447 15
4 14.137347 45
5 10.678489 25

一通り調べ終わって debug モードを終了したい場合は、console からコマンド Q を実行すればよい。

最小二乗法によるモデル推定の実装

一通り実装の辿り方がわかったところで、lm 関数の中核となるモデル推定をどのように実装しているか見ていく。

lm 関数による線形回帰では最小二乗法によって係数・定数項などのパラメータ推定を行うが、その際に最適化手法として行列の QR 分解というものを用いている。

以降では理論について軽く触れたうえで、コード上ではどこで最小二乗法による係数推定を実行しているかを見ていく。

理論の復習

実装を辿る前に、線形回帰の理論の復習を簡単に行う。

以下のような $p-1$ 個の説明変数を用いた線形回帰モデルを考える：

y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... +\beta_{p-1} x_{p-1} + \varepsilon.

$\varepsilon$ は期待値が0で独立同時な分布に従う誤差項である。

ここで、 $n$ 組の目的変数 $y_i$ と説明変数 $x_{1i}, x_{2i}, ..., x_{p-1,i}$ からなるデータセット( $i=1, 2, ..., n$ )が与えられたとき、定数項と回帰係数の推定量 $\hat{\boldsymbol{\beta}} = (\hat{\beta}_0, \hat{\beta}_1, \hat{\beta}_2, ..., \hat{\beta}_{p-1})^\top$ は最小二乗法により以下のように求められる：^[2]

\hat{\boldsymbol{\beta}} = \left( \boldsymbol{X}^\top \boldsymbol{X}\right)^{-1} \boldsymbol{X}^\top \boldsymbol{Y}.

$\boldsymbol{Y} = (y_1, y_2, ..., y_n)^\top$ は目的変数ベクトル、 $\boldsymbol{X}$ は計画行列といって定数項に対応するダミー変数と説明変数から成る $n \times p$ 行列である：

\boldsymbol{X} = \begin{pmatrix} 1 & x_{11} & x_{21} & ... & x_{p-1,1}\\ 1 & x_{12} & x_{22} & ... & x_{p-1,2}\\ \vdots & \vdots & \vdots & ... & \vdots\\ 1 & x_{1n} & x_{2n} & ... & x_{p-1,n}\\ \end{pmatrix}.

ただし、数値計算上では上記の推定式を直接使わず、計画行列 $\boldsymbol{X}$ を直交行列 $\boldsymbol{Q}$ と上三角行列 $\boldsymbol{R}$ の積の形 $\boldsymbol{X} = \boldsymbol{Q}\boldsymbol{R}$ に分解してから計算を行う。
このような行列分解を、QR 分解と呼ぶ。

これを用いると、先ほどの回帰推定量は以下の形に書き表すことができる：

\hat{\boldsymbol{\beta}} = \boldsymbol{R}^{-1} \boldsymbol{Q}^\top \boldsymbol{Y}

詳細は以下の記事を参照いただきたい：

どこで係数推定をおこなっているか？

上記で述べたような最小二乗法による係数推定は、lm 関数では FORTRAN で実装されたプログラムにより実行している。

といっても、lm 関数から直接 FORTRAN によるプログラムを呼び出しているわけではなく、まずは推定部分を切り出した lm.fit 関数を呼び、その中から Cdqrls という C言語で書かれたプログラムを呼び出していた。これは、FORTRAN で書かれた dqrls をラップしたものである。

さらに dqrls は内部で2つの別の FORTRAN で書かれたサブルーチンを呼び出していた。
一つ目は dqrdc2 といって入力された計画行列を QR 分解するものである。
二つ目は dqrsl といって dqrdc2 で QR 分解した結果を用いて最小二乗法による回帰係数を算出するものである。

以下、順に見ていく。

`lm.fit`

lm 関数内で入力引数のパース・整形を行ったあと、以下の箇所で lm.fit 関数ないし lm.wfit 関数が呼ばれ、推定の処理が行われる：

変数 w はサンプル重みであり、引数 weights から来ている。
サンプル重みを指定していれば lm.wfit が呼ばれ、指定していなければ lm.fit が呼ばれる。

以下、lm.fit に絞って見ていく。

lm.fit 関数は先ほどまでと同じファイルの以下の箇所で定義されている：

引数のうち x は計画行列、y は目的変数である。

中核となる係数推定の処理は以下の部分で C 言語で書かれたプログラムを呼び出して行っている：^[3]

`Cdqrls`

上記で呼び出される関数は、以下のファイルにて Cdqrls として定義されている：
[https://github.com/wch/r-source/blob/master/src/library/stats/src/lm.c]

コメントによるとこちらの関数は FORTRAN の dqrls というサブルーチンをメモリ消費を抑えるために C 言語でラップしたものとのことで、以下の箇所で dqrls を呼び出している：

`dqrls`, `dqrdc2`, `dqrsl`

dqrls は線形一次方程式の最小二乗解を求めるためのサブルーチンである。
[https://github.com/wch/r-source/blob/master/src/appl/dqrls.f]

この中で、dqrdc2 と dqrsl という別のサブルーチンを呼び出している：

dqrdc2 は、Householder reflection というアルゴリズムを用いて入力された計画行列を QR 分解するものである:
[https://github.com/wch/r-source/blob/master/src/appl/dqrdc2.f]

そして dqrsl では、QR分解を行った結果をもとに回帰係数を求める:
[https://github.com/wch/r-source/blob/master/src/appl/dqrsl.f]

dqrdc2, dqrsl は LINPACK という線形方程式を解くためのパッケージ群に由来している：

https://www.netlib.org/linpack/

なお、 dqrdc2 は LINPACK にある dqrdc そのものではなく、lm 関数用に改変したものである。
以前の記事でも触れたように、lm 関数では入力した説明変数に完全な多重共線性がある場合は説明変数を一部除外して推定を行う。
この変数選択の処理は dqrdc2 で計画行列を QR 分解する際に行われており、この部分が LINPACK オリジナルの dqrdc と大きく異なる箇所の一つである。
こちらのアルゴリズムの詳細については、また別の機会に触れたい。

推定結果表示の実装

次に、推定結果を表示するための機能がどのように実装されているかを見ていく。

冒頭の基本的な使い方の箇所に書いたように、lm 関数で作成したモデル内容を print 関数や summary 関数で表示することができた。

一方で print 関数は lm の推定結果に限らず汎用的に変数などの中身を表示する際に使われる：

a <- "hoge"
print(a)

[1] "hoge"

また、 summary についても lm 以外にも例えば glm（一般化線形モデル）のようにほかのモデル推定結果を表示することができる。

このような処理を行うためには、Python であれば class を定義したうえで特殊 method __str__() や .summary() のような要約用 method を実装することになるだろう。
では、R の場合にはどのように実装しているのだろうか？

結論から言ってしまうと、print 関数と summary 関数を実行した際の表示内容は、それぞれ以下の print.lm 関数と summary.lm 関数で定義されている：

ではなぜ print 関数・summary 関数でこれらが呼び出されるかというと、lm 関数による推定結果オブジェクトには lm という class 属性が付与されており、それにより print.lm と summary.lm がそれぞれ print 関数と summary 関数の lm に対応する method として認識されるからである。

R における class と method

上記を理解するために、R において class と method がどのような仕様になっているかを簡単に紹介する。

R のオブジェクトには class 属性を付与することができる。
例えば、以下では変数 x に sample_class という class 属性を付与してみた：

x <- c(1,2,3)
# print(x)
# > [1] 1 2 3
class(x) <- "sample_class"
print(x)

[1] 1 2 3
attr(,"class")
[1] "sample_class"

一方で、上記で見た print 関数や summary 関数は、引数の class 属性に応じて振る舞いを変えることができる。
このような関数をジェネリック関数と呼ぶ。

ジェネリック関数を実行する際に呼び出される関数は、<ジェネリック関数名>.<class 属性名> という名前で定義されている。

例えば、以下のように print.sample_class という関数を定義して、そのうえで先ほどクラス属性 sample_class を付与した変数 x を print 関数に渡してみる：

print.sample_class <- function(x) {
  print("Hi! print.sample_class() is called!")
  invisible(x)
}

class(x) <- "sample_class"
print(x)

[1] "Hi! print.sample_class() is called!"

すると、上記のように print(x) を実行した際に新たに定義した print.sample_class が呼び出されることがわかる。

この時、このようにジェネリック関数から入力引数の class 属性に応じて呼び出される関数のことを method といい、例えば print.sample_class は、「print 関数の method」と呼ばれる。

なお、ジェネリック関数は、UseMethod という特別な関数を使って定義されている。
例えば、 print 関数の実装を見てみると、以下のようになっている。

上記に述べた R のクラスシステムは S3 と呼ばれている（AWSのオブジェクトストレージとは無関係）。
S3 システムの入門については、以下の書籍の第8章が優れており大変参考になった：
Garrett Grolemund 著（大橋真也監訳、長尾高弘訳）「RStudioではじめるRプログラミング入門」（O'REILLY, 2015）
https://amzn.asia/d/aE3Iuje

`print.lm` と `summary.lm`

S3 システムにおける class と method の扱いがわかったところで、 lm 関数における print 関数と summary 関数の話に戻る。

lm 関数の返り値には、以下の箇所で lm という class 属性が付与されている：

なお、変数 mlm は目的変数 y が複数カラムからなる場合に TRUE となるフラグであり、多値予測モデルの場合は mlm という class 属性になる。

これにより、以下のように lm 関数の出力に対して print 関数を作用させると、lm class に対応した print 関数の method である print.lm が呼び出される：

model <- lm(y ~ x, data=df)
print(model)

一方で summary 関数も print と同様に以下のようにジェネリック関数として定義されている：

したがって、summary(model) のように lm 関数の出力を引数に取って実行すると、summary.lm が呼び出されるわけである。

なお、print と summary 以外にも lm class に対応した method を持つジェネリック関数がいくつかあり、以下のように確認することができる：

methods(class="lm")

 [1] add1           alias          anova          case.names     coerce         confint
 [7] cooks.distance deviance       dfbeta         dfbetas        drop1          dummy.coef
[13] effects        extractAIC     family         formula        fortify        hatvalues
[19] influence      initialize     kappa          labels         logLik         model.frame
[25] model.matrix   nobs           plot           predict        print          proj
[31] qqnorm         qr             residuals      rstandard      rstudent       show
[37] simulate       slotsFromS3    summary        variable.names vcov
see '?methods' for accessing help and source code

例えば plot は以下のように予測残差プロットを描画することができる：

plot(model)

こちらについてもコードを読んでみると例によって以下のような plot.lm 関数が定義されている：

まとめにかえて

以上で lm 関数の実装をざっと眺めつつ、R のパッケージがどのように実装されているのかを見てみた。

最小二乗法による係数推定部分については C 言語や FORTRAN で実装したプログラムを呼び出していた。
FORTRAN（≠ Fortran）で書かれたコードは go to 文など昔ながらの記法が含まれていて馴染みが薄いと読みにくいが、以下の記事が参考になった：

推定結果の表示などは、S3 システムを活用して class と method を定義することにより print や summary といったジェネリック関数を利用できるように実装されていた。
S3 をはじめとした R におけるオブジェクト指向プログラミングのより詳細な話題については、以下の書籍の第12~16章が詳しい：
Hadley Wickham, Advanced R, Second Edition
First Edition の日本語訳：石田基広・市川太祐・高柳慎一・福島真太朗翻訳、R言語徹底解説（共立出版、2016）
https://adv-r.hadley.nz/

また、将来自分で R のパッケージを実装する際には、上述の著者らによる下記の書籍が参考になりそうである：
Hadley Wickham and Jennifer Bryan, R Packages, Second Edition
https://r-pkgs.org/

脚注

もう一つ別の R ソースコードのミラーが存在する: https://github.com/SurajGupta/r-source 本文中に挙げた repository は1時間間隔で Subversion の repository と同期しているのに対して、こちらはリリースごとに1コミットされるためバージョンごとの違いを比較する際に読みやすいとのことである（参考）。 ↩︎
説明変数間に完全な多重共線性がある場合はここでは考えない。多重共線性についてはこちらを、完全な多重共線性がある場合の最小二乗推定量についてはこちらを参照。 ↩︎
ここで用いている .Call 関数の使い方と C 言語による外部関数の書き方については、次の記事の解説が分かりやすかった： https://www.r-bloggers.com/2012/11/using-r-callhello/ また、参考文献として挙げた Advanced R の First Edition にも解説が存在していたようである http://adv-r.had.co.nz/C-interface.html ↩︎