ETC3580_tutorial

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(visreg)
library(broom)
library(faraway)

fit2 = lm(Fertility ~ (.)^2, data = swiss) %>% 
  step(trace = FALSE)
aug = augment(fit2)
aug %>% 
  ggplot(aes(x = .fitted, y = .resid)) +
  geom_point() + 
  geom_smooth()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

aug %>% 
  ggplot(aes(x = Agriculture, y = .resid)) + 
  geom_point() + 
  geom_smooth()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

aug %>% 
  ggplot(aes(x = Education, y = .resid)) + 
  geom_point() + 
  geom_smooth()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

このコードは、R言語を使って線形回帰モデルの構築、モデル選択、そして診断プロットの作成を行っています。各ステップを詳しく解説すると、以下のようになります。

1. モデルの構築と変数選択

モデル式の指定
lm(Fertility ~ (.)^2, data = swiss)
ここでは、スイスのデータセット（swiss）を用いて、目的変数 Fertility（出生率や不妊指数など）を、残りのすべての説明変数とその2次の交互作用項（(.)^2）を使って説明する線形回帰モデルを作成しています。
ステップワイズ法による変数選択
%>% step(trace = FALSE)
作成した大規模なモデルに対して、step() 関数を適用し、AIC（赤池情報量基準）などを基に不要な項を削除して、よりシンプルで適切なモデルへと絞り込んでいます。trace = FALSEは、選択過程の詳細な出力を抑制するためのオプションです。

2. モデルの診断データの生成

augment() 関数の利用
aug = augment(fit2)
ここでは、broom パッケージの augment() 関数を使って、モデルオブジェクト fit2 に対して、各観測値に対応するフィッティング値（.fitted）や残差（.resid）などの診断情報を元のデータに追加しています。これにより、各データ点の予測値と誤差を簡単にプロットできるようになります。

3. 残差のプロットによるモデル診断

フィッティング値 vs. 残差のプロット
```
aug %>% 
  ggplot(aes(x = .fitted, y = .resid)) +
  geom_point() + 
  geom_smooth()
```
このプロットでは、モデルのフィッティング値（予測値）と残差の関係を視覚化しています。散布図で各点をプロットし、geom_smooth() により平滑化曲線を追加することで、残差にパターンやシステマティックな偏りがないかを確認します。もしパターンが見られると、モデルに非線形性や他の問題（例：異分散性）が存在する可能性があります。
Agriculture変数 vs. 残差のプロット
```
aug %>% 
  ggplot(aes(x = Agriculture, y = .resid)) + 
  geom_point() + 
  geom_smooth()
```
このプロットでは、説明変数の一つである Agriculture と残差の関係をプロットしています。目的は、特定の説明変数がモデルの残差にどのような影響を与えているかをチェックすることです。ここでも平滑化曲線を用いることで、非線形な関係やモデルが捉え切れていない変動がないかを視覚的に評価できます。

まとめ

全体として、このコードは以下のプロセスを実行しています：

全ての1次項と2次交互作用を含むモデルの構築し、ステップワイズ法を使って最適なモデルに絞り込む。
augment() 関数を用いて、各観測値のフィッティング値や残差などの診断情報を追加。
モデル診断のためのプロット（フィッティング値 vs. 残差、特定の変数 vs. 残差）を作成し、モデルの適合性や潜在的な問題（非線形性、異分散性など）を視覚的に確認する。

このような手法は、線形回帰モデルの評価や改善のための一般的なプロセスであり、モデルがデータにどれだけうまく適合しているか、またはどの部分に問題があるかを見極めるのに役立ちます。

USER

以下の疑問点に関してさらに詳しく解説してくれ。

フィッティング値ってなんだ？なんでそれと残差をプロットしているの？
この図をどう見たらモデルの適合性や潜在的な問題が確認できるのか？

ASSISTANT

以下、それぞれの疑問点について詳細に解説します。

フィッティング値とは？

定義
フィッティング値（または予測値）は、回帰モデルによって各観測点に対して計算された値です。
例えば、線形回帰モデルの場合、各観測点の説明変数の値を用いて、回帰係数をかけ合わせた結果がフィッティング値となります。
意味合い
モデルがデータにどれだけうまく当てはまっているか、またはどの程度の値を予測しているかを示す指標です。
実際の観測値との差（残差）が小さいほど、モデルの予測が実データに近いことを意味します。

なぜフィッティング値と残差をプロットするのか？

目的
フィッティング値と残差をプロットすることで、モデルの適合性や仮定の妥当性を評価できます。
このプロットは「残差プロット」と呼ばれ、回帰分析における重要な診断ツールです。
チェックする主な点
- ランダム性
  理想的な線形モデルでは、残差はフィッティング値に対してランダムに散らばっているはずです。
  特定のパターン（例えば、曲線状の傾向）が見られる場合、モデルに捉えきれていない非線形性が存在する可能性があります。
- 平均ゼロ
  残差の平均は0であるべきです。これにより、モデルがバイアスなくデータを捉えていると判断できます。
- 等分散性（ホモセダスティシティ）
  残差の分散がフィッティング値の大小にかかわらず一定であるかどうかを確認します。
  もし残差の散らばりが、フィッティング値の大小によって大きく変わる（例えば、広がっていく）場合は、異分散性の問題がある可能性があります。

残差プロットの見方とモデル適合性の評価方法

パターンが見られない場合（理想的な状態）
- 点が無作為に散らばっていれば、モデルはデータの変動をうまく捉えていると考えられます。
- 残差が平均0で、等分散であるならば、線形回帰モデルの基本的な前提（線形性、独立性、等分散性）が満たされていると判断できます。
パターンがある場合（問題の兆候）
- 曲線状のパターン
  - もし平滑化曲線が曲線を描いている場合、モデルがデータの非線形な関係をうまく捉えられていない可能性があります。
    → この場合、非線形項や別のモデル（例えば、ポリノミアル回帰など）の導入が検討されます。
- ファンネル（漏斗）状のパターン
  - フィッティング値が大きくなるにつれて残差の散らばりが広がる（または狭まる）場合、異分散性（ホモセダスティシティの違反）が疑われます。
    → その場合、重み付け回帰や誤差分散の調整が必要かもしれません。
- 外れ値や影響の大きい点
  - 一部の点が他と大きく離れて散布している場合、外れ値や影響の大きい観測値が存在する可能性があります。
    → これらはモデルの推定結果に大きな影響を与えるため、詳細な検討が必要です。
Agriculture変数とのプロットの場合
- 特定の説明変数（この場合はAgriculture）と残差の関係を見ることで、その変数がモデルの予測誤差にどのような影響を及ぼしているかがわかります。
- もしAgricultureの値と残差の間に特定のパターン（例えば、増加・減少の傾向）が見られる場合、その変数と目的変数の関係が単純な線形ではない可能性があるため、モデルの再検討が必要かもしれません。

まとめ

フィッティング値はモデルが各観測値に対して予測する値であり、実際の観測値との差分（残差）によってモデルの精度を評価できます。
残差プロットを作成することで、モデルの前提条件（線形性、等分散性、独立性）が満たされているかを確認し、非線形性や異分散性、外れ値などの問題点を視覚的に検出することができます。

このように、フィッティング値と残差の関係を見ることで、モデルの信頼性や問題点を評価し、必要に応じてモデルの改善に役立てることができるのです。

aug %>% 
  ggplot(aes(sample = .resid)) + 
  geom_qq() + 
  geom_qq_line()

このコードは、QQプロット（Quantile-Quantile Plot）を作成しており、モデルの残差（.resid）の分布が理論上の正規分布にどれだけ近いかを視覚的に評価するためのものです。

コードの詳細

ggplot(aes(sample = .resid))
- ここでは、aug（augment関数で得られた診断情報を含むデータセット）から残差（.resid）をサンプルとして指定しています。
geom_qq()
- この関数は、指定されたサンプル（ここでは残差）の分位点と、理論上の正規分布の分位点をプロットし、QQプロットを作成します。
geom_qq_line()
- 理論上の正規分布に基づいた参考線（45度線に相当する）をプロットに追加します。この線は、データが正規分布に従っていれば、点がこの線上に並ぶことを期待するための目安となります。

QQプロットの目的と解釈

正規性の検定
- 回帰モデルでは、残差が正規分布に従うという仮定があります。QQプロットはこの仮定が成り立っているかどうかを視覚的にチェックするためのツールです。
見方
- 点が参考線に沿って並んでいる場合：
  残差が理論上の正規分布に近いと判断でき、モデルの正規性の仮定が妥当であるとみなせます。
- 点が大きく参考線から逸脱している場合：
  残差が正規分布から逸脱している可能性があり、例えば、裾の重い分布（尖度が高い）、非対称性（歪度がある）などが疑われます。この場合、推定結果の信頼性や統計検定の結果に影響が出る可能性があります。

まとめ

目的：
QQプロットは、残差が正規分布に従っているかどうかを確認するためのグラフです。
評価方法：
点がgeom_qq_line() で追加された参考線上にほぼ沿っていれば、残差の正規性が担保されていると判断できます。一方で、体系的な逸脱が見られる場合は、モデルの前提条件に問題がある可能性があり、さらなる検討や改善（例えば、変数変換や別の分布を仮定したモデル）が必要となることが示唆されます。このプロットは、モデルの適合性や正規性の仮定の確認に非常に有用な診断ツールです。

aug %>% 
  ggplot(aes(x = .hat)) + 
  geom_dotplot()

## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

aug %>% 
  ggplot(aes( x = .cooksd, y = .hat)) +
  geom_point()

aug %>% 
  filter(.cooksd > 0.6) %>% 
  glimpse()

## Rows: 1
## Columns: 13
## $ .rownames        <chr> "V. De Geneve"
## $ Fertility        <dbl> 35
## $ Agriculture      <dbl> 1.2
## $ Examination      <int> 37
## $ Education        <int> 53
## $ Catholic         <dbl> 42.34
## $ Infant.Mortality <dbl> 18
## $ .fitted          <dbl> 32.60391
## $ .resid           <dbl> 2.396086
## $ .hat             <dbl> 0.8741296
## $ .sigma           <dbl> 6.208944
## $ .cooksd          <dbl> 0.6811614
## $ .std.resid       <dbl> 1.084899

このコードは、augment() で得られた診断情報を含むデータセット（aug）から、各観測値の Cook’s Distance (.cooksd) と レバレッジ (.hat) の値をプロットして、どのデータ点が回帰モデルに大きな影響を与えているかを視覚的に確認するためのものです。

コードの解説

ggplot(aes( x = .cooksd, y = .hat)) + geom_point()
- x軸：.cooksd
  - 各観測点の Cook’s Distance をプロットしています。これは、その点が回帰係数にどれだけ影響を与えるかの指標です。
- y軸：.hat
  - 各観測点のレバレッジ（hat値）をプロットしています。これは、その点の説明変数の値が、他のデータ点と比べてどれだけ「極端」かを示す指標です。
- geom_point()
  - 各観測点を散布図としてプロットしています。

Cook’s Distance (.cooksd) とは？

定義
Cook’s Distance は、ある観測点を除外したときにモデル全体の推定結果（回帰係数）がどれだけ変化するかを測る指標です。
意味
- 値が大きいほど、その観測点はモデルの結果に強い影響を与えていることを意味します。
- 高い Cook’s Distance を持つ点は、外れ値や影響の大きい点である可能性があるため、詳細に検討する必要があります。

Hat 値 (.hat) とは？

定義
Hat 値は、回帰分析における「レバレッジ」と呼ばれる指標で、各観測点が回帰モデルのフィッティングにどれだけ影響を与えているかを示す値です。
- 数学的には、ハット行列（H）の対角成分として計算されます。
意味
- 観測点の説明変数の値が全体の平均からどれだけ離れているかが反映されます。
- 高い hat 値を持つ観測点は、他の点と比べて説明変数の値が「極端」であるため、フィッティングに大きな影響を及ぼす可能性があります。

このプロットの解釈

影響の大きい観測点の発見
- 散布図上で、Cook’s Distance とレバレッジの両方が高い観測点は、モデルの推定に対して強い影響を持つ可能性があるため、特に注意が必要です。
診断目的
- モデル全体の信頼性を検証するため、影響の大きいデータ点（外れ値やレバレッジの高い点）が存在するかを確認します。
- もしこれらの点が存在する場合、モデルの再評価や、必要に応じた対処（例えば、外れ値の除去、ロバスト回帰の検討など）が求められることがあります。

このように、Cook’s Distance とレバレッジの関係をプロットすることで、どの観測点がモデルのフィッティングに大きな影響を与えているかを視覚的に把握でき、モデル診断や改善のための重要な情報を得ることができます。

ETC3580_tutorial_2

2025-03-17