日々のつれづれ

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

モザイクプロットを活用したクロス集計の可視化

R Python

モザイクプロットとは

モザイクプロットは、データの分布や変数間の関係性を視覚的に表現するグラフの一種です。データを行と列のセルに分割し、各セルの面積を値の大きさに応じて調整することで、データの全体像を一目で把握できるようになっています。このグラフは、複雑なデータ構造を簡潔に表現することができるため、データ分析の初期段階で有効活用されます。

グラフの見方

モザイクプロットでは、各セルの面積が値の大きさを表しています。大きなセルは値が大きいことを、小さなセルは値が小さいことを示しています。

利用できる状況

モザイクプロットは、主に以下のような状況で活用されます:

  • 複数の変数間の関係性を一度に確認したい場合
  • データの全体的な分布傾向を把握したい場合
  • 2 つの変数間の相関関係を把握したい場合

これらの状況では、モザイクプロットが有効な可視化手段となります。データの構造が複雑な場合でも、一目で全体像を把握できるため、分析の初期段階で活用されることが多いです。

デモデータを使った作図

では、実際にデモデータを使ってモザイクプロットを作成してみましょう。

データの準備

まずは、サンプルデータを用意します。今回は、ある企業の従業員データとして、従業員の性別、年齢、所属部門の情報が含むデータを作成します。

set.seed(123)
# データの準備
dat <- data.frame(
  "DEPT." = sample(
    c("Sales", "General Affairs", "Development"),
    size = 20, replace = TRUE, prob=c(3,1,6)/10
  ),
  "SEX" = sample(c("M", "F"), size = 20, replace = TRUE, prob=c(1,4)/5),
  "AGE" = sample(25:30, size = 20, replace = TRUE)
)

knitr::kable(dat)

クロス集計表の作成

このとき、部門ごとの男女構成を集計します。クロス集計表は、行と列の変数の組み合わせによる度数を表したものです。

tbl = table(dat[,1:2])
knitr::kable(tbl)
F M
Development 10 2
General Affairs 2 1
Sales 4 1

この集計表を視覚的な理解を助けるグラフがモザイクプロットです。

モザイクプロットの作成

上記のデータを使って、モザイクプロットを作成します。モザイクプロットは、クロス集計表の情報を視覚的に表現したグラフです。各セルの面積が値の大きさを表します。 つまり、人数が一番多い Development の Female の面積が一番大きく、人数が一番少ない General Affair の Male と Sales の Male の面積が小さくなります。

X 軸の幅が各水準の合計比率を表していて、

$$Developmentの幅 : General Affairの幅 : Salesの幅 = Developmentの人数 : General Affairの人数 : Salesの人数$$

になります。

# モザイクプロットの作成
mosaicplot(x = tbl)

このように、モザイクプロットを使うことで、従業員の性別と部門の関係性を一目で把握することができました。このグラフを見ることで、例えば開発部門には男性が多く、総務部門には女性が多いといった傾向を簡単に読み取ることができます。

plot.tableの機能（補足）

モザイクプロットはtable(x)に対するグラフになるので、R ではplot.table(x)でも作図できます。

plot(x=tbl) # 引数がtableオブジェクトであれば、plotはplot.tableとして機能する

Python の例（補足）

Python では、pandasのpd.crosstab関数を使ってクロス集計表を作成し、statsmodels.graphics.mosaicplot.mosaic関数を使ってモザイクプロットを描画することができます。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.graphics.mosaicplot import mosaic

# 乱数シードの設定
np.random.seed(123)

# データの準備
dept_choices = ["Sales", "General Affairs", "Development"]
dept_probs = [0.3, 0.1, 0.6]
sex_choices = ["M", "F"]
sex_probs = [0.2, 0.8]
age_range = range(25, 31)

dept = np.random.choice(dept_choices, size=20, p=dept_probs)
sex = np.random.choice(sex_choices, size=20, p=sex_probs)
age = np.random.choice(age_range, size=20)

df = pd.DataFrame({"DEPT.": dept, "SEX": sex, "AGE": age})

# クロス集計表の作成と表示
tbl = pd.crosstab(df["SEX"], df["DEPT."])
print("\nクロス集計表:")
print(tbl)

# モザイクプロットの描画
plt.figure(figsize=(8, 6))
mosaic(df, ['DEPT.', 'SEX'], title='Gender composition in each department')
plt.show()

pythonのモザイクプロットはデフォルトで色分けになります。このようにセルを色分けすると、データの傾向をより視覚的に表現できます。

例えば、赤色のセルは値が高く、青色のセルは値が低いといった具合にすると、グラフはより直感的な表現になります。

回帰モデルの係数の分散分析のまとめ・その３

R 勉強会

R Advent Calender 2023として投稿することにしました

共分散分析で回帰分析の分散分析が出てきたので、少し整理します。

回帰分析を分散分析する目的

回帰分析の分散分析は、データを意思決定に利用するための方法です。

これを使えば、回帰モデルに用いた異なる要因（説明変数、因子など）が目的変数に与える影響を統計的に評価できます。

一般的な回帰モデルは以下の式で表されます。

$$ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_n x_n + \varepsilon $$

  • $y$ : 目的変数
  • $x_1, x_2, \ldots, x_n$ : 説明変数、要因
  • $\beta_0, \beta_1, \beta_2, \ldots, \beta_n$ : 回帰係数（各説明変数の影響度）
  • $\varepsilon$ : 誤差項

この回帰モデルにおいて、異なる説明変数が目的変数に与える影響を評価するために分散分析を利用します。そして、回帰モデルや説明変数により目的変数の総変動（ばらつき）が減少（説明）する量を見ることで、その影響が分かります。

これには、以下の 2 つのアプローチがあります。

  • 回帰モデル全体の有意性評価：回帰モデルが目的変数をどれだけ説明できるかを見るために、回帰モデルと単なる平均値（切片だけのモデル）を比較する。
  • 特定の説明変数や切片に対する検定：各説明変数や切片が目的変数にどれだけ影響を与えるかを見るために、その係数が 0 であるかどうかを検定する。

重回帰モデルの場合

重回帰モデルは説明変数が複数あるため、「回帰モデル全体の有意性評価」と「特定の説明変数や切片に対する検定」の 2 つのアプローチは異なります。

回帰モデル全体の有意性評価の計算手順

ここでは、回帰モデル全体の有意性評価を整理します。この評価では、回帰モデルが目的変数を説明できる程度を見ています。

次の計算手順は、サンプル数が $m$、説明変数の総数が $n$の例になります。

  1. 目的変数の総変動を回帰変動と残差変動に分解する
    • $SS_y$ : 目的変数の総変動、目的変数のばらつきの程度
      • $y_i$ : $i$ 番目の実測値
      • $\bar{y}$ : 目的変数の平均値 $$ SS_y = \sum_{i = 1}^{n}(y_i - \bar{y})2 $$
    • $SS_R$ : 回帰変動、回帰モデルが目的変数を説明できる程度
      • $\hat{y}i$ : $i$ 番目の予測値 $$ SS_R = \sum{i = 1}^{n}(\hat{y}_i - \bar{y})2 $$
    • $SS_E$ : 残差変動、回帰モデルが目的変数を説明できない程度 $$ SS_E = \sum{i = 1}^{n}(y_i - \hat{y}i)2 $$
  2. それぞれの変動を自由度で割って平均二乗を求める
    • $df_y$ : 全体の自由度、観測値の数 - 1 $$ df_y = m - 1 $$
    • $df_R$ : 回帰の自由度、説明変数の個数 $$ df_R = n $$
    • $df_E$ : 残差の自由度、全体の自由度 - 回帰の自由度 $$ df_E = df_y - df_R = m - n - 1 $$
    • $MS_y$ : 全体の平均二乗、総変動を全体の自由度で割った値 $$ MS_y = \frac{SS_y}{df_y} $$
    • $MS_R$ : 回帰の平均二乗、回帰変動を回帰の自由度で割った値 $$ MS_R = \frac{SS_R}{df_R} $$
    • $MS_E$ : 残差の平均二乗、残差変動を残差の自由度で割った値 $$ MS_E = \frac{SS_E}{df_E} $$
  3. 回帰の平均二乗を残差の平均二乗で割って、F 統計量$F_R$を求める $$ F_R = \frac{MS_R}{MS_E} $$
  4. 回帰と残差の自由度の F 分布から p 値$p_R$を求めて有意水準で判断する $$ p_R = P(F > F_R) = 1 - P(F \leq F_R) $$ なお、このときの帰無仮説は
    • 帰無仮説 ( $H_0$ ): $\beta_1 = 0$（説明変数は目的変数に影響しない）
    • 対立仮説 ( $H_1$ ): $\beta_1 \neq 0$（説明変数は目的変数に影響する）

R で計算する

回帰分析で分散分析すると、異なる要因（説明変数、因子など）が目的変数に与える影響を統計的に評価できます。そして、これはデータに基づいた合理的な意思決定に有効です。

説明変数が 2 つの重回帰モデルを例に確認します。式は次のようになります。

できるだけスクラッチ実装をしてみます。

$$ y{mpg} = \beta_0 + \beta_1 x{wt} + \beta{2} x{hp} + \varepsilon $$

これは、mtcarsデータセットを使った回帰モデルで、目的変数が mpg（ガロンあたりのマイル数）、説明変数が wt（車両重量）と hp（馬力）です。そして、wt と hp に交互作用はないと考えています。

# データの読み込み
data(mtcars)

# 回帰モデルの作成
model <- lm(mpg ~ wt + hp, data = mtcars)
model

Call:
lm(formula = mpg ~ wt + hp, data = mtcars)

Coefficients:
(Intercept)           wt           hp
   37.22727     -3.87783     -0.03177

係数と切片の推定値が、それぞれ $\hat{\beta}0 = 37.22727$ 、$\hat{\beta}1 = -3.87783$ 、$\hat{\beta}_2 = -0.03177$ になりました。

続いて、回帰モデル全体の有意性を評価します。計算の流れは単回帰モデルと同じです。

  1. 目的変数の総変動 $SS_y$ を回帰変動 $SS_R$ と残差変動 $SS_E$ に分解する 回帰変動は、回帰モデルの目的変数の平均値からのズレの程度です。良い回帰モデルほど値は大きくなります。
# 目的変数とその推定値
y <- mtcars$mpg # 目的変数
bar_y <- mean(y) # 平均値
hat_y <- fitted(model) # 推定値
# hat_y <- predict(model) # 推定値

# 総変動と回帰変動と残差変動
SS_y <- sum((y - bar_y)^2) # 総変動
SS_R <- sum((hat_y - bar_y)^2) # 回帰変動
SS_E <- sum((y - hat_y)^2) # 残差変動
  1. それぞれの自由度で割って、それぞれの平均二乗 $MS_R$ と $MS_E$ を求める
# 自由度
df_y <- nrow(mtcars) - 1 # 全体の自由度
df_R <- length(coef(model)) - 1 # 回帰の自由度（説明変数の数）
df_E <- df_y - df_R # 残差の自由度

# 平均二乗、分散
MS_y <- SS_y / df_y # 全体の平均二乗
MS_R <- SS_R / df_R # 回帰の平均二乗
MS_E <- SS_E / df_E # 残差の平均二乗
  1. 回帰の平均二乗を残差の平均二乗で割って、F 統計量を計算する F 統計量は、回帰モデルが目的変数をどれだけ説明できるかを表す指標です。回帰モデルの当てはまりが良ければ、F 統計量は大きくなります。
# F統計量の計算
F_R <- MS_R / MS_E
  1. 自由度が（$df_R$, $df_E$）の F 分布から p 値を求めて有意水準で判断する p 値は、F 統計量が得られる確率を表す指標です。回帰モデルが目的変数に有意な影響を与えていれば、p 値が小さくなります。
# p値の計算
p_R <- pf(F_R, df_R, df_E, lower.tail = FALSE)
# p_R <- 1 - pf(F_R, df_R, df_E)
  1. 分散分析表を作成する
# 分散分析表の作成
res <- data.frame(
   Df = c(df_R, df_E, df_y),
   `Sum Sq` = c(SS_R, SS_E, SS_y),
   `Mean Sq` = c(MS_R, MS_E, MS_y),
   `F value` = c(F_R, NA, NA),
   `Pr(>F)` = c(p_R, NA, NA),
   row.names = c("Regression", "Residuals", "Total"),
   check.names = FALSE
)
res

           Df    Sum Sq    Mean Sq  F value       Pr(>F)
Regression  2  930.9994 465.499716 69.21121 9.109054e-12
Residuals  29  195.0478   6.725785       NA           NA
Total      31 1126.0472  36.324103       NA           NA

なお、これはanova(model)では取得できません。

summary(model)の返り値を使って計算します。

f_stat = summary(model)$fstatistic
p_value = pf(f_stat[1], f_stat[2], f_stat[3], lower.tail=FALSE)
cat("F statistic: ", f_stat[1], "\n")
cat("p value: ", p_value, "\n")

F statistic:  69.21121
p value:  9.109054e-12

回帰モデルの係数の分散分析のまとめ・その２

R

共分散分析で回帰分析の分散分析が出てきたので、少し整理します。

回帰分析を分散分析する目的

回帰分析の分散分析は、データを意思決定に利用するための方法です。

これを使えば、回帰モデルに用いた異なる要因（説明変数、因子など）が目的変数に与える影響を統計的に評価できます。

一般的な回帰モデルは以下の式で表されます。

$$ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_n x_n + \varepsilon $$

  • $y$ : 目的変数
  • $x_1, x_2, \ldots, x_n$ : 説明変数、要因
  • $\beta_0, \beta_1, \beta_2, \ldots, \beta_n$ : 回帰係数（各説明変数の影響度）
  • $\varepsilon$ : 誤差項

この回帰モデルにおいて、異なる説明変数が目的変数に与える影響を評価するために分散分析を利用します。そして、全ての$x_n$が$y$に与える影響を見る「回帰モデル全体の有意性評価」と、個々の$x_n$が$y$に与える影響を個別にみる「特定の説明変数や切片に対する検定」の 2 つのアプローチがあります。

R で計算する

まずは、説明変数が 1 種類の単回帰モデルを考えます。単回帰モデルの分散分析は、説明変数が目的変数をどれだけ説明できるかを見ています。

$$y = \beta_0 + \beta_1 x + \varepsilon$$

なお、単回帰モデルの場合、説明変数が 1 つしかないので、「回帰モデル全体の有意性評価」と「特定の説明変数や切片に対する検定」は同じです。

ここからは、mtcarsデータセットを使い、目的変数を mpg（ガロンあたりのマイル数）、説明変数を wt（車両重量）にした回帰モデルで進めます。

# データセットの読み込み
data(mtcars)

# モデルの適合
model <- lm(mpg ~ wt, data = mtcars)

model

Call:
lm(formula = mpg ~ wt, data = mtcars)

Coefficients:
(Intercept)           wt
     37.285       -5.344

係数と切片の推定値が、$\hat{\beta}0 = 37.285$ 、$\hat{\beta}1 = -5.344$になりました。

続いて、回帰モデル全体の有意性を評価します。計算の流れは以下の通りです。

  1. 目的変数の総変動 $SS_y$ を回帰変動 $SS_R$ と残差変動 $SS_E$ に分解する
  2. それぞれの自由度で割って、それぞれの平均二乗 $MS_R$ と $MS_E$ を求める
  3. 回帰の平均二乗を残差の平均二乗で割って、F 統計量を計算する
  4. 自由度が（$df_R$, $df_E$）の F 分布から p 値を求めて有意水準で判断する

目的変数の総変動 $SS_y$

目的変数の総変動（$SS_y$）は、目的変数 $y$ の値がその平均値 $\bar{y}$ からどれだけずれているかを表す指標です。

目的変数とその平均値の差の二乗和を求めます。

$$SS_y = \sum_{i=1}^{n}{(y_i - \bar{y})2}$$

# 目的変数とその平均値
y <- mtcars$mpg # 目的変数
bar_y <- mean(y) # 平均値

# 総変動
SS_y <- sum((y - bar_y)^2)

回帰変動 $SS_R$

回帰変動（$SS_R$）は、回帰モデルが目的変数の平均値からどれだけずれているかを表す指標です。

回帰モデルによる目的変数の予測値 $\hat{y}$ とその平均値 $\bar{y}$ の差の二乗和を求めます。

$$SS_R = \sum{i=1}^{n}{(\hat{y}i - \bar{y})2}$$

# モデルの予測値
hat_y <- fitted(model)
# hat_y <- predict(model) # こちらでも可

# 回帰変動
SS_R <- sum((hat_y - bar_y)^2)

残差変動 $SS_E$

残差変動（$SS_E$）は、回帰モデルで説明できない目的変数の残りの変動を表す指標です。

目的変数 $y$ と回帰モデルによる予測値 $\hat{y}$ の差の二乗和を求めます。また、残差変動は総変動から回帰変動を引いた値になっています。

$$SS_E = SS_y - SS_R = \sum{i=1}^{n}{(y_i - \hat{y}i)2}$$

# 残差変動
SS_E <- sum((y - hat_y)^2)
# SS_E <- sum(residuals(model)^2) # こちらでも可

目的変数の総変動 $SS_y$

目的変数の総変動は、目的変数 $y$ の値がその平均値 $\bar{y}$ からどれだけずれているかを表す指標です。総変動は、目的変数とその平均値の差の二乗和で求められます。

$$SS_y = \sum_{i=1}^{n}{(y_i - \bar{y})2}$$

# 目的変数とその平均値
y <- mtcars$mpg # 目的変数
bar_y <- mean(y) # 平均値

# 総変動
SS_y <- sum((y - bar_y)^2)

回帰変動の自由度 $df_R$、残差変動の自由度 $df_E$

回帰変動と残差変動の自由度は、それぞれの平方和を求める際に使ったデータの個数や制約の数によって決まります。

  • 全体の自由度はサンプル数 - 1（$df_y = n - 1$）です。
  • 回帰変動の自由度は回帰係数の個数から 1 を引いた値です。単回帰モデルは回帰係数は $\beta_0$ と $\beta_1$ の 2 つなので、$df_R = 2 - 1 = 1$ です。
  • 残差変動の自由度は観測値の個数から回帰係数の個数を引いた値で、全体の自由度から回帰変動の自由度を引いた値です。今回、観測値は 32 個、回帰係数は 2 つなので、$df_E = df_y - df_R = 32 - 2 = 30$ です。
# 自由度
df_y <- nrow(mtcars) - 1 # 全体の自由度
df_R <- length(coefficients(model)) - 1 # 回帰変動の自由度
df_E <- df_y - df_R # 残差変動の自由度

回帰変動の平均平方 $MS_R$、残差変動の平均平方 $MS_E$

回帰変動と残差変動の平均平方は、それぞれの変動がどれだけ大きいかを表す指標です。

それぞれの平方和を対応する自由度で割って求めます。

$$MS_R = \frac{SS_R}{df_R}$$

$$MS_E = \frac{SS_E}{df_E}$$

# 平均平方
MS_y <- SS_y / df_y # 総変動の平均平方
MS_R <- SS_R / df_R # 回帰変動の平均平方
MS_E <- SS_E / df_E # 残差変動の平均平方

F 統計量と p 値

F 統計量は、回帰モデルが目的変数をどれだけ説明できるかを表す指標です。回帰モデルの当てはまりが良ければ、F 統計量は大きくなります。

回帰変動と残差変動の平均平方の比を求めます。

$$F = \frac{MS_R}{MS_E}$$

p 値は、F 統計量が得られる確率を表す指標です。回帰モデルが目的変数に有意な影響を与えていれば、p 値が小さくなります。

なお、このときの帰無仮説は次になります。

  • 帰無仮説 ( $H_0$ ): $\beta_1 = 0$（説明変数は目的変数に影響しない）
  • 対立仮説 ( $H_1$ ): $\beta_1 \neq 0$（説明変数は目的変数に影響する）
# F 統計量と p 値
F_R <- MS_R / MS_E # F統計量
p_R <- pf(F, df_R, df_E, lower.tail = FALSE) # p 値

分散分析表

回帰変動と誤差変動の平方和（SS）、自由度（DF）、平均平方（MS）、F 統計量（F_R）、および p 値（p_R）を分散分析表にまとて、単回帰モデルを評価します。

res <- data.frame(
  ID = c("Regression", "Residual", "Total"),
  DF = c(df_R, df_E, df_y),
  SS = c(SS_R, SS_E, SS_y),
  MS = c(MS_R, MS_E, MS_y),
  F = c(F_R, NA, NA),
  p = c(p_R, NA, NA)
)
cat("ANOVA table:")
res

ANOVA table:          ID DF        SS         MS        F            p
1 Regression  1  847.7252 847.725250 91.37533 1.293959e-10
2   Residual 30  278.3219   9.277398       NA           NA
3      Total 31 1126.0472  36.324103       NA           NA

F 統計量が 91.37533、p 値が 1.293959e-10 になりました。

p 値は十分に小さいので、この回帰モデルは単なる平均値よりも優れた予測力を持つとなります。

なお、これはanova(model)で同じ結果を取得できます。

anova(model)

Analysis of Variance Table

Response: mpg
          Df Sum Sq Mean Sq F value    Pr(>F)
wt         1 847.73  847.73  91.375 1.294e-10 ***
Residuals 30 278.32    9.28
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

回帰モデルの係数の分散分析と二元配置分散分析の違い

R

回帰モデルの係数の分散分析は、二元配置分散分析と計算方法が似ているので、分散分析の知識が役に立ちます。ですが、この 2 つは目的や仮定が違っています。注意が必要です。

目的 説明変数 説明変数間、因子間の関係 差の検定
回帰モデルの係数の分散分析 回帰モデルの精度、有意性を見る 連続値 相関、非線形性を考慮 係数がゼロと異なるかを t 検定
二元配置分散分析 因子が値に与える影響を見る カテゴリー 交互作用を考慮 因子水準間の平均値を F 検定

二元配置分散分析

カテゴリー変数の因子が値に与える影響を見るために、各因子水準間の平均値に差を F 検定で評価する。

平方和 自由度 平均平方 F 値 p 値
因子 A $SS_A = nb \sum{i = 1}^a (\bar{y}{i..} − \bar{y})2$ $a-1$ $MS_A = SS_A/(a-1)$ $F_A = MS_A/MS_E$ $p_A$
因子 B $SS_B = na \sum{j = 1}^b (\bar{y}{.j.} − \bar{y})2$ $b-1$ $MS_B = SS_B/(b-1)$ $F_B = MS_B/MS_E$ $p_B$
交互作用 $SS{AB} = n \sum{i = 1}^a \sum{j = 1}^b (\bar{y}{ij.} − \bar{y}{i..} − \bar{y}{.j.} + \bar{y})2$ $(a-1)(b-1)$ $MS{AB} = SS{AB}/*1$ $F{AB} = MS{AB}/MS_E$ $p_{AB}$
誤差 $SS_E = \sum{i = 1}^a \sum{j = 1}^b \sum{k = 1}^n ( y{ijk} − \bar{y}_{ij.})2$ $ab(n-1)$ $MS_E = SS_E/(ab(n-1))$ - -
総計 $SS_T = SS_A + SS_B + SS{AB} + SS_E = \sum{i = 1}^a \sum{j = 1}^b \sum{k = 1}^n ( y_{i j k} − \bar{y})2$ $abn-1$ - - -
  • $\bar{y}_{i..}$ :因子 A の水準 i の平均値
  • $\bar{y}_{.j.}$ :因子 B の水準 j の平均値
  • $\bar{y}_{ij.}$ :因子 A の水準 i と因子 B の水準 j の平均値
  • $\bar{y}$ :全体の平均値
  • $y_{ijk}$ :因子 A の水準 i、因子 B の水準 j、および実験単位 k の値
  • $a$ :因子 A の水準数
  • $b$ :因子 B の水準数
  • $n$ :各因子水準の組み合わせに対する実験単位の数

回帰モデルの係数の分散分析

回帰モデルの各説明変数が目的変数に与える影響、回帰モデルの精度や有意性を見るために、係数がゼロと異なるかどうかを t 検定で評価する。

変数 平方和 自由度 平均平方 F 値 p 値
説明変数 R1 $SS{R1} = \sum{i = 1}^n (\hat{y}_1 − \bar{y})2$ $ 1 $ $MS{R1} = SS{R1}/1$ $F_1 = MS_{R1}/MS_E$ $p_1$
説明変数 R2 $SS{R2} = \sum{i = 1}^n (\hat{y}2 − \hat{y}1)2$ $ 1 $ $MS{R2} = SS{R2}/1$ $F_2 = MS_{R2}/MS_E$ $p_2$
誤差 $SS_E = \sum{i = 1}^n (y_i − \hat{y}2)2$ $n-2$ $MS_E = SS_E/(n-2)$ - -
総計 $SS_T = SS{R1} + SS{R2} + SS_E = \sum_{i = 1}^n (y_i − \bar{y})2$ $n-1$ - - -

ここで、

  • $\hat{y}_1$ :説明変数 R1 のみを使った回帰モデルで予測された値
  • $\hat{y}_2$ :説明変数 R1 と R2 を使った回帰モデルで予測された値
  • $\bar{y}$ :目的変数の平均値
  • $y_i$ :目的変数の実測値

*1:a-1)(b-1

平方和のまとめ・その１

R

2 つの要因による効果とその交互作用を検討する統計手法が二元配置分散分析でした。

この方法は、ある値に対して要因 A と要因 B が紐づいたデータあるとき、要因 A の主効果、要因 B の主効果、要因 A と要因 B の交互作用の 3 つの効果を評価します。

二元配置分散分析の平方和

この評価に偏差平方和を使います。偏差平方和は集団の平均値からの乖離の程度を示しており、この差の大小でバラツキの程度が分かります。

分散分析は、ある値を要因と誤差で説明できると考えており、ある値のバラツキ（総変動, $SS_T$）は要素のバラツキ（要因の変動, $SS_A, SS_B$）と誤差のバラツキ（残差変動, $SS_E$）の合計と考えます。更に二元配置分散分析は、要因が組み合わさった時に生じる変動（交互作用の変動, $SS_{AB}$）も加わります。

$SS_T = SS_A + SS_B + SS_{AB} + SS_E$

これをまとめると次のようになります。

平方和の種類 名称 乖離の程度 数式
全体の平方和 総変動 各値と全体平均値 $$SS_T = \sum{i=1}^a \sum{j=1}^b \sum{k=1}^{n{ij}} (x_{ijk} - \bar{x})2$$
要因 A の平方和 要因 A の変動 要因 A の水準ごとの平均値と全体平均値 $$SS_A = b \sum{i=1}^a (\bar{x}{i\cdot\cdot} - \bar{x})2$$
要因 B の平方和 要因 B の変動 要因 B の水準ごとの平均値と全体平均値 $$SS_B = a \sum{j=1}^b (\bar{x}{\cdot j\cdot} - \bar{x})2$$
交互作用の平方和 交互作用の変動 要因 A と要因 B の水準ごとの組み合わせの平均値と要因 A と要因 B の水準ごとの平均値 $$SS{AB} = n \sum{i=1}^a \sum{j=1}^b (\bar{x}{ij\cdot} - \bar{x}{i\cdot\cdot} - \bar{x}{\cdot j\cdot} + \bar{x})2$$
残差（誤差）の平方和 残差（誤差）変動 各値と水準組み合わせごとの平均値 $$SS_E = \sum{i=1}^a \sum{j=1}^b \sum{k=1}^{n{ij}} (x{ijk} - \bar{x}{ij\cdot})2$$

不揃いデータの場合

ただし、データの偏りや欠損について注意が必要です。

二元配置分散分析に用いるデータは、すべての要因と水準でデータ数が等しい「等例数」、かつ欠損がない「完全配置」が理想的です。しかし、現実には実現できないこともあり、結果に影響を与える可能性が出てきます。

要因・水準の例数が違う場合（不等例数）

種類 方法 利点 欠点
タッカー法 水準の組み合わせ別の平均値を使う 各水準の組み合わせが 1 つのデータになるので、等例数と同じように扱える 各水準組み合わせのばらつき（残差変動）を無視するので、F 比が過大になりやすい
ウェルチ法 水準の組み合わせ別の分散で重み付けした平方和を使う 各水準の組み合わせのばらつき（残差変動）を考慮するので、F 比が適切になる 自由度や F 比の計算が複雑になる
ハインツ法 水準の組み合わせ別の標準偏差で標準化した値を使う 各水準の組み合わせのばらつき（残差変動）を均一化するので、F 比が適切になる 標準化により単位や尺度を失うので、解釈が難しくなる

要因・水準が欠損値を持つ場合（欠損配置）

種類 方法 利点 欠点
除外 欠損値が含まれる水準組み合わせを分析から除外する 欠損値の影響を受けないので、F 比が適切になる データ数が減るので、分析の精度や検出力が低下する
補完 欠損値を何らかの方法で推定して補完する データ数が減らないので、分析の精度や検出力が低下しない 補完方法によって分析結果が変わる可能性があり、補完の妥当性や信頼性の評価が必要になる
混合モデル 固定効果とランダム効果を考慮した混合モデルを使う 欠損値のパターンやメカニズムに依存しないので、F 比が適切になる モデルの選択や推定に専門的な知識や技術が必要になる

このように、不均等数や欠損配置には慎重な対応が必要になります。

そして、この場合、平方和の計算方法も重要になります。これは、上にある総変動、要因 A、要因 B、交互作用、残差変動の平方和とは異なる平方和です。

この平方和については、別でまとめます。

回帰モデルの係数の分散分析のまとめ・その１

R

共分散分析で回帰分析の分散分析が出てきたので、少し整理します。

回帰分析を分散分析する目的

回帰分析の分散分析は、データを意思決定に利用するための方法です。

これを使えば、回帰モデルに用いた異なる要因（説明変数、因子など）が目的変数に与える影響を統計的に評価できます。

一般的な回帰モデルは以下の式で表されます。

$$ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_n x_n + \varepsilon $$

  • $y$ : 目的変数
  • $x_1, x_2, \ldots, x_n$ : 説明変数、要因
  • $\beta_0, \beta_1, \beta_2, \ldots, \beta_n$ : 回帰係数（各説明変数の影響度）
  • $\varepsilon$ : 誤差項

この回帰モデルにおいて、異なる説明変数が目的変数に与える影響を評価するために分散分析を利用します。そして、回帰モデルや説明変数により目的変数の総変動（ばらつき）が減少（説明）する量を見ることで、その影響が分かります。

これには、以下の 2 つのアプローチがあります。

  • 回帰モデル全体の有意性評価：回帰モデルが目的変数をどれだけ説明できるかを見るために、回帰モデルと単なる平均値（切片だけのモデル）を比較する。
  • 特定の説明変数や切片に対する検定：各説明変数や切片が目的変数にどれだけ影響を与えるかを見るために、その係数が 0 であるかどうかを検定する。

評価の流れ

評価は大きく 2 つに分かれます。

  • 回帰モデルに対して目的変数の総変動、説明変数の変動（回帰変動）、残差の変動（誤差変動）を使った分散分析
  • 回帰係数に対する t 検定で、回帰係数が 0 でないことを統計的に確認する

そして、モデルの適合度、説明力、および各説明変数の重要性を評価します。

評価対象
モデルの適合度 回帰モデルの予測の信頼性を確認する。適合していれば、正しく予測できる可能性が高い。
説明力 説明変数が目的変数（応答変数）を説明する程度を評価する。重要な変数の特定に使えて、現象やプロセスの理解が深まる。
変数の重要性 目的変数に最も影響する変数を選択する（特徴選択:feature selection）。不要な変数を除き、モデルの複雑性を下げて、モデルの解釈性が上げる。

モデルの適合度の評価

誤差変動に対する各説明変数の統計的有意差を F 検定して、観測データに対するモデルの適合度を評価します。

説明力の評価

各説明変数の係数（$\beta_1, \beta_2, \ldots, \beta_n$）に対して t 検定を行い、説明変数が目的変数に与える影響を評価します。

目的変数に対して有意な影響を持つ説明変数は、この p 値が小さい傾向があります。

変数の重要性の評価

各説明変数の係数の大きさ（絶対値）を評価します。

目的変数に与える影響が大きい説明変数は、係数が大きい傾向があります。

そのため、係数の絶対値が小さい説明変数の除外は、回帰モデルの解釈に有効です。

また、他の様々な情報を使って目的変数の要否を判断することも有効です。

分散分析表について

分散分析表は分散分析の結果の解釈に重要です。

回帰モデルの分散分析は、要因がカテゴリデータから連続値に代わるので、自由度の求め方が変わりますが、それ以外の要因に対しては変わりません。

モデル全体、説明変数、残差について、各要素の平方和（Sum of Squares, SS）、自由度（Degrees of Freedom, DF）、および平均平方（Mean Square, MS）を計算します。

項目 意味 記号 式
モデル全体 データ全体のばらつき $SS{y}$, $SS{T}$ $$SS{y} = \sum{i=1}^{n} (y_i - \bar{y})2$$
^ ^ $df{y}$, $df{T}$ $$df_{y} = n - 1$$
^ ^ $MS{y}$, $MS{T}$ $$MS{y} = \frac{SS{y}}{df_{y}}$$
説明変数 回帰モデルが説明できる目的変数の変動 $SS_{R}$ $$SS{R} = \sum{i=1}^{n} (\hat{y}_i - \bar{y})2$$
^ ^ $df_{R}$ $$df_{R} = \text{説明変数の数}$$
^ ^ $MS_{R}$ $$MS{R} = \frac{SS{R}}{df_{R}}$$
残差 モデルが説明できなかった目的変数の変動 $SS_{E}$ $$SS{E} = SS{y} - SS{R} = \sum{i=1}^{n} (y_i - \hat{y}_i)2$$
^ ^ $df_{E}$ $$df{E} = df{y} - df_{R} = n - 1 - \text{説明変数の数}$$
^ ^ $MS_{E}$ $$MS{E} = \frac{SS{E}}{df_{E}}$$
  • $n$: サンプルサイズ
  • $Y_i$: 個々の観測値、個々のデータ
  • $\bar{Y}$: 目的変数の平均
  • $\hat{Y}_i$: 回帰モデルの予測値

回帰係数の検定

分散分析で有意になれば、個々の回帰係数（$\beta_i$）を t 検定して、目的変数に与える影響を確認します。

ここでの帰無仮説と対立仮説は、係数が 0 かどうかです。

  • $H_0: \beta_j = 0$
  • $H_1: \beta_j \neq 0$

そして、t 値は次の式になります。

$$ t_j = \frac{\hat{\beta}j}{\hat{\sigma}{\hat{\beta}_j}} $$

このとき、$\hat{\beta}j$ は $x_j$ の係数の推定値で、$\hat{\sigma}{\hat{\beta}j}$ は $\hat{\beta}j$ の標準誤差です。

そして、t 値は自由度 $n - (\text{説明変数の数}) - 1$ の t 分布に従うので、t 値を棄却域と比較すれば、回帰係数の目的変数に対する影響が分かります。

つまり、仮に単回帰分析で有意水準を $\alpha = 0.05$ にするなら、棄却域は $t_j < -t(0.025, n - 2)$ または $t_j > t(0.025, n - 2)$ です。

このt検定の部分は、別でまとめようと思います。

共分散分析のメモ・その１

R

交互作用が見つかる場合、二元配置分散分析をして、ad-hoc 解析をしてきました。そのため、共分散分析（ANCOVA）の必要性は高くなく、使う機会に恵まれなかった。

とはいえ、使わないから知らなくて良い、ってこともないので、まとめることにしました。

共分散分析(ANCOVA: Analysis of co-variance)とは

共分散分析（ANCOVA）は、分散分析（ANOVA: Analysis of Variance）と回帰分析を組み合わせた統計解析手法です。分散分析は、要因の水準が観測値の平均値に与える影響を検定するための方法で、要因が 1 つなら一元配置分散分析、2 つなら二元配置分散分析でした。

共分散分析は、回帰分析の特性を取り入れつつ、比較対象の要因（独立変数）以外の因子（共変量、Covariate）の影響を観測値（従属変数）から取り除き、要因の水準が測定値の平均値の差にどの程度影響を与えるのかを検討する方法です。なお、共分散分析は、回帰分析の性質を持つため、共変量が従属変数と線形の関係にあり、回帰直線が平行である場合に適しています。

共分散分析のモデル式： \displaystyle{ Y{ij} = \beta _ 0 + \beta _ 1 X{ij} + \beta _ 2 C _ i + \varepsilon _ {ij} }

  • $Y_{ij}$: i 番目の条件（グループ）と j 番目の観測値（従属変数）
  • $X_{ij}$: i 番目の条件と j 番目の観測における要因の値（独立変数）
  • $C_i$: i 番目の条件に関連する共変量の値
  • $\beta_0$: 切片、$\beta_1$: 独立変数の係数、$\beta_2$: 共変量の係数
  • $\varepsilon_{ij}$: 誤差項（平均 0 の正規分布で互いに独立）

分散分析、重回帰分析との違い

共分散分析は分散分析と回帰分析を組み合わせた手法なので、数式の形も重回帰分析と類似しています。

方法 調査対象 説明変数 目的変数 式の形 特徴
共分散分析 群間の平均値の差 名義尺度、計量尺度 計量尺度 $y=a_1x_1+a_2x_2+...+b+e$ 共変量が平均値に与える効果を調整できる
分散分析 群間の平均値の差 名義尺度 計量尺度 $y=a_1x_1+a_2x_2+...+b+e$ 3 群以上でも比較できる
重回帰分析 説明変数と目的変数の関係 計量尺度 計量尺度 $y=a_1x_1+a_2x_2+...+b+e$ 回帰係数や決定係数などで関係性を評価できる

また、共分散分析を行う際に注意が必要な確認事項もあります。

対象 内容 対応
独立変数と共変量 相関、多重共線性 強ければ、共変量を除外する
独立変数と共変量 交互作用 有意ならば、交互作用項を除外する
従属変数 正規性、等分散性 保証できなければ、観測値を変換する

計算の流れ

ANCOVA の計算手順は一般的に以下の通りです：

  1. データの収集および準備：従属変数、独立変数、共変量のデータを収集、整理する。
  2. モデルの設定：数式を用いて ANCOVA のモデルを設定する。
  3. パラメータ推定：最小二乗法などを用いて、モデルのパラメータ（係数）を推定する。
  4. 仮説検定：推定されたパラメータを用いて、独立変数の水準が従属変数に与える影響を検定する。
  5. 共変量の調整：共変量の効果を考慮し、要因の水準間を比較する。
  6. 結果の解釈：統計的有意性の評価結果を解釈する。