第15章 無相関の検定

参考サイト.

• 当サイト 看護応用「相関係数」Advanced.

• 相関係数は意外と難しい.
  

• 相関係数をカンゼンニリカイする.

• 相関係数の検定の仕方.

• 猫も杓子も構造化 様「相関係数の標本分布」.
  （数理的根拠やシミュレーションも示されています）

標準誤差.
標本から測定された統計量の標準偏差であり、標本統計量の精度を表す。 その中でも特に標本平均の標準偏差を指すことが多い。

1 はじめに.

平均は確率変数である.
1変量母集団から抽出された標本の標本平均は確率変数である。
標本平均は正規分布に分布収束する。
これにより、正規分布に基づく推定や検定が可能となる。

相関係数は確率変数である.
調査によって得られる2変量データは、2変量母集団から抽出された標本である。
標本相関係数は確率変数である。
では、

• 標本相関係数はどのような分布に従うのであろうか？
  
• 標本相関係数は正規分布で近似できるだろうか？

1.1 類比.

標本平均と標本相関係数を類比によって考えてみよう。

1.1.1 標本平均.

　　　\(\bar{X}\sim N\left(\mu,\dfrac{\sigma^2}{n}\right)\)　　　\(Z=\dfrac{\bar{X}-\mu}{\frac{\sigma}{\sqrt{n}}}\sim N(0,1)\).

　母分散が未知のとき、
　\(\bar{X}\)の母標準偏差\(\dfrac{\sigma}{\sqrt{n}}\)を推定値である標準誤差\(\dfrac{s}{\sqrt{n}}\) で代用すると.
　　　\(t=\dfrac{\overline{X}-\mu}{\left. s \right/ \sqrt{n}} \sim t(n-1)\).

1.1.2 標本相関係数.

南風原朝和先生「心理統計学の基礎」p114では、
母相関係数 \(\rho\) の2変量正規分布から大きさNの標本を抽出するとき、
標本相関係数 \(r\) の期待値と標準誤差は近似的に次の値になる.

• 期待値：\(\displaystyle{ \mu_r = \rho - \frac{ \rho (1- \rho^2)}{2N} }\).
  （第二項はバイアスで、標本相関係数の値は少し小さめに出る）.
• 標準誤差：\(\displaystyle{ \sigma_r = \frac{ 1- \rho^2}{\sqrt{N}} }\).

帰無仮説 \(\rho=0\)の仮定のもとでは、

• 期待値：\(\displaystyle{ \mu_r = 0}\).
• 標準誤差：\(\displaystyle{ s_r = \sqrt{\frac{ 1- r^2}{N-2} }}\)　（自由度は\(N-2\)）.

従って、t値は. 　　　\(t=\dfrac{r-0}{\frac{ \sqrt{1- r^2}}{\sqrt{N-2}}}=\dfrac{r\sqrt{N-2}}{ \sqrt{1- r^2}}=r\sqrt{\dfrac{N-2}{1-r^2}}\)

数理的な導出は上記の参考サイトに説明があるが、かなり難しい。
シミュレーションで体験的に理解する「ひょっこり方式」の出番である。

• 100個の標本相関係数 r の期待値 .
• 100個のr の標準誤差の真値　\(\displaystyle{ \sigma_r = \frac{ 1- \rho^2}{\sqrt{N}} }\).
• 100個のr の標準誤差の推定値　\(\displaystyle{ s_r = \sqrt{\frac{ 1- r^2}{N-2} }}\).

これらをシミュレーションで確かめる。

母集団は二変量正規分布に従い、大きさは10000とする。.

2 シミュレーション

2.1 母相関係数０のとき.

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])
cor(data$X,data$Y)

## [1] 0.02179983

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

R=c()
K=c()
TT=c()
for (i in 1:1000){
  S <- sample_n(tbl = data, size = 10)
  r <- cor(S$X,S$Y)
  R <- c(R,r)
  k <- sqrt((1-r^2)/(10-2))
  K <- c(K,k)
  t <- r/k
  TT <- c(TT,t)
}
hist(R)

norm<-fitdist(R,"norm")
plot(norm)

summary(norm)

## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters : 
##        estimate  Std. Error
## mean 0.01615226 0.010219242
## sd   0.32316079 0.007225784
## Loglikelihood:  -289.3333   AIC:  582.6665   BIC:  592.482 
## Correlation matrix:
##              mean           sd
## mean 1.000000e+00 1.678973e-11
## sd   1.678973e-11 1.000000e+00

公式による値との違い.

#大きさ10の標本による標本標準偏差
mean(R)

## [1] 0.01615226

sd(R)#相関係数の標準誤差

## [1] 0.3233225

#相関係数の標準誤差の公式による値
(1-0^2)/sqrt((10))

## [1] 0.3162278

#sqrt((1-0^2)/(10-2))
#(1-0^2)/sqrt((10-2))
mean(K)

## [1] 0.333522

hist(K,prob=T)

\(\displaystyle{ s_r = \sqrt{\frac{ 1- r^2}{N-2} }}\) の値は、よく近似できているといえる。.

t分布への当てはまり.

hist(TT,prob=T)
curve(dt(x,8), col="blue",add=T)
curve(dnorm(x,0,1), col="red",add=T)

t分布のほうが、当てはまりが良いようである。

参考.
　　　\({f(r)=\dfrac{(1-r^2)^\frac{n-4}{2}}{B(\frac{1}{2},\frac{1}{2}(n-2))} }\)

2.2 母相関係数 0.75 のとき.

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])
cor(data$X,data$Y)

## [1] 0.7433995

R=c()
K=c()
TT=c()
for (i in 1:1000){
  S <- sample_n(tbl = data, size = 10)
  r <- cor(S$X,S$Y)
  R <- c(R,r)
  k <- sqrt((1-r^2)/(10-2))
  K <- c(K,k)
  t <- r/k
  TT <- c(TT,t)
}
hist(R)

norm<-fitdist(R,"norm")
plot(norm)

summary(norm)

## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters : 
##       estimate  Std. Error
## mean 0.7209155 0.005506537
## sd   0.1741320 0.003893132
## Loglikelihood:  329.0031   AIC:  -654.0062   BIC:  -644.1907 
## Correlation matrix:
##               mean            sd
## mean  1.000000e+00 -5.483659e-12
## sd   -5.483659e-12  1.000000e+00

#大きさ10の標本による標本標準偏差
mean(R)

## [1] 0.7209155

sd(R)#相関係数の標準誤差

## [1] 0.1742191

#相関係数の標準誤差の公式による値
(1-0.75^2)/sqrt((10))

## [1] 0.1383496

mean(K)

## [1] 0.2302323

hist(K,prob=T)

母相関係数が０でない場合は、あ当てはまりもよくないし、近似もよくない。

2.2.1 確率密度関数.

母相関係数が０でない場合、標本相関係数は正規分布にならず、少し右に偏った分布になる。
その場合の確率密度関数については、猫も杓子も構造化 様「相関係数の標本分布」に説明があるので、コードを使用させていただくと.

library(hypergeo)
rho=0.75;size=10
d.samplecor <- function(r, rho, size){
  shisuu1 <- (size-1)/2
  shisuu2 <- (size-4)/2
  shisuu3 <- size - 3/2

  bunshi <- (size-2) * gamma(size-1) * ((1 - rho^2)^shisuu1) * ((1-r^2)^shisuu2)
  bunbo <- sqrt(2*pi)* gamma(size - 1/2)*((1-rho*r)^shisuu3)
  choukika <- hypergeo(1/2,1/2,(2*size-1)/2,(rho*r+1)/2)

  value <- (bunshi/bunbo) * choukika
  return(Re(value))
}
#d.samplecor(0.75,0.75, 10)

hist((R),prob=T)
curve(d.samplecor(x,rho, size), col="red",add=T)

確かによく当てはまっている。
参考までに、相関係数をカンゼンニリカイするによれば、確率密度関数は次の式になる。

\({f(r)=\frac{(n-2)\Gamma(n-1)(1-\rho^2)^\frac{n-1}{2}(1-r^2)^\frac{n-4}{2}}{\sqrt{2\pi}\Gamma(n-\frac{1}{2})(1-\rho r)^{n-\frac{3}{2}}} {}_2 F_1\Bigl(\frac{1}{2},\frac{1}{2};\frac{1}{2}(2n-1);\frac{1}{2}(\rho r+1)\Bigr) }\)

2.3 フィッシャーの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 正規確率紙.

木暮仁様「正規確率紙」によれば、 横軸に一様な目盛りで縦軸が正規確率目盛りを持つグラフ用紙のことです。 測定した標本データを正規確率紙上にプロットするだけで、正規分布になっているかどうか、正規分布になっているならば、標準偏差や上位２０％点はいくつになるかなどを、簡単に知ることができます.

このほか 横浜国立大学酒井信介先生が詳しく説明しておられます。

• 「確率紙の考え方」には、数理まで含めた詳しい説明.
  
• 「確率紙で学ぶ確率・統計」及び、.
• 「R入門 第2回 グラフィクス・データセット」 には、Rによる実装例があります.

参考サイト.
QCプラネット様「正規確率紙がよく分かる」

以下、

群馬大学 青木繁先生「正規確率紙」. に従ってRで行ってみます。

https://bellcurve.jp/statistics/course/12931.html

正規確率紙による方法 https://www.jstage.jst.go.jp/article/jve1977/1978/2/1978_2_1/_pdf/-char/ja

3.1 実行例.

青木先生の関数を用いて、

• \(N(0,1)\) に従う正規乱数100個（左図）.
• \(B(100,0.5)\) に従う二項乱数50個（右図）.

正規確率紙にプロットする。

source("http://aoki2.si.gunma-u.ac.jp/R/src/npp2.R", encoding="euc-jp")
par(mfrow = c(1,2))
npp2(rnorm(100))
npp2(rbinom(50,100,1/2))

直線状に点が並んでいれば正規性があると考えられる。
二項分布も回数が大きければ正規分布に近いことがわかる。

3.2 150個の正規乱数.

probs <- c(0.01, 0.1, 1, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 99.9, 99.99)/100
par(mfrow = c(1,2))
x <- rnorm(150,5,2)
n <- length(x)               # データの個数
x <- sort(x)         # 昇順に並べ替える
y <- (1:n-0.5)/n     # 累積相対確率

• 左図は分布関数の逆関数.
• 右図は 観測値x → 経験確率分布値 → 逆関数値x.

par(mfrow = c(1,2))
plot(y,qnorm(y,0,1),ylim=c(-3,3),yaxt="n",xlab="prob",main="prob -> x")
axis(2, qnorm(probs), probs*100)#重要
plot(x, qnorm(y),ylim=c(-3,3),yaxt="n",main="x -> prob -> x")
axis(2, qnorm(probs), probs*100)#重要
abline(v=5,h=qnorm(0.5),col="lightblue")

右図について、
縦軸の50%のときのxを見ると、xの平均が読み取れる 。

参考サイト.
Rjpwiki様「確率紙プロット.
当サイト番外編 正規確率紙.
福島大学.
　標準正規分布の分布関数G(k)について.
　縦軸の目盛kのところにG(k)の目盛をとり、横軸は等間隔にとったグラフ用紙.

4 例題15.

参考サイト. Rを使ったデータ処理「相関係数の検定」.

元のデータが (x,y) として与えられていれば、

関数 cor.test(x,y) で相関について検定することができる.
ここでは関数を使わず、標本相関係数から値を構成する。

r=0.25
t <- r*sqrt((30-2)/(1-r^2))
t

## [1] 1.36626

1-pt(t,28,lower.tail = TRUE)

## [1] 0.09136525

P値は 0.09>0.05 であるので、帰無仮説は棄却されない。
即ち、無相関は棄却されない。

参考までに、N=60 と標本数が増えれば、

r=0.25
t <- r*sqrt((60-2)/(1-r^2))
t

## [1] 1.966384

1-pt(t,58,lower.tail = TRUE)

## [1] 0.02702299

P値は 0.027<0.05 であるので、無相関仮説は棄却される。
サンプル数によって検定結果は変わってくる。

5 演習15.

5.1 (1).

r=0.21
t <- r*sqrt((50-2)/(1-r^2))
t

## [1] 1.488105

1-pt(t,29,lower.tail = TRUE)

## [1] 0.07375872

有意ではなく、無相関は棄却されない.

5.2 (2).

度数分布表から、正規確率紙を用いて、

• 正規性の判定.
• 母平均、母標準偏差の推定.

を行う。

• データから経験累積分布関数の値を計算し.
• 経験累積分布関数の分位点から理論分布関数によりxを逆算する.

本来は度数分布表から母平均と母標準偏差の推定値を求めるべきであるが、簡単のため、mean=53, sd=13　を用いて理論分布関数を描くと、ほぼ一致していることがわかる（左図）。

• xの逆算で qnorm(y1,50,10) を用いると、(x1,x2) 直線上に位置する.
  
• 回帰直線を描くと、50%点に x=53 が対応し、これが平均の予測値となる（右図）.

y <- c(3,12,18,35,15,5)
x1 <- c(29,39,49,59,69,79)
y1 <- cumsum(c(3,12,18,35,15,5))/90#probsに相当
par(mfrow = c(1,2))
plot(x1,y1)
curve(pnorm(x, mean=53, sd=13),add=T,col="red")
#予測点
x2 <- qnorm(y1)
data <- data.frame(x1=x1,x2=x2)
model <- lm(x2 ~ x1, data = data)

plot(x1,x2,yaxt="n")
axis(2, qnorm(probs),round(100*probs,0))#座標目盛
abline(model,col="lightblue")
abline(h=qnorm(0.5),col="red",lty = 3)
abline(v=53,col="red",lty = 3)

注意. もし、qnorm(y1,53,13) を用いれば、回帰直線は y=x になる。 平均・標準偏差を変えても直線性は保たれる。
（一次変換で移るからと考えられる）.

右図を大きく描いてみよう.

probs <- c(0.01, 0.1, 1, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 99.9, 99.99)/100
x3 <- qnorm(probs,50,10)
#x3
x2 <- qnorm(y1)
data <- data.frame(x1=x1,x2=x2)
model <- lm(x2 ~ x1, data = data)

plot(x1,qnorm(y1),yaxt="n")
abline(model,col="blue")
axis(2, qnorm(y1), round(100*y1,0))#座標目盛
abline(h=qnorm(0.5),col="red",lty = 3)
abline(v=53,col="red",lty = 3)
abline(h=qnorm(0.841),col="blue",lty = 3)
abline(v=66,col="blue",lty = 3)

縦軸50に対応するx1は 53 に近い。
標準偏差については、縦軸 84.1 に対応するx1は 66 に近いので、
標準偏差は　66-53=13 と推定される。

お詫び.
縦軸の目盛の取り方を試行錯誤しているうちに、分かりにくいコードになってしまいました。不十分な点はご容赦ください。

5.3 平均0、標準偏差1.

\(N(0,1)\)に従う100個の乱数を出力し、正規確率紙にプロットする。

x <- rnorm(100,0,1)
x <- sort(x)         # 昇順に並べ替える
n <- length(x) 

y <- pnorm(x)#理論値
y1 <- (1:n-0.5)/n　　# 経験値

par(mfrow = c(1,1))
plot(x,qnorm(y1),yaxt="n")
axis(2, qnorm(y1),round(100*y1,0))#座標目盛
abline(h=qnorm(0.5),col="red",lty = 3)
abline(v=0,col="red",lty = 3)
abline(h=qnorm(0.841),col="blue",lty = 3)
abline(v=1,col="blue",lty = 3)

正規確率紙が妥当な結果を与えてくれることが分かる。

正規確率紙の特徴は、縦軸の座標の取り方を工夫することで、プロットを容易にしていることにある。
演習15(2)の縦軸目盛をフル目盛にして出力しておこう。

y <- c(3,12,18,35,15,5)
x1 <- c(29,39,49,59,69,79)
y1 <- cumsum(c(3,12,18,35,15,5))/90

x <- c(29,39,49,59,69,79)        # 昇順に並べ替える
y <- cumsum(c(3,12,18,35,15,5))/90

data <- data.frame(x1=x1,x2=qnorm(y1))
model <- lm(x2 ~ x1, data = data)
#y1 <- (1:n-0.5)/n　　# 経験値

par(mfrow = c(1,1))
plot(x,qnorm(y1),yaxt="n")
abline(model)
axis(2, qnorm(probs),round(100*probs,0))#座標目盛
abline(h=qnorm(0.5),col="red",lty = 3)
abline(v=53,col="red",lty = 3)
abline(h=qnorm(0.841),col="blue",lty = 3)
abline(v=66,col="blue",lty = 3)

6 終わりにかえて.

小寺平治先生の「はじめての統計15講」は、統計を初めて学ぶ人向けに、数式などをできるだけ使わずに書かれています。
しかし、講義の教材として読み解いていくと、このテキストが堅固な数理的基盤のもとで書かれているのがわかります。確率変数など理解の難しいところは、乱数によるシミュレーションなどで説明すると良いと思います。

最後の章に「無相関の検定」を持ってきているところは意味深長です。
相関係数の分析は、「多変量解析」への入り口であるからです。
多変量解析はデータサイエンスの中で重要な役割を果たしており、その有用性から高校情報などにも取り入れられています。

数学的に高度であるため、統計学は高等学校で教えられてきませんでした。
今日、RやPythonなどの無料ソフトを活用し体験的に学ぶ道が開かれています。 当サイトで提唱する「ひょっこり方式」は、その1つの例と考えてください。

RやRstudioを利用することは、慣れるまでは大変と思います。

数値シミュレーションを活用した統計の学習教材は、 専修大学教授 小杉考司先生の. 「数値シミュレーションで読み解く統計のしくみ～Rでためしてわかる心理統計」.
統計ソフトRについても丁寧に説明されており、Rを使って統計を学びたい方にお勧めです。

最後に、
これから教材化したいテーマについて自分のメモとして列記し、自分への宿題としたいと思います。

• 相関分析.

• 多変量解析.

• パス係数.

• 因子分析.

観測変数間の相関.
Bartlettの球面性検定(共分散行列が球形性を持たないという帰無仮説の検定).
有意であれば変数間に相関があり因子分析を行うには妥当(帰無仮説は無相間).

Bartlett の球形検定と Mauchly の球形検定の違い. Bartlettは相関行列の無相関性の検定. Mauchlyは被験者内要因の分散が等しい検定. （被験者間計画の分散分析でいう分散の等質性の仮定に相当する）.
http://riseki.php.xdomain.jp/index.php?ANOVA君%2F球面性検定の出力.

相関が弱ければ共通因子がないことになり、因子分析は意味をなさなくなる。（相関0.3以上の項目がなければ因子分析には適さない）.
相関が0.9以上のものは多重共線性（同一内容測定）が考えられ、どちらかの変数を除外する。

• 信頼性係数.

• 主成分分析.

変数間の相関関係を用いて無相関な合成変数を求めることで、多くの変数を少ない変数に縮約する手法です。

http://hyoccori-toukei.la.coocan.jp/id-8/id-3/id-9.html

http://hyoccori-toukei.la.coocan.jp/id-5/id-5/id-3.html.
分割表の分析

Bartlett (1950) の相関行列の検定