Rで多項式回帰するときに使うpoly関数のオプションrawの働きと多項式の直交性

コメントはありません

Rで多項式回帰をするとき、モデル式を作るためにpolyという関数を使いました。
poly関数は変数と次数を与えると、その次数の多項式モデルを返してくれる関数なのですが、我々が普通にイメージする多項式回帰を行うためにはオプションraw=TRUEを指定しなければなりません。
ではrawは何を指定するオプションなのか、とヘルプを引くと

raw: if true, use raw and not orthogonal polynomials.

とorthogonal polynomials、すなわち直交多項式などというパワーワードが出てきます。
ヘルプの中には直交多項式が何なのかに関するヒントは全く記されていないのでググって調べるとチェビシェフ多項式、ルジャンドル多項式とかいう新しい用語がヒットします。
それらは微積分が混じった恐ろしげな数式で定義されており、その定義式を解読しするだけで中々しんどい。
そしてようやくそれらを解読し計算できるようになったところでふと気づくのは、結局それが多項式回帰とどんな関係があるのかは分からない。
ここらで心が折れてきます。

そこをなんとか私なりに理解をしたつもりになったので、その内容をここにまとめておきます。
一言で結論を述べますと、ここでいう直交多項式とはn次多項式の各項の観測ベクトルをGram-Schmidt直交化によって直交化したものになっています。

【統計モデル】

n次元ベクトルxをpoly関数に与え、オプションをdegree=3, raw=FALSEとして目的変数yの線形回帰をするとします。
すなわち

#R
lm(y ~ poly(x, degree = 3, raw = FALSE))

このとき、実際に行われる回帰のモデル式は

\beta_0 + \beta_1 x' + \beta_2 x'' + \beta_3 x'''

になります。
ただしx'は正規化されたxであり
\displaystyle x' = \frac{1}{\sigma'}(x - \bar{x'})

\displaystyle \bar{x'} = \frac{1}{n}\sum_{i=1}^nx_i,\ \sigma' = \sqrt{\frac{1}{n}\sum_{i=1}^n (x_i-\bar{x'})^2}

です。
またx''はx^2をGram-Schmidt直交化によってx'と直交にしたベクトル
\displaystyle b_2 = x^2 - x'(x^2\cdot x')

を正規化したベクトル
\displaystyle x'' = \frac{1}{\sigma''}(b_2-\bar{b_2})

です。
同様にしてx'''もx^3をGram-Schmidt直交化によってx',x''と直交にしたベクトルを正規化したベクトルとなります。

degreeが3以上のときにもGram-Schmidt直交可は一次独立なベクトルがあるかぎり手続きを継続できるので、観測ベクトルのベキ乗の系列が一次独立であるかぎり目的の次数まで観測ベクトルの変形ができます。
逆に観測ベクトルのベキ乗系列が一次従属になれば関数polyはエラーを返します。

【計算】

具体的にpoly関数の働き、すなわちどういう計算をしているかを示します。

ベクトルxをpoly関数に与え、degree=3, raw=FALSEとしたとき

#R
poly(x, degree = 2, raw = FALSE)

の戻り値は行列になり、その行列は一列目が

#R
b1 <- (x - mean(x)) / sqrt(sum((x - mean(x))^2))

で計算されます。なおこれは平均と分散とで正規化されたxです。
また二列目は

#R
b2 <- x^2 - b1 * sum(x * b1)/sum(b1^2)
b2 <- (b2 - mean(b2)) / sqrt(sum((b2 - mean(b2))^2))

と上節で述べたようにx^2をGram-Schmidt直交化によってxと直交なベクトルを正規化したベクトルになります。

degreeが2よりも大きくなったとしても同じ要領で計算されます。

【余談】

上記の結論はpoly関数の内部定義を読み解くことで見出しました。
内部定義を見ると分かるのですが、実際にはqrという関数を使って多項式を直交化しています。
関数qrは行列をQR分解するR関数で、QR分解とは行列を高速に上三角化する処理です。
ベクトルの直交化がqrでできるというテクニックは使い道があるのかもしれません。

関連記事

コメントを残す