R の線形回帰: lm()、summary()、predict()

考え方: 最小二乗直線

線形回帰は点の雲の中に直線を当てはめます: y = intercept + slope × x。ありうるすべての直線のうち、lm() は 残差の平方和 を最小にするものを選びます。残差とは点と直線の間の垂直方向の隔たりです。2 乗することで大きな外しが不釣り合いに重く数えられ、それが 1 つの極端な外れ値が当てはめ全体を傾けうる理由です。

相関が「どれだけ密に連動しているか」を単位のない 1 つの数値で与えるのに対し、回帰は 方程式 を与えます。単位があり、解釈できる傾きがあり、予測のための仕組みが備わっています。

fit <- lm(mpg ~ wt, data = mtcars)
print(coef(fit))

formula は「mpg を wt の関数としてモデル化する」と読みます。2 つの係数が当てはめられた直線です: mpg ≈ 37.3 − 5.3 × 重量。傾きには実際の単位があります。車が 1000 ポンド重くなる（wt は 1000 ポンド単位です）ごとに、およそ 5.3 マイル/ガロンの燃費を失います。切片（重量ゼロで 37.3 mpg）は直線がゼロを横切る位置にすぎません。重さのない車は存在しないので、読み込みすぎないでください。

summary() を順に読む

summary(fit) は、どの統計学の授業でも解釈を求められる出力です。実行して、ブロックごとに見ていきましょう。

fit <- lm(mpg ~ wt, data = mtcars)
summary(fit)

Call — 当てはめたモデルを反復表示します。今は些細ですが、6 つのモデルオブジェクトをやりくりしているときには命綱になります。

Residuals — 残り（実測値 − 予測値）の五数要約です。中央値がゼロ付近で、Min/Max と 1Q/3Q がおおむね対称であることが望ましい状態です。強い非対称性は、直線モデルが何かを見落としていることを示唆します。

Coefficients — 出力の心臓部で、項ごとに 1 行です。

• Estimate — 当てはめられた値です。wt では −5.34、つまり 1000 ポンド増えるごとに、燃費がおよそ 5.3 mpg 下がることと関連している ということです。傾きは必ず単位付きの文に翻訳しましょう。その一文がモデルの実用的な内容のすべてです。
• Std. Error — 標本を取り直したときに推定値がどれだけ揺らぐかを表します。ゼロから標準誤差 2 個分以内にある推定値はぐらついています。
• t value — Estimate ÷ Std. Error。係数がゼロから標準誤差何個分離れているかです（ここでは −9.56）。
• Pr(>|t|) — 「この係数は本当にゼロでありうるか?」に対する p 値です。wt では約 1.3e-10 です。重量が mpg と線形の関係を本当に持たないなら、32 個の標本でこれほど急な傾きが現れることは事実上ありません。p 値が小さい = 関連が存在する証拠であって、モデルが 正しい ことの証明ではなく、重要性の尺度でもありません（精確に推定されたごく小さな効果にも小さな p 値が付きます）。
• Signif. codes / 星印 — p 値の列を視覚的に略記したものです。便利ですが、情報は増えません。

Residual standard error: 3.05 on 30 degrees of freedom — 予測の外し方の典型的な大きさを、応答変数自身の単位で表したものです。予測はおおむね 3 mpg ほどずれます。mpg のスケール（およそ 10〜34 の範囲）と照らして判断しましょう。

Multiple R-squared: 0.75 — 重量が mpg の分散のおよそ 75% を説明します。Adjusted R-squared (0.74) は説明変数ごとの罰則を加えて同じものを計算し直したものです。生の版は変数を追加すれば — たとえ純粋なノイズでも — 増えることしかできないからです。説明変数の数が異なるモデルを比較するときは、調整済みの方が誠実です。そして「良いモデル = 高い R²」という反射は抑えましょう。本当に有用な効果が低 R² のモデルの中にあることもあります（応答がノイジー、多数の要因のひとつ、など）し、高い R² が過適合や漏洩した変数から来ていることもあります。

F-statistic: 91.4 ... p-value: 1.29e-10 — モデル全体の検定です。このモデルは「全員に平均を予測するだけ」に勝るか? 説明変数が 1 つのときは傾きの t 検定と重複します（9.56² ≈ 91.4 に注目）。複数の説明変数があるときは、少なくとも 1 つ の係数がゼロでないという同時検定になります。その仕組みは分散分析と同じ分散分解です。

重回帰: 他を一定に保つ

説明変数は + で追加します。

fit2 <- lm(mpg ~ wt + hp, data = mtcars)
summary(fit2)

解釈が 1 つの決定的な点で変わります。各 Estimate は、いまや 他を一定に保ったうえでの その説明変数の効果です。wt の係数（−5.3 から下がって約 −3.9）は、同じ馬力の車どうしを比べたときの 重量増加による燃費の代償です。単回帰の −5.3 は、重い車ほど馬力も高い傾向があるという事実を黙って抱き合わせていました。重回帰はそれを解きほぐします。これはまた、変数を追加すると係数がずれる理由でもあります。新しい説明変数が既存のものと相関していれば、既存のものの職務内容が変わるのです。決定係数は約 0.83 に上がり、ここでは 調整済み の版が単一説明変数モデルとの公正な比較になります。

予測: predict()

当てはめたモデルは関数であり、predict() がそれを評価します。列名が説明変数と正確に一致する newdata のデータフレームを作ります。

fit <- lm(mpg ~ wt, data = mtcars)

new_cars <- data.frame(wt = c(2.5, 3.5))
print(predict(fit, newdata = new_cars))

print(predict(fit, newdata = new_cars, interval = "confidence"))
print(predict(fit, newdata = new_cars, interval = "prediction"))

2 種類の区間は異なる問いに答えるもので、これを取り違えるのは定番の試験でのミスです。

• interval = "confidence" — 平均 に関する不確実性です。「2500 ポンドの車すべてについて、平均の mpg はどこにあるか?」狭く、データが増えるほど縮みます。
• interval = "prediction" — その重量の 個々の新しい車 が収まりそうな範囲です。はるかに広くなります。1 台の車は直線のまわりに自分自身のばらつきを持ち、そのばらつきはどれだけデータを増やしても平均化して消えないからです。

1 つの新しい観測値についての問いに対して信頼区間を報告すると、精度を大幅に誇張することになります。

診断と外挿の罠

summary() はモデルが何を推定したかを教えてくれます。それを信じてよいかを教えてくれるのは残差プロットです。対話的なセッションでは次のようにします。

par(mfrow = c(2, 2))
plot(fit)   # four diagnostic plots

見るべきポイントはこうです。Residuals vs Fitted は形のない雲であるべきです。曲線が見えるなら関係は直線的ではありません。漏斗状（予測値が大きくなるほど散らばりが広がる）なら分散が一定でなく、標準誤差が狂っています。Q-Q プロット の点は直線に沿うべきです。裾が重ければ、外れ値が当てはめを歪めています。Scale-Location は再び漏斗の確認です。Residuals vs Leverage は影響力の大きい点 — それ単独で係数を引きずる観測値 — を示します（mtcars では Chrysler Imperial のような特殊な車がここに現れがちです）。当てはめる前に生データの散布図を手早く描けば、この多くは早い段階で捕まえられます。

最後に、どんな診断でも捕まえられない罠があります。外挿 です。モデルはおよそ 1500〜5400 ポンドの車から学習しました。predict() に wt として 8 を渡せば、平然と 負の mpg を返してきます。数学は直線をどこまでも延長しますが、証拠はデータの端で止まっています。予測は、当てはめに使った範囲の内側（またはその近く）でのみ行いましょう。

この章のまとめ

• fit <- lm(y ~ x, data = df) が最小二乗直線を当てはめます。coef(fit) が方程式、summary(fit) が完全な報告です。
• Estimate は単位付きの文として読みましょう。Pr(>|t|) は「これはゼロでありうるか?」を問うのであって、「これは重要か?」ではありません。
• 残差標準誤差は実際の単位での典型的な外し幅、調整済み決定係数がモデル比較のための公正な数値です。
• 重回帰ではすべての係数が「他を一定に保ったうえで」を意味し、相関する説明変数が加わると係数はずれます。
• predict(fit, newdata, interval = ...): 平均には "confidence"、1 件の新しいケースには "prediction"（広い方）を使います。
• plot(fit) で曲がり、漏斗、影響力の大きい点を確認しましょう。データの範囲外の予測は決して信用してはいけません。

次は、結果が数値ではなく yes/no のとき — glm() によるロジスティック回帰です。