前回のおさらいと今回のイントロ

  前回は平均、分散、標準偏差、中央値など、データの特徴を捉える指標について学んだ。 前回の資料では、数学のテストの得点のみ、身長のみ、など1種類のデータの扱いについて学んだ。 今回の資料では、数学のテストと英語のテストの関係、身長と体重の関係、など2種類のデータの関係性を記述する指標について学ぶ。

前回の演習課題の主な回答

3: データの数が\(N\)(自然数)、\(i\)番目のデータの値が\(x_i\)(\(i = 1, ..., N\))、その平均を\(\bar{x}\)、分散を\(\mathrm{Var}(x)\)とするとき、以下の等式を証明せよ。

\(\mathrm{Var}(x) = \frac{1}{N}\sum_{i=1}^Nx_i^2 - \bar{x}^2\)

分散は2乗の平均から平均の2乗を引いたものになる。これが必ず0より大きくなることも確認せよ。


回答:

\[\begin{align} \mathrm{Var}(x) &= \frac{1}{N}\sum_{i=1}^N\left(x_i - \frac{1}{N}\sum_{i=1}^Nx_i\right) \nonumber\\ &= \frac{1}{N}\sum_{i=1}^Nx_i^2 -2*\frac{1}{N}\sum_{i=1}^Nx_i\frac{1}{N}\sum_{i=1}^Nx_i + \left( \frac{1}{N}\sum_{i=1}^Nx_i\right)^2 \nonumber\\ &= \frac{1}{N}\sum_{i=1}^Nx_i^2 - \bar{x}^2 \end{align}\]


(コーシーシュワルツの不等式): ベクトル\(\boldsymbol{x} = (x_1, ..., x_N)\)、\(\boldsymbol{u} = (1/N, 1/N, ..., 1/N)\)を考える。内積の定義式より、\(|\boldsymbol{x}|^2|\boldsymbol{u}|^2 \ge (\boldsymbol{x}\cdot\boldsymbol{u})^2\)。\((\boldsymbol{x}\cdot\boldsymbol{u})^2 = \bar{x}^2\)、\(|\boldsymbol{x}|^2|\boldsymbol{u}|^2 = \frac{1}{N}\sum_{i=1}^Nx_i^2\)、従って、 \[\begin{equation} \frac{1}{N}\sum_{i=1}^Nx_i^2 \ge \bar{x}^2 \end{equation}\]

データを見る → 散布図

  1種類のデータについて見るためには、ヒストグラムが有効であることは前回の資料にて説明した。2種類のデータの関係性をまず見るために有効な手段が、散布図である。以下では、横軸に数学のテストの得点、縦軸に物理のテストの得点を示した散布図を提示する。

library(readxl) # R特有。便利なライブラリを読み込む。ここでは.xlsxの読み込み関数を含むパッケージをインストール
sampledata <- read_excel("sampledata.xlsx") # データ読み込み
phys = sampledata$phys  # データ読み込み
math = sampledata$math  # データ読み込み
# phys
# hist(phys) # 物理の得点のヒストグラム
par(pty="s")
plot(math, phys) # scatter plot + line, scatter plotのみであれば、普通のplot(math, phys)で十分

上の散布図を見ると、数学の得点が高い人ほど物理の得点が高いように見える。これだけでは定性的な議論にとどまるが、このデータを見た直感や定性的な議論に基づき、定量的な解析を進めていくことは良いアイデアである。

相関係数

  2種類のデータの関係性を記述する指標として、相関係数が有効である。相関係数\(r\)は、\(-1\le r \le 1\)を満たす。-1であれば負の相関、0であれば無相関、1であれば正の相関をもつという。各々のケースは以下の図。

library(MASS)
par(pty="s", mfrow=c(1,3))

# 平均
mu = c(0, 0)

# 分散共分散行列(相関~0)
Sigma = rbind(
  c(1, 0), 
  c(0, 1)
)
dat = mvrnorm(1000,mu,Sigma)
r   = cor(dat)
plot(dat, main = bquote("r = " ~ .(r[1,2])), asp = 1)

# 分散共分散行列(相関 < 0)
Sigma = rbind(
  c(1, -0.5), 
  c(-0.5, 1)
)
dat = mvrnorm(1000,mu,Sigma)
r   = cor(dat)
plot(dat, main = bquote("r = " ~ .(r[1,2])), asp = 1)

# 分散共分散行列(相関 > 0)
Sigma = rbind(
  c(1, 0.5), 
  c(0.5, 1)
)
dat = mvrnorm(1000,mu,Sigma)
r   = cor(dat)
plot(dat, main = bquote("r = " ~ .(r[1,2])), asp = 1)

左図のように、データ1とデータ2の間に特に傾向がなさそうな場合、相関係数は0に近くなる。中図のように、データ1が大きいときにはデータ２は小さくなり、データ1が小さいときにはデータ2が大きくなるとき、相関係数は0より小さくなる。右図のように、データ1が大きいときにはデータ２は大きくなり、データ1が小さいときにはデータ2が小さくなるとき、相関係数は0より大きくなる。

  データ1を\(\boldsymbol{x} = (x_1, ..., x_N)\)、データ2を\(\boldsymbol{y} = (y_1, ..., y_N)\)とする。このとき、\(\boldsymbol{x}\)と\(\boldsymbol{y}\)の相関係数の計算方法は、 \[\begin{equation} r = \frac{\frac{1}{N}\sum_{i=1}^N(x_i - m_x)(y_i - m_y)}{\sigma_x\sigma_y} \end{equation}\] である。ただし、\(m_x\)、\(\sigma_x\)は\(\boldsymbol{x}\)の平均と標準偏差。\(\boldsymbol{y}\)も同様。

  分子の \[\begin{equation} \frac{1}{N}\sum_{i=1}^N(x_i - m_x)(y_i - m_y) \end{equation}\] は共分散と呼ばれ、２変数の関係性を定量化するためにしばしば利用される。ただし、共分散は単位があり(例えばxがcm、yはkgという単位をもつとき、共分散はcm*kgという単位をもつ)、相関係数が単位なし。そのため、共分散はmかcmかに依存して100倍される。相関係数は標準偏差で割り算することで無単位化しているため、cmだろうがmだろうが値は同じ。例えば、下記に、1) もとの相関係数(cor)、2) もとの共分散(cov)、3) データを100倍した相関係数、4) データを100倍した共分散を示す。

# 平均
mu = c(0, 0)

# 分散共分散行列(相関~0)
Sigma = rbind(
  c(1, 0), 
  c(0, 1)
)
dat = mvrnorm(1000,mu,Sigma)

cor(dat[,1], dat[,2])
## [1] 0.02087745
cov(dat[,1], dat[,2])
## [1] 0.02152662
cor(dat[,1]*100, dat[,2]*100)
## [1] 0.02087745
cov(dat[,1]*100, dat[,2]*100)
## [1] 215.2662

共分散についての考察

  共分散 \[\begin{equation} \mathrm{Cov}(x,y) = \frac{1}{N}\sum_{i=1}^N(x_i - m_x)(y_i - m_y) \end{equation}\] について議論を進めてみよう。特に、上式をよく眺めると、 \[\begin{equation} \mathrm{Cov}(x,x) = \mathrm{Var}(x) \end{equation}\] を満たすことがわかる。さらに、演習問題で計算したとおり、 \[\begin{equation} \mathrm{Var}(x) = \frac{1}{N}\sum_{i = 1}^Nx_i^2 - \left( \frac{1}{N}\sum_{i = 1}^Nx_i \right)^2 \end{equation}\] であることから、共分散についても同様の計算が可能であることが期待される。

  共分散の式を展開していくと、 \[\begin{align} \mathrm{Cov}(x,y) &= \frac{1}{N}\sum_{i=1}^Nx_iy_i - 2\frac{1}{N}\sum_{i=1}^Nx_i\frac{1}{N}\sum_{i=1}^Ny_i + \frac{1}{N}\sum_{i=1}^Nx_i\frac{1}{N}\sum_{i=1}^Ny_i \nonumber\\ &= \frac{1}{N}\sum_{i=1}^Nx_iy_i - m_xm_y \end{align}\] を満たす。すなわち、 \[\begin{equation} \mathrm{Cov}(x,y) = \frac{1}{N}\sum_{i=1}^Nx_iy_i - m_xm_y \end{equation}\] を得る。事実、上式から\(\mathrm{Cov}(x,x) = \frac{1}{N}\sum_{i=1}^Nx_i^2 - m_x^2\)として分散に関する式を得られることも留意。

相関係数をさらに議論する。

  この講義の後半で出てくる範囲 + 少しはみ出した範囲の内容であるものの、統計に関する重要なポイントを議論する。 それは、「信頼区間」という考えである。相関係数の信頼区間とは、得られたデータから真の相関係数(この意味も講義の後半で議論する)を推定するとき、とある信頼性をもって相関係数が収まるであろう範囲を意味する。下記ではしばしば議論される95%信頼区間を算出する。すなわち、95%の信頼性をもって真の相関係数が入るであろう区間を算出する。

# 平均
mu = c(0, 0)

# 分散共分散行列(相関~0)
Sigma = rbind(
  c(1, 0.5), 
  c(0.5, 1)
)
dat = mvrnorm(1000,mu,Sigma)

cor.test(dat[,1], dat[,2])
## 
##  Pearson's product-moment correlation
## 
## data:  dat[, 1] and dat[, 2]
## t = 19.718, df = 998, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4833670 0.5726882
## sample estimates:
##       cor 
## 0.5294936

上記が示している中で重要な情報は (データをサンプルし直せば値は変わる)、

95% confidence interval:

sample estimates: cor

という部分。特に、95%の信頼区間に0が入っていないため、この相関係数は0より大きい可能性が高いことが示唆される。 データに最も即している値としては、相関係数はsample estimatesにより得られる値がよい。

  相関係数の信頼区間の計算は下記がわかりやすそう(正規分布の信頼区間に関する知識が必要) https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/PASS/Confidence_Intervals_for_Pearsons_Correlation.pdf

相関係数の気をつけるべき点。

  上記で導出したように、相関係数はあくまで二変数間の角度のみに着目する。 要は、直線の傾きに着目する。すなわち、直線以外のものは定量化できない。

x = seq(0, 2, by = 0.01)
y = 2*x - x^2 + rnorm(length(x), 0, 0.2)
r = cor.test(x, y)
plot(x, y, main = bquote("95% confidence interval: " ~ .(r$conf.int[1]) ~ "~ r  ~" ~ .(r$conf.int[2])), asp = 1)

上記は、\((x-1)^2+y^2 = 1 \ (y \ge 0)\)に従うx, yを多数生成した後、わずかにノイズを加えたものである。すなわち、この2つのデータの間には関係がある。しかし、推定した相関係数の信頼区間は0を含むため、この2つのデータの間には関係がない可能性が高いと推定されてしまっている。

  上記のように直線以外の関係性では相関係数は信頼にたりうる値にならない。ただし、下記のように２変数間の関係性が線形になるような適切な変数を取ることで、相関係数は適切な値になりうる。

x = seq(0, 2, by = 0.01)
z = -(x-1)^2
y = 2*x - x^2 + rnorm(length(x), 0, 0.2)
r = cor.test(z, y)
plot(z, y, main = bquote("95% confidence interval: " ~ .(r$conf.int[1]) ~ "~ r  ~" ~ .(r$conf.int[2])), asp = 1)

他にも相関係数の種類がある

  上記は全てピアソンの相関係数。これらは連続値や離散値を扱うことを想定している。 順位を扱うときは、スピアマンの相関係数やケンドールの相関係数がある。

演習課題

1: 相関係数が\(-1 \le r \le 1\)となることを証明せよ。


2: \(\boldsymbol{x} = (1, 2, 3, 4, 5)\)、\(\boldsymbol{y} = (1, 2, 3, 4, 5)\)とする。


3: 2種類のデータ\(\boldsymbol{x} = (x_1, ..., x_N)\)と\(\boldsymbol{y} = (y_1, ..., y_N)\)を計測したとき、

\[\begin{equation} y_i = ax_i + b \end{equation}\]

として傾きa、切片bの直線をフィットすることを線形回帰と呼ぶ。そして、

\[\begin{equation} E = \frac{1}{2}\sum_{i=1}^N(y_i - ax_i -b)^2 \end{equation}\]

として、データとフィットした直線の間の2乗誤差を最小化する手法を最小二乗法と呼ぶ。

\(\boldsymbol{x}\)と\(\boldsymbol{y}\)の共分散を\(\sigma_{\rm xy}\)、相関係数を\(r\)、\(\boldsymbol{x}\)の分散を\(\sigma_{\rm x}\)^2とする。最小二乗法を利用して推定した\(a\)、\(b\)は、\(a = \frac{\sigma_{\rm xy}}{\sigma_{\rm x}} = r\frac{\sigma_{\rm y}}{\sigma_{\rm x}}\)となることを示せ。

すなわち、相関係数\(r\)は線形回帰した直線の傾きを意味する。