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言語を使って線形回帰モデルの構築、モデル選択、そして診断プロットの作成を行っています。各ステップを詳しく解説すると、以下のようになります。
lm(Fertility ~ (.)^2, data = swiss)
swiss
)を用いて、目的変数
Fertility
(出生率や不妊指数など)を、残りのすべての説明変数とその2次の交互作用項((.)^2
)を使って説明する線形回帰モデルを作成しています。%>% step(trace = FALSE)
step()
関数を適用し、AIC(赤池情報量基準)などを基に不要な項を削除して、よりシンプルで適切なモデルへと絞り込んでいます。trace = FALSE
は、選択過程の詳細な出力を抑制するためのオプションです。aug = augment(fit2)
broom
パッケージの augment()
関数を使って、モデルオブジェクト fit2
に対して、各観測値に対応するフィッティング値(.fitted
)や残差(.resid
)などの診断情報を元のデータに追加しています。これにより、各データ点の予測値と誤差を簡単にプロットできるようになります。フィッティング値 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
と残差の関係をプロットしています。目的は、特定の説明変数がモデルの残差にどのような影響を与えているかをチェックすることです。ここでも平滑化曲線を用いることで、非線形な関係やモデルが捉え切れていない変動がないかを視覚的に評価できます。
全体として、このコードは以下のプロセスを実行しています:
augment()
関数を用いて、各観測値のフィッティング値や残差などの診断情報を追加。このような手法は、線形回帰モデルの評価や改善のための一般的なプロセスであり、モデルがデータにどれだけうまく適合しているか、またはどの部分に問題があるかを見極めるのに役立ちます。
以下の疑問点に関してさらに詳しく解説してくれ。
以下、それぞれの疑問点について詳細に解説します。
定義
フィッティング値(または予測値)は、回帰モデルによって各観測点に対して計算された値です。
例えば、線形回帰モデルの場合、各観測点の説明変数の値を用いて、回帰係数をかけ合わせた結果がフィッティング値となります。
意味合い
モデルがデータにどれだけうまく当てはまっているか、またはどの程度の値を予測しているかを示す指標です。
実際の観測値との差(残差)が小さいほど、モデルの予測が実データに近いことを意味します。
目的
フィッティング値と残差をプロットすることで、モデルの適合性や仮定の妥当性を評価できます。
このプロットは「残差プロット」と呼ばれ、回帰分析における重要な診断ツールです。
チェックする主な点
ランダム性
理想的な線形モデルでは、残差はフィッティング値に対してランダムに散らばっているはずです。
特定のパターン(例えば、曲線状の傾向)が見られる場合、モデルに捉えきれていない非線形性が存在する可能性があります。
平均ゼロ
残差の平均は0であるべきです。これにより、モデルがバイアスなくデータを捉えていると判断できます。
等分散性(ホモセダスティシティ)
残差の分散がフィッティング値の大小にかかわらず一定であるかどうかを確認します。
もし残差の散らばりが、フィッティング値の大小によって大きく変わる(例えば、広がっていく)場合は、異分散性の問題がある可能性があります。
aug %>%
ggplot(aes(sample = .resid)) +
geom_qq() +
geom_qq_line()
このコードは、QQプロット(Quantile-Quantile Plot)を作成しており、モデルの残差(.resid)の分布が理論上の正規分布にどれだけ近いかを視覚的に評価するためのものです。
ggplot(aes(sample = .resid))
geom_qq()
geom_qq_line()
点が参考線に沿って並んでいる場合:
残差が理論上の正規分布に近いと判断でき、モデルの正規性の仮定が妥当であるとみなせます。
点が大きく参考線から逸脱している場合:
残差が正規分布から逸脱している可能性があり、例えば、裾の重い分布(尖度が高い)、非対称性(歪度がある)などが疑われます。この場合、推定結果の信頼性や統計検定の結果に影響が出る可能性があります。
目的:
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()