Ma.05.Rで簡単な「予測」を試みる

Rを使って「予測」を試みて、その結果をどう解釈すればよいかを復習します。
相関分析から線形回帰、ロジスティック回帰まで、ステップバイステップで進めます。
「モデルを作って終わり」ではなく、そのモデルが「本当に使えるのか？」を診断し、考察するプロセスを重視します。

Step 1. 関連性の確認（相関分析）

予測モデルを作る前に、そもそも変数同士に関係があるのか（相関）を調べます。

• 2つの変数がどの程度「直線的に」関連しているか（ピアソン）？
• または、順序として関連しているか（スピアマン、ケンドール）？
• 多数の変数間の関連性を一覧（ヒートマップ）で視覚化したい。

コピー用コード

# --- 1. 2変数間の相関検定 ---
# ピアソンの積率相関係数 (2変数の正規性・独立性を仮定)
cor.test(df$var_A, df$var_B, method = "pearson")

# スピアマンの順位相関係数 (外れ値に頑健。順序関係)
cor.test(df$var_A, df$var_B, method = "spearman")

# ケンドールのτ (スピアマンよりサンプルサイズの影響を受けにくい)
cor.test(df$var_A, df$var_B, method = "kendall")


# --- 2. 多数の変数間の相関マトリックス (視覚化) ---
library(corrplot)
cor_mat <- cor(df13, method = "pearson", use = "complete.obs") # NAを除外

corrplot.mixed(cor_mat,
               upper = "color",  # 上半分は色
               lower = "number", # 下半分は数値
               tl.pos = "lt")    # 変数名は左上

出力のどこを見るか？ (cor.test)

# cor.test(…) の出力例
data: df$var_A and df$var_B
t = 10.5, df = 1000, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.29 0.35
sample estimates:
cor
0.32

• cor (相関係数): 0.32
  • メインの結果です。「+1」に近いほど強い正の相関、「-1」に近いほど強い負の相関、「0」は無相関を示します。
• 95 percent confidence interval: [0.29, 0.35]
  • 真の相関係数が含まれる範囲を示します。

Step 2. 線形回帰 (OLS)

相関が確認できたら、次に「X が 1 増えたら、Y がどれだけ増えるか？」を予測するモデル（数式）に関心が向きます。
まずは線形回帰です。

• Y = b0 + b1 * X という最も単純な予測モデルを作成する。
• モデルの当てはまり（フィット）が「良いか・悪いか」を視覚的に診断する。

コピー用コード

# --- 1. 単回帰モデルの作成 →　結果確認---
# lm( 目的変数Y ~ 説明変数X, data = データ )
model1 <- lm(jinko ~ syotoku, data = df)

summary(model1)

# --- 3. 残差診断 (必須！) ---
par(mfrow = c(2, 2))
plot(model1)
par(mfrow = c(1, 1)) # プロット設定を元に戻す

出力のどこを見るか？ (summary(model1))

# summary(model1) の出力例
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.7e+06 1.1e+06 -5.12 1.2e-06 ***
syotoku 2.9e+03 3.4e+02 8.50 3.3e-12 ***
—
Residual standard error: 1800000 on 98 degrees of freedom
Multiple R-squared: 0.424, Adjusted R-squared: 0.418
F-statistic: 72.2 on 1 and 98 DF, p-value: 3.3e-12

• Coefficients:
  • Estimate: モデルの「係数（切片と傾き）」です。
    • (Intercept) (切片): syotoku が0の時の jinko の予測値。
    • syotoku (傾き): syotoku が1単位増えると jinko が 2.9e+03（2900）増える、と予測されます。
  • Pr(>|t|) (p値): 係数が「今回以上に極端な統計量が出る」確率。
• Adjusted R-squared (決定係数 R2): 0.418
  • モデルの「説明力」です。jinko のバラツキのうち、syotoku によって 42.4% を説明できたことを示します。
  • 説明変数が一個のときはMultiple R-squaredで問題ないのですが、調整済みのほうを見るクセに害はないと思います。
• F-statistic の p-value: 3.3e-12
  • モデル「全体」が統計的に有意かどうか。このモデルは（ランダムに作るよりは）予測に役立つか？を示す。
    （仮説に基づいているので、ランダムにYを決めるよりは役立つように見える場合が多いため、残差診断プロットを優先する）

出力のどこを見るか？ (plot(model1))

• Residuals vs Fitted (残差 vs 予測値):
  • 理想: 点が水平な赤い線の周りに均一に散らばっている。
  • 注目点: 「フィットが悪い」「線が折れ曲がっている」。
  • もし点が均一でなく、U字型やラッパ型になっていたり、赤い線が水平でなく曲がっていたりする場合、このモデルの前提（線形性）が崩れています。
• Normal Q-Q (正規Q-Qプロット):
  • 理想: 点が点線の上にきれいに並んでいる。
  • 注目点: 点が点線から大きく外れている場合、残差が正規分布していない（＝予測が特定の方向に偏っている）ことを示します。

【その場での自問】このモデルは使えるか？
　→　Residuals vs Fitted が折れ曲がっている（均一でない）場合
　　　summary() の結果が良く見えても、線形性/等分散の仮定に懸念があるため、結論で活用はしない。（このモデルは信頼できないと考える）

対策:

1. 他の重要な変数（例: 東京ダミー）が抜けているのでは？ →　 Step 3: 重回帰へ
2. そもそも syotoku と jinko の関係は直線ではない？（例: log(jinko) を使うべき？）
3. 「そもそもY（jinko）の設定は適切か？」と問い直す勇気が重要。

Step 3. 重回帰分析 – モデルの改善と評価

複数の説明変数（X1, X2…）を投入して、より現実に即したモデルを目指します。

• 複数の要因（例: 所得＋大都市ダミー）を同時に考慮したモデルを作成する。
• モデルが満たすべき前提条件（多重共線性など）をチェックする。
• どの変数が「最もYに影響を与えているか」を比較する。

コピー用コード

# --- 1. 重回帰モデルの作成 →　結果---
# 説明変数を + でつなげて投入
qol_model <- lm(q1_1 ~ q3s1_1 + q3s2_1 + q3s3_1 + q3s4_1 + 
                       q3s5_1 + q3s6_1 + q3s7_1 + q3s8_1 + 
                       q3s9_1 + q3s10_1 + q3s11_1 + q3s12_1 + 
                       q3s13_1 + age_num,
                data = naiman_c12)

summary(qol_model)

# --- 2. 標準化偏回帰係数 (β) の算出 ---
# 係数(Estimate)の単位がバラバラだと比較できないため、
# 係数を標準化して「影響度」を比較できるようにする

library(lm.beta)
lm.beta(qol_model)

# --- 3. 多重共線性 (VIF) のチェック ---
# 説明変数同士が似すぎていないか？ (例: "数学Aの点数" と "数学Ⅰの点数")
# 似すぎているとモデルが不安定になる
library(car)
vif(qol_model)

# --- 4. 残差診断 (必須！) ---
par(mfrow = c(2, 2))
plot(qol_model)
par(mfrow = c(1, 1))

出力のどこを見るか？ (重回帰の評価ステップ)

1. モデルの「前提条件」は満たしているか？

• vif(qol_model) の出力:
  • 注目点: 各変数のVIF値。
  • 判断基準: この値が 10（厳しく見るなら5）を超えている変数があるか？
    あれば、その変数は他の変数と似すぎている（多重共線性）ため、モデルから除外するか、別の変数と統合することを検討します。
• plot(qol_model) の出力:
  • 注目点: Residuals vs Fitted が（ある程度）均一に散らばっているか？ Normal Q-Q が（ある程度）直線に乗っているか？
  • 判断基準: Step 2 と同じ。ここがひどく崩れていると、summary の結果は信頼できません。

2. モデルは「予測に役立っていない」わけではないか？

• summary(qol_model) の F-statistic の p-value:
  • 注目点: モデル全体のp値（一応）。
  • 判断基準: p < 0.05 か？ 万が一にも有意でなければ、そのモデルは単なるノイズを拾っているだけかもしれません。

3. Yを「どの程度」予測できるか？

• summary(qol_model) の Adjusted R-squared (調整済みR2):
  • 注目点: 重回帰では 変数の数を考慮した Adjusted R-squared を見ます。
  • 判断基準: この値が「どの程度か」（例: 0.42ならYの変動の42%を説明）。

4. 予測の「誤差」はどの程度か？

• summary(qol_model) の Residual standard error (RSE):
  • 注目点: 予測値の「平均的な誤差の大きさ」です。
  • 判断基準: この値（例: 10.5）を、目的変数Yのスケール（例: q3の範囲が0-130, 平均69）と比較します。
    「範囲が130で、平均69点のものに対し、予測誤差が平均10.5点もある」のは許容できるか？を考察します。

5. 「どの変数」が予測に重要か？

• lm.beta(qol_model) の Standardized Coefficients (Estimate):
  • 注目点: 標準化された係数（標準化偏回帰係数）です。
  • 判断基準: この絶対値が最も大きい変数が、Yの予測に最も強く（プラスまたはマイナスに）影響している因子だと解釈できます。
  • 回帰式に当てはめるものではないですが、一体どの説明変数、因子が強く影響を与えているのだろうか？という疑問に応えてくれる。

Step 4. ロジスティック回帰 – Yが「0か1か」の場合（Glmでも活用しやすい）

目的変数Yが jinko（連続値）ではなく、k6_h（0=low, 1=high）のような2値の場合、線形回帰 (lm) は使えません。代わりに glm を使います。

• 結果が0/1の変数（例: ハイリスク群かどうか）の「発生確率」を予測する。
• 各説明変数が、結果を「1」にしやすくする（オッズを上げる）か、「0」にしやすくする（オッズを下げる）かを評価する。

コピー用コード

# --- 1. ロジスティック回帰モデルの作成 ---
# family = binomial (二項分布) が「ロジスティック回帰」を指定するおまじない
logi_model1 <- glm(k6_h ~ q1_1 + q4s13_1 + q21 + q31_1 + age_num + sex,
                   data = naiman_c12, 
                   family = binomial)

summary(logi_model1)

# --- 2. オッズ比 (OR) の算出 ---
# 係数(Estimate)は「対数オッズ」で解釈しにくいため、
# exp() で「オッズ比」に変換する
exp(coef(logi_model1))

# --- 3. オッズ比の信頼区間 (CI) の算出 ---
# この信頼区間が「1」をまたぐかが最重要
exp(confint.default(logi_model1))

出力のどこを見るか？ (summary(logi_model1))

# summary(logi_model1) の出力例
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.500 0.200 7.50 < 2e-16 ***
q1_1 -0.085 0.010 -8.50 < 2e-16 ***
q4s13_1 0.107 0.020 5.35 8.0e-08 ***
age_num -0.042 0.003 -14.00 < 2e-16 ***
—
Null deviance: 8500 on 10000 degrees of freedom
Residual deviance: 7800 on 9994 degrees of freedom
AIC: 7814

• Coefficients (Estimate, Pr(>|z|):
  • Pr(>|z|) (p値) が 0.05 より小さいかを確認し、有意な変数（q1_1, q4s13_1, age_num）を特定します。
• Null Deviance (ヌル逸脱度): 8500
  • 「ベースラインの不確実性」。説明変数を一切使わなかった場合のモデルの当てはまりの悪さ（不確実性）です。
• Residual Deviance (残差逸脱度): 7800
  • 「モデルの不確実性」。q1_1 などの説明変数を使った後でも、なお残っている当てはまりの悪さ（不確実性）です。
• AIC (赤池情報量基準): 7814
  • モデルの「悪さ」を示す指標（小さいほど良い）。これ単体では意味がなく、モデルを修正（変数を追加・削除）した後に、前のモデルとAICを比較するために使います。

モデルの評価と解釈 (重要)

1. モデルは役に立っているか？ (逸脱度の比較)

• 注目点: Null Deviance (8500) と Residual Deviance (7800) の差。
• 解釈: モデルをあてはめることで、不確実性は 8500 - 7800 = 700 減少しました。
• (考察): 不確実性は (700 / 8500) * 100 約 8.2% 減少しました。これは十分か？
  残りの約92%の不確実性を説明する、他の重要な変数（例: 経済状況、社会関係など）が抜けていないか？と考え続けます。

2. どの変数が強く関連しているか？ (オッズ比の解釈)

# exp(coef(logi_model1)) の出力例
(Intercept)     q1_1   q4s13_1   age_num 
     4.481     0.918     1.113     0.959 

# exp(confint.default(logi_model1)) の出力例 (95% CI)
                2.5 %    97.5 %
(Intercept)    3.031     6.621
q1_1           0.899     0.938  <- 1をまたいでいない
q4s13_1        1.070     1.158  <- 1をまたいでいない
age_num        0.953     0.965  <- 1をまたいでいない

• 注目点: exp(coef) (オッズ比, OR) と、その 95% CI が 「1」をまたいでいるか。
• 解釈:
  • q1_1: OR = 0.918 (CI: [0.899, 0.938])
    • CIが1をまたいでいないので、信頼できる結果です。
    • ORが1より小さい (0.918) ので、k6_h=1 (ハイリスク) に「なりにくい」因子（保護因子）です。
    • → 「q1_1 が1増えるごとに、ハイリスクになるオッズは 0.918 倍 (約 8.2% 減少) する」
  • q4s13_1: OR = 1.113 (CI: [1.070, 1.158])
    • CIが1をまたいでいない。
    • ORが1より大きい (1.113) ので、k6_h=1 に「なりやすい」因子（危険因子）です。
    • → 「q4s13_1 が1増えるごとに、ハイリスクになるオッズは 1.113 倍 (約 11.3% 増加) する」

【考察】モデルをどう改善するか？

1. summary で AIC (例: 7814) をメモしておく。
2. 仮説に基づき、新しい変数（例: income）をモデルに追加する (logi_model2)。
3. summary(logi_model2) を実行し、Residual Deviance が logi_model1 (7800) より減っているか？ AIC が (7814) より下がっているか？ を比較します。
4. AICが下がり、説明力（逸脱度の減少率）が上がれば、logi_model2 の方が良いモデルだと判断できます。