（2変量を扱うので少し難しくなります）.

相関係数は、2変量 x,y の直線関係の強さの指標である。
では、相関の検定とはどういうことであろうか？
前提として 2変量正規分布をする母集団があり、得られた相関係数は標本から算出したものである。
p228のKey Point によれば、母集団の相関係数が０で無相関のとき、標本相関係数は平均0の正規分布に従う。
つまり、標本相関係数は 0に近い値を変動する。
従って、標本相関係数が正だからといって、母集団に正の相関があると結論するわけにはいかない。
標本相関係数から母相関係数を推定し、本当に相関があるのかどうかを検定する必要がある。

1 標本数100.

1.1 母相関係数 0 のとき.

p225によれば、母相関係数 \(\rho=0\) のもとでは、

n個の標本相関係数 rの標準誤差は　\(s_r=\sqrt{\dfrac{1-r^2}{n-2}}\).
標準誤差で標準化したt値は　　　　\(t=\dfrac{r}{s_r} \sim t(n-2)\).

これらのことが成立するという。
（母相関係数が 0 でない場合は成立しないので、後述のFisher Z 変換が必要になる）

母集団：相関０の2変量正規分布.

sigma <- matrix(c(1,0,0,1), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)  
data <- data.frame(X=x[,1],Y=x[,2])

大きさ100の標本を抽出し相関係数を算出.

R=c()
for (i in 1:1000){
  S <- sample_n(tbl = data, size = 100)
  r <- cor(S$X,S$Y)
  R <- c(R,r)
}
sd(R)#相関係数の標準誤差

## [1] 0.09882158

sqrt((1-0^2)/(100-2))#相関係数の標準誤差の理論値

## [1] 0.1010153

シミュレーションによる標準偏差と \(s_r=\sqrt{\dfrac{1-r^2}{n-2}}\)　は一致している。

hist(R,breaks="sturges",prob=T)
curve(dnorm(x,0,sqrt((1-0^2)/(100-2))), col="red",add=T)

ヒストグラムも理論曲線と一致している（理論曲線は正規分布）。

1.1.1 正規QQプロット.

標本相関係数が正規分布に従うかどうか、正規QQプロットで調べることもできる。
点が直線上に分布していれば正規分布と判断できる。

qqnorm(R,  main = "QQ Plot by qqnorm function")
qqline(R, col="red")

点は直線上に配置されており、正規分布と判断できる。

1.2 母相関係数 0.75のとき

母集団：相関.75の2変量正規分布.
母相関係数が 0 でないので相関係数は正規分布にならない。
シミュレーションではどうだろうか。

sigma <- matrix(c(1,0.75,0.75,1), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)  
data <- data.frame(X=x[,1],Y=x[,2])

R=c()
for (i in 1:1000){
  S <- sample_n(tbl = data, size = 100)
  r <- cor(S$X,S$Y)
  R <- c(R,r)
}
sd(R)

## [1] 0.04562831

sqrt((1-0.75^2)/(100-2))

## [1] 0.06681531

標準偏差について、
　（理論値）＞（実現値）　　となっているようである。
　

hist(R,breaks="sturges",prob=T)
curve(dnorm(x,0.75,sqrt((1-0.75^2)/(100-2))), col="red",add=T)

ヒストグラムも一致せず、偏りがあるように見える。

qqnorm(R,  main = "QQ Plot by qqnorm function")
qqline(R, col="red")

正規QQプロットも直線からズレがあり、正規分布とは言えないようである。

2 標本数10.

2.1 母相関係数 0 のとき.

母集団：相関０の2変量正規分布.

sigma <- matrix(c(1,0,0,1), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)  
data <- data.frame(X=x[,1],Y=x[,2])

標本相関係数とt値を１万個出力し、分布をヒストグラムに表示する。

R=c()
ST=c()
for (i in 1:1000){
  S <- sample_n(tbl = data, size = 10)
  r <- cor(S$X,S$Y)
  st <- r/(sqrt((1-r^2)/(10-2)))
  R <- c(R,r)
  ST <- c(ST,st)
}
sd(R)

## [1] 0.3262019

sqrt((1-0^2)/(10-2))

## [1] 0.3535534

シミュレーションの標準偏差は、理論値とほぼ一致しているようである。

hist(R,breaks="sturges",prob=T)

hist(ST,breaks="sturges",prob=T)
curve(dt(x, 8), col="red",add=T)

自由度８のt分布に近似的に一致している。

qqnorm(R,  main = "QQ Plot by qqnorm function")
qqline(R, col="red")

母相関係数が 0 なので相関係数は正規分布に従うはずであるが、正規QQプロットは直線から少しズレがあるように見える。（原因は不明）

2.2 母相関係数 0.75

母集団：相関係数 0.75 の2変量正規分布.

sigma <- matrix(c(1,0.75,0.75,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)

R=c()
ST=c()
for (i in 1:1000){
  S <- sample_n(tbl = data, size = 10)
  r <- cor(S$X,S$Y)
  st <- (r-0.75)/(sqrt((1-r^2)/(10-2)))
  R <- c(R,r)
  ST <- c(ST,st)
}
mean(R)

## [1] 0.7303601

sd(R)

## [1] 0.1712686

sqrt((1-0.75^2)/(10-2))

## [1] 0.2338536

標準偏差もズレがあるようである。
ヒストグラムを描くと、分布に偏りがあることが分かる。

hist(R,breaks="sturges",prob=T)

hist(ST,breaks="sturges",prob=T)
curve(dt(x, 8), col="red",add=T)

t値の分布も理論通りとはならないようである。

qqnorm(R,  main = "QQ Plot by qqnorm function")
qqline(R, col="red")

正規QQプロットは直線から少しズレがあるように見える。
母相関係数が 0.75 のとき、標本相関係数は正規分布に従うとは言えないようである。

3 FisherのZ変換.

母相関係数が 0 でない場合は. 　　\(Z=\dfrac{1}{2}\log\left(\dfrac{1+r}{1-r}\right)\)　　によって変換すると正規分布に近くなる。

この場合、母相関係数を \(\rho\) として.　

平均値　\(\mu_Z=\dfrac{1}{2}\log\left(\dfrac{1+\rho}{1-\rho}\right)\).
標準偏差　\(\dfrac{1}{\sqrt{n-3}}\)　　　となる。

母集団：相関係数 0.75 の2変量正規分布.

ZR=c()
for (i in 1:1000){
  S <- sample_n(tbl = data, size = 10)
  r <- cor(S$X,S$Y)
  zr <- 1/2*log((1+r)/(1-r))
  ZR <- c(ZR,zr)
}
hist(ZR,breaks="sturges",prob=T)
curve(dnorm(x,1/2*log((1+0.75)/(1-0.75)),1/sqrt(10-3)), col="red",add=T)

qqnorm(ZR,  main = "QQ Plot by qqnorm function")
qqline(ZR, col="red")

正規分布と言って良さそうである。

3.1 例7_4(p226)

dat<- data.frame(X=c(8,7,6,6,6,5,4,4,3,2),Y=c(6,5,7,6,4,5,6,3,4,3))

fit = lm( formula = Y ~ X, data = dat )
summary(fit)

## 
## Call:
## lm(formula = Y ~ X, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.39806 -0.69660 -0.08252  0.55340  1.68932 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   2.5728     1.1100   2.318   0.0491 *
## X             0.4563     0.2058   2.218   0.0574 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.144 on 8 degrees of freedom
## Multiple R-squared:  0.3807, Adjusted R-squared:  0.3033 
## F-statistic: 4.918 on 1 and 8 DF,  p-value: 0.05739

相関係数の検定の手計算はテキストに詳しく書かれているので割愛する。
ここではRの関数を用いた算出法を示す。
回帰と相関の両方について整理しておく。

3.1.1 回帰係数の信頼区間.

関数confint().

confint(fit)

##                   2.5 %    97.5 %
## (Intercept)  0.01320958 5.1324215
## X           -0.01817883 0.9308002

回帰係数の信頼区間は、confint()関数で得られる。

サンプル数は10なので、Z変換を適用する。

p227に従い計算する.

r <- cor(dat$X,dat$Y);r

## [1] 0.6170161

Zr <- 1/2*log((1+r)/(1-r));Zr

## [1] 0.7201725

Rの関数を紹介しよう。

参考サイト　Rを使ったデータ処理.

3.1.2 相関係数の信頼区間.

cor.test関数を用いると、相関係数の信頼区間が出力される。

result <- cor.test(dat$X,dat$Y)
result

## 
##  Pearson's product-moment correlation
## 
## data:  dat$X and dat$Y
## t = 2.2177, df = 8, p-value = 0.05739
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.02062134  0.89784068
## sample estimates:
##       cor 
## 0.6170161

result$conf.int#相関係数の95%信頼区間

## [1] -0.02062134  0.89784068
## attr(,"conf.level")
## [1] 0.95

4 Advanced.

4.1 回帰の残差分散.

2変量正規分布の分散との関係.

回帰の残差分散は、\(s^2=(1-\rho^2)\sigma_y^2\).
残差分散はxに依存しません。
\(\displaystyle \beta=\frac{\sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) \left( y_{i}-\overline{y}\right) }{ \sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) ^{2}}=\dfrac{s_{xy}}{s_x^2}=\dfrac{s_{xy}}{s_xs_y}\cdot\dfrac{s_y}{s_x}=r_{xy}\cdot\dfrac{s_y}{s_x}\).
\(V(\beta)=\dfrac{\sigma^2}{T_{xx}}=\dfrac{\sigma^2}{ns_{x}^2}=\dfrac{(1-\rho^2)\sigma_y^2}{ns_{x}^2}\).

回帰係数の分散は、回帰の残差分散とxの偏差平方和で決まる。
たった一回のサンプリングで、回帰係数の推定値だけでなく、その分散までわかってしまうということになる。

相関係数についても同様のことが成立している。

4.2 相関係数の分散.

母相関係数\(\rho=0\) のもとでは、
n個の標本相関係数 rの標準誤差は　\(s_r=\sqrt{\dfrac{1-r^2}{n-2}}\)
標準誤差で標準化したt値は　　　　\(t=\dfrac{r}{s_r} \sim t(n-2)\).
たった一回のサンプリングで、相関係数の推定値だけでなく、その分散までわかってしまう。
数理的根拠については、
Znn様「標本相関係数の密度関数の導出」にありましたが難しいです。 .
猫も杓子も構造化様「相関係数の標本分布」にはシミュレーションが示されています。

数理的導出に突っ込みたくなりますが、シミュレーションで我慢しておいた方が良さそうです。

4.2.1 Z変換を用いた場合.

p229によれば.　\(Z=\dfrac{1}{2}\log\left(\dfrac{1+r}{1-r}\right)\).

平均値　\(\mu_Z=\dfrac{1}{2}\log\left(\dfrac{1+\rho}{1-\rho}\right)\).
標準偏差　\(\dfrac{1}{\sqrt{n-3}}\).

5 相関係数の差の検定.

２標本が共通の2変量正規母集団から得られたと仮定（母相関係数：\(\rho_1=\rho_2\) ）.

標本相関係数の差：\(r_1-r_2\) を統計量として検定する。

母相関係数を０としていないので、
r を Z変換したうえで、 \(Z_1-Z_2\)　を検定する。
　　　\((Z_1-Z_2) \sim N \left(0~,~~\dfrac{1}{\sqrt{\dfrac{1}{n_1-3}+\dfrac{1}{n_2-3}}}\right)\).

検定は標準正規分布を用いた検定であり難しくはない。

Rのパッケージpsychのr.test関数を用いて出力しておく。

5.1 例7-4(p226).

library(psych)
r.test(n=10, r12=0.617)

## Correlation tests 
## Call:r.test(n = 10, r12 = 0.617)
## Test of significance of a  correlation 
##  t value 2.22    with probability < 0.057
##  and confidence interval  -0.02 0.9

5.2 例7-5(p233).

r.test(n=35, r12=0.68, r34=0.52, n2=40)

## Correlation tests 
## Call:r.test(n = 35, r12 = 0.68, r34 = 0.52, n2 = 40)
## Test of difference between two independent correlations 
##  z value 1.05    with probability  0.3