相関係数と散布図
Ryota Mugiyama
March 9th , 2018
目的
よく教科書に載っている「相関係数がのときの散布図」をRで作る。
*この記事はこちらのウェブサイトに置かれています。Rのコードに関して前嶋直樹さん（東京大学大学院）よりアドバイスをいただきました。記して感謝いたします。
パッケージの準備
#以下のパッケージを使用する。入っていない場合はインストール。
#install.packages(tidyverse)
#install.packages(knitr)
#install.packages(purrr)
library(tidyverse)
library(knitr)
library(purrr) #別解のやり方をするときに使用。手順
2次元正規乱数を発生させる関数を定義する（http://cse.naro.affrc.go.jp/takezawa/r-tips/r/60.html）。
r2norm <- function(n, mu, sigma, rho) {
tmp <- rnorm(n)
x <- mu+sigma*tmp
y <- rho*x + sqrt(1-rho^2)*rnorm(n)
return(data.frame(x=x,y=y))
}この関数を用いて、相関係数-1から+1まで合計21個のサンプルを生成する。ここではそれぞれサンプルサイズは5000とした。
time <- seq(-1, 1, length=21)
size <- 5000
cordata <- data.frame()
for (i in time){
cor <- i
cordata_mini <- r2norm(size, 0, 1, i)
cordata_mini$cor <- cor
cordata <- rbind(cordata, cordata_mini)
}なおこの部分はpurrrというパッケージを使って書くこともできる。今回はコメントアウトしておく。
#cor = seq(-1,1,length = 21)
#data_list <- purrr::map(.x = cor,.f = r2norm,n = 5000,mu = 0,sigma = 1)
#data_list <- purrr::map2(.f=cbind,.x=data_list,.y=cor)
#cordata <- dplyr::bind_rows(data_list) %>% magrittr::set_colnames(c("x","y","cor"))生成したサンプルの相関係数が真の値に近いかどうかを確認しておく。生成したサンプルのサイズが大きいほど、真の値に近い値を得られる。今回はそれぞれ5000ケースとたくさんのケースを抽出したので、いずれもそれなりに真の値に近い値が得られると思われる。
corsample <- data.frame()
for (i in time){
corsub <- cor(cordata$x[cordata$cor == i], cordata$y[cordata$cor == i])
cortrue <- i
corsub <- cbind(cortrue, corsub)
corsample <- rbind(corsample, corsub)
}
## 95%信頼区間も確認しておく
corsample$lowerCI <- corsample$corsub - 1.96*sqrt((1 - corsample$corsub^2)/(size - 2))
corsample$upperCI <- corsample$corsub + 1.96*sqrt((1 - corsample$corsub^2)/(size - 2))
colnames(corsample) <- c("True", "Sample", "lower 95% CI", "upper 95% CI")
kable(corsample)| True | Sample | lower 95% CI | upper 95% CI |
|---|---|---|---|
| -1.0 | -1.0000000 | -1.0000000 | -1.0000000 |
| -0.9 | -0.9030635 | -0.9149712 | -0.8911558 |
| -0.8 | -0.7984840 | -0.8151744 | -0.7817936 |
| -0.7 | -0.7015711 | -0.7213273 | -0.6818149 |
| -0.6 | -0.5972078 | -0.6194450 | -0.5749706 |
| -0.5 | -0.5063585 | -0.5302656 | -0.4824513 |
| -0.4 | -0.4073308 | -0.4326508 | -0.3820109 |
| -0.3 | -0.2838236 | -0.3104077 | -0.2572396 |
| -0.2 | -0.2338958 | -0.2608509 | -0.2069407 |
| -0.1 | -0.0840388 | -0.1116648 | -0.0564127 |
| 0.0 | 0.0228504 | -0.0048665 | 0.0505672 |
| 0.1 | 0.1119374 | 0.0843875 | 0.1394873 |
| 0.2 | 0.1823351 | 0.1550757 | 0.2095945 |
| 0.3 | 0.2971524 | 0.2706806 | 0.3236243 |
| 0.4 | 0.3827103 | 0.3570969 | 0.4083238 |
| 0.5 | 0.4877326 | 0.4635296 | 0.5119356 |
| 0.6 | 0.6006996 | 0.5785348 | 0.6228643 |
| 0.7 | 0.7074617 | 0.6878676 | 0.7270557 |
| 0.8 | 0.7908676 | 0.7739008 | 0.8078344 |
| 0.9 | 0.9010751 | 0.8890521 | 0.9130980 |
| 1.0 | 1.0000000 | 1.0000000 | 1.0000000 |
サンプルサイズ（size）が大きければ大きいほど生成したサンプルから得られる相関係数は真の値に近づくが、その分あとで作成する散布図がつぶれてしまうので難しいところ。
相関係数の値ごとに、散布図をプロット。相関係数-1の場合は自明なので省略。
cordata1 <- subset(cordata, cor > -1)
p1 <- ggplot(cordata1, aes(x = x, y = y)) +
geom_point() +
facet_wrap(~cor) +
xlim(-4,4) +
ylim(-4,4)
p1## Warning: Removed 3 rows containing missing values (geom_point).
いくつかピックアップして抜き出すと次のようになる。
cordata$cor <- round(cordata$cor, digits = 1)
cordata2 <- subset(cordata, abs(cor) == 0.9 | abs(cor) == 0.7 | abs(cor) == 0.5 | abs(cor) == 0.3 | abs(cor) == 0)
p2 <- ggplot(cordata2, aes(x = x, y = y)) +
geom_point() +
facet_wrap(~cor) +
xlim(-4,4) +
ylim(-4,4)
p2## Warning: Removed 2 rows containing missing values (geom_point).
ただこれだとちょっと点が大きすぎ、相関0.7、0.5、0.3の違いがやや分かりにくい。そこで少し点を小さくしてみる。
cordata$cor <- round(cordata$cor, digits = 1)
cordata2 <- subset(cordata, abs(cor) == 0.9 | abs(cor) == 0.7 | abs(cor) == 0.5 | abs(cor) == 0.3 | abs(cor) == 0)
p3 <- ggplot(cordata2, aes(x = x, y = y)) +
geom_point(aes(x = x, y = y), size = 0.5) +
facet_wrap(~cor) +
xlim(-4,4) +
ylim(-4,4)
p3## Warning: Removed 2 rows containing missing values (geom_point).