日々のつれづれ

不惑をむかえ戸惑いを隠せない男性の独り言

分散分析のメモ・その3

分散分析と F 検定の関係が分からなくなることが多いのでメモ

二元配置分散分析(Two-Way ANOVA)

分散分析は、3 つ以上の群の平均値に差があるかどうかを評価する統計手法で、異なる条件や処理が平均値に与える影響を評価できる。な お、二元配置分散分析は、要因数が 2 種類の分散分析。

各々の群について、処理の影響を示す群間分散と処理と関係のない群内分散を求める。次に、群内分散に対する群間分散の比を求める。そして、この比を F 値として、F 分布から p 値を求める。このとき、二元配置分散分析は、異なる条件や処理が同時に作用したときの効果(交互作用)を評価することが特徴になる。なお、このときの帰無仮説は「すべての群の平均値が等しい」になる。

対応のない二元配置分散分析

一元配置分散分析は 1 要因の各水準からなる群に対して要因分散と誤差分散を求めて、誤差分散に対する要因分散の比が F 値になりました。

二元配置分散分析は、要因が 2 に増えますが、計算手順は一元配置分散分析と同じく誤差分散に対する要因分散の比を求めます。ただし、2 つの要因に加えて、その交互作用を計算する点が異なります。 そして、$F_A$、$F_B$、$F_{AB}$のそれぞれについて F 検定を行い、有意差を確認します。

F 値の式

対応のない二元配置分散分析の場合、個体間の変動を考慮しない。

交互作用と誤差の自由度の計算に注意して、要因 A、要因 B、交互作用、誤差について、それぞれの変動、自由度から分散を計算できる。

項目 変動 ($SS$) 自由度 ($\nu$, $df$) 分散 F 値
全体 $$SS_T = \sum{i=1}^a \sum{j=1}^b \sum{k=1}^{n{ij}} (X{ijk} - \overline{X}{..})2$$ $$\nu_T = N - 1$$ $$-$$ $$-$$
要因 A $$SS_A = \sum{i=1}^a n{i} (\overline{X}{i.} - \overline{X}{..})2$$ $$\nu_A = a - 1$$ $$MS_A = \frac{SS_A}{\nu_A}$$ $$F_A = \frac{MS_A}{MS_E}$$
要因 B $$SS_B = \sum{j=1}^b n{j} (\overline{X}{.j} - \overline{X}{..})2$$ $$\nu_B = b - 1$$ $$MS_B = \frac{SS_B}{\nu_B}$$ $$F_B = \frac{MS_B}{MS_E}$$
交互作用 $$SS{AB} = \sum{i=1}^a \sum{j=1}^b n{ij} (\overline{X}{ij} - \overline{X}{i.} - \overline{X}{.j} + \overline{X}{..})2$$ $$\nu_{AB} = ab - 1 - \nu_A - \nu_B = \nu_A \times \nu_B$$ $$MS{AB} = \frac{SS{AB}}{\nu_{AB}}$$ $$F{AB} = \frac{MS{AB}}{MS_E}$$
誤差 $$SS_E = SS_T - SS_A -SS_B -SS{AB} \= \sum{i=1}^a \sum{j=1}^b \sum{k=1}^{n{ij}} (X{ijk} - \overline{X}_{ij})2$$ $$\nu_E = N - ab = \nu_T - \nu_A - \nu_B - \nu_{AB}$$ $$MS_E = \frac{SS_E}{\nu_E}$$ $$-$$
  • サンプル数
    • $N$ $\cdots$ 全観測値の数
    • $n_i$ $\cdots$ 要因 A の第 $i$ 水準に属する観測値の数
    • $n_j$ $\cdots$ 要因 B の第 $j$ 水準に属する観測値の数
    • $n_{ij}$ $\cdots$ 要因 A の第 $i$ 水準と要因 B の第 $j$ 水準に属する観測値の数 (繰り返し測定数)
  • 水準数
    • $a$ $\cdots$ 要因 A の水準数
    • $b$ $\cdots$ 要因 B の水準数
  • 観測値
    • $X_{ijk}$ $\cdots$ 要因 A の第 $i$ 水準と要因 B の第 $j$ 水準に属する $k$ 番目の観測値
    • $\overline{X}_{i.}$ $\cdots$ 要因 A の第 $i$ 水準に属する観測値の平均
    • $\overline{X}_{.j}$ $\cdots$ 要因 B の第 $j$ 水準に属する観測値の平均
    • $\overline{X}_{ij}$ $\cdots$ 要因 A の第 $i$ 水準と要因 B の第 $j$ 水準に属する観測値の平均
    • $\overline{X}_{..}$ $\cdots$ 全観測値の平均

R で実装する

F 値と p 値は、対応無しの場合と同様に計算します。

デモデータの作成

デモ用のデータを作成します。このデータは、3 つの水準(A, B, C)を持つ要因 A(gA)と、2 つの水準(a, b)を持つ要因 B (gB)によって分類された測定値(value)からなります。

rm(list=ls(all=TRUE))
set.seed(123) # 乱数のシードを固定
data <- data.frame(
  gA = factor(rep(LETTERS[1:3], each = 4)), # 要因A、3水準
  gB = factor(rep(letters[1:2], each = 2)), # 要因B、2水準
  value = unlist(lapply(c(1, 4) * 40, sample, size = 6, replace = TRUE)) # 測定値
)

knitr::kable(sapply(levels(data$gA), function(i) {
  sapply(levels(data$gB), function(j) {
    paste(subset(data, subset = gA == i & gB == j)$value, collapse = ",")
  })
}))
A B C
a 31,15 37,14 90,91
b 14,3 118,153 91,92

表の形は対応のある一元配置分散分析と同じ。

でも、対応のある一元配置分散分析が a と b は互いに影響しない(独立になる)ことに対して、対応のない二元配置分散分析は a と b は互いに影響しあう可能性、交互作用を考えることになる。

ここが大きく違う点で、要因A、要因Bに加えて、交互作用に対してもF検定をすることになる。

分散分析の関数を使わない

# 各水準の観測個数
nij <- table(data[,-3])

# 総変動
# 平均
mean_T <- mean(data$value)
# 変動、偏差平方和
ss_T <- sum((data$value - mean_T)^2)
# 自由度、総データ数 - 1
df_T <- nrow(data) - 1

# 要因A
# 平均
mean_A <- as.vector(tapply(data$value, data$gA, mean))
# 変動、偏差平方和
ss_A <- sum((mean_A[data$gA] - mean_T)^2)
# 自由度、水準数 - 1
df_A <- length(levels(data$gA)) - 1
# 分散
ms_A <- ss_A / df_A

# 要因B
# 平均
mean_B <- as.vector(tapply(data$value, data$gB, mean))
# 変動、偏差平方和
ss_B <- sum((mean_B[data$gB] - mean_T)^2)
# 自由度、水準数 - 1
df_B <- length(levels(data$gB)) - 1
# 分散
ms_B <- ss_B / df_B

# 交互作用
# 平均
mean_AB <- tapply(data$value, list(data$gA, data$gB), mean)
# 変動、偏差平方和
ss_AB <- sum(nij*t(t(mean_AB - mean_A) - mean_B + mean_T)^2)
# 自由度、(要因A水準数 - 1)(要因B水準数 - 1)
df_AB <- df_A*df_B
# 分散
ms_AB <- ss_AB / df_AB

# 誤差
# 変動、偏差平方和
ss_E <- ss_T - ss_A - ss_B - ss_AB
# ss_E <- sum((data$value - mapply(rep, t(mean_AB), t(nij)))^2)
# 自由度、(要因A水準数 - 1)(要因B水準数 - 1)
df_E <- df_T - df_A - df_B - df_AB
# 分散
ms_E <- ss_E / df_E

# F値、要因分散/誤差分散
f_A <- ms_A / ms_E
f_B <- ms_B / ms_E
f_AB <- ms_AB / ms_E

# F値からp値を求める
p_A <- pf(f_A, df_A, df_E, lower.tail = FALSE)
p_B <- pf(f_B, df_B, df_E, lower.tail = FALSE)
p_AB <- pf(f_AB, df_AB, df_E, lower.tail = FALSE)

# 結果の表示
cat("factor A - F-value:", round(f_A, 2), "/ p-value:", round(p_A, 4), "\n")
cat("factor B - F-value:", round(f_B, 2), "/ p-value:", round(p_B, 4), "\n")
cat("factor AB - F-value:", round(f_AB, 2), "/ p-value:", round(p_AB, 4), "\n")
factor A - F-value: 37.38 / p-value: 4e-04
factor B - F-value: 17.46 / p-value: 0.0058
factor AB - F-value: 25.9 / p-value: 0.0011

分散分析の関数を使う

二元配置分散分析は 2 要因を$\ast$で繋ぐ。

# 対応のある一元配置分散分析
res <- aov(value ~ gA * gB, data = data)
summary(res)
            Df Sum Sq Mean Sq F value  Pr(>F)
gA           2  13287    6644   37.38 0.00041 ***
gB           1   3104    3104   17.46 0.00582 **
gA:gB        2   9207    4604   25.90 0.00112 **
Residuals    6   1067     178
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
要因 F 値 p 値 $\alpha$ = 0.05 $\alpha$ = 0.01
要因 A 37.38 0.00041 有意差あり 有意差あり
要因 B 17.46 0.00582 有意差あり 有意差あり
交互作用 25.90 0.00112 有意差あり 有意差あり

つまり、いずれの要因も有意差があって、交互作用もある、という結果になる。

ただし、交互作用に有意差があっても、要因 A と要因 B のどの水準の組み合わせに有意差があるのかは分からない。

この先の解析は多重比較になる。

分散分析のメモ・その2

分散分析と F 検定の関係が分からなくなることが多いのでメモ

一元配置分散分析

分散分析は、3 つ以上の群の平均値に差があるかどうかを評価する統計手法で、異なる条件や処理が平均値に与える影響を評価できる。な お、一元配置分散分析は、要因数が 1 種類の分散分析。

各々の群について、処理の影響を示す群間分散と処理と関係のない群内分散を求める。次に、群内分散に対する群間分散の比を求める。そして、この比を F 値として、F 分布から p 値を求める。なお、このときの帰無仮説は「すべての群の平均値が等しい」になる。

対応がある場合の一元配置分散分析

1 つの要因につき同じ個体や対象を複数回測定したとき「対応がある」と言う。

各群のデータが同じ個体や対象ごとに対応関係があるので、同じ人に異なる処置を施してその効果を測定した場合や、同じ物体に異なる条件で実験した場合などが該当する。

F 値の計算では、この個体や対象の違いによる変動(個体間変動)がばらつきに加わるので、群内変動からこれを除いた調整済み誤差分散を使う。

F 値の式

\begin{eqnarray}

F_{(\nu_1,\nu_2)} &=& \frac{要因分散}{調整済み誤差分散}
&=& \frac{\frac{群間変動}{群間変動の自由度}}{\frac{調整済み群内変動}{調整済み群内変動の自由度}}\\
&=& \frac{\frac{群間変動}{群間変動の自由度}}{\frac{郡内変動-個体間変動}{群内変動の自由度-個体間変動の自由度}}
&=& \frac
{\frac{反復数 \times \sum(群平均 - 総平均)^2}{処理数 - 1}}
{\frac
{\sum\sum(測定値 - 群平均)^2 - \sum\sum(個体別平均 - 総平均)^2}
{(データ総数 - 処理数)-(個体数 - 1)}
}\\
&=& \frac
{\frac{n\sum(\bar{x_i}-\bar{\bar{x}})^2}{\nu_1}}
{\frac{\sum\sum(x-\bar{x_i})^2 - \sum\sum(\bar{x_j}-\bar{\bar{x}})^2}{\nu_2}}

\end{eqnarray}
  • 調整済み群内変動は、群内変動から個体間変動を引いた値。要因とは別の影響によるデータのばらつきを示す。
  • $\nu_1$ は、群間変動の自由度で、処理数(群の数)から 1 を引いた値。
  • $\nu_2$ は、調整済み群内変動の自由度で、群内変動の自由度から個体間変動の自由度を引いた値。
  • $n$ は、各群のデータ数(反復数)。
  • $\bar{x_i}$ は、i 番目の群の平均値。
  • $\bar{\bar{x}}$ は、すべてのデータの平均値(総平均)。
  • $x$ は、各データの測定値。
  • $\bar{x_j}$ は、j 番目の個体や対象の平均値。

計算の流れ

  1. 群間変動を求める。各グループの平均と全体の平均の差の二乗和を計算する。
  2. 群間変動の自由度を求める。処理数(グループ数)から 1 を引く。
  3. 群内変動を求める。各データと、そのグループの平均の差の二乗和を計算する。
  4. 群内変動の自由度を求める。総データ数から処理数を引く。
  5. 個体間変動を求める。各個体の平均と全体の平均の差の二乗和を計算する。 $個体間変動 = \sum\sum(個体別平均 - 総平均)2 = \sum\sum(\bar{x_j}-\bar{\bar{x}})2$
  6. 調整済み郡内変動を求める。郡内変動から個体間変動を引く。 $調整済み郡内変動 = 郡内変動 - 個体間変動$
  7. 個体間変動の自由度を求める。個体数から 1 を引く。 $個体間変動の自由度 = 個体数 - 1$
  8. 調整済み郡内変動の自由度を求める。郡内変動の自由度から個体間変動の自由度を引く。 $調整済み郡内変動の自由度 = 群内変動の自由度-個体間変動の自由度 = \nu_{2}$
  9. 調整済み誤差分散を求める。調整済み郡内変動を調整済み郡内変動の自由度で割る。 $調整済み誤差分散 = \frac{調整済み群内変動}{調整済み群内変動の自由度}$
  10. F 値を求める。要因分散を調整済み誤差分散で割る。 $F_{(\nu_1,\nu_2)} = \frac{要因分散}{調整済み誤差分散}$
  11. p 値を求めて、判定する。F 値と自由度を用いて、p 値を計算する。

R で実装する

F 値と p 値は、対応無しの場合と同様に計算します。

デモデータの作成

デモ用のデータを作成します。このデータは、3 つの水準(A, B, C)を持つ要因(group)と、2 つの水準(a, b)を持つ個体差(individual)によって分類された測定値(value)からなります。

set.seed(123) # 乱数のシードを固定
data <- data.frame(
  group = factor(rep(LETTERS[1:3], each = 4)), # 要因、3水準
  individual = factor(rep(letters[1:2], each = 2)), # 個体差、2個体
  value = unlist(lapply(c(1, 4) * 40, sample, size = 6, replace = TRUE)) # 測定値
)

knitr::kable(sapply(levels(data$group), function(i) {
  sapply(levels(data$individual), function(j) {
    paste(subset(data, subset = group == i & individual == j)$value, collapse = ",")
  })
}))
A B C
a 31,15 37,14 90,91
b 14,3 118,153 91,92

それぞれの個体 (a と b) が、それぞれ要因の各水準 (A、B、C) の値を持っている。このように、1 つの個体が複数の値を持つデータを「対応のある」と言う。

分散分析の関数を使わない

# 総平均
overall_mean <- mean(data$value)
# 個体間平均
individual_mean <- tapply(data$value, data$individual, mean)
# 個体間変動、「個体別平均-総平均」の偏差平方和
ss_individual = sum((individual_mean[data$individual] - overall_mean)^2)
# 個体間変動の自由度、水準数(個体数) - 1
df_individual <- length(levels(data$individual)) - 1

# 個体差調整済み郡内変動、郡内変動-個体間変動
ss_within_new <- ss_within - ss_individual
# 個体差調整済み郡内変動の自由度、郡内変動の自由度-個体間変動の自由度=(総データ数-処理数)-(個体数-1)
df_within_new <- df_within - df_individual
# 個体差調整済み誤差分散
ms_within_new <- ss_within_new / df_within_new

# F値、要因分散/誤差分散
f_value <- ms_between / ms_within_new

# F値からp値を求める
p_value <- pf(f_value, df_between, df_within_new, lower.tail = FALSE)

# 結果の表示
cat("F-value:", round(f_value, 2), "\n")
cat("p-value:", round(p_value, 4), "\n")
F - value: 5.17
p - value: 0.0362

分散分析の関数を使う

対応のある場合は、要因と個体を$+$で繋ぐ。

# 対応のある一元配置分散分析
res <- aov(value ~ group + individual, data = data)
summary(res)
            Df Sum Sq Mean Sq F value Pr(>F)
group        2  13287    6644   5.173 0.0362 *
individual   1   3104    3104   2.417 0.1586
Residuals    8  10274    1284
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F 値は 5.17 で、p 値は 0.0362。 有意水準が 5%なら0.05 より小さいので、有意差あり、平均値に差がある、となる。

分散分析のメモ・その1

分散分析と F 検定の関係が分からなくなることが多いのでメモ

一元配置分散分析

分散分析は、3 つ以上の群の平均値に差があるかどうかを評価する統計手法で、異なる条件や処理が平均値に与える影響を評価できる。な お、一元配置分散分析は、要因数が 1 種類の分散分析。

各々の群について、処理の影響を示す群間分散と処理と関係のない群内分散を求める。次に、群内分散に対する群間分散の比を求める。そして、この比を F 値として、F 分布から p 値を求める。なお、このときの帰無仮説は「すべての群の平均値が等しい」になる。

対応がない場合の一元配置分散分析

同じ個体や対象から複数回測定されたものではないとき、「対応がない」になる。

F 値の式

\begin{eqnarray}

F_{(\nu_1,\nu_2)} &=& \frac{要因分散}{誤差分散}
&=& \frac{\frac{群間変動}{群間変動の自由度}}{\frac{群内変動}{群内変動の自由度}}
&=& \frac{\frac{反復数 \times \sum(群平均 - 総平均)^2}{処理数 - 1}}{\frac{\sum\sum(測定値 - 群平均)^2}{データ総数 - 処理数}}
&=& \frac{\frac{n\sum(\bar{x_i}-\bar{\bar{x}})^2}{\nu_1}}{\frac{\sum\sum(x-\bar{x_i})^2}{\nu_2}}

\end{eqnarray}
  • $F_{(\nu_1,\nu_2)}$F は、群間の平均値のばらつき(要因分散)が、群内のデータのばらつき(誤差分散)に比べてどれだけ大きいかを表す指標。
  • $\nu_1$ は、群間変動の自由度で、処理数(群の数)から 1 を引いた値。
  • $\nu_2$ は、群内変動の自由度で、データ総数から処理数を引いた値。
  • $n$ は、各群のデータ数(反復数)。
  • $x$ は、各データの測定値。
  • $\bar{x_i}$ は、i 番目の群の平均値。
  • $\bar{\bar{x}}$ は、すべてのデータの平均値(総平均)。

計算の流れ

  1. 群間変動を求める。各グループの平均と全体の平均の差の二乗和を計算する。

$群間変動 = 反復数 \times \sum(群平均 - 総平均) = n\sum(\bar{x_i}-\bar{\bar{x}})2$

  1. 群間変動の自由度を求める。処理数(グループ数)から 1 を引く。

$群間変動の自由度 = 処理数 - 1$

  1. 群内変動を求める。各データと、そのグループの平均の差の二乗和を計算する。

$群内変動\sum\sum(測定値 - 群平均)2$

  1. 群内変動の自由度を求める。総データ数から処理数を引く。

$群間変動の自由度 = データ総数 - 処理数$

  1. 要因分散を求める。群間変動を群間変動の自由度で割る。

$要因分散 = \frac{群間変動}{群間変動の自由度}$

  1. 誤差分散を求める。群内変動を群内変動の自由度で割る。

$誤差分散 = \frac{群内変動}{群内変動の自由度}$

  1. F 値を求める。要因分散を誤差分散で割る。

$F値 = \frac{要因分散}{誤差分散}$

  1. p 値を求めて、判定する。F 値と自由度を用いて、p 値を計算する。

R で実装する

対応のない場合を実装してみる。

デモデータの作成

3 つの水準(A, B, C)を持つ要因(group)と、2 つの水準(a, b)を持つ個体差(individual)によって分類された測定値(value)からなるデータを作る。

set.seed(123) # 乱数のシードを固定
data <- data.frame(
  group = factor(rep(LETTERS[1:3], each = 4)), # 要因、3水準
  individual = factor(rep(letters[1:2], each = 2)), # 個体差、2個体
  value = unlist(lapply(c(1,4)*40, sample, size = 6, replace = TRUE)) # 測定値
)

knitr::kable(t(sapply(levels(data$group), function(i){
  paste(subset(data, subset=group==i)$value,collapse=",")
})))
A B C
31,15,14,3 37,14,118,153 90,91,91,92

この表は、各群と各個体に対応する測定値を持つが、対応がない場合は、個体差(a, b)がいらない。 このように、各水準が複数の値を持ち、水準以外の区別がないデータを「対応のある」と言う。

分散分析の関数を使わない

# 総平均
overall_mean <- mean(data$value)
# 群平均
group_means <- tapply(data$value, data$group, mean)

# 群間変動、「値-群平均」の偏差平方和
ss_between <- sum((group_means[data$group] - overall_mean)^2)
# 群間変動の自由度、水準数(処理数) - 1
df_between <- length(levels(data$group)) - 1
# 要因分散
ms_between <- ss_between / df_between

# 群内変動、誤差、「群平均-総平均」の偏差平方和
ss_within <- sum((data$value - group_means[data$group])^2)
# 群内変動の自由度、総データ数-水準数(処理数)
df_within <- length(data$value) - length(levels(data$group))
# 誤差分散
ms_within <- ss_within / df_within

# F値、要因分散/誤差分散
f_value <- ms_between / ms_within

# F値からp値を求める
p_value <- pf(f_value, df_between, df_within, lower.tail = FALSE)

# 結果の表示
cat("F-value:", round(f_value, 2), "\n")
cat("p-value:", round(p_value, 4), "\n")
F-value: 4.47
p-value: 0.0449

分散分析の関数を使う

# 対応のない一元配置分散分析
res <- aov(value ~ group, data = data)
summary(res)
            Df Sum Sq Mean Sq F value Pr(>F)
group        2  13287    6644    4.47 0.0449 *
Residuals    9  13378    1486
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F 値は 4.47 で、p 値は 0.0449。

有意水準が 5%なら0.05 より小さいので、有意差あり、平均値に差がある、となる。

LARSの整理・その3

LASSOをまとめていると、LARS(最小角回帰法:Least Angle Regression)が出てきたので、これを機会に手順をちゃんと理解しようと思う。

行列と角度、高校や大学で習ったものの、いざ実データで、となるとどうしてもピンと来ない。実装してようやく理解できたので、その記録。

LARS をまとめる

LARS を理解する上で、引っかかったのは「最小二乗法」と「プロジェクションの方向」でした。

最小二乗法

高校までは二次元だった最小二乗法は、大学で多次元になります。行列で処理するわけですが、行列に苦手意識があって、理解が曖昧になりがちです。改めて、最小二乗法を整理してみると、意外と理解できました。

実装したコードだと、最小二乗法は、G <- solve(t(X_active) %*% X_active) %*% t(X_active)で使っていて、いわゆる疑似逆行列を求めて係数$\beta$を推定しています。

  1. 残差の二乗和(RSS)を、目的変数を(y)、予測値を(\hat{y})、その残差を(\varepsilon)は、次の関係になる。 [ \varepsilon = y - \hat{y} ] [ \text{RSS} = \varepsilonT \varepsilon = (y - \hat{y})T (y - \hat{y}) ]

  2. また、説明変数行列(X)と係数ベクトル(\beta)の積が予測値(\hat{y})になる。 [ \hat{y} = X \beta ]

  3. (\beta)を微分して RSS を 0 にすると、式は$y$、$X$、$\beta$の式になる。 [ \frac{\partial \text{RSS}}{\partial \beta} = \frac{\partial}{\partial \beta} \left( (y - X \beta)T (y - X \beta) \right) = 0 ] [ \frac{\partial \text{RSS}}{\partial \beta} = -2 XT (y - X \beta) = 0 ] [ -2 XT X \beta + 2 XT y = 0 ] [ XT X \beta = XT y ]

  4. (X^{T} X)が正則なので、疑似逆行列を使って式を変形する。 [ (XT X)^{-1}(XT X) \beta = (XT X)^{-1}XT y ] [ I \beta = (XT X)^{-1}XT y ] [ \beta = (XT X)^{-1} XT y ]

プロジェクションの方向

LARS は、残差ベクトルと絶対値が最も高い説明変数を選び、アクティブセットと呼ばれる変数の集合に加えます。説明変数の係数を更新する際に、残差ベクトルとアクティブセット内の説明変数の符号付き線形結合ベクトルとの内積が最大になる(相関が等しくなる)ようにします。なお、符号付き線形結合ベクトルとは、説明変数の係数にその相関係数の符号(正か負か)を掛けたものです。残差ベクトルと符号付き線形結合ベクトルの内積が最大になると(相関が等しくなると)、なす角が 0 に近づきます。これが「最少角」の名前の由来です。

このとき、符号付き線形結合ベクトルがプロジェクションの方向を示します。例えば、アクティブセット $X$ が $x_1$ と $x_2$ からなる行列で、係数ベクトルが $\beta$ とすると、残差ベクトルは $y - X\beta$ になります。LARS はこの $y - X\beta$ と $x_1$ と $x_2$ との内積が最大に(相関が等しく)なるようにします。つまり、

$$ XT(y - X\beta) = (c_1, c_2)T $$

$$ XTy - XTX\beta = (c_1, c_2)T $$

$$ \beta = (XTX)^{-1}XTy - (XTX)^{-1}(c_1, c_2)T $$

ここで、$c_1$ と $c_2$ はそれぞれ $y - X\beta$ と $x_1$ と $x_2$ との内積です。この式は、残差ベクトルとプロジェクション方向(符号付き線形結合ベクトル)のなす角が最小(0 に近い)ことを意味します。

この $(XTX)^{-1}XTy$ がプロジェクション項、$(XTX)^{-1}(c_1, c_2)T$ が補正項になります。プロジェクション項は、教師データをアクティブセット内の説明変数で近似したものです。これは、教師データをアクティブセット内の説明変数で張られた部分空間に射影したものに相当します。補正項は、プロジェクション項から係数ベクトルをずらすためのものです。これは、残差ベクトルとプロジェクション方向が同じ角度になるように設定されます。

このとき、$s$ を補正項の大きさを決めるインデックス、相関ベクトルのノルムを $C = \sqrt{c_12 + c_22}$ と置くと、

$$ (XTX)^{-1}(c_1, c_2)T = (XTX)^{-1}C(c_s/C)T $$

と変換できます。なお、この相関ベクトルをノルムで割った $(c_s/C)$ は、その $s$ 番目の要素以外がゼロであるベクトルになります。

LARS では、残差ベクトルとアクティブセット外の説明変数との内積(相関)が等しくなるように$s$を選びます。これにより、補正項が最も大きくなります。補正項が大きくなれば、係数ベクトルはプロジェクション項から離れます。つまり、プロジェクション方向は残差ベクトルと同じ角度になるように決められます。

LASSO回帰の復習・その4

LASSO回帰を実装する

Rのmtcarsを使った実装例を作った。

回帰係数(w)と正則化パラメータ($\lambda$)の初期値は適当に決めている。

回帰係数の初期値

ゼロベクトルや小さなランダムノイズを含んだベクトルなどがある。

  1. ゼロベクトル: 回帰係数をすべてゼロでスタートすると、最小二乗法の推定と同じくらいの大きさになる。
  2. ランダムノイズ: 小さなランダムノイズでスタートすると、初期値からバラバラになるので、局所最適解に陥る可能性が軽減する。
  3. データの特性に基づく初期値: データの特性などに基づいて選択する。例えば、特定の変数が目的変数と関連性が高いときにその係数を初期値にする、など。

正則化パラメータ

$\lambda$の大きさは、残る変数の数に影響する。$\lambda$が大きいと0になる係数が増え、変数は残らない。

$\lambda$の値は実際のデータやモデルの複雑さに依存する、つまり明確な設定基準はなさそう。

モデルの予測性能を探索して、その結果から$\lambda$を選択するようだ。

  1. クロスバリデーション(Cross-validation): データを訓練データとテストデータに分割して、モデルを評価する。複数の$\lambda$に対して評価を繰り返して、予測性能が最も高い$\lambda$を選ぶ。
  2. 情報量基準(Information Criterion): AIC赤池情報量基準)やBICベイズ情報量基準)などを使用して、モデルの適合度とモデルの複雑さのバランスが良い$\lambda$を選ぶ。
  3. パス探索(Path Searching): 一連の$\lambda$の値でのでのLasso回帰の推定結果を求めて、変数の選択状況や回帰係数の変化を視覚的に確認して、適切な$\lambda$を選ぶ。

回帰係数wの求め方

i番目の回帰係数は、説明変数と残差の内積と$\lambda$の大きさできまる。

この内積が$\pm \frac{\lambda}{2}$に入れば、$\lambda$と説明変数の大きさで調整した値、入らなければ0になる。

もちろん、説明変数が全て0なら、回帰係数は0になる。

説明変数X[,i]と残差r内積 回帰係数w[i]
説明変数と残差の内積 < $\frac{\lambda}{2}$ $\lambda$ の半分を足して説明変数の二乗和で割る
説明変数と残差の内積 < -$\frac{\lambda}{2}$ $\lambda$ の半分を引いて説明変数の二乗和で割る
その外側 0
説明変数が全て0 0
lasso <- function(X, y, w, lambda) {
  # 座標降下法(Coordinate Descent)でLasso回帰の係数を計算する関数
  ## X: 説明変数のマトリクス
  ## y: 目的変数のベクトル
  ## lambda: L1正則化のλ

  # ラムダの半分
  lambda2 <- lambda / 2
  for (i in 1:ncol(X)) {
    # 残差の二乗和とL1正則化項の和をrにする
    r <- sum((y - X %*% w)^2) + lambda * sum(abs(w))
    if (is.na(r) | is.infinite(r)) r <- 0 # rを計算できないときのエラー処理

    # 説明変数の列と残差の内積
    sumXr <- sum(X[, i] * r)
    # 説明変数の列の二乗和
    sumX2 <- sum(X[, i]^2)

    if (sumX2 == 0) {
      w[j] <- 0
    } else {
      if (sumXr > lambda2) {
        w[i] <- (sumXr - lambda2) / sumX2
      } else if (sumXr < -lambda2) {
        w[i] <- (sumXr + lambda2) / sumX2
      } else {
        w[i] <- 0
      }
    }
  }
  return(w)
}

data(mtcars)
X <- as.matrix(mtcars[, -1])
y <- mtcars$mpg
w = rep(0, ncol(X))
names(w) <- colnames(X)
lambda <- .0000000001

lasso(apply(X,2,scale), y, w, lambda)

これで、L1正則化を利用した回帰係数が分かる。

でも実際は、目的変数と予測値の残差が小さくなる(収束する)まで繰り返すことになる。

実際は、最大繰り返し数(max_iter)と収束値(tol)を組み込むことになる。

lasso_iter <- function(X, y, lambda, max_iter = 1000, tol = lambda/1000) {
  # 座標降下法(Coordinate Descent)でLasso回帰の係数を計算する関数
  ## X: 説明変数のマトリクス
  ## y: 目的変数のベクトル
  ## lambda: L1正則化のλ
  # 収束判定するための引数
  ## max_iter: 繰り返し回数の上限
  ## tol: 収束判定の下限値

  # 回帰係数の初期値
  w <- rep(0, ncol(X))
  names(w) = colnames(X)

  # 収束の確認フラグ
  converged <- FALSE

  # 反復回数のカウント
  iter <- 0

  # 収束するか最大反復回数に達するまで繰り返す
  while (!converged && iter < max_iter) {
    # 前の回帰係数のベクトルを保存する
    w_old <- w

    # 各説明変数に対して回帰係数を更新する
    w = lasso(X, y, w, lambda)

    # 反復回数を増やす
    iter <- iter + 1

    # 前の回帰係数と現在の回帰係数の差分ベクトルのノルムを計算する
    diff_norm <- sqrt(sum((w - w_old)^2))

    # 差分ベクトルのノルムが許容誤差より小さくなったら収束したと判断する
    if (diff_norm < tol) {
      converged <- TRUE
    }
  }

  # 収束判定、反復回数、回帰係数を返す
  return(list(converged = converged, iter = iter, w = w))
}

data(mtcars)
X <- as.matrix(mtcars[, -1])
y <- mtcars$mpg
w = rep(0, ncol(X))
names(w) <- colnames(X)
lambda <- .0000000001

lasso_iter(apply(X, 2, scale), y, lambda)

LASSO回帰の復習・その3

Lasso回帰の最適化アルゴリズム

最適化アルゴリズムは座標降下法が有名ですが、調べてみると他にもいくつかありました。 分かる範囲でリストアップしてみた。 それぞれの理解をまとめておく。

名称 特徴 メリット デメリット
座標降下法 1変数ずつ重みの更新を繰り返して
最適解に近似させる
計算が簡単で高速 収束性が保証されない
サイクリック座標降下法 座標降下法の一種
重みを固定した順番で更新する
座標降下法よりも収束が速い 更新順序によっては収束しない
シューティングアルゴリズム 座標降下法の一種
重みをランダムに選んで更新する
重みの更新順序に依存しない 収束までに多くの反復が必要
LARS(最小角回帰法) 重みを最も急速に増加させる
方向に沿って更新する
変数選択と正則化を同時に行える 計算量が大きい
プロキシミティ演算子 重みを現在の値から一定量だけ動かす 非線形回帰や非凸最適化問題にも適用できる 収束速度が遅い

関数にする

数式を見ていても、今ひとつピンとこないので、動作の理解も含めてコーディングしてみた。

座標降下法(Coordinate Descent)

  1. 残差の二乗和とL1正則化項を計算し、rとして保存する。
  2. もし説明変数の列がすべて0の場合は、回帰係数を0にする。
  3. そうでなければ、回帰係数を以下の条件で計算する。
    • もし「説明変数の列と残差の内積」が「ラムダの半分」より大きければ、回帰係数を「説明変数の列と残差の内積 + lambda / 2」を「説明変数の列の二乗和」で割った値にする。
    • もし「説明変数の列と残差の内積」が「ラムダの半分の負値」より小さければ、回帰係数を「説明変数の列と残差の内積 - lambda / 2」を「説明変数の列の二乗和」で割った値にする。
    • それ以外は、回帰係数を0にする。

ゼロを中心にした$\lambda$の範囲外なら回帰係数はゼロ、範囲内なら非ゼロという感じです。

説明変数X[,i]と残差r内積 回帰係数w[i]
説明変数と残差の内積 < $\lambda$の半分 $\lambda$ の半分を足して説明変数の二乗和で割る
説明変数と残差の内積 < -$\lambda$の半分 $\lambda$ の半分を引いて説明変数の二乗和で割る
その外側 0
説明変数が全て0 0
lasso <- function(X, y, w, lambda) {
  # ラムダの半分
  lambda2 <- lambda / 2
  for (i in 1:ncol(X)) {
    # 残差の二乗和とL1正則化項の和をrにする
    r <- sum((y - X %*% w)^2) + lambda * sum(abs(w))
    if (is.na(r) | is.infinite(r)) r <- 0 # rを計算できないときのエラー処理

    # 説明変数の列と残差の内積
    sumXr <- sum(X[, i] * r)
    # 説明変数の列の二乗和
    sumX2 <- sum(X[, i]^2)

    if (sumX2 == 0) {
      w[j] <- 0
    } else {
      if (sumXr > lambda2) {
        w[i] <- (sumXr - lambda2) / sumX2
      } else if (sumXr < -lambda2) {
        w[i] <- (sumXr + lambda2) / sumX2
      } else {
        w[i] <- 0
      }
    }
  }
  return(w)
}

サイクリック座標降下法(Cyclic Coordinate Descent)

  1. 残差をrにする。
  2. 説明変数の列ごとに次の手順を繰り返す。
    • 「説明変数の列と残差の内積」をzにする。
    • 既存の回帰係数w[-i]の符号とlambda/nrow(X)を掛けた値をaにする。
    • もし、zの絶対値がlambda/nrow(X)より大きければ、回帰係数をaにする。
    • それ以外なら、回帰係数を0にする。
lasso <- function(X, y, w, lambda) {
  for (i in 1:ncol(X)) {
    r <- y - X %*% w
    z <- X[, i] %*% r / nrow(X)
    a <- z + sign(w[-i]) * lambda / nrow(X)
    w[j] <- ifelse(abs(z) > lambda / nrow(X), a, 0)
  }
  return(w)
}

シューティングアルゴリズム(Shooting Algorithm)

  1. 「説明変数の列と目的変数の内積」の符号から回帰係数の符号を決める。
  2. abs(x * y) - lambda / 20以下に制約した結果を、max(abs(x * y) - lambda / 2, 0) にする。
  3. この値を説明変数の列の二乗和sum(x^2)で割り、回帰係数にする。
lasso <- function(X, y, w, lambda) {
  for (i in 1:ncol(X)) {
    pm <- sign(sum(X[, i] * y))
    w[i] <- pm * max(abs(X[, i] * y) - lambda / 2, 0) / sum(X[, i]^2)
  }
  return(w)
}

最小角回帰法(LARS: Least Angle Regression)

  1. 最も相関の高い説明変数を見つける。
  2. アクティブな説明変数の部分行列を取得し、最小二乗法で回帰係数を計算する。
  3. アクティブな説明変数の係数と残差を更新する。
  4. アクティブな説明変数の係数のノルムを計算し、プロジェクション方向を決定する。

LARSは理解があいまいなので、別でまとめる。

lars <- function(X, y) {
  n <- nrow(X) # サンプル数
  p <- ncol(X) # 説明変数の数
  w <- matrix(0, nrow = p, ncol = p) # 回帰係数の初期値
  mu <- rep(0, n) # 予測値の初期値
  active_set <- c() # アクティブな説明変数のインデックス
  residual <- y # 残差、初期値は目的変数

  for (k in 1:p) {
    # 最も相関の高い説明変数を見つける
    c <- t(X) %*% residual # 相関ベクトル
    c_max <- which.max(abs(c)) # 相関の絶対値が最大のインデックス
    active_set <- c(active_set, c_max) # アクティブな説明変数のインデックスに追加

    # アクティブな説明変数の部分行列を取得し、最小二乗法で回帰係数を推定
    X_active <- X[, active_set, drop = FALSE] # アクティブな説明変数の部分行列
    G <- solve(t(X_active) %*% X_active) %*% t(X_active) # Gは疑似逆行列
    a <- drop(G %*% y) # 回帰係数の推定値
    w[active_set, k] <- a # 係数行列に係数を追加(k列目に格納)

    # 予測値の更新
    mu <- X_active %*% a # アクティブな説明変数のみを用いて予測値を計算
    residual <- y - mu # 残差の更新
  }

  # アクティブセットのインデックス active_set、回帰係数 w 、予測値 mu、残差 residual を返す
  list(active_set = active_set, w = w, mu = mu, residual = residual)
}

プロキシミティ演算子

  1. 入力ベクトル v正則化項の重みパラメータ lambda を受け取る。
  2. L1正則化に基づくプロキシミティ演算子を計算する。
  3. 要素の絶対値から正則化項の重みを引く。
  4. 入力ベクトルの要素ごとに、ゼロ以下にする。
lasso <- function(v, lambda) {
  sign(v) * pmax(abs(v) - lambda, 0)
}

LARSの整理・その2

LASSOをまとめていると、LARS(最小角回帰法:Least Angle Regression)が出てきたので、これを機会に手順をちゃんと理解しようと思う。

行列と角度、高校や大学で習ったものの、いざ実データで、となるとどうしてもピンと来ない。実装してようやく理解できたので、その記録。

LARSとは

LARSは、LASSOを効率的に計算できるスパースモデリングの一種です。

回帰分析において、変数選択と回帰係数の推定を同時に行うアルゴリズムであり、特徴量が非常に多い高次元データセットに対して効率的に特徴選択できる強みがある。

LARS は最小二乗法をベースに次の流れで進みます。

  1. 初期値として全ての回帰係数を 0 にする。
  2. 目的変数となす角度が最も小さい(目的変数との相関係数の絶対値が最も大きい)説明変数を求め、この選択した説明変数をアクティブセットにする。
  3. アクティブセットの線形結合で表される方向に沿って、目的変数との残差が小さくなるように回帰係数を更新する。
  4. 更新した後の残差となす角度が最も小さい説明変数を求め、この選択した説明変数をアクティブセットに加える。
  5. アクティブセットの線形結合で表される方向(プロジェクションの方向)に沿って、残差を減らすように回帰係数を更新する。
  6. この残差とアクティブセット以外の列ベクトルとの角度が同じになるまで、4,5 を繰り返す。
  7. 残差とアクティブセット以外の列ベクトルとの角度が同じになったら、その列ベクトルをアクティブセットに加える。また、残差ベクトルとアクティブセット内の列ベクトルとの角度が 180 度になったら、その列ベクトルをアクティブセットから除く。
  8. 全ての説明変数がアクティブセットになるか、残差が 0 になったら終了する。

実行手順

この手順を書き下してみると、こんな流れになる。

lars <- function(X, y) {
# 説明変数行列 `X` と目的変数ベクトル `y` が引数。
# 変数 `n` と `p` は、それぞれサンプル数と説明変数の数。
# `w` は回帰係数の初期値を保持する行列で、0が初期値。
# `mu` は予測値の初期値を保持するベクトルで、0 が初期値。
# `active_set` はアクティブな説明変数のインデックスを保持するベクトルで、空のベクトルが初期値。
# `residual` は残差を保持するベクトルで、初期値は目的変数 `y`。
  n <- nrow(X) # サンプル数
  p <- ncol(X) # 説明変数の数
  w <- matrix(0, nrow = p, ncol = p) # 回帰係数の初期値
  mu <- rep(0, n) # 予測値の初期値
  active_set <- c() # アクティブな説明変数のインデックス
  residual <- y # 残差、初期値は目的変数

# 変数 `k` のループは、説明変数の数 `p` に対して実行されます。
  for (k in 1:p) {
    # `k` が `p` より小さければ、アクティブな説明変数の係数と残差を更新
    c <- t(X) %*% residual # 相関ベクトル
    c_max <- which.max(abs(c)) # 相関の絶対値が最大のインデックス
    active_set <- c(active_set, c_max) # `active_set` に `c_max` を追加

    # アクティブな説明変数の部分行列を取得し、最小二乗法で回帰係数を推定
    X_active <- X[, active_set, drop = FALSE] # アクティブな説明変数の部分行列
    G <- solve(t(X_active) %*% X_active) %*% t(X_active) # Gは疑似逆行列
    a <- drop(G %*% y) # 回帰係数の推定値
    w[active_set, k] <- a # 係数行列に係数を追加(k列目に格納)

    # 予測値の更新
    mu <- X_active %*% a # アクティブな説明変数のみを用いて予測値を計算
    residual <- y - mu # 残差の更新

  # アクティブセットのインデックス active_set、回帰係数 w 、予測値 mu、残差 residual を返す
  list(active_set = active_set, w = w, mu = mu, residual = residual)
}

「最小二乗法」を振り返ったこと、「プロジェクションの方向」を幾何的に理解できたこと、この2つの理解がLARSの理解に大きくかかわりました。