Rの関数:poly {stats}

Rの関数から poly {stats} を確認します。

関数 poly とは

poly は、指定された数値ベクトルに基づいて「Orthogonal Polynomials(以降、直交多項式)」または「Raw Polynomials(以降、本投稿では通常多項式)」の行列を生成するための関数です。

非線形な関係性をモデル化する際、線形回帰モデルの中に説明変数の2乗や3乗の項を追加する手法（多項式回帰）が用いられます。

しかし、\(X\) と \(X^2\)、\(X^3\) といった「生の累乗値」をそのまま回帰モデルに投入すると、変数間に強い相関が生じます。

当該現象は多重共線性（Multicollinearity）と呼ばれ、推計された係数の分散を著しく増大させ、各項の統計的有意性を誤って評価する原因となります。

poly 関数（デフォルト設定）は、QR分解を用いて、互いに完全に無相関（相関係数がゼロ）となるように変換された新しい多項式変数のセットを生成します。

当該機能を利用することで、モデルの全体的な予測精度を損なうことなく、1次、2次、3次といった各次数が持つ純粋で独立した説明力を評価することが可能となります。

関数 poly の引数

args(poly)
function (x, ..., degree = 1, coefs = NULL, raw = FALSE, simple = FALSE) 
NULL
  1. x

    • 基礎となる説明変数を表す数値ベクトル、または行列です。欠損値（NA）を含めることはできず、事前に除外する必要があります。
  2. ...

    • 追加の数値ベクトルを指定します。複数のベクトルを与えた場合、内部で polym 関数が呼び出され、変数の交互作用を含む多変量の多項式が生成されます。
  3. degree

    • 生成する多項式の最高次数を指定する整数です。
    • デフォルト: 1。
    • 1以上であり、かつデータに含まれる「一意な値の数（ユニークな観測値の種類数）」よりも小さい値を指定する必要があります。
  4. coefs

    • 直交化の基準として用いる係数（alpha および norm2）のリストです。
    • デフォルト: NULL（入力データ x から新しく係数を計算します）。
    • モデルの学習時に生成した直交多項式の基準を、新しい未知のデータ（テストデータ）に対して適用し、同一のスケールで直交変換を行うために使用されます。
  5. raw

    • どのような多項式を生成するかを決定する論理値です。
    • デフォルト: FALSE。
    • FALSE を指定すると、無相関化およびスケーリングが行われた「直交多項式」を出力します。TRUE を指定すると、単純なべき乗計算（\(x, x^2, x^3 \dots\)）による「通常多項式」を出力します。
  6. simple

    • 戻り値のオブジェクト形式を簡略化するかどうかの論理値です。
    • デフォルト: FALSE。
    • TRUE を指定すると、係数リストなどの付加的な属性（Attributes）を一切持たない、純粋な数値行列のみを返します。

シミュレーションコード

以下に、poly 関数の機能を確認するためのサンプルデータを用いたシミュレーションコードを示します。

本シミュレーションでは、3次関数に従う非線形なデータを生成し、raw = TRUE（通常多項式）と raw = FALSE（直交多項式）の双方が生成する変数の性質（相関構造）の違いを比較します。

さらに、双方を回帰モデルに組み込んだ際、予測結果は完全に一致する一方で、各変数の有意性の評価がどのように変化するかを検証します。

なお、有意水準は5%とします。

サンプルデータの生成

# パッケージの読み込み
library(ggplot2)
library(gridExtra)

# 乱数シードの固定
seed <- 20260506
set.seed(seed)

# 1. サンプルデータの生成
# 説明変数 X (等間隔なデータ)
n_obs <- 100
var_X <- seq(2, 12, length.out = n_obs)

# 応答変数 Y (3次関数に従う非線形データ + 正規ノイズ)
# 真のモデル: Y = 10 + 2*X - 1.5*X^2 + 0.5*X^3
error_term <- rnorm(n_obs, mean = 0, sd = 10)
var_Y <- 10 + 2 * var_X - 1.5 * var_X^2 + 0.5 * var_X^3 + error_term

df_sim <- data.frame(X = var_X, Y = var_Y)

cat('=== 生成したサンプルデータの一部を確認 ===\n')
str(df_sim)
=== 生成したサンプルデータの一部を確認 ===
'data.frame':   100 obs. of  2 variables:
 $ X: num  2 2.1 2.2 2.3 2.4 ...
 $ Y: num  12.5 18.3 14.4 26.1 13.8 ...

poly関数による多項式行列の生成

# 2. poly関数による多項式行列の生成
degree_val <- 3

# 通常多項式 (raw = TRUE)
poly_raw <- poly(df_sim$X, degree = degree_val, raw = TRUE)

# 直交多項式 (raw = FALSE, デフォルト)
poly_ortho <- poly(df_sim$X, degree = degree_val, raw = FALSE)

cat("=== 生成された行列の相関構造の比較 ===\n")
cat("【通常多項式 (raw = TRUE) の列間相関係数】\n")
print(round(cor(poly_raw), 3))

cat("\n【直交多項式 (raw = FALSE) の列間相関係数】\n")
print(round(cor(poly_ortho), 3))
=== 生成された行列の相関構造の比較 ===
【通常多項式 (raw = TRUE) の列間相関係数】
      1     2     3
1 1.000 0.983 0.947
2 0.983 1.000 0.989
3 0.947 0.989 1.000

【直交多項式 (raw = FALSE) の列間相関係数】
  1 2 3
1 1 0 0
2 0 1 0
3 0 0 1

通常多項式では、各項の間に強い相関が存在します。

一方、直交多項式では、すべての列間の相関係数が完全にゼロ（無相関）となっていることが確認できます。

回帰モデルへの適用と係数の比較

# 3. 回帰モデルへの適用と係数の比較
# それぞれの多項式行列を説明変数として線形回帰を実行
model_raw <- lm(Y ~ poly_raw, data = df_sim)
model_ortho <- lm(Y ~ poly_ortho, data = df_sim)

cat("=== 回帰モデルのサマリー比較 ===\n")
cat("【通常多項式モデルの係数推定結果】\n")
print(round(coef(summary(model_raw)), 4))

cat("\n【直交多項式モデルの係数推定結果】\n")
print(round(coef(summary(model_ortho)), 4))
=== 回帰モデルのサマリー比較 ===
【通常多項式モデルの係数推定結果】
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   5.5837    13.0216  0.4288   0.6690
poly_raw1     5.9872     6.7763  0.8835   0.3791
poly_raw2    -2.1616     1.0491 -2.0605   0.0421
poly_raw3     0.5293     0.0496 10.6733   0.0000

【直交多項式モデルの係数推定結果】
             Estimate Std. Error  t value Pr(>|t|)
(Intercept)  199.2313     0.9651 206.4272        0
poly_ortho1 1796.8664     9.6514 186.1767        0
poly_ortho2  680.7039     9.6514  70.5290        0
poly_ortho3  103.0127     9.6514  10.6733        0
最高次項（3次項）のt値の一致

両モデルの3次項の t value（t値）が一致していることを確認できます。

  • poly_raw3 の t値: 10.6733
  • poly_ortho3 の t値: 10.6733

最高次項がもたらす「それ以下の次数（1次および2次）では説明できない純粋な追加情報」の価値は、変数の表現方法をどのように変更しても不変であり、上記の一致は、両モデルが根本的に同一の現象を捉えている証左と言えます。

1次項・2次項における多重共線性の影響

下位の次数（1次および2次）の推定結果には、多重共線性の有無による差異が表れています。

  • 通常多項式 (model_raw)

    • 1次項 (poly_raw1) のP値が 0.3791 となり、設定した有意水準を超えています。データは \(X\) が 2 から 12 へ増加する過程で明確な上昇トレンド（1次の効果）を含んでいますが、通常多項式では \(X\)、\(X^2\)、\(X^3\) が同時に単調増加して強い相関を持つため、モデルでは「どの項のおかげで \(Y\) が増えているのか」を判別できなくなっていることを示唆しています。
  • 直交多項式 (model_ortho)

    • すべての項のP値が 0と設定した有意水準を下回っています。各次数が完全に無相関化されているため、「直線的な上昇（1次）」「U字の曲がり（2次）」のそれぞれの効果が、互いに干渉することなく純粋に評価されています。
切片と標準誤差の性質
  • 標準誤差の均一化:

    • 直交多項式では、1次から3次までのすべての標準誤差（Std. Error）が 9.6514 に統一されています。ベクトルの長さが揃えられているため、係数（Estimate）の大きさをそのまま「変数ごとの影響力の強さ」として比較することが可能です。
  • 切片の意味の違い:

    • 通常多項式の切片（5.5837）は「\(X=0\) のときの \(Y\) の値」を意味します。今回のデータは \(X\) が 2 から 12 であるため、\(X=0\) は遠くの外挿となり、推定精度が低く、P値が設定した有意水準を上回っています（P=0.6690）。
    • 対照的に、直交多項式の切片（199.2313）は、中心化処理の恩恵により「\(Y\) の平均値」を表します。それゆえ、高いt値と設定した有意水準を下回る有意性が付与されています。

モデルの予測値の比較

# 4. モデルの予測値の比較
# 二つのモデルの最終的な予測性能が等価であることを確認
df_sim$Pred_Raw <- predict(model_raw)
df_sim$Pred_Ortho <- predict(model_ortho)

cat(paste0("両モデルの予測値は一致しているか: ", all.equal(df_sim$Pred_Raw,df_sim$Pred_Ortho)))
両モデルの予測値は一致しているか: TRUE

変数の内部的な表現方法（スケールや相関）が異なるだけで、モデル全体としての予測精度は完全に一致しています。

それゆえ、両モデルの予測値は全て一致しています。

生成された変数の視覚化

# 5. 生成された変数の視覚化
# 各変数の形状をプロット用のデータフレームに格納
df_poly <- data.frame(
  X = var_X,
  Raw_1 = poly_raw[, 1],
  Raw_2 = poly_raw[, 2],
  Raw_3 = poly_raw[, 3],
  Ortho_1 = poly_ortho[, 1],
  Ortho_2 = poly_ortho[, 2],
  Ortho_3 = poly_ortho[, 3]
)

# 通常多項式の形状プロット
p1 <- ggplot(df_poly, aes(x = X)) +
  geom_line(aes(y = Raw_1, color = "1次"), linewidth = 1) +
  geom_line(aes(y = Raw_2, color = "2次"), linewidth = 1) +
  geom_line(aes(y = Raw_3, color = "3次"), linewidth = 1) +
  labs(
    title = "図1: 通常多項式の形状 (raw = TRUE)",
    y = "変数の値",
    color = "次数"
  ) +
  scale_color_manual(values = c("1次" = "blue", "2次" = "green", "3次" = "red")) +
  theme_minimal()

# 直交多項式の形状プロット
p2 <- ggplot(df_poly, aes(x = X)) +
  geom_line(aes(y = Ortho_1, color = "1次"), linewidth = 1) +
  geom_line(aes(y = Ortho_2, color = "2次"), linewidth = 1) +
  geom_line(aes(y = Ortho_3, color = "3次"), linewidth = 1) +
  labs(
    title = "図2: 直交多項式の形状 (raw = FALSE)",
    y = "変数の値",
    color = "次数"
  ) +
  scale_color_manual(values = c("1次" = "blue", "2次" = "green", "3次" = "red")) +
  theme_minimal()

grid.arrange(p1, p2, ncol = 1)
Figure 1
図1（上段）：通常多項式の形状

Figure 1 の図1は、単なる \(X, X^2, X^3\) の計算結果を描画したものです。

1次項（青線）は値が小さいためほぼ平坦に見えますが、3次項（赤線）は \(X=12\) において \(12^3 = 1728\) に達し、縦軸のスケールを押し上げています。このように変数の桁数が大きく異なると、回帰分析の内部計算において精度の低下を招く原因となります。

図2（下段）：直交多項式の形状

Figure 1 の図2は、poly 関数がQR分解を用いて無相関化とスケーリングを行った後の変数の形状を描画したものです。

縦軸のスケールを確認しますと、すべての曲線が概ね \(\pm 0.2\) の範囲に収められています。全項の標準誤差が 9.6514 に統一されていたのは、このように各ベクトルの長さ（分散）が等しくなるよう内部で調整されているためです。

1次（青線） は、データの中心付近（\(X=7\)付近）でゼロと交差する純粋な直線トレンド。2次（緑線） は、1次の直線的要素を排除した、純粋なU字型（または逆U字型）の曲率。3次（赤線）は、 1次と2次の要素を排除した、純粋なS字型のうねり、と中心が原点からズレた非対称なデータであっても、poly 関数は各曲線を「互いに全く異なる動きをする（相関ゼロの）波形」へと変換しています。

以上です。