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

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

USER

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

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

ASSISTANT

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


フィッティング値とは?

  • 定義
    フィッティング値(または予測値)は、回帰モデルによって各観測点に対して計算された値です。
    例えば、線形回帰モデルの場合、各観測点の説明変数の値を用いて、回帰係数をかけ合わせた結果がフィッティング値となります。

  • 意味合い
    モデルがデータにどれだけうまく当てはまっているか、またはどの程度の値を予測しているかを示す指標です。
    実際の観測値との差(残差)が小さいほど、モデルの予測が実データに近いことを意味します。


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

  • 目的
    フィッティング値と残差をプロットすることで、モデルの適合性仮定の妥当性を評価できます。
    このプロットは「残差プロット」と呼ばれ、回帰分析における重要な診断ツールです。

  • チェックする主な点

    • ランダム性
      理想的な線形モデルでは、残差はフィッティング値に対してランダムに散らばっているはずです。
      特定のパターン(例えば、曲線状の傾向)が見られる場合、モデルに捉えきれていない非線形性が存在する可能性があります。

    • 平均ゼロ
      残差の平均は0であるべきです。これにより、モデルがバイアスなくデータを捉えていると判断できます。

    • 等分散性(ホモセダスティシティ)
      残差の分散がフィッティング値の大小にかかわらず一定であるかどうかを確認します。
      もし残差の散らばりが、フィッティング値の大小によって大きく変わる(例えば、広がっていく)場合は、異分散性の問題がある可能性があります。


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

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

まとめ

  1. フィッティング値はモデルが各観測値に対して予測する値であり、実際の観測値との差分(残差)によってモデルの精度を評価できます。
  2. 残差プロットを作成することで、モデルの前提条件(線形性、等分散性、独立性)が満たされているかを確認し、非線形性や異分散性、外れ値などの問題点を視覚的に検出することができます。
このように、フィッティング値と残差の関係を見ることで、モデルの信頼性や問題点を評価し、必要に応じてモデルの改善に役立てることができるのです。
aug %>% 
  ggplot(aes(sample = .resid)) + 
  geom_qq() + 
  geom_qq_line()

何をしているのか?

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

コードの詳細

QQプロットの目的と解釈

まとめ

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) の値をプロットして、どのデータ点が回帰モデルに大きな影響を与えているかを視覚的に確認するためのものです。


コードの解説


Cook’s Distance (.cooksd) とは?


Hat 値 (.hat) とは?


このプロットの解釈

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