Principal axis regression (PAR)は、主成分分析 (PCA)と似たような考え方で、2変量データの主要な変動方向を直線として捉えるために使用されます。 文献では、関連する2変量のline-fitting方法として、major axis (MA) や standardised major axis (SMA) がよく扱われます (Warton et al. 2006)。本記事ではこのうち、共分散行列の第1固有ベクトルを主軸として使うMAに近い考え方を扱います。
本記事では、共分散行列の第1固有ベクトルを主軸として使う、major axis regressionに近い実装を扱います。 以下では、RでPrincipal axis regressionを実装してみます。
Principal axis regressionは、データの共分散行列を計算し、その固有値と固有ベクトルを求めることで、データの主要な変動方向を特定します。 第1固有ベクトルはもっとも分散が大きい方向、第2固有ベクトルはそれに直交する方向を表します。
よく用いられるordinary least squares regression (OLS; 最小二乗回帰)では、目的変数と説明変数を分けて、目的変数方向の誤差が小さくなるように直線を当てはめます。 一方でPARでは、どちらか一方を目的変数として固定するのではなく、2変量の点群全体の変動方向を表す直線を求めます。
RでPrincipal axis regressionを実装するためには、以下の手順を踏む必要があります。
なお、ここでは\(x\)と\(y\)を同じスケールの模擬データとして扱います。 実データで2つの変数の単位やスケールが大きく異なる場合は、log変換や標準化を行うか、共分散行列ではなく相関行列を使うかを検討する必要があります。 共分散行列を使う方法と相関行列を使う方法では、推定される軸の向きが変わることがあります。
Rの実装のセクションでは、詳しい数式の説明は省略し、コードを中心に説明します。 数式を用いたメカニズムの説明は、次のセクションで行います。
まず、適当なデータセットを用意します。 ここでは、rnorm()関数を使用して、ランダムなデータを生成します。 主軸を \(y = \alpha + \beta x\) として、ノイズを加えたデータを生成します。
set.seed(123)
n <- 100
beta <- 1.5
alpha <- 1
# 主軸方向のスコア
t <- rnorm(n, mean = 0, sd = 2)
# 主軸に直交するノイズ
e <- rnorm(n, mean = 0, sd = 0.5)
# 主軸方向ベクトルと直交方向ベクトル
v1 <- c(1, beta) / sqrt(1 + beta^2)
v2 <- c(-beta, 1) / sqrt(1 + beta^2)
x <- t * v1[1] + e * v2[1]
y <- alpha + t * v1[2] + e * v2[2]
plot(x, y, asp = 1, pch = 16)
共分散行列は、以下のように表されます。
\[ \Sigma = \begin{bmatrix}\sigma_{xx} & \sigma_{xy} \\ \sigma_{yx} & \sigma_{yy}\end{bmatrix} \]
ここで、\(\sigma_{xx}\)は\(x\)の分散、\(\sigma_{yy}\)は\(y\)の分散、\(\sigma_{xy}\)は\(x\)と\(y\)の共分散を表します。 なお、共分散とは、2つの変数がどの程度一緒に変動するかを示す指標です。 以下のように計算されます。
\[ \sigma_{xy} = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y}) \]
Rでは、cov()関数を使用して共分散行列を簡単に計算することができます。
共分散行列の固有値と固有ベクトルを計算するためには、eigen()関数を使用します。 ちなみに、固有値はデータの変動の大きさを表し、固有ベクトルはその変動の方向を表します。 Rのeigen()は固有値を大きい順に返すため、vectorsの1列目が第1主軸、2列目が第2主軸に対応します。
eigen_result <- eigen(cov_matrix, symmetric = TRUE)
eigen_resulteigen() decomposition
$values
[1] 3.3335480 0.2331491
$vectors
[,1] [,2]
[1,] 0.5663796 -0.8241445
[2,] 0.8241445 0.5663796eigen()関数は、固有値と固有ベクトルをリスト形式で返します。 固有値はvalues、固有ベクトルはvectorsに格納されます。 \(x\)と\(y\)の2次元のデータの場合、2つの固有値と対応する固有ベクトルが得られます。
[1] "第1主軸の固有値: 3.33354798737934"[1] "第2主軸の固有値: 0.233149107828597"[1] "第1主軸の固有ベクトル: (0.566379566968877, 0.824144517739545)"[1] "第2主軸の固有ベクトル: (-0.824144517739545, 0.566379566968877)"固有ベクトルの符号は任意です。例えば、\((u_x, u_y)\) と \((-u_x, -u_y)\) は同じ軸方向を表します。 そのため、直線の傾き \(u_y/u_x\) は変わりませんが、次の記事で扱う「軸どうしの角度」を計算する場合には、この符号の任意性に注意する必要があります。
固有ベクトル、すなわち主要な変動方向をデータポイントとともにプロットしてみます。 データの中心をベクトルが通るようにします。 つまり、直線が\(y=a + bx\)のとき、\(a\)と\(b\)は以下のように計算されます。
\[ a = \bar{y} - b \bar{x} \]
ここで、\(\bar{x}\)と\(\bar{y}\)はそれぞれ\(x\)と\(y\)の平均値です。 また、\(b\)は固有ベクトルの比率で計算されます。 第1主軸を赤、第2主軸を青で描きます。
# 第1主軸と第2主軸の傾き
major_axis_slope <- eigen_vectors[2, 1] / eigen_vectors[1, 1]
minor_axis_slope <- eigen_vectors[2, 2] / eigen_vectors[1, 2]
# データの中心を通るように切片を計算
major_axis_intercept <- mean(y) - major_axis_slope * mean(x)
minor_axis_intercept <- mean(y) - minor_axis_slope * mean(x)
# データポイントと主要な変動方向をプロット
plot(x, y, asp = 1, pch = 16, main = "Principal Axis Regression")
# 第1主軸の線を追加
abline(
a = major_axis_intercept,
b = major_axis_slope,
col = 2,
lwd = 3
)
# 第2主軸の線を追加
abline(
a = minor_axis_intercept,
b = minor_axis_slope,
col = 4,
lwd = 3
)
ここで行った計算は、2変量データに対するPCAと同じです。 PCAでは、共分散行列の固有値・固有ベクトルを求めることで、データのばらつきが大きい方向から順に新しい軸を定義します。
今回の例では、元の軸は \(x\) と \(y\) ですが、固有値分解によって新しく得られる軸は第1主軸と第2主軸です。
したがって、principal axis regression や major axis regression は、2変量PCAで得られる第1主軸を、2変量データを要約する直線として解釈したものと考えることができます。
今回は、共分散行列の固有値と固有ベクトルを用いて、データの主要な変動方向を特定しました。 このメカニズムについて、本文では詳しく説明していませんでした。 どうして共分散行列や固有ベクトルを使うと、主軸が求められるのでしょうか。
直感的には、主軸とは散布図の点群の真ん中を通り、データのばらつきを最もよく表す直線です。
たとえば、散布図に適当な直線を1本引いたとします。 このとき、それぞれの点から直線に向かって垂線を下ろすと、点がその直線からどれだけ離れているかを測ることができます。
主軸としてよい直線とは、各点から直線までの距離が全体として小さくなるような直線です。 つまり、点群の流れに沿っている直線ほど、主軸としてよい直線だと考えられます。
一方で、これは別の見方もできます。 点群にうまく沿った直線では、各点はその直線方向に大きくばらついています。 逆に、点群に合っていない直線では、直線方向のばらつきは小さく、直線から外れた方向のずれが大きくなります。
したがって、主軸はデータをその方向に見たとき、ばらつきが最も大きくなる方向と考えることができます。
この「ある方向に見たときのばらつき」を計算するために使うのが、共分散行列です。
2変量データの各点を、平均を中心にした形で
\[ x_i - \bar{x}, \quad y_i - \bar{y} \]
とします。 ここで、ある方向を表すベクトルを\(\mathbf{u}\)を以下のように定義します。
\[ \mathbf{u} = \begin{bmatrix} u_x \\ u_y \end{bmatrix} \]
このベクトルは、直線の向きを表します。
各点がこの方向にどれくらい進んでいるかは、次のように計算できます。
\[ s_i = u_x (x_i - \bar{x}) + u_y (y_i - \bar{y}) \]
ここで、\(s_i\)は点\(i\)が方向\(\mathbf{u}\)に沿ってどれだけ進んでいるかを表すスカラー値です。 つまり、それぞれの点の「直線方向の位置」のようなものです。 この \(s_i\) のばらつきが大きいほど、その直線方向にデータが大きく広がっていることになります。
このばらつきを計算するために、共分散行列を使います。 共分散行列\(\Sigma\)は、データのばらつきと変動の方向を表す行列です。 ある方向\(\mathbf{u}\)に沿ったばらつきは、次のように計算されます。
\[ \text{Var}(s) = \mathbf{u}^T \Sigma \mathbf{u} \]
また、共分散行列\(\Sigma\)は次のように表されます。
\[ \Sigma = \begin{bmatrix}\sigma_{xx} & \sigma_{xy} \\ \sigma_{yx} & \sigma_{yy}\end{bmatrix} \]
ここで、\(\sigma_{xx}\)は\(x\)の分散、\(\sigma_{yy}\)は\(y\)の分散、\(\sigma_{xy}\)は\(x\)と\(y\)の共分散を表します。 分散は、データが平均からどれだけ広がっているかを表す指標であり、共分散は2つの変数がどの程度一緒に変動するかを示す指標です。
分散は以下のように計算されます。
\[ \sigma_s^2 = \frac{1}{n-1} \sum_{i=1}^{n} (s_i - \bar{s})^2 \]
また、共分散は以下のように計算されます。
\[ \sigma_{xy} = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y}) \]
つまり、共分散行列は、任意の方向 \(\mathbf{u}\) に対して、その方向のばらつきがどれくらいかを計算するための情報を持っていると考えることができます。
では、どの方向 \(\mathbf{u}\) を選べば、ばらつきが最大になるのでしょうか。
ある方向を表す単位ベクトルを
\[ \mathbf{u} = \begin{bmatrix} u_x \\ u_y \end{bmatrix} \]
とします。
この方向に沿ったデータのばらつきは、共分散行列 \(\Sigma\) を使って
\[ \mathbf{u}^{\mathrm{T}} \Sigma \mathbf{u} \]
と表すことができます。
ここで、\(\mathbf{u}^{\mathrm{T}} \Sigma \mathbf{u}\) は1つの値になります。 この値が大きいほど、データは方向 \(\mathbf{u}\) に沿って大きくばらついていることを意味します。
したがって、主軸を求めることは、
\[ \mathbf{u}^{\mathrm{T}} \Sigma \mathbf{u} \]
が最も大きくなる方向 \(\mathbf{u}\) を探すことに対応します。 ただし、方向だけを考えるため、\(\mathbf{u}\) の長さは1に固定します。
\[ \mathbf{u}^{\mathrm{T}}\mathbf{u} = 1 \]
この条件のもとで \(\mathbf{u}^{\mathrm{T}} \Sigma \mathbf{u}\) を最大にする問題を解くと、
\[ \Sigma \mathbf{u} = \lambda \mathbf{u} \]
という形の式が得られます。
これは、共分散行列 \(\Sigma\) の固有値・固有ベクトルの式です。 ここで、\(\mathbf{u}\) が固有ベクトル、\(\lambda\) が固有値です。
この式は、線形代数の分野で「固有値問題」と呼ばれるものです。 この関係を導くには、制約付き最大化問題としてラグランジュの未定乗数法を使うのが一般的です。 ただし、実際の計算はRの eigen() が行ってくれるため、ここでは詳しい導出は省略します。
つまり、共分散行列の固有ベクトルは「ばらつきが大きくなる方向」を表し、固有値は「その方向にどれだけばらついているか」を表します。
特に、最も大きい固有値に対応する固有ベクトルが、第1主軸の方向になります。
簡単にまとめると、以下のようになります。
2変量データの場合、共分散行列の固有ベクトルは2つ得られます。
このうち、第1固有ベクトルの方向を直線として描いたものが、今回のprincipal axis regressionにおける主軸です。
今回の仮想データを利用して、視覚的に理解を深めてみます。 まずは、冒頭の散布図を再度表示します。
次に、共分散行列を計算し、固有値と固有ベクトルを求めます。 固有ベクトルの符号は任意であるため、固有ベクトルの符号を統一するために、例えば第1主軸の\(x\)成分が正になるように調整します。
固有値や固有ベクトルは英語ではeigenvalueやeigenvectorと呼ばれます。 このため、Rの関数もeigen()という名前になっています。
次に、比較用の適当な方向を用意し、それぞれの方向に対して、各点から直線への足を計算する関数を定義します。 ベクトルは、方向を表すために正規化して単位ベクトルにしておきます。
# データの中心
center <- c(mean(x), mean(y))
u0 <- c(1, 0.4) # 適当な方向のベクトル
u0 <- u0 / sqrt(sum(u0^2)) # 正規化して単位ベクトルにする
# 直線からの距離を計算する関数
get_foot_points <- function(x, y, u, center) {
X <- cbind(x - center[1], y - center[2])
score <- as.vector(X %*% u) # 各点のスコアを計算
foot_centered <- outer(score, u) # スコアに基づいて足の位置を計算
foot <- sweep(foot_centered, 2, center, "+") # 中心を足す
foot
}
# 各点から直線への距離を計算
foot0 <- get_foot_points(x, y, u0, center)
foot1 <- get_foot_points(x, y, u1, center)データの上に、適当に引いた直線と各点からその直線への足をプロットしてみます。
次に、ここに第1主軸の線と、各点から第1主軸への足を追加してみます。
plot(x, y, asp = 1, pch = 16, xlab = "x", ylab = "y")
# 適当な方向の線と距離
abline(
a = center[2] - (u0[2] / u0[1]) * center[1],
b = u0[2] / u0[1],
col = "gray30",
lwd = 2
)
segments(x, y, foot0[, 1], foot0[, 2], col = "gray70")
# 第1主軸の線と距離
abline(
a = center[2] - (u1[2] / u1[1]) * center[1],
b = u1[2] / u1[1],
col = 2,
lwd = 3
)
segments(x, y, foot1[, 1], foot1[, 2], col = 3)
この時の各点から直線までの距離の二乗和を計算して比較します。
[1] 85.83191[1] 37.87125第1主軸の方が、適当な方向よりも点から直線への距離の二乗和が小さくなっていることがわかります。 つまり、データのばらつきにうまく沿った直線を選ぶと、点から直線への距離が全体として小さくなることがわかります。 第1主軸は、この距離の合計が最小になるような方向を表していることがわかります。
第2主軸も同様にプロットしてみます。
u2 <- eigen_vectors[, 2]
if (u2[1] < 0) {
u2 <- -u2
}
plot(x, y, asp = 1, pch = 16, xlab = "x", ylab = "y")
# 第1主軸の線
abline(
a = center[2] - (u1[2] / u1[1]) * center[1],
b = u1[2] / u1[1],
col = 2,
lwd = 3
)
# 第2主軸の線
abline(
a = center[2] - (u2[2] / u2[1]) * center[1],
b = u2[2] / u2[1],
col = 4,
lwd = 3
)
第1主軸はデータのばらつきに沿った方向を表し、第2主軸は第1主軸に直交する方向を表していることがわかります。
cov()関数で共分散行列を計算し、eigen()関数で固有値と固有ベクトルを求めることができます。この記事では、2変量データから主軸を1本求める基本的な考え方を確認しました。 次の記事では、この第1固有ベクトルを使って、種ごとの代表値から求めた主軸と、種内変異を含めたときの主軸のずれを調べます。