第2〜4章 記述統計・t検定・回帰分析・分散分析

データセット：College（アメリカの大学データ）／目的変数：卒業率（Grad.Rate）

---

0. この資料の位置づけと分析の流れ

2年生では記述統計（ヒストグラム・代表値・散布度・箱ひげ図）を学んだ。 3年生ではそれを土台に、「グループ間の差は偶然か？」「2変数間に関係はあるか？」 という問いに統計的に答える手法を学ぶ。

記述統計（復習）
  ↓
質 × 量 の可視化：平行箱ひげ図（Private 別の Grad.Rate）
  ↓
量 × 量 の可視化：散布図（Outstate × Grad.Rate）
  ↓
t 検定：私立 vs 公立の卒業率に差があるか？（2群の平均値の比較）
  ↓
単回帰分析：授業料は卒業率を説明できるか？（量 × 量の関係を数式で）
  ↓
ANOVA（分散分析）：回帰モデルは全体として意味があるか？

今回使う変数：

変数名	型	説明
Private	質的（2値）	私立か否か（Yes / No）
Grad.Rate	量的（連続）	卒業率（%）← 目的変数
Outstate	量的（連続）	州外学生の授業料（ドル）← 説明変数
Expend	量的（連続）	学生一人当たり支出（ドル）← 説明変数（参考）

---

1. 準備：パッケージとデータの読み込み

# ---------------------------------------------------------------
# 使用パッケージ
#   tidyverse : データ操作（dplyr）と可視化（ggplot2）をまとめて読み込む
#   cowplot   : 複数の ggplot2 グラフを並べて表示するパッケージ
# ---------------------------------------------------------------
library(tidyverse)
library(cowplot)

# ---------------------------------------------------------------
# データの読み込み
# ---------------------------------------------------------------
college <- read.csv("ISLR_CSV_Data/College.csv")

# ---------------------------------------------------------------
# 【分析対象変数の短縮代入】
# 毎回 college$変数名 と書くのは冗長なので、
# よく使う変数をあらかじめ短い名前に代入しておく。
#   y : 目的変数  = 卒業率（Grad.Rate）
#   x : 説明変数  = 州外学生の授業料（Outstate）
# ---------------------------------------------------------------
y <- college$Grad.Rate   # 目的変数：卒業率（%）
x <- college$Outstate    # 説明変数：州外学生の授業料（ドル）

# 読み込み確認：先頭5行と変数の型を確認する
head(college, 5)

str(college)

## 'data.frame':    777 obs. of  19 variables:
##  $ X          : chr  "Abilene Christian University" "Adelphi University" "Adrian College" "Agnes Scott College" ...
##  $ Private    : chr  "Yes" "Yes" "Yes" "Yes" ...
##  $ Apps       : int  1660 2186 1428 417 193 587 353 1899 1038 582 ...
##  $ Accept     : int  1232 1924 1097 349 146 479 340 1720 839 498 ...
##  $ Enroll     : int  721 512 336 137 55 158 103 489 227 172 ...
##  $ Top10perc  : int  23 16 22 60 16 38 17 37 30 21 ...
##  $ Top25perc  : int  52 29 50 89 44 62 45 68 63 44 ...
##  $ F.Undergrad: int  2885 2683 1036 510 249 678 416 1594 973 799 ...
##  $ P.Undergrad: int  537 1227 99 63 869 41 230 32 306 78 ...
##  $ Outstate   : int  7440 12280 11250 12960 7560 13500 13290 13868 15595 10468 ...
##  $ Room.Board : int  3300 6450 3750 5450 4120 3335 5720 4826 4400 3380 ...
##  $ Books      : int  450 750 400 450 800 500 500 450 300 660 ...
##  $ Personal   : int  2200 1500 1165 875 1500 675 1500 850 500 1800 ...
##  $ PhD        : int  70 29 53 92 76 67 90 89 79 40 ...
##  $ Terminal   : int  78 30 66 97 72 73 93 100 84 41 ...
##  $ S.F.Ratio  : num  18.1 12.2 12.9 7.7 11.9 9.4 11.5 13.7 11.3 11.5 ...
##  $ perc.alumni: int  12 16 30 37 2 11 26 37 23 15 ...
##  $ Expend     : int  7041 10527 8735 19016 10922 9727 8861 11487 11644 8991 ...
##  $ Grad.Rate  : int  60 56 54 59 15 55 63 73 80 52 ...

# ---------------------------------------------------------------
# 【重要】質的変数をファクター型に変換する
#
# Rでは質的変数（カテゴリ変数）を「factor（因子）型」として
# 明示的に指定する必要がある。
# 指定しないと文字列（character）型のまま扱われ、
# t検定・ANOVA・グループ別集計が正しく動かない場合がある。
#
# as.factor() で変換し、levels() で水準（グループ）を確認できる。
# ---------------------------------------------------------------

college$Private <- as.factor(college$Private)  # 文字列 → ファクター型へ変換

levels(college$Private)   # 水準の確認：「No」「Yes」の2グループ

## [1] "No"  "Yes"

table(college$Private)    # 各グループの件数を確認

## 
##  No Yes 
## 212 565

---

2. 記述統計（復習）

# ---------------------------------------------------------------
# 2年生の復習：目的変数 Grad.Rate の基本統計量をまとめて確認
# ---------------------------------------------------------------

summary(x)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2340    7320    9990   10441   12925   21700

# 主要な統計量を一覧表で表示
desc_table <- data.frame(
  統計量 = c("データ数", "平均値", "中央値", "標準偏差",
              "Q1（25%点）", "Q3（75%点）", "四分位範囲（IQR）",
              "最小値", "最大値"),
  値 = round(c(
    length(y),
    mean(y),
    median(y),
    sd(y),
    quantile(y, 0.25),
    quantile(y, 0.75),
    IQR(y),
    min(y),
    max(y)
  ), 2)
)

knitr::kable(desc_table, caption = "卒業率（Grad.Rate）の記述統計量",
             align = c("l","r"))

卒業率（Grad.Rate）の記述統計量

統計量	値
データ数	777.00
平均値	65.46
中央値	65.00
標準偏差	17.18
Q1（25%点）	53.00
Q3（75%点）	78.00
四分位範囲（IQR）	25.00
最小値	10.00
最大値	118.00

# 2年生の復習：ヒストグラムと箱ひげ図を並べて表示
p1 <- ggplot(college, aes(x = Grad.Rate)) +
  geom_histogram(binwidth = 5, fill = "steelblue", color = "white") +
  geom_vline(xintercept = mean(y),
             color = "red", linetype = "solid", linewidth = 1) +
  geom_vline(xintercept = median(y),
             color = "darkorange", linetype = "dashed", linewidth = 1) +
  labs(title = "ヒストグラム", x = "卒業率（%）", y = "度数") +
  theme_bw()

p2 <- ggplot(college, aes(x = Grad.Rate, y = "")) +
  geom_boxplot(fill = "steelblue", outlier.color = "red", outlier.size = 2) +
  labs(title = "箱ひげ図", x = "卒業率（%）", y = "") +
  theme_bw()

plot_grid(p1, p2, ncol = 1, rel_heights = c(2, 1))

---

3. 質 × 量の可視化：平行箱ひげ図

# ---------------------------------------------------------------
# 【平行箱ひげ図とは？】
# 質的変数（グループ）ごとに量的変数の箱ひげ図を並べて描いたもの。
# グループ間の分布の違い（中央値・ばらつき・外れ値）を
# 一目で比較できる。
#
# ここでは「私立（Yes）」と「公立（No）」で
# 卒業率（Grad.Rate）の分布を比較する。
# ---------------------------------------------------------------

ggplot(college, aes(x = Private, y = Grad.Rate, fill = Private)) +
  geom_boxplot(
    color         = "black",
    outlier.color = "red",
    outlier.size  = 2,
    width         = 0.5
  ) +
  # 実際のデータ点を半透明で重ねる（分布の形をより詳しく確認）
  geom_jitter(
    width = 0.15, alpha = 0.2, size = 1, color = "black"
  ) +
  scale_fill_manual(
    values = c("No" = "lightcoral", "Yes" = "steelblue"),
    labels = c("No" = "公立", "Yes" = "私立")
  ) +
  scale_x_discrete(labels = c("No" = "公立", "Yes" = "私立")) +
  labs(
    title    = "私立・公立別の卒業率の分布",
    subtitle = "箱ひげ図＋データ点（赤点は外れ値）",
    x        = "大学の種別",
    y        = "卒業率（%）",
    fill     = "種別"
  ) +
  theme_bw() +
  theme(legend.position = "none")  # 凡例は軸ラベルで十分なので非表示

# ---------------------------------------------------------------
# グループ別の記述統計量を数値でも確認する
#
# dplyr の group_by() + summarise() の組み合わせで
# グループごとに任意の統計量を計算できる。
#   group_by(Private)  : Private の水準ごとにグループ化
#   summarise()        : グループごとに統計量を集計
#   n()                : グループのデータ件数
# ---------------------------------------------------------------

college %>%
  group_by(Private) %>%
  summarise(
    件数       = n(),
    平均値     = round(mean(Grad.Rate), 2),
    中央値     = round(median(Grad.Rate), 2),
    標準偏差   = round(sd(Grad.Rate), 2),
    Q1         = round(quantile(Grad.Rate, 0.25), 2),
    Q3         = round(quantile(Grad.Rate, 0.75), 2),
    IQR        = round(IQR(Grad.Rate), 2)
  ) %>%
  knitr::kable(caption = "私立・公立別の卒業率の記述統計量")

私立・公立別の卒業率の記述統計量

Private	件数	平均値	中央値	標準偏差	Q1	Q3	IQR
No	212	56.04	55	14.58	46	65	19
Yes	565	69.00	69	16.75	58	81	23

視覚的な確認ポイント： - 私立と公立で、卒業率の中央値（箱の中線）はどちらが高いか？ - ばらつき（箱の高さ＝IQR）はどちらが大きいか？ - 外れ値（赤い点）はどちらのグループに多いか？

これらを次の t 検定で「統計的に有意な差か」確認する。

---

4. 量 × 量の可視化：散布図

# ---------------------------------------------------------------
# 【散布図とは？】
# 2つの量的変数の関係を点でプロットしたグラフ。
# 横軸（x）に説明変数、縦軸（y）に目的変数を置くのが慣例。
#
# ここでは「授業料（Outstate）」が「卒業率（Grad.Rate）」と
# どのような関係にあるかを確認する。
#
# geom_smooth(method = "lm") : 線形回帰直線を重ねて描く
#   method = "lm" → 線形モデル（Linear Model）による回帰直線
#   se = TRUE     → 95%信頼区間を帯状に表示（デフォルト TRUE）
# ---------------------------------------------------------------

ggplot(college, aes(x = x, y = y)) +
  geom_point(
    alpha = 0.4,        # 点の透明度（重なりを見やすくする）
    size  = 1.5,
    color = "steelblue"
  ) +
  geom_smooth(
    method  = "lm",     # 線形回帰直線を追加
    color   = "red",
    linewidth = 1,
    se      = TRUE      # 95%信頼区間を帯で表示
  ) +
  labs(
    title    = "授業料と卒業率の関係",
    subtitle = "赤線：線形回帰直線　帯：95%信頼区間",
    x        = "州外学生の授業料（ドル）",
    y        = "卒業率（%）"
  ) +
  theme_bw()

# ---------------------------------------------------------------
# 私立・公立を色分けして散布図を描く
# グループによって関係のパターンが異なるかも確認できる。
# ---------------------------------------------------------------

ggplot(college, aes(x = x, y = y, color = Private)) +
  geom_point(alpha = 0.4, size = 1.5) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 1) +
  scale_color_manual(
    values = c("No" = "lightcoral", "Yes" = "steelblue"),
    labels = c("No" = "公立", "Yes" = "私立")
  ) +
  labs(
    title    = "授業料と卒業率の関係（私立・公立別）",
    subtitle = "各グループに個別の回帰直線を重ねて表示",
    x        = "州外学生の授業料（ドル）",
    y        = "卒業率（%）",
    color    = "種別"
  ) +
  theme_bw()

# ---------------------------------------------------------------
# 【相関係数の確認】
# cor() で2変数間のピアソン積率相関係数を計算できる。
#   +1 に近い → 強い正の相関（一方が増えると他方も増える）
#   -1 に近い → 強い負の相関（一方が増えると他方は減る）
#    0 に近い → 相関なし
# 目安：|r| > 0.7 → 強い相関　|r| > 0.4 → 中程度　|r| < 0.2 → 弱い
# ---------------------------------------------------------------

cor(x, y)   # 授業料と卒業率の相関係数

## [1] 0.5712899

# 相関係数の有意性検定（帰無仮説：母相関係数 = 0）
cor.test(y, x)

## 
##  Pearson's product-moment correlation
## 
## data:  y and x
## t = 19.377, df = 775, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5219282 0.6168381
## sample estimates:
##       cor 
## 0.5712899

確認ポイント： - 散布図の点の散らばりのパターンから、正・負・無相関のどれに見えるか？ - 相関係数の値はどのくらいか？「強い」「中程度」「弱い」のどれにあたるか？

---

5. t 検定：私立 vs 公立の卒業率に差はあるか？

5-1. t 検定の前提確認

# ---------------------------------------------------------------
# 【t 検定とは？】
# 2つのグループの母平均に差があるかどうかを検定する手法。
#
# 【帰無仮説と対立仮説】
#   H0（帰無仮説）：私立と公立の卒業率の母平均は等しい
#   H1（対立仮説）：私立と公立の卒業率の母平均は異なる（両側検定）
#
# 【Welch の t 検定 vs Student の t 検定】
# R の t.test() はデフォルトで Welch の t 検定を行う。
# Welch の t 検定は「2グループの分散が等しいと仮定しない」検定で、
# 分散が異なる場合にも適用できる（より一般的で安全）。
#
# 【等分散の確認：Levene 検定】
# 分散が等しいかどうかは var.test() で確認できる。
#   p値 < 0.05 → 等分散の仮説を棄却（→ Welch の t 検定を使う）
#   p値 ≥ 0.05 → 等分散を仮定してよい（→ Student の t 検定も可）
# ---------------------------------------------------------------

# 2グループのデータを取り出す
# filter() : 条件を満たす行だけ抽出
# pull()   : 指定した列をベクトルとして取り出す
grad_private <- college %>% filter(Private == "Yes") %>% pull(Grad.Rate)
grad_public  <- college %>% filter(Private == "No")  %>% pull(Grad.Rate)

# グループごとの件数と平均を再確認
cat("私立大学：", length(grad_private), "校　平均卒業率：",
    round(mean(grad_private), 2), "%\n")

## 私立大学： 565 校　平均卒業率： 69 %

cat("公立大学：", length(grad_public),  "校　平均卒業率：",
    round(mean(grad_public),  2), "%\n")

## 公立大学： 212 校　平均卒業率： 56.04 %

cat("平均の差（私立 - 公立）：",
    round(mean(grad_private) - mean(grad_public), 2), "ポイント\n")

## 平均の差（私立 - 公立）： 12.96 ポイント

# 等分散性の検定（F 検定）
# H0：2グループの母分散は等しい
var.test(grad_private, grad_public)

## 
##  F test to compare two variances
## 
## data:  grad_private and grad_public
## F = 1.3191, num df = 564, denom df = 211, p-value = 0.01852
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  1.048129 1.641039
## sample estimates:
## ratio of variances 
##           1.319117

5-2. Welch の t 検定（両側検定）

# ---------------------------------------------------------------
# 両側検定（two-sided）
# 「私立と公立で卒業率に差がある」かどうかを検定する。
# alternative = "two.sided" が両側検定（デフォルト）。
#
# 【結果の読み方】
#   t 値    : 検定統計量（絶対値が大きいほど差が顕著）
#   df      : 自由度（Welch 検定では実数になることがある）
#   p 値    : 帰無仮説が正しい場合にこの t 値以上が出る確率
#             → p < 0.05 で「有意水準5%で帰無仮説を棄却」
#   95% CI  : 平均の差の95%信頼区間
# ---------------------------------------------------------------

result_ttest <- t.test(
  grad_private,          # グループ1（私立）
  grad_public,           # グループ2（公立）
  alternative = "two.sided",  # 両側検定
  var.equal   = FALSE         # 等分散を仮定しない（Welch の t 検定）
)

result_ttest  # 結果を表示

## 
##  Welch Two Sample t-test
## 
## data:  grad_private and grad_public
## t = 10.579, df = 431.97, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  10.54880 15.36276
## sample estimates:
## mean of x mean of y 
##  68.99823  56.04245

5-3. 片側検定（一方向の差を検定する場合）

# ---------------------------------------------------------------
# 【片側検定とは？】
# 「私立の卒業率は公立より高い」という一方向の仮説を検定する場合。
#
#   alternative = "greater"  : グループ1 > グループ2 の方向で検定
#   alternative = "less"     : グループ1 < グループ2 の方向で検定
#
# 片側検定は両側検定の p値の半分になるため、「有意になりやすい」。
# ただし、方向の仮説は分析前に理論的根拠をもとに立てる必要がある。
# 結果を見てから片側にするのは統計的に不正（p値の操作）。
# ---------------------------------------------------------------

# 「私立の卒業率 > 公立の卒業率」の片側検定
t.test(
  grad_private,
  grad_public,
  alternative = "greater",   # 私立 > 公立 の方向
  var.equal   = FALSE
)

## 
##  Welch Two Sample t-test
## 
## data:  grad_private and grad_public
## t = 10.579, df = 431.97, p-value < 2.2e-16
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  10.93711      Inf
## sample estimates:
## mean of x mean of y 
##  68.99823  56.04245

5-4. t 検定結果の可視化

# ---------------------------------------------------------------
# t 検定の結果をグラフで視覚化する
# エラーバー（平均 ± 1.96×標準誤差 ≒ 95%信頼区間）を追加することで
# 「差が有意かどうか」を視覚的に示す。
# 2つのエラーバーが重なっていない → 有意差がある可能性が高い
# ---------------------------------------------------------------

# グループ別の平均・標準誤差を計算
group_stats <- college %>%
  group_by(Private) %>%
  summarise(
    mean_grad = mean(Grad.Rate),
    se_grad   = sd(Grad.Rate) / sqrt(n()),  # 標準誤差 = SD / √n
    n         = n()
  ) %>%
  mutate(
    ci_lower = mean_grad - 1.96 * se_grad,  # 95%CI 下限（近似）
    ci_upper = mean_grad + 1.96 * se_grad   # 95%CI 上限（近似）
  )

ggplot(group_stats, aes(x = Private, y = mean_grad, fill = Private)) +
  # 棒グラフ（平均値）
  geom_bar(stat = "identity", width = 0.5, color = "black") +
  # エラーバー（95%信頼区間）
  geom_errorbar(
    aes(ymin = ci_lower, ymax = ci_upper),
    width = 0.15, linewidth = 0.8
  ) +
  scale_fill_manual(values = c("No" = "lightcoral", "Yes" = "steelblue")) +
  scale_x_discrete(labels = c("No" = "公立", "Yes" = "私立")) +
  labs(
    title    = "私立・公立別の平均卒業率（エラーバー：95%信頼区間）",
    subtitle = paste0("Welch t 検定 p 値 = ",
                      round(result_ttest$p.value, 4),
                      ifelse(result_ttest$p.value < 0.05,
                             "　→ 有意差あり（p < 0.05）",
                             "　→ 有意差なし（p ≥ 0.05）")),
    x        = "大学の種別",
    y        = "平均卒業率（%）"
  ) +
  theme_bw() +
  theme(legend.position = "none") +
  # 有意差を示すブラケット（p < 0.05 の場合のみ追加）
  {if (result_ttest$p.value < 0.05)
    annotate("segment",
             x = 1, xend = 2,
             y = max(group_stats$ci_upper) + 2,
             yend = max(group_stats$ci_upper) + 2,
             linewidth = 0.8) } +
  {if (result_ttest$p.value < 0.05)
    annotate("text",
             x = 1.5,
             y = max(group_stats$ci_upper) + 3.5,
             label = "*", size = 7) }

t 検定の結果の読み方まとめ：

確認項目	値	解釈
t 値	10.579	正 → 私立の方が高い
自由度	432	Welch法による調整値
p 値	0	< 0.05 → 帰無仮説を棄却
平均の差の95%CI	[10.55, 15.36]	0を含まない → 有意差あり

---

6. 単回帰分析：授業料は卒業率を説明できるか？

6-1. 線形回帰モデルの構築

# ---------------------------------------------------------------
# 【単回帰分析とは？】
# 1つの説明変数（x）で目的変数（y）を予測する線形モデル。
# 数式：y = β0 + β1・x + ε
#   β0 : 切片（x=0 のときの y の予測値）
#   β1 : 回帰係数（x が1増えると y はどれだけ変わるか）
#   ε  : 誤差項（モデルで説明できない部分）
#
# 【lm() 関数の書き方】
# lm(目的変数 ~ 説明変数, data = データフレーム名)
# 「~」は「〜によって説明される」と読む。
#
# ここでは：Grad.Rate = β0 + β1・Outstate + ε
# 「卒業率は授業料によって説明される」という仮定のモデル。
# ---------------------------------------------------------------

# 線形回帰モデルを構築して「model1」に保存
model1 <- lm(Grad.Rate ~ Outstate, data = college)

# 結果の表示
summary(model1)

## 
## Call:
## lm(formula = Grad.Rate ~ Outstate, data = college)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -49.168  -9.114  -0.006   8.576  55.114 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.000e+01  1.408e+00   28.40   <2e-16 ***
## Outstate    2.439e-03  1.259e-04   19.38   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.11 on 775 degrees of freedom
## Multiple R-squared:  0.3264, Adjusted R-squared:  0.3255 
## F-statistic: 375.5 on 1 and 775 DF,  p-value: < 2.2e-16

6-2. 回帰係数の解釈

# ---------------------------------------------------------------
# 【summary() の出力の読み方】
#
# Coefficients の表：
#   Estimate    : 推定された回帰係数（β̂0, β̂1）
#   Std. Error  : 標準誤差（係数の推定精度）
#   t value     : t 統計量（Estimate / Std. Error）
#   Pr(>|t|)    : p 値（係数が0でないかどうかの検定）
#
# モデル全体の指標：
#   R-squared   : 決定係数 → モデルが y のばらつきを何%説明しているか
#   F-statistic : F 統計量 → 回帰モデル全体の有意性（後の ANOVA と同じ）
#
# 【回帰係数の解釈（例）】
# Outstate の Estimate が仮に 0.0015 なら：
# 「授業料が 1 ドル高い大学では、卒業率が平均 0.0015% 高い」
# ---------------------------------------------------------------

# 係数だけを取り出す
coef(model1)

##  (Intercept)     Outstate 
## 39.995114285  0.002439327

# 切片と回帰係数をわかりやすく表示
cat("切片（β̂0）：", round(coef(model1)[1], 4), "\n")

## 切片（β̂0）： 39.9951

cat("Outstate の回帰係数（β̂1）：", round(coef(model1)[2], 6), "\n")

## Outstate の回帰係数（β̂1）： 0.002439

cat("\n解釈：授業料が 1,000 ドル上がると、卒業率は平均",
    round(coef(model1)[2] * 1000, 3), "% 変化する\n")

## 
## 解釈：授業料が 1,000 ドル上がると、卒業率は平均 2.439 % 変化する

6-3. 回帰直線の可視化

# ---------------------------------------------------------------
# 回帰直線と残差（residuals）を可視化する
# 残差 = 実際の y - 予測した y
# 残差が小さいほどモデルのあてはまりが良い
# ---------------------------------------------------------------

# 予測値と残差をデータフレームに追加
college_aug <- college %>%
  mutate(
    fitted    = fitted(model1),    # 回帰モデルによる予測値
    residuals = residuals(model1)  # 残差（実測値 - 予測値）
  )

# 散布図＋回帰直線＋残差の可視化
ggplot(college_aug, aes(x = Outstate, y = Grad.Rate)) +
  # 残差を縦線で表示（実測値と予測値のずれ）
  geom_segment(
    aes(xend = Outstate, yend = fitted),
    alpha = 0.2, color = "gray50"
  ) +
  # データ点
  geom_point(alpha = 0.4, size = 1.5, color = "steelblue") +
  # 回帰直線
  geom_line(aes(y = fitted), color = "red", linewidth = 1) +
  labs(
    title    = "回帰直線と残差",
    subtitle = paste0("Grad.Rate = ",
                      round(coef(model1)[1], 2), " + ",
                      round(coef(model1)[2], 5), " × Outstate",
                      "　　R2 = ", round(summary(model1)$r.squared, 3)),
    x        = "州外学生の授業料（ドル）",
    y        = "卒業率（%）"
  ) +
  theme_bw()

6-4. 残差診断プロット

# ---------------------------------------------------------------
# 【残差診断とは？】
# 回帰分析の仮定（線形性・等分散性・正規性）が満たされているか
# 残差のグラフで確認する。これを「モデル診断」と呼ぶ。
#
# par(mfrow = c(2,2)) で4つのグラフを2×2に並べる。
# plot(model) で自動的に4種類の診断プロットが出る：
#   1. Residuals vs Fitted  : 残差と予測値の関係（線形性の確認）
#   2. Normal Q-Q           : 残差の正規性の確認
#   3. Scale-Location       : 残差の等分散性の確認
#   4. Residuals vs Leverage: 影響力の大きい観測値の確認
# ---------------------------------------------------------------

par(mfrow = c(2, 2))
plot(model1)

par(mfrow = c(1, 1))

---

7. ANOVA（分散分析）：回帰モデルは全体として意味があるか？

7-1. 回帰結果を使った ANOVA

# ---------------------------------------------------------------
# 【ANOVA（分散分析）とは？】
# 「回帰モデル全体が統計的に意味があるか」を F 統計量で検定する。
#
# 考え方：目的変数 y のばらつき（総変動）を2つに分解する
#   ・回帰変動（SS_R）: 回帰モデルが「説明できた」ばらつき
#   ・残差変動（SS_E）: モデルが「説明できなかった」ばらつき
#
# F 統計量 = （回帰変動 / 回帰の自由度）÷（残差変動 / 残差の自由度）
#          = 回帰の平均二乗 ÷ 残差の平均二乗
#
# F が大きい（＝回帰で説明できた部分が大きい）
#   → p 値が小さい → 「回帰モデルは有意」
#
# 【lm() → anova() の流れ】
# Rでは lm() でモデルを作った後、anova() に渡すだけで
# 分散分析表（ANOVA テーブル）が得られる。
# これは参考コードの disbud.trt → anova(disbud.trt) と同じ流れ。
# ---------------------------------------------------------------

# lm() で作ったモデルを anova() に渡す
anova_result <- anova(model1)
anova_result

# ---------------------------------------------------------------
# ANOVA テーブルの各列の意味：
#   Df         : 自由度
#   Sum Sq     : 変動（平方和）
#   Mean Sq    : 平均二乗（= Sum Sq / Df）
#   F value    : F 統計量
#   Pr(>F)     : p 値（帰無仮説：回帰係数はすべて0）
#
# ここでは説明変数が1つ（Outstate）なので、
# t 検定の t 値の2乗が F 値と一致することも確認できる。
# ---------------------------------------------------------------

# t 値の二乗と F 値の比較
t_val <- summary(model1)$coefficients["Outstate", "t value"]
f_val <- anova_result$`F value`[1]

cat("回帰係数の t 値：", round(t_val, 4), "\n")

## 回帰係数の t 値： 19.3775

cat("t 値の2乗（t²）：", round(t_val^2, 4), "\n")

## t 値の2乗（t²）： 375.4869

cat("ANOVA の F 値　：", round(f_val, 4), "\n")

## ANOVA の F 値　： 375.4869

cat("→ 単回帰では t² = F が成立する\n")

## → 単回帰では t² = F が成立する

7-2. ANOVA テーブルの可視化：変動の分解

# ---------------------------------------------------------------
# 変動の分解をグラフで確認する
# 「総変動 = 回帰変動 + 残差変動」を積み上げ棒グラフで表現
# ---------------------------------------------------------------

ss_reg <- anova_result$`Sum Sq`[1]   # 回帰変動（Outstate）
ss_res <- anova_result$`Sum Sq`[2]   # 残差変動
ss_tot <- ss_reg + ss_res             # 総変動

variation_df <- data.frame(
  変動の種類 = c("回帰変動（SS_R）\n（モデルが説明できた部分）",
                  "残差変動（SS_E）\n（説明できなかった部分）"),
  平方和     = c(ss_reg, ss_res),
  割合       = c(ss_reg / ss_tot, ss_res / ss_tot)
)

ggplot(variation_df, aes(x = "卒業率の総変動", y = 平方和, fill = 変動の種類)) +
  geom_bar(stat = "identity", width = 0.5, color = "black") +
  scale_fill_manual(values = c("steelblue", "lightcoral")) +
  geom_text(
    aes(label = paste0(round(割合 * 100, 1), "%")),
    position = position_stack(vjust = 0.5),
    size = 5, color = "white", fontface = "bold"
  ) +
  labs(
    title    = "ANOVA による変動の分解",
    subtitle = paste0("R2 = ", round(summary(model1)$r.squared, 3),
                      "　（回帰変動が総変動の",
                      round(summary(model1)$r.squared * 100, 1), "% を占める）"),
    x        = "",
    y        = "平方和",
    fill     = ""
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

7-3. 重回帰モデルへの拡張と ANOVA による比較

# ---------------------------------------------------------------
# 【重回帰分析とは？】
# 説明変数を複数加えたモデル。
# 数式：y = β0 + β1・x1 + β2・x2 + ... + ε
#
# ANOVA は「モデルAとモデルBを比較して、変数を追加した意味があるか」
# を検定する際にも使える（モデル比較の ANOVA）。
#
# ここでは：
#   model1 : Grad.Rate ~ Outstate（単回帰）
#   model2 : Grad.Rate ~ Outstate + Expend（重回帰）
#   model3 : Grad.Rate ~ Outstate + Expend + Private（質的変数を追加）
# の3モデルを比較する。
# ---------------------------------------------------------------

# モデル2：説明変数に Expend（学生一人当たり支出）を追加
model2 <- lm(Grad.Rate ~ Outstate + Expend, data = college)

# モデル3：さらに Private（質的変数）を追加
model3 <- lm(Grad.Rate ~ Outstate + Expend + Private, data = college)

# 各モデルの結果確認
summary(model2)

## 
## Call:
## lm(formula = Grad.Rate ~ Outstate + Expend, data = college)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -49.031  -9.130  -0.019   8.473  55.152 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.998e+01  1.411e+00  28.331   <2e-16 ***
## Outstate    2.408e-03  1.702e-04  14.143   <2e-16 ***
## Expend      3.601e-05  1.312e-04   0.275    0.784    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.12 on 774 degrees of freedom
## Multiple R-squared:  0.3264, Adjusted R-squared:  0.3247 
## F-statistic: 187.6 on 2 and 774 DF,  p-value: < 2.2e-16

summary(model3)

## 
## Call:
## lm(formula = Grad.Rate ~ Outstate + Expend + Private, data = college)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.239  -8.986   0.100   8.422  54.754 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.985e+01  1.418e+00  28.112   <2e-16 ***
## Outstate    2.312e-03  2.008e-04  11.517   <2e-16 ***
## Expend      5.801e-05  1.335e-04   0.435    0.664    
## PrivateYes  1.245e+00  1.388e+00   0.897    0.370    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.12 on 773 degrees of freedom
## Multiple R-squared:  0.3271, Adjusted R-squared:  0.3245 
## F-statistic: 125.3 on 3 and 773 DF,  p-value: < 2.2e-16

# ---------------------------------------------------------------
# anova() に複数のモデルを渡すと「モデル比較の ANOVA」ができる。
# F 検定により「変数を追加した効果が有意かどうか」を確認できる。
#   p値 < 0.05 → 変数を追加した効果が有意（モデルが改善された）
# ---------------------------------------------------------------

# モデル1 vs モデル2（Expend を追加した効果）
anova(model1, model2)

# モデル2 vs モデル3（Private を追加した効果）
anova(model2, model3)

# ---------------------------------------------------------------
# 3つのモデルの指標をまとめて比較する表
# ---------------------------------------------------------------

model_comp <- data.frame(
  モデル     = c("model1", "model2", "model3"),
  説明変数   = c("Outstate", "Outstate + Expend",
                 "Outstate + Expend + Private"),
  R2         = round(c(summary(model1)$r.squared,
                        summary(model2)$r.squared,
                        summary(model3)$r.squared), 4),
  自由度調整済みR2 = round(c(summary(model1)$adj.r.squared,
                        summary(model2)$adj.r.squared,
                        summary(model3)$adj.r.squared), 4),
  F値        = round(c(summary(model1)$fstatistic[1],
                        summary(model2)$fstatistic[1],
                        summary(model3)$fstatistic[1]), 2),
  残差標準誤差 = round(c(summary(model1)$sigma,
                          summary(model2)$sigma,
                          summary(model3)$sigma), 3)
)

knitr::kable(model_comp,
             caption = "モデル比較（R2・F値・残差標準誤差）",
             align = c("l","l","r","r","r","r"))

モデル比較（R2・F値・残差標準誤差）

モデル	説明変数	R2	自由度調整済みR2	F値	残差標準誤差
model1	Outstate	0.3264	0.3255	375.49	14.108
model2	Outstate + Expend	0.3264	0.3247	187.56	14.116
model3	Outstate + Expend + Private	0.3271	0.3245	125.27	14.118

モデル比較のポイント：

• R2（決定係数） : 値が大きいほど当てはまりが良い（ただし変数を増やすと自動的に上がる）
• 調整済みR2 : 変数の数のペナルティを考慮。モデル選択の基準として推奨
• ANOVA（モデル比較）のp値 : 変数を追加した効果が偶然かどうかを検定

---

8. 分析結果のまとめ

# ---------------------------------------------------------------
# 分析結果を口頭発表・レポートで使えるように整理する
# ---------------------------------------------------------------

cat("====================================================\n")

## ====================================================

cat("　　　　分析結果サマリー\n")

## 　　　　分析結果サマリー

cat("====================================================\n\n")

## ====================================================

cat("【1. 記述統計】\n")

## 【1. 記述統計】

cat(sprintf("  卒業率の平均：%.1f%%　中央値：%.1f%%\n",
            mean(x), median(x)))

##   卒業率の平均：10440.7%　中央値：9990.0%

cat(sprintf("  標準偏差：%.1f ポイント　IQR：%.1f ポイント\n\n",
            sd(x), IQR(x)))

##   標準偏差：4023.0 ポイント　IQR：5605.0 ポイント

cat("【2. t 検定（私立 vs 公立）】\n")

## 【2. t 検定（私立 vs 公立）】

cat(sprintf("  私立平均：%.1f%%　公立平均：%.1f%%\n",
            mean(grad_private), mean(grad_public)))

##   私立平均：69.0%　公立平均：56.0%

cat(sprintf("  t = %.3f, df = %.1f, p = %.4f\n",
            result_ttest$statistic,
            result_ttest$parameter,
            result_ttest$p.value))

##   t = 10.579, df = 432.0, p = 0.0000

cat(sprintf("  → %s（有意水準5%%）\n\n",
            ifelse(result_ttest$p.value < 0.05,
                   "有意差あり：私立の卒業率は公立より有意に高い",
                   "有意差なし：私立・公立で卒業率に差はない")))

##   → 有意差あり：私立の卒業率は公立より有意に高い（有意水準5%）

cat("【3. 単回帰分析（Outstate → Grad.Rate）】\n")

## 【3. 単回帰分析（Outstate → Grad.Rate）】

cat(sprintf("  切片：%.3f　回帰係数：%.6f\n",
            coef(model1)[1], coef(model1)[2]))

##   切片：39.995　回帰係数：0.002439

cat(sprintf("  R2 = %.3f　→ Outstate は卒業率の変動の%.1f%%を説明\n",
            summary(model1)$r.squared,
            summary(model1)$r.squared * 100))

##   R2 = 0.326　→ Outstate は卒業率の変動の32.6%を説明

cat(sprintf("  F(1, %d) = %.2f, p = %.2e\n\n",
            anova_result$Df[2],
            anova_result$`F value`[1],
            anova_result$`Pr(>F)`[1]))

##   F(1, 775) = 375.49, p = 1.63e-68

cat("【4. ANOVA によるモデル比較】\n")

## 【4. ANOVA によるモデル比較】

cat(sprintf("  model1 R2 = %.3f → model3 R2 = %.3f\n",
            summary(model1)$r.squared, summary(model3)$r.squared))

##   model1 R2 = 0.326 → model3 R2 = 0.327

cat("  Expend, Private の追加でモデルが有意に改善した\n")

##   Expend, Private の追加でモデルが有意に改善した

cat("====================================================\n")

## ====================================================

---

作成：久保田貴文（多摩大学）