R の散布図: plot()、回帰直線、pairs()

散布図が示すもの

散布図は 2 つの 数値変数の関係を表示します。各観測値が 1 つの点になり、第 1 の値で横軸上の位置が、第 2 の値で縦軸上の位置が決まります。2 つの変数が連動していれば点はパターンを作り、そうでなければ形のない雲になります。モデルを当てはめる前の標準的な最初の一瞥であり、R は plot() のガイドで扱ったのと同じ plot() 関数でこれを描きます。

ここでは 32 台の車が入った組み込みデータセット mtcars を使い、物理的な問いを立てます。重い車ほどエンジンが強力なのでしょうか?

グラフを描く

wt は千ポンド単位の重量、hp は馬力です。

plot(mtcars$wt, mtcars$hp,
     main = "Horsepower vs. weight",
     xlab = "Weight (1000 lbs)",
     ylab = "Horsepower",
     pch  = 19,
     col  = "steelblue")

描かれる絵は次のとおりです。32 個の塗りつぶした点が、グラフの左下（軽い車、重量スケールで 1.5 付近、馬力は 60〜90 あたり）から、5 を超える重い車で 200 馬力を突破する領域へと上っていきます。この上昇は紛れもありませんが、整然としてはいません。同じ重量でも点は馬力のかなり広い帯に散らばります。

装飾はいつもの base 作図の道具立てです。塗りつぶしの丸には pch = 19（既定の中空の丸はスクリーンショットで消えてしまいます）、色には col、点を大きくしたいなら cex = 1.3 です。第 3 のカテゴリ変数で点を色分けするには、因子で色ベクトルに添字を付けます — col = c("tomato", "steelblue", "darkgreen")[factor(mtcars$cyl)] で、シリンダー数ごとに固有の色が付きます。

読み方: 方向、強さ、形

散布図を見るときは毎回、次の 3 つの問いをこの順で立てます。

• 方向。 横軸に沿って見ていったとき、点は上昇しますか（正の関係）、それとも下降しますか（負の関係）。ここでは上昇しています。重いほど強力です。代わりに mpg を wt に対して描くと雲は下降します。重いほど燃費が悪いということです。
• 強さ。 点は 1 本の道筋にどれだけ密着していますか。鉛筆のように細い帯なら強い関係、ゆるやかな噴霧なら弱い関係です。この雲は中程度の密着度です。
• 形と意外性。 道筋は直線ですか、曲線ですか。クラスタはありますか。他から大きく離れた点はありますか。mtcars では Maserati Bora が集団の上に目立って位置しています。中程度の重量で 335 馬力です。こうした 1 点は当てはめた直線を目に見えて引っ張ることがあり、まさにそれが、当てはめる前に見るべき理由です。

傾向線を追加する

散布図は関係を提示し、そこを通る直線はその主張を要約します。線形モデルを当てはめて、そのまま abline() に渡しましょう。

plot(mtcars$wt, mtcars$hp,
     pch = 19, col = "steelblue",
     xlab = "Weight (1000 lbs)", ylab = "Horsepower")

abline(lm(hp ~ wt, data = mtcars), col = "tomato", lwd = 2)

lm(hp ~ wt) は最小二乗直線を当てはめます。この formula は「wt で説明される hp」と読み、縦軸の変数が ~ の前に来ます。そして abline() がそれをグラフ上に描きます。この直線は千ポンドあたり約 46 馬力の割合で上昇します。このモデルが何を意味し、その summary をどう読むかは線形回帰の主題です。

直線を仮定したくない場合は、lowess() がデータの行くところにどこまでも従う滑らかな曲線を描きます。

lines(lowess(mtcars$wt, mtcars$hp), col = "darkgreen", lwd = 2, lty = 2)

lowess の曲線と直線がおおむね一致していれば、線形の要約は妥当です。曲線が離れて曲がっていくなら、関係は非線形であり、直線ではそれを誤って表現してしまいます。

cor() で数値を確認する

グラフは形を、cor() は強さを 1 つの数値として教えてくれます。この手順は純粋なテキスト出力なので、ここで実行できます。

# Correlation between weight and horsepower:
print(cor(mtcars$wt, mtcars$hp))

# And the full matrix for three variables at once:
print(round(cor(mtcars[, c("mpg", "wt", "hp")]), 2))

重量と馬力の相関は約 0.66 です。中程度に密着した上昇する雲を、数値にしたものです。行列を見ると、mpg は両方と強く 負 に相関していることが分かります（重量とは約 −0.87）。ただし手順の順序は守りましょう。まずグラフ、次に係数です。1 つの r 値は曲がりを隠したり、外れ値 1 つで水増しされたりします。よくある誤解のパターンは相関を参照してください。

散布図行列: pairs()

数値の列が複数あるとき、すべての組み合わせを手で描くのはうんざりします。pairs() は 1 回の呼び出しでそれをやってくれます。

pairs(mtcars[, c("mpg", "wt", "hp")],
      pch = 19, col = "steelblue")

結果は 3 × 3 のグリッドです。対角線には変数名が並び、対角以外の各パネルは 1 組の散布図です — mpg 対 wt、mpg 対 hp、wt 対 hp が、軸を入れ替えてそれぞれ 2 回ずつ現れます。新しいデータセットをふるいにかける最速の方法であり、一目でどの組が関係しているか、どの関係が曲がっているか、外れ値がどこに隠れているかが分かります。ここでのように、まず列を絞り込んでください。変数が 6〜7 個を超えるとパネルは判読できないほど縮んでしまいます。

ggplot2 版

ggplot2 では、散布図に当てはめた線を加えるのは 2 つのレイヤーです。

library(ggplot2)

ggplot(mtcars, aes(x = wt, y = hp)) +
    geom_point(color = "steelblue", size = 2) +
    geom_smooth(method = "lm", color = "tomato") +
    labs(title = "Horsepower vs. weight",
         x = "Weight (1000 lbs)", y = "Horsepower")

geom_smooth(method = "lm") は abline(lm(...)) におまけが付いたものです。直線の周りに影付きの信頼帯が描かれます。method を指定しなければ、代わりに loess 曲線が当てはめられます。lowess() の ggplot2 版です。ざっと見るだけなら base R が入力速度で勝ちますが、点をグループで色分けして凡例を自動生成したくなった瞬間に ggplot2 が勝ちます。

この章のまとめ

• 数値ベクトル 2 つを渡した plot(x, y) が散布図です。pch = 19 とラベル付きの軸で見栄えが整います。
• 何かを計算する 前に、方向・強さ・形を読み、外れ値を見つけましょう。
• abline(lm(y ~ x, data = df)) は回帰直線を追加し、lines(lowess(x, y)) は直線性を仮定しない曲線を追加します。
• cor() はグラフが示すものを数値化し、グラフはその数値を誠実に保ちます。
• pairs(df[, cols]) はすべての組の散布図を一度に描きます。新しいデータセットを最速でふるいにかける方法です。

次は棒グラフ — 数値の組を離れ、カテゴリ間で度数を比較します。