第４章 相関係数
#tidyr はデータフレームの並べ方を展開したり集約したりするパッケージ
library(tidyr)
#データフレームの操作に特化したパッケージ
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union相関係数は２つの変量 x,y の関係を表します。難しく感じますが、具体例で確かめながら進めましょう。
1 復習.
＜定義＞.
偏差は平均からのズレです。
あなたはテストの30点を持って帰ったらお家の人に叱られました。ところが平均が20点のテストだったとしたらどうでしょう。テストが難しかったのであって、自分はそんなに悪くないと主張することができるでしょう。平均を基準として、平均からのズレを計算したものが偏差です。偏差平方和は偏差の二乗の合計.
平均から離れたデータが多いと大きくなります。分散は偏差平方の平均で、平均からの離れ具合の平均です。二乗しているのでルートにしたものが標準偏差です。偏差積は２つのデータがどれくらい似通っているかを表します.
似通っていると (x,y) は直線 y=x に近づきます。
x+yが一定であれば、xyが大きいほど２つの値x,yは似た値になります。
（数学的根拠は \((x+y)/2\geq \sqrt{xy}\) ですが立ち入る必要はありません）偏差積和は偏差積の合計. ２つのデータが直線y=xにどれくらい近いかを表します。偏差積の平均が共分散で、(x,y)がどれくらいy=xに近く分布しているかを表します。
＜関数＞.
関数covは共分散であるが、n-1で割った値であることに注意.
（cov( )にデータフレームを代入すると分散共分散行列を出力する）．関数varは分散であるが、n-1で割った値であることに注意.
記述統計の分散・共分散は n で割って求めるが、推測統計の分散・共分散は n-1 で割る。そのほうが推測の精度が良いからである。Rやexcelの関数は推測統計であり n-1 で割っている。
2 データが２つの点.
2.1
共分散（偏差積の平均）を計算しましょう。
data <- data.frame(x=c(-2,2),y=c(-2,2))
mutate(data,偏差積=x*y)| x | y | 偏差積 |
|---|---|---|
| -2 | -2 | 4 |
| 2 | 2 | 4 |
attach(data)
plot(x,y,asp = 1,xlim=c(-3,3),ylim=c(-3,3))
abline(h=0,v=0,col="lightblue")
abline(a=0,b=1,col="pink")
偏差積の計算.
(x-mean(x))*(y-mean(y))#偏差積## [1] 4 44,4の平均なので共分散は４です。
共分散はプラスに大きくなり、データの位置が y=x
に近いことが分かります。
2.1.1 参考.
関数cov()は推測統計の関数なので、使う場合は注意が必要です。
cov()は n-1=1 で割っているので、修正して２で割る計算をします。
cov(x,y)/2#関数## [1] 4定義から計算した値と一致しています。
2.2
共分散（偏差積の平均）を計算しましょう。
data <- data.frame(x=c(-2,2),y=c(2,-2))
mutate(data,偏差積=x*y)| x | y | 偏差積 |
|---|---|---|
| -2 | 2 | -4 |
| 2 | -2 | -4 |
attach(data)## The following objects are masked from data (pos = 3):
##
## x, yplot(x,y,asp = 1,xlim=c(-3,3),ylim=c(-3,3))
abline(h=0,v=0,col="lightblue")
abline(a=0,b=1,col="pink")
偏差積の計算.
(x-mean(x))*(y-mean(y))## [1] -4 -4-4,-4 の平均なので共分散は -4 です。
共分散はマイナスになり、データの位置が y=x
から離れていることが分かります。
3 データが４つの点.
共分散（偏差積の平均）を計算しましょう。
x <- c(-2,2)
y <- c(-2,2)
data <- crossing(x, y)
mutate(data,偏差積=x*y)| x | y | 偏差積 |
|---|---|---|
| -2 | -2 | 4 |
| -2 | 2 | -4 |
| 2 | -2 | -4 |
| 2 | 2 | 4 |
X <- data$x
Y <- data$y
plot(X,Y,asp = 1,xlim=c(-3,3),ylim=c(-3,3))
abline(h=0,v=0,col="lightblue")
abline(a=0,b=1,col="pink")
(X-mean(X))*(Y-mean(Y))## [1] 4 -4 -4 44 -4 -4 4 の平均なので共分散は０です。まんべんなく散らばっていると共分散は０に近くなります。
4 様々なケース.
4.1 値の大きいデータ.
データの値が大きいと、共分散の値も大きくなります.
x <- c(-5,5)
y <- c(-5,5)
data <- data.frame(x=c(-5,5),y=c(-5,5))
mutate(data,偏差積=x*y)| x | y | 偏差積 |
|---|---|---|
| -5 | -5 | 25 |
| 5 | 5 | 25 |
plot(x,y,asp = 1,xlim=c(-6,6),ylim=c(-6,6))
abline(h=0,v=0,col="lightblue")
abline(a=0,b=1,col="pink")
(x-mean(x))*(y-mean(y))## [1] 25 25値は大きくなりましたが、x,yの見た目の関係は同じです。
値に左右されず関係だけを知りたい場合は、標準化という操作をします。
標準化は平均を引き、標準偏差で割ることです。
すでに平均を引いてあれば、標準偏差で割るだけです。
4.2 標準化して計算.
4.2.1 標準化の例.
x <- c(-2,2)
y <- c(-2,2)
data <- data.frame(x=c(-2,2),y=c(-2,2))
mutate(data,偏差積=x*y)| x | y | 偏差積 |
|---|---|---|
| -2 | -2 | 4 |
| 2 | 2 | 4 |
x,y の標準偏差は 2 なので 2 で割ります。
X <- (x-mean(x))/2
Y <- (x-mean(y))/2
data <- data.frame(x=X,y=Y)
mutate(data,偏差積=x*y)| x | y | 偏差積 |
|---|---|---|
| -1 | -1 | 1 |
| 1 | 1 | 1 |
標準化した偏差積.
(X-mean(X))*(Y-mean(Y))#偏差積## [1] 1 1標準化すると偏差積の平均（共分散）は１になり、x,y
の関係を判断する共通尺度として利用することができるようになります。
標準化して求めた共分散を相関係数といいます。（p21のポイント参照）.
4.2.2
x <- c(-5,5)
y <- c(-5,5)
data <- data.frame(x=c(-5,5),y=c(-5,5))
mutate(data,偏差積=x*y)| x | y | 偏差積 |
|---|---|---|
| -5 | -5 | 25 |
| 5 | 5 | 25 |
X <- (x-mean(x))/5
Y <- (x-mean(y))/5
data <- data.frame(x=X,y=Y)
mutate(data,偏差積=x*y)| x | y | 偏差積 |
|---|---|---|
| -1 | -1 | 1 |
| 1 | 1 | 1 |
標準化した偏差積.
(X-mean(X))*(Y-mean(Y))#偏差積## [1] 1 1標準化して求めた共分散は１になっています。
標準化して求めることで、値の大きさに左右されず判断できるようになりました。
標準化した共分散が１なら、データは y=x
上に並んでいると判断できます。
5 相関係数.
相関係数は、標準化して求めた偏差積の平均（共分散）です。
標準化は、平均を引いて標準偏差で割るので、
xの標準偏差を\(s_x\)
、yの標準偏差を\(s_y\)とすると、
標準化して求めた偏差積は \(\left(\dfrac{x
- \bar{x}}{s_x}-0\right) \left(\dfrac{y -
\bar{y}}{s_y}-0\right)=\dfrac{x - \bar{x}}{s_x}\cdot\dfrac{y -
\bar{y}}{s_y}\)
その平均が相関係数\(r_{xy}\)なので、
\(r_{xy} =\displaystyle
\dfrac{1}{n}\sum_{i=1}^n \dfrac{x_i - \bar{x}}{s_x} \cdot \frac{y_i -
\bar{y}}{s_y} =\left(\dfrac{1}{n}\sum_{i=1}^n (x_i - \bar{x})(y_i -
\bar{y})\right)\dfrac{1}{s_xs_y}\
\).
共分散 \(\displaystyle Cov(x,y)=\dfrac{1}{n}\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})\) とおくと、
- \(r_{xy}=\dfrac{Cov(x,y)}{s_x~s_y}\).
（共分散を標準偏差の積で割った値：p22参考の公式）.
計算しやすい式に書き直します。
\(s_x=\sqrt{\dfrac{1}{n}\sum_{i = 1}^n (x_i - \bar{x})^2}~,~~s_y=\sqrt{\dfrac{1}{n}\sum_{i = 1}^n (y_i - \bar{y})^2}\) なので、
\(s_x~s_y=\dfrac{1}{n}\sqrt{\sum_{i = 1}^n (x_i - \bar{x})^2}\sqrt{\sum_{i = 1}^n (y_i - \bar{y})^2}\).
- \(r_{xy} = \dfrac{\sum_{i = 1}^n (x_i -
\bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i = 1}^n (x_i -
\bar{x})^2}\sqrt{\sum_{i = 1}^n (y_i - \bar{y})^2}}\).
（1/nが約分で消え簡単になっています）.
公式の表現が２通りあり、覚えようとすると大変に感じます。
表による計算で、どちらも同じ値になることを確認しておきましょう。
5.1 確認.
参考サイト.
Acaric 様「文書作成ツール
LaTeXって何？Wordとの違いや便利なポイントを解説！.
LaTeXについて.
- LaTeX（ラテフ）は数式表現に優れた文書ツールであり、理工系の分野を中心に論文執筆の標準となっています。
- 画像に依存しないためホームページの数式の多くはLaTeXで作成されています。（右クリックなどでTeX Commands が表示される）.
ところがプログラムコード内ではLaTeXは使用できないため、以下の表の列名は次のような表現になっています。
\(x-\bar{x}\) =
x-mx \(y-\bar{y}\) =
y-my \((x-\bar{x})(y-\bar{y})\) =
(x-mx)(y-my).
\((x-\bar{x})^2\) =
(x-mx)^2 \((y-\bar{y})^2\) =
(y-my)^2.
もとの表.
| \(x\) | \(y\) | \(x-\bar{x}\) | \(y-\bar{y}\) | \((x-\bar{x})(y-\bar{y})\) | \((x-\bar{x})^2\) | \((y-\bar{y})^2\) |
|---|---|---|---|---|---|---|
| -2 | -2 | |||||
| 2 | 2 | |||||
| 合計 | ||||||
| 平均 | ||||||
共分散/標準偏差：表の平均で計算
（最下列に列平均を表示します）.
x=c(-2,2);y=c(-2,2)
mx=mean(x);my=mean(y)
data <- data.frame(x=x,y=y,A=x-mx,B=y-my,C= (x-mx)*(y-my),D=(x-mx)^2,E=(y-my)^2 )
colnames(data) <- c("x","y","x-mx","y-my", "(x-mx)(y-my)","(x-mx)^2","(y-my)^2)")
rbind(data, apply(data,2,mean))| x | y | x-mx | y-my | (x-mx)(y-my) | (x-mx)^2 | (y-my)^2) |
|---|---|---|---|---|---|---|
| -2 | -2 | -2 | -2 | 4 | 4 | 4 |
| 2 | 2 | 2 | 2 | 4 | 4 | 4 |
| 0 | 0 | 0 | 0 | 4 | 4 | 4 |
\(r_{xy}=\dfrac{Cov(x,y)}{s_x~s_y}=\dfrac{ \dfrac{1}{n}\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\dfrac{1}{n}\sum_{i = 1}^n (x_i - \bar{x})^2}\sqrt{\dfrac{1}{n}\sum_{i = 1}^n (y_i - \bar{y})^2}}=\dfrac{4}{\sqrt{4}\sqrt{4}}=1\).
偏差積和/偏差平方和：表の合計で計算
（最下列に列合計を表示します）.
rbind(data, apply(data,2,sum))| x | y | x-mx | y-my | (x-mx)(y-my) | (x-mx)^2 | (y-my)^2) |
|---|---|---|---|---|---|---|
| -2 | -2 | -2 | -2 | 4 | 4 | 4 |
| 2 | 2 | 2 | 2 | 4 | 4 | 4 |
| 0 | 0 | 0 | 0 | 8 | 8 | 8 |
\(r_{xy} = \dfrac{\sum_{i = 1}^n (x_i -
\bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i = 1}^n (x_i -
\bar{x})^2}\sqrt{\sum_{i = 1}^n (y_i -
\bar{y})^2}}=\dfrac{8}{\sqrt{8}\sqrt{8}}=1\)
［注意］
分母分子の\(\dfrac{1}{n}\)が約分で消え、個数nの影響を受けなくなっていることが分かります。
推測統計では\(\dfrac{1}{n-1}\)で割った値を用いますが、約分で消え影響を受けないので、
記述統計の相関係数＝推測統計の相関係数.
が成立します。
6 まとめ.
共分散は、偏差積の平均.
分散は、偏差の二乗の平均.
相関係数は2変量の関係の度合いを表す.
分布を表すものとしてヒストグラムが有用であるが、分布を簡潔に表現する指標として、平均・分散が有効.
2変量の関係を視覚的に捉えるには散布図が有用であるが、 散布図の状況を簡潔に表現する指標が相関係数.
相関係数は、標準化した偏差積の平均.
小平先生のテキストは意味を大切にする立場から、nで割って共分散、分散を求めてから計算しています. 実用上は、偏差積和と偏差平方和から計算することができます。
練習です。
6.1 例.
p24の下のほうにある例です。
| x | y |
|---|---|
| 2 | 4 |
| 7 | 3 |
| 9 | 8 |
\(r_{xy} = \dfrac{\sum_{i = 1}^n (x_i -
\bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i = 1}^n (x_i -
\bar{x})^2}\sqrt{\sum_{i = 1}^n (y_i -
\bar{y})^2}}\) を用いて相関係数を計算せよ。
（偏差積和/偏差平方和：表の合計で計算）
もとの表.
| \(x\) | \(y\) | \(x-\bar{x}\) | \(y-\bar{y}\) | \((x-\bar{x})(y-\bar{y})\) | \((x-\bar{x})^2\) | \((y-\bar{y})^2\) |
|---|---|---|---|---|---|---|
| 2 | 4 | |||||
| 7 | 3 | |||||
| 9 | 8 | |||||
| 合計 | ||||||
| 平均 | ||||||
解答.
（最下列に列合計を表示します）.
x=c(2,7,9);y=c(4,3,8)
mx=mean(x);my=mean(y)
data <- data.frame(x=x,y=y,A=x-mx,B=y-my,C= (x-mx)*(y-my),D=(x-mx)^2,E=(y-my)^2 )
colnames(data) <- c("x","y","x-mx","y-my", "(x-mx)(y-my)","(x-mx)^2","(y-my)^2)")
rbind(data, apply(data,2,sum))| x | y | x-mx | y-my | (x-mx)(y-my) | (x-mx)^2 | (y-my)^2) |
|---|---|---|---|---|---|---|
| 2 | 4 | -4 | -1 | 4 | 16 | 1 |
| 7 | 3 | 1 | -2 | -2 | 1 | 4 |
| 9 | 8 | 3 | 3 | 9 | 9 | 9 |
| 18 | 15 | 0 | 0 | 11 | 26 | 14 |
\(r_{xy} = \dfrac{\sum_{i = 1}^n (x_i - \bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i = 1}^n (x_i - \bar{x})^2}\sqrt{\sum_{i = 1}^n (y_i - \bar{y})^2}}=\dfrac{11}{\sqrt{26}\sqrt{14}}=0.576\)
別解.
テキストp25のⅠ では、
\(r_{xy}=\dfrac{Cov(x,y)}{s_x~s_y}=\dfrac{
\frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})(y_i -
\bar{y})}{\sqrt{\frac{1}{n}\sum_{i = 1}^n (x_i -
\bar{x})^2}\sqrt{\frac{1}{n}\sum_{i = 1}^n (y_i -
\bar{y})^2}}\) で計算している.
即ち、共分散を標準偏差の積で割った値で計算している.
同じ値になることを確かめておこう。
（最下列に列平均を表示します）.
x=c(2,7,9);y=c(4,3,8)
mx=mean(x);my=mean(y)
data <- data.frame(x=x,y=y,A=x-mx,B=y-my,C= (x-mx)*(y-my),D=(x-mx)^2,E=(y-my)^2 )
colnames(data) <- c("x","y","x-mx","y-my", "(x-mx)(y-my)","(x-mx)^2","(y-my)^2)")
rbind(data, apply(data,2,mean))| x | y | x-mx | y-my | (x-mx)(y-my) | (x-mx)^2 | (y-my)^2) |
|---|---|---|---|---|---|---|
| 2 | 4 | -4 | -1 | 4.000000 | 16.000000 | 1.000000 |
| 7 | 3 | 1 | -2 | -2.000000 | 1.000000 | 4.000000 |
| 9 | 8 | 3 | 3 | 9.000000 | 9.000000 | 9.000000 |
| 6 | 5 | 0 | 0 | 3.666667 | 8.666667 | 4.666667 |
\(r_{xy}=\dfrac{Cov(x,y)}{s_x~s_y}=\dfrac{
\dfrac{1}{n}\sum_{i=1}^n (x_i - \bar{x})(y_i -
\bar{y})}{\sqrt{\dfrac{1}{n}\sum_{i = 1}^n (x_i -
\bar{x})^2}\sqrt{\dfrac{1}{n}\sum_{i = 1}^n (y_i -
\bar{y})^2}}=\dfrac{3.667}{\sqrt{8.667}\sqrt{4.667}}=0.576\).
共分散/標準偏差：表の平均で計算すると小数計算が面倒になるので、偏差積和/偏差平方和：表の合計で計算する方法がよく使われる訳です。
6.1.1 注意.
Rの関数は推測統計なので
(n-1)で割っています。記述統計の分散にするには、\(\dfrac{n-1}{n}\)を乗ずる️必要があります。
var(data$x,data$y)*2/3## [1] 3.666667しかし、相関係数は修正の必要はありません。
記述統計の相関係数も、推測統計の相関係数も同じです。
なぜでしょうか？
相関係数の商の式において、nもn-1も約分によって消えているからです。
cor(data$x,data$y)## [1] 0.5765567コンピュータは n-1
で割って計算していますが、修正しなくても記述統計の相関係数と一致しています。
7 例題４.
もとの表.
| \(x\) | \(y\) | \(x-\bar{x}\) | \(y-\bar{y}\) | \((x-\bar{x})(y-\bar{y})\) | \((x-\bar{x})^2\) | \((y-\bar{y})^2\) |
|---|---|---|---|---|---|---|
| 165 | 160 | |||||
| 170 | 162 | |||||
| 168 | 166 | |||||
| 172 | 164 | |||||
| 169 | 160 | |||||
| 164 | 156 | |||||
| 162 | 164 | |||||
| 172 | 162 | |||||
| 170 | 158 | |||||
| 168 | 158 | |||||
| 合計 | ||||||
| 平均 | ||||||
散布図.
x <- c(165,170,168,172,169,164,162,172,170,168)
y <- c(160,162,166,164,160,156,164,162,158,158)
data <- data.frame(x=x,y=y)
plot(x,y,asp = 1,xlim=c(160,175),ylim=c(150,170))
abline(h=0,v=0,col="lightblue")
\(r_{xy} = \dfrac{\sum_{i = 1}^n (x_i -
\bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i = 1}^n (x_i -
\bar{x})^2}\sqrt{\sum_{i = 1}^n (y_i -
\bar{y})^2}}\) を用いて相関係数を計算せよ。
（最下列に列合計を表示します）.
mx=mean(x);my=mean(y)
data <- data.frame(x=x,y=y,A=x-mx,B=y-my,C= (x-mx)*(y-my),D=(x-mx)^2,E=(y-my)^2 )
colnames(data) <- c("x","y","x-mx","y-my", "(x-mx)(y-my)","(x-mx)^2","(y-my)^2)")
rbind(data, apply(data,2,sum))| x | y | x-mx | y-my | (x-mx)(y-my) | (x-mx)^2 | (y-my)^2) |
|---|---|---|---|---|---|---|
| 165 | 160 | -3 | -1 | 3 | 9 | 1 |
| 170 | 162 | 2 | 1 | 2 | 4 | 1 |
| 168 | 166 | 0 | 5 | 0 | 0 | 25 |
| 172 | 164 | 4 | 3 | 12 | 16 | 9 |
| 169 | 160 | 1 | -1 | -1 | 1 | 1 |
| 164 | 156 | -4 | -5 | 20 | 16 | 25 |
| 162 | 164 | -6 | 3 | -18 | 36 | 9 |
| 172 | 162 | 4 | 1 | 4 | 16 | 1 |
| 170 | 158 | 2 | -3 | -6 | 4 | 9 |
| 168 | 158 | 0 | -3 | 0 | 0 | 9 |
| 1680 | 1610 | 0 | 0 | 16 | 102 | 90 |
\(r_{xy} = \dfrac{\sum_{i = 1}^n (x_i - \bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i = 1}^n (x_i - \bar{x})^2}\sqrt{\sum_{i = 1}^n (y_i - \bar{y})^2}}=\dfrac{16}{\sqrt{102}\sqrt{90}}=0.167\)
別解.
\(r_{xy}=\dfrac{Cov(x,y)}{s_x~s_y}=\dfrac{
\frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})(y_i -
\bar{y})}{\sqrt{\frac{1}{n}\sum_{i = 1}^n (x_i -
\bar{x})^2}\sqrt{\frac{1}{n}\sum_{i = 1}^n (y_i -
\bar{y})^2}}\)（共分散を標準偏差の積で割って計算せよ）.
（最下列に列平均を表示します）.
mx=mean(x);my=mean(y)
data <- data.frame(x=x,y=y,A=x-mx,B=y-my,C= (x-mx)*(y-my),D=(x-mx)^2,E=(y-my)^2 )
colnames(data) <- c("x","y","x-mx","y-my", "(x-mx)(y-my)","(x-mx)^2","(y-my)^2)")
rbind(data, apply(data,2,mean))| x | y | x-mx | y-my | (x-mx)(y-my) | (x-mx)^2 | (y-my)^2) |
|---|---|---|---|---|---|---|
| 165 | 160 | -3 | -1 | 3.0 | 9.0 | 1 |
| 170 | 162 | 2 | 1 | 2.0 | 4.0 | 1 |
| 168 | 166 | 0 | 5 | 0.0 | 0.0 | 25 |
| 172 | 164 | 4 | 3 | 12.0 | 16.0 | 9 |
| 169 | 160 | 1 | -1 | -1.0 | 1.0 | 1 |
| 164 | 156 | -4 | -5 | 20.0 | 16.0 | 25 |
| 162 | 164 | -6 | 3 | -18.0 | 36.0 | 9 |
| 172 | 162 | 4 | 1 | 4.0 | 16.0 | 1 |
| 170 | 158 | 2 | -3 | -6.0 | 4.0 | 9 |
| 168 | 158 | 0 | -3 | 0.0 | 0.0 | 9 |
| 168 | 161 | 0 | 0 | 1.6 | 10.2 | 9 |
\(r_{xy}=\dfrac{Cov(x,y)}{s_x~s_y}=\dfrac{ \dfrac{1}{n}\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\dfrac{1}{n}\sum_{i = 1}^n (x_i - \bar{x})^2}\sqrt{\dfrac{1}{n}\sum_{i = 1}^n (y_i - \bar{y})^2}}=\dfrac{1.6}{\sqrt{10.2}\sqrt{9}}=0.167\).
8 演習４.
| \(x\) | \(y\) | \(x-\bar{x}\) | \(y-\bar{y}\) | \((x-\bar{x})(y-\bar{y})\) | \((x-\bar{x})^2\) | \((y-\bar{y})^2\) |
|---|---|---|---|---|---|---|
| 2 | 7 | |||||
| 4 | 8 | |||||
| 9 | 10 | |||||
| 1 | 3 | |||||
| 6 | 11 | |||||
| 8 | 9 | |||||
| 合計 | ||||||
| 平均 | ||||||
散布図.
x <- c(2,4,9,1,6,8)
y <- c(7,8,10,3,11,9)
data <- data.frame(x=x,y=y)plot(x,y,asp = 1,xlim=c(0,12),ylim=c(0,12))
abline(h=0,v=0,col="lightblue")
\(r_{xy} = \dfrac{\sum_{i = 1}^n (x_i -
\bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i = 1}^n (x_i -
\bar{x})^2}\sqrt{\sum_{i = 1}^n (y_i -
\bar{y})^2}}\) を用いて相関係数を計算せよ。
（最下列に列合計を表示します）.
mx=mean(x);my=mean(y)
data <- data.frame(x=x,y=y,A=x-mx,B=y-my,C= (x-mx)*(y-my),D=(x-mx)^2,E=(y-my)^2 )
colnames(data) <- c("x","y","x-mx","y-my", "(x-mx)(y-my)","(x-mx)^2","(y-my)^2)")
rbind(data, apply(data,2,sum))| x | y | x-mx | y-my | (x-mx)(y-my) | (x-mx)^2 | (y-my)^2) |
|---|---|---|---|---|---|---|
| 2 | 7 | -3 | -1 | 3 | 9 | 1 |
| 4 | 8 | -1 | 0 | 0 | 1 | 0 |
| 9 | 10 | 4 | 2 | 8 | 16 | 4 |
| 1 | 3 | -4 | -5 | 20 | 16 | 25 |
| 6 | 11 | 1 | 3 | 3 | 1 | 9 |
| 8 | 9 | 3 | 1 | 3 | 9 | 1 |
| 30 | 48 | 0 | 0 | 37 | 52 | 40 |
\(r_{xy} = \dfrac{\sum_{i = 1}^n (x_i - \bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i = 1}^n (x_i - \bar{x})^2}\sqrt{\sum_{i = 1}^n (y_i - \bar{y})^2}}=\dfrac{37}{\sqrt{52}\sqrt{40}}=0.811\)
別解. 共分散を標準偏差の積で割った値で計算しなさい.
（最下列に列平均を表示します）.
mx=mean(x);my=mean(y)
data <- data.frame(x=x,y=y,A=x-mx,B=y-my,C= (x-mx)*(y-my),D=(x-mx)^2,E=(y-my)^2 )
colnames(data) <- c("x","y","x-mx","y-my", "(x-mx)(y-my)","(x-mx)^2","(y-my)^2)")
rbind(data, apply(data,2,mean))| x | y | x-mx | y-my | (x-mx)(y-my) | (x-mx)^2 | (y-my)^2) |
|---|---|---|---|---|---|---|
| 2 | 7 | -3 | -1 | 3.000000 | 9.000000 | 1.000000 |
| 4 | 8 | -1 | 0 | 0.000000 | 1.000000 | 0.000000 |
| 9 | 10 | 4 | 2 | 8.000000 | 16.000000 | 4.000000 |
| 1 | 3 | -4 | -5 | 20.000000 | 16.000000 | 25.000000 |
| 6 | 11 | 1 | 3 | 3.000000 | 1.000000 | 9.000000 |
| 8 | 9 | 3 | 1 | 3.000000 | 9.000000 | 1.000000 |
| 5 | 8 | 0 | 0 | 6.166667 | 8.666667 | 6.666667 |
\(r_{xy}=\dfrac{Cov(x,y)}{s_x~s_y}=\dfrac{
\dfrac{1}{n}\sum_{i=1}^n (x_i - \bar{x})(y_i -
\bar{y})}{\sqrt{\dfrac{1}{n}\sum_{i = 1}^n (x_i -
\bar{x})^2}\sqrt{\dfrac{1}{n}\sum_{i = 1}^n (y_i -
\bar{y})^2}}=\dfrac{6.166}{\sqrt{8.66}\sqrt{6.66}}=0.811\).
9 研究.
相関のある2変量乱数は2変量正規乱数rmvnormで出力できる。
回帰と相関の違いについては市原清志先生の「バイオサイエンスの統計学」第
７章に詳しい。
以下の説明は考査の出題には入りませんが、図でイメージを見ておくとよいでしょう。
| 回帰 | 相関 | |
|---|---|---|
| 誤 差 | y方向のみ | x,y両方が確率変数 |
| 信頼域 | yの信頼区間 | x,yの分布の中心からの等確率距離 |
| 目 的 | どのように直線に関係づけられるか | どの程度直線的か |
| 単 位 | 単位に依存 | 単位に依存しない |
ここでは、相関係数によって散布図がどのようになるかを Rの乱数でシミュレーションする。イメージを掴んでおこう。
黒線は回帰直線、青の縦線は x=1.5 の直線を表します。
どんなことが読み取れるでしょうか？
library(mvtnorm)
library(scatterplot3d)9.1 標準偏差 \(s_x=1,s_y=1\).
9.1.1 相関係数 r=0.
sigma.zero <- matrix(c(1,0,0,1), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma.zero)
par(mfrow = c(1,1))
plot(x,xlim=c(-4,4),ylim=c(-4,4),asp=1,col="grey") 
9.1.2 相関係数 r=0.5.
sigma <- matrix(c(1,0.5,0.5,1), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)
data <- data.frame(X=x[,1],Y=x[,2])
model <- lm(Y ~ X, data = data)
plot(x,asp=1,col="grey")
abline(model)
abline(v=1.5,col="blue")
model$coef## (Intercept) X
## 0.0014491 0.4905962回帰直線は、\(y=0.5x\)
標準化されていれば、回帰係数と相関係数は一致します。
回帰直線はxを固定したときの推定値を与えます。
例えていうと、相関係数はサツマイモの太さのようなもので、右斜めにしてサツマイモを切ったときに、切り口の中点を繋いだ線が回帰直線です。
9.2 標準偏差 \(s_x=1,s_y=2\).
9.2.1 相関係数 r=0.
sigma <- matrix(c(1,0,0,4), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)
plot(x,asp=1,col="grey")
#var(x)
#分散共分散行列は、matrix(c(1,0,0,4), ncol=2) になっています。 9.2.2 相関係数 r=05.
sigma <- matrix(c(1,1,1,4), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)
data <- data.frame(X=x[,1],Y=x[,2])
model <- lm(Y ~ X, data = data)
plot(x,asp=1,col="grey")
abline(model)
abline(v=2,col="blue")
model$coef## (Intercept) X
## 2.847874e-05 1.016359e+00回帰直線は、\(y=1\cdot x\)
一般に回帰係数は、
標準偏差 \(s_x,s_y\)、相関係数 \(r_{xy}\) とすると \(\displaystyle
b=\dfrac{s_y}{s_x}r_{xy}\) であるので、
本問では、 \(b=2\cdot 0.5=1\) となっています。
9.2.3 相関係数 r=1.
sigma <- matrix(c(1,2,2,4), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)
data <- data.frame(X=x[,1],Y=x[,2])
model <- lm(Y ~ X, data = data)
plot(x,asp=1,col="grey")
abline(model)
abline(v=2,col="blue")
model$coef## (Intercept) X
## -1.90586e-17 2.00000e+00回帰直線は、\(y=2\cdot x\)
9.2.4 相関係数 r=0.4.
sigma <- matrix(c(1,0.8,0.8,4), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)
data <- data.frame(X=x[,1],Y=x[,2])
model <- lm(Y ~ X, data = data)
plot(x,asp=1,col="grey")
abline(model)
abline(v=2,col="blue")
model$coef## (Intercept) X
## -0.005835981 0.783000493回帰直線は、\(y=0.8\cdot x\).
回帰分析については、看護統計の看護応用
の＜回帰と相関＞の項目に関連記事があります。
（2変量を扱うので少し難しくなります）.