2変量の関連(3)

2つの連続変数同士の関連を調べる方法として、相関と回帰についてまとめておきます。

このページでも、データはすでに読み込んでいるという前提でコードを書いております。

相関

2つの連続変数同士の相関は、cor.test()関数を用いて調べることができます。

cor.test(1つめの連続変数, 2つめの連続変数, method = "pearsonもしくは指定なしか、spearman、kendall")

身長と体重に相関関係が見られるかどうかを検討するためのコードは下記の通りになります。最初に図示して、分布を確かめることはもちろん必須です。

## 相関の検定(身長と体重)

2つの変数間に相関関係があるか確かめる場合に行う。

```{r cortest}

plot(dat$height, dat$weight)

cor.test(dat$height, dat$weight)

```

> cor.test(dat$height, dat$weight)

Pearson's product-moment correlation

data: dat$height and dat$weight
t = 9.3422, df = 93, p-value = 5.041e-15
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.5748646 0.7869766
sample estimates:
cor
0.6957931

結果から、身長と体重の間には、有意な相関関係があり、相関係数は0.70であることがわかります。また、このとき、相関係数の95%信頼区間は、[0.57, 0.79]となります。

相関関係をみる場合、集中楕円を描くことがある。利用する関数はdataEllipse()関数で、carパッケージに入っている。

## 集中楕円の描画

集中楕円は、欠損値があると描画されないので、最初にサブセットを作成してから実行します。

```{r corell}

library("car")

dat.E = dat[, c("height", "weight")]

dat.E = subset(dat.E, complete.cases(dat.E))

dataEllipse(dat.E$height, dat.E$weight, levels = 0.95)

```

回帰

回帰は、一方の変数(応答変数)のばらつきを他方の変数(説明変数)のばらつきでどの程度説明できるのかをモデルを当てはめて考える手法です。説明力が高ければ、説明変数の値から、対応する応答変数の値が推定することが期待できます。

回帰分析は、lm()関数summary()関数を用いて行うことが一般的です。

結果を入れる変数 = lm(応答変数 ~ 説明変数, data = データフレーム)
summary(結果を入れる変数)

収入のばらつきが体重のばらつきをどの程度説明できるかを検討するコードは、下記の通りです。この場合も図示は必須です。

## 回帰分析(収入と体重)

説明変数のばらつきがどの程度応答変数のばらつきを説明するのか調べる。

```{r lmtest}

plot(weight ~ income, data = dat)

res.lm = lm(weight ~ income, data = dat)

summary(res.lm)

```

> res.lm = lm(weight ~ income, data = dat)
> summary(res.lm)

Call:
lm(formula = weight ~ income, data = dat)

Residuals:
 Min   1Q  Median  3Q   Max
-25.0568 -8.9324 0.7626 7.5220 24.6571

Coefficients:
     Estimate Std. Error t value Pr(>|t|)
(Intercept) 66.9190 2.6152  25.59 < 2e-16 ***
income   0.8977  0.2250  3.99 0.000134 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11.3 on 90 degrees of freedom
(8 observations deleted due to missingness)
Multiple R-squared: 0.1503, Adjusted R-squared: 0.1409
F-statistic: 15.92 on 1 and 90 DF, p-value: 0.0001343

結果から、モデルが有意で、説明変数である収入のばらつきが応答変数である体重のばらつきを説明していることが言えそうです。分散分析表の中の変数が有意であるということは、モデル式における回帰係数が有意に0とは異なることを示しています。この場合、収入の回帰係数であるincomeが0.90程度で、有意に0ではないということです。下のコードを用いることで、信頼区間と予測区間を図示することもできます。

### 信頼区間と予測区間の図示

信頼区間と予測区間は、predict()関数を用いて図示できます。

```{r pred}

income_new = data.frame(income = 1:33) # 予測したい収入のデータを作成

conf_95 = predict(res.lm, income_new, interval = 'confidence', level = 0.95) # 95%信頼区間

pred_95 = predict(res.lm, income_new, interval = 'prediction', level = 0.95) # 95%予測区間

# 信頼区間シンプル版プロット

matplot(income_new, conf_95[, 1:3], type = "l", lty = c(1, 2, 2), col = c(1, 4, 4), lwd = c(3, 1, 1))

points(weight ~ income, data = dat)

# 予測区間シンプル版プロット

matplot(income_new, pred_95[, 1:3], type = "l", lty = c(1, 2, 2), col = c(1, 4, 4), lwd = c(3, 1, 1))

points(weight ~ income, data = dat)

```

新しい関数がいくつか出てくるので面食らうかも知れませんが、ぜひヘルプを参照して使えるようになってください。predict()関数を用いて、信頼区間(conf_95)と予測区間(pred_95)を作成し、matplot()関数で描画し、最後にpoint()関数で散布図を重ねています。慣れれば、図の重ね描きも問題なくできるようになると思います。