3.4_任意の相関係数を持つ変数が従う確率分布

数値シミュレーションで読み解く統計の仕組み・サポートサイト.

1 相関係数と2変量正規分布.

1.1 相関係数の例.

set.seed(123) # 乱数のシードを指定
n <- 100 # 標本サイズ
# 数学の点数が従う正規分布の平均
math_mu <- 50
# 数学の点数が従う正規分布の標準偏差
math_sigma <- 10
# 正規分布に従う、数学テストの得点の実現値
math <- rnorm(n, math_mu, math_sigma)

# 身長が従う正規分布の平均
height_mu <- 160
# 身長が従う正規分布の標準偏差
height_sigma <- 15
# 正規分布に従う、身長の実現値
height <- rnorm(n, height_mu, height_sigma)

cor(math, height) # 相関係数の実現値

## [1] -0.04953215

plot(math, height) # 散布図

1.2 2変量正規分布.

dmnorm: Multivariate Normal Density Function.

• dmnorm(X, Mu, Sigma).

library(mnormt)
x <- y <- seq(from = -5, to = 5, by = 0.1)
mu_vec <- c(0, 0)
rho <- 0.4
covmatrix <- matrix(c(1, rho, rho, 1), nrow = 2, ncol = 2)
density2d <- function(x, y) {
  dmnorm(cbind(x, y), mu_vec, covmatrix)
}
z <- outer(x, y, density2d)

par(family = "HiraKakuProN-W3") #日本語フォントの指定
# par(mfrow = c(2, 1))
persp(x, y, z, theta = 0, phi = 25, expand = 0.6, ticktype = "detailed", main = "X軸側から見た正規分布")

persp(x, y, z, theta = 90, phi = 25, expand = 0.6, ticktype = "detailed", main = "Y軸側から見た正規分布")

1.3 2変量正規乱数の生成.

当サイト.

• 2変量正規分布と線形回帰分析.
  
• 多変量正規分布 参照.

\((x,y)\) が2変量正規分布に従うとき、確率密度関数は 平均と分散共分散行列によって決まる。
\(\begin{pmatrix} x \\ y \end{pmatrix} \sim N_2\left( \begin{pmatrix} \mu_x \\ \mu_y \end{pmatrix},\begin{pmatrix} \sigma_x^2 & \sigma_{xy}\\ \sigma_{xy} & \sigma_y^2 \end{pmatrix}\right) = N_2\left( \begin{pmatrix} \mu_x \\ \mu_y \end{pmatrix},\begin{pmatrix} \sigma_x^2 & \rho\sigma_x \sigma_y \\ \rho\sigma_x \sigma_y & \sigma_y^2 \end{pmatrix}\right)\) のとき、
\(\Sigma^{-1}= \dfrac{1}{\sigma_x^2\sigma_y^2(1-\rho^2)} \begin{pmatrix} \sigma_y^2 & -\sigma_{xy}\\ -\sigma_{xy} & \sigma_x^2 \end{pmatrix}\)　であり、

• \(f(x,y)=\dfrac{1}{\left(\sqrt{2\pi}\right)^2\sqrt{|\Sigma|}}\exp\left[{-\dfrac{1}{2}(X-\mu)'\ \Sigma^{-1}(X-\mu)}\right]\).

• \(x\sim N(\mu_x,\sigma_x^2)\)、　\(y\sim N(\mu_y,\sigma_y^2)\).

• \(E(y|x)=\mu_y+\rho \sigma_y\left(\dfrac{x-\mu_x}{\sigma_x}\right)\)、　\(V(y|x)=(1-\rho^2)\sigma_y^2\).
  \(x\) を固定したときの条件付き期待値は、 \(y\) の \(x\) による回帰直線に等しい.

1.3.1 \((x,y)=\)（数学テスト得点,身長）.

共分散は　\(\sigma_{xy}=\rho\sigma_x \sigma_y\) によって計算される（相関係数によって共分散の値は定まる）。

# それぞれ、数学テストの得点と身長が従う正規分布の平均
mu_vec <- c(math_mu, height_mu)
mu_vec

## [1]  50 160

rho <- 0.4 # 数学と身長が従う2変量正規分布の相関パラメータ
tau <- rho * math_sigma * height_sigma # 数学と身長の共分散
tau

## [1] 60

2 ３変量正規乱数.

age_mu <- 15 # 年齢の平均
age_sigma <- 2 # 年齢の標準偏差
# それぞれ、数学テストの得点、身長、年齢が従う正規分布の平均
mu_vec <- c(math_mu, height_mu, age_mu)
# それぞれ、数学テストの得点、身長、年齢が従う正規分布の標準偏差
sigma_vec <- c(math_sigma, height_sigma, age_sigma)

2.0.1 相関行列.

rho_matrix <- matrix(c(1, 0.4, 0.5, 0.4, 1, 0.8, 0.5, 0.8, 1),
  nrow = length(sigma_vec)
)
# 行名を付ける。説明のために命名しているので、実用上は必須ではない
rownames(rho_matrix) <- c("math", "height", "age")
# 列名を付ける。説明のために命名しているので、実用上は必須ではない
colnames(rho_matrix) <- c("math", "height", "age")
rho_matrix

##        math height age
## math    1.0    0.4 0.5
## height  0.4    1.0 0.8
## age     0.5    0.8 1.0

2.0.2 分散共分散行列の算出.

分散ベクトルと相関行列から、分散共分散行列を計算する.

cov_matrix <- function(sigma_vec, rho_mx) {
  return(
    # 線形代数による演算
    diag(sigma_vec) %*% rho_mx %*% diag(sigma_vec)
  )
}

## ----------------------------------------------------
diag(sigma_vec)

##      [,1] [,2] [,3]
## [1,]   10    0    0
## [2,]    0   15    0
## [3,]    0    0    2

rho_matrix

##        math height age
## math    1.0    0.4 0.5
## height  0.4    1.0 0.8
## age     0.5    0.8 1.0

cov_matrix(sigma_vec = sigma_vec, rho_mx = rho_matrix)

##      [,1] [,2] [,3]
## [1,]  100   60   10
## [2,]   60  225   24
## [3,]   10   24    4

2.1 ３変量正規乱数の生成.

set.seed(123)
dat_3norm <- MASS::mvrnorm(
  n = 2000,
  mu = mu_vec,
  Sigma = cov_matrix(sigma_vec = sigma_vec, rho_mx = rho_matrix)
)

colnames(dat_3norm) <- c("math", "height", "age")

head(as.data.frame(dat_3norm), 5) # 冒頭5行を表示

math	height	age
49.16093	169.8511	16.12952
53.23926	162.5290	16.11891
36.40559	138.8688	13.15220
59.49592	155.1357	13.48781
50.74273	157.7575	12.51995

cor(dat_3norm)

##             math    height       age
## math   1.0000000 0.4167597 0.5183609
## height 0.4167597 1.0000000 0.8019000
## age    0.5183609 0.8019000 1.0000000

pairs(dat_3norm)

2.2 偏相関係数の算出.

pcor3 <- function(x, y, covariate) {
  # 偏相関係数の計算における分子
  pcor_top <- cor(x, y) - (cor(x, covariate) * cor(y, covariate))
  # 偏相関係数の計算における分母
  pcor_bottom <- sqrt(1 - cor(x, covariate)^2) * sqrt(1 - (cor(y, covariate)^2))
  return(pcor_top / pcor_bottom) # 偏相関係数
}

dat_3norm <- as.data.frame(dat_3norm)
pcor3(x = dat_3norm$math, y = dat_3norm$height, covariate = dat_3norm$age)

## [1] 0.002125782

#パッケージによる計算
library(ppcor)

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

result <- pcor(dat_3norm[,1:3])
result$estimate

##               math      height       age
## math   1.000000000 0.002125782 0.3390930
## height 0.002125782 1.000000000 0.7536671
## age    0.339093043 0.753667137 1.0000000

2.3 相関係数と偏相関係数.

cor(dat_3norm)

##             math    height       age
## math   1.0000000 0.4167597 0.5183609
## height 0.4167597 1.0000000 0.8019000
## age    0.5183609 0.8019000 1.0000000

result$estimate

##               math      height       age
## math   1.000000000 0.002125782 0.3390930
## height 0.002125782 1.000000000 0.7536671
## age    0.339093043 0.753667137 1.0000000

2.4 標準化単回帰係数.

• 標準化単回帰係数は相関係数に等しい.
  
• 双方向の標準化単回帰係数は同じである.
  （逆方向の回帰で、標準化単回帰係が互いに逆数になる訳ではない）.

math_age <- lm(math ~ age, data = as.data.frame(scale(dat_3norm)))
round(math_age$coefficients,3)

## (Intercept)         age 
##       0.000       0.518

age_math <- lm(age ~ math, data = as.data.frame(scale(dat_3norm)))
round(age_math$coefficients,3)

## (Intercept)        math 
##       0.000       0.518

2.5 標準化偏回帰係数.

• 偏相関係数.
  　　　\(\displaystyle r [(y|x_1), (x_2|x_1)] = \frac{r[y, x_2] - r[y, x_1] r[x_1,x_2]}{\sqrt{1-r^2[y, x_1]} \sqrt{1-r^2[x_1, x_2]}}\).

• 標準偏回帰係数.
  　　　\(\displaystyle b^{*}[ y, (x_2|x_1) ] =b^{*}[ (y|x_1), (x_2|x_1) ]= \frac{r[y, x_2] - r[y, x_1] r[x_1,x_2]}{ 1-r^2[x_1, x_2] }\).

• 標準偏回帰係数と偏相関係数の関係.

　　　　　　\(\displaystyle b^{*}[ y, (x_2|x_1) ] = r [(y|x_1), (x_2|x_1)] \cdot\sqrt{\frac{1-r^2[y, x_1]} {1-r^2[x_2, x_1]} }\)

いずれも、相関係数だけから算出できることに注意しておこう.

math_age_height <- lm(math ~ age+height, data = as.data.frame(scale(dat_3norm)))
round(math_age_height$coefficients,3)

## (Intercept)         age      height 
##       0.000       0.516       0.003

age_math_height <- lm(age ~ math+height, data = as.data.frame(scale(dat_3norm)))
round(age_math_height$coefficients,3)

## (Intercept)        math      height 
##       0.000       0.223       0.709

height_math_age <- lm(height ~ math+age, data = as.data.frame(scale(dat_3norm)))
round(height_math_age$coefficients,3)

## (Intercept)        math         age 
##       0.000       0.001       0.801

• 標準化偏回帰係数は相関係数と若干異なる.
  
• 双方向の標準化偏回帰係数も互いに異なっている.
  （標準化偏回帰係数の相乗平均は、標準化偏相関係数に等しい）.

実際、

• 標準化偏回帰係数 math-age：0.516.
  （相関係数は0.518：heightの影響は小さいと考えられる）.

• 標準化偏回帰係数 math-height：0.003.
  （相関係数0.417：ageの影響は大きいと考えられる）.

• 双方向偏回帰係数.
  　math-age：0.516　　　age-math：0.223.
  （相乗平均は、偏相関係数 0.339 に等しい）.

まとめ.

• 相関係数は双方向の回帰係数の相乗平均.

• 偏相関係数も双方向の標準化偏回帰係数の相乗平均である.

2.6 構造方程式.

dat = as.data.frame(scale(dat_3norm))

2.6.1 　

attach(dat)

## The following objects are masked _by_ .GlobalEnv:
## 
##     height, math

model.RM <- '

  height ~ age
  math ~ age
'

2.6.2 　

attach(dat)

## The following objects are masked _by_ .GlobalEnv:
## 
##     height, math

## The following objects are masked from dat (pos = 3):
## 
##     age, height, math

model.RM <- '
  math ~ height
  height ~ age
  math ~ age
'

2.6.3 　

attach(dat)

## The following objects are masked _by_ .GlobalEnv:
## 
##     height, math

## The following objects are masked from dat (pos = 3):
## 
##     age, height, math

## The following objects are masked from dat (pos = 4):
## 
##     age, height, math

model.RM <- '
  math ~ height + age


'

2.7 影響力. 　

単純な相関は、本来の学歴の作用である総合効果と、性別によるニセの作用である擬似相関に分解. https://www.ritsumei.ac.jp/~tsutsui/misc/regression.pdf

# 数学が目的変数、年齢が説明変数の回帰分析
residual_math <- lm(math ~ age, data = dat_3norm)
residual_math$coefficients

## (Intercept)         age 
##   11.268043    2.570916

  # 残差の取り出し
#  residuals()


# 身長が目的変数、年齢が説明変数の回帰分析
residual_height <- lm(height ~ age, data = dat_3norm) |>
  # 残差の取り出し
  residuals()

#cor(residual_math, residual_height) # 残差同士の相関係数


## -------------------------------------------------------------------------------------------
# 年齢と身長の相関係数（0 ~ 0.9）
rho_age_height <- seq(from = 0, to = 0.9, by = 0.3)

set.seed(123)
for (i in 1:length(rho_age_height)) {
  rho_matrix <- matrix(
    c(
      1, 0.4, 0.5,
      0.4, 1, rho_age_height[i],
      0.5, rho_age_height[i], 1
    ),
    nrow = length(sigma_vec)
  )

  dat_3norm <- MASS::mvrnorm(
    n = 2000,
    mu = mu_vec,
    Sigma = cov_matrix(sigma_vec, rho_matrix)
  )

  print(
    paste(
      "年齢(共変量)と身長のρ = ",
      rho_age_height[i],
      ", 数学と身長の偏相関係数 = ",
      round(pcor3(dat_3norm[, 1], dat_3norm[, 2], dat_3norm[, 3]), 4)
    ) # paste()は、複数の引数を文字列に変換したうえで結合する関数
  )
}

## [1] "年齢(共変量)と身長のρ =  0 , 数学と身長の偏相関係数 =  0.4744"
## [1] "年齢(共変量)と身長のρ =  0.3 , 数学と身長の偏相関係数 =  0.2536"
## [1] "年齢(共変量)と身長のρ =  0.6 , 数学と身長の偏相関係数 =  0.1657"
## [1] "年齢(共変量)と身長のρ =  0.9 , 数学と身長の偏相関係数 =  -0.1259"

3 統計検定人文科学2018③.

rho_matrix <- matrix(c(1, -0.4, -0.2, 0.34,
                       -0.4, 1, 0.5, 0.2,
                       -0.2,0.85, 1,0.4,
                       0.34,0.2,0.4,1),
  nrow = 4
)
# 行名を付ける。説明のために命名しているので、実用上は必須ではない
rownames(rho_matrix) <- c("X1", "X2", "X3","X4")
# 列名を付ける。説明のために命名しているので、実用上は必須ではない
colnames(rho_matrix) <- c("X1", "X2", "X3","X4")
rho_matrix

##       X1   X2    X3   X4
## X1  1.00 -0.4 -0.20 0.34
## X2 -0.40  1.0  0.85 0.20
## X3 -0.20  0.5  1.00 0.40
## X4  0.34  0.2  0.40 1.00

　　　\(\displaystyle r [(y|x_1), (x_2|x_1)] = \frac{r[y, x_2] - r[y, x_1] r[x_1,x_2]}{\sqrt{1-r^2[y, x_1]} \sqrt{1-r^2[x_1, x_2]}}\).

　　　\(\displaystyle b^{*}[ y, (x_2|x_1) ] = r [(y|x_1), (x_2|x_1)] \cdot\sqrt{\frac{1-r^2[y, x_1]} {1-r^2[x_2, x_1]} }\).

# 偏相関 (partial correlation) 
r42 <- rho_matrix[4,2]
r41 <- rho_matrix[4,1]
r21 <- rho_matrix[2,1]

pry42 <- (r42-r41*r21)/((sqrt(1-r41^2))*(sqrt(1-r21^2)))
pry42*sqrt((1-r41^2)/(1-r21^2))

## [1] 0.4

3.0.1 \(\alpha_{32}\)

# 偏相関 (partial correlation) 
r32 <- rho_matrix[3,2]
r31 <- rho_matrix[3,1]
r21 <- rho_matrix[2,1]

pry32 <- (r32-r31*r21)/((sqrt(1-r31^2))*(sqrt(1-r21^2)))
pry32*sqrt((1-r31^2)/(1-r21^2))

## [1] 0.5

3.0.2 \(\alpha_{41}\)

# 偏相関 (partial correlation) 
r41 <- rho_matrix[4,1]
r42 <- rho_matrix[4,2]
r21 <- rho_matrix[2,1]

pry41 <- (r41-r42*r21)/((sqrt(1-r42^2))*(sqrt(1-r21^2)))
pry41*sqrt((1-r42^2)/(1-r21^2))

## [1] 0.5

3.0.3 \(\alpha_{43}\)

# 偏相関 (partial correlation) 
r43 <- rho_matrix[4,1]
r42 <- rho_matrix[4,2]
r23 <- rho_matrix[2,3]

pry43 <- (r43-r42*r23)/((sqrt(1-r42^2))*(sqrt(1-r23^2)))
pry43*sqrt((1-r42^2)/(1-r23^2))

## [1] 0.6126126

（\(X_4\)）については、共変数が２種類（\(X_3,X_2,X_1\)）とあるので、偏相関係数から計算するのは無理である. 乱数を生成し、lm関数で求めてみよう.

3.1 乱数シミュレーション.

set.seed(123)
dat_4norm <- MASS::mvrnorm(
  n = 2000,
  mu = c(0,0,0,0),
  Sigma = matrix(c(1, -0.4, -0.2, 0.34,
                       -0.4, 1, 0.5, 0.2,
                       -0.2,0.85, 1,0.4,
                       0.34,0.2,0.4,1),nrow = 4))

colnames(dat_4norm) <- c("X1", "X2", "X3","X4")

head(as.data.frame(dat_4norm), 5) # 冒頭5行を表示

X1	X2	X3	X4
-0.3824782	0.7931178	0.4226029	-0.4068081
-0.2089962	0.2614478	-0.2564340	0.7303176
0.5722769	-0.7712809	-1.5941763	-1.3105192
0.9240242	-0.8605498	0.7813284	0.7854991
-0.2380381	-1.1998006	0.7489498	0.0948816

3.1.1 重回帰.

result <- lm(X4 ~ X1 + X2 +X3,data=as.data.frame(dat_4norm))
result$coefficients

## (Intercept)          X1          X2          X3 
## -0.01870683  0.49508646  0.20338667  0.40423172

result1 <- lm(X2 ~ X1,data=as.data.frame(dat_4norm))
result1$coefficients

## (Intercept)          X1 
## -0.01644867 -0.39144704

result2 <- lm(X3 ~ X2
,data=as.data.frame(dat_4norm))
result2$coefficients

## (Intercept)          X2 
##  -0.0157985   0.5043357

dat <- as.data.frame(dat_4norm)
attach(dat)
model.RM <- '
  X2 ~ X1
  X3 ~ X2
  X4 ~ X1 + X2 +X3
'  

相関行列とモデルが決まれば、パス図は一意に定まる。
パス係数は、問題の解答と一致している

3.2 モデルから生成.

e2 <- rnorm(100,0,0.92)
e3 <- rnorm(100,0,0.86)
e4 <- rnorm(100,0,0.79)
X1 <- rnorm(100,0,1)
X2 <- -0.4*X1+e2
X3 <- 0.5*X2+e3
X4 <-0.5*X1 + 0.2*X2 +0.4*X3+e4
data <- data.frame(X1=X1,X2=X2,X3=X3,X4=X4)
head(data)

X1	X2	X3	X4
-0.1092598	0.5138243	-0.6169785	-0.5898894
-0.4735071	1.8527562	0.4791925	-0.4897068
0.1995025	-1.6462067	-0.7623974	-1.0825718
0.6273856	0.0134986	0.3590428	-0.5407261
1.6919625	-0.9243700	-1.8853274	0.7943057
1.6634383	-1.0145381	-0.8449257	0.5917286

round(cor(data),2)

##       X1    X2    X3   X4
## X1  1.00 -0.39 -0.26 0.41
## X2 -0.39  1.00  0.57 0.11
## X3 -0.26  0.57  1.00 0.37
## X4  0.41  0.11  0.37 1.00

e <- rnorm(100,0,1)
X1 <- rnorm(100,0,1)
X2 <- -0.4*X1+e
cor(X1,X2)

## [1] -0.3652245

4 統計検定人文科学2019④.

e10 <- rnorm(1000,0,sqrt(0.59))
e11 <- rnorm(1000,0,sqrt(0.998))
f1 <- rnorm(1000,0,1)
f2 <- 0.64*f1+e10
f3 <- 0.014*f1+0.032*f2+e11

e1 <- rnorm(1000,0,sqrt(0.572))
e2 <- rnorm(1000,0,sqrt(0.687))
e3 <- rnorm(1000,0,sqrt(0.725))
y1 <- 0.654*f1+e1
y2 <- 0.559*f1+e2
y3 <- 0.524*f1+e3

e4 <- rnorm(1000,0,sqrt(0.6))
e5 <- rnorm(1000,0,sqrt(0.512))
e6 <- rnorm(1000,0,sqrt(0.668))
y4 <- 0.632*f2+e4
y5 <- 0.699*f2+e5
y6 <- 0.576*f2+e6

e7 <- rnorm(1000,0,sqrt(0.605))
e8 <- rnorm(1000,0,sqrt(0.723))
e9 <- rnorm(1000,0,sqrt(0.723))
y7 <- 0.628*f3+e7
y8 <- 0.527*f3+e8
y9 <- 0.526*f3+e9

data <- data.frame(y1=y1,y2=y2,y3=y3,y4=y4,y5=y5,y6=y6,y7=y7,y8=y8,y9=y9)
str(data)

## 'data.frame':    1000 obs. of  9 variables:
##  $ y1: num  0.6518 0.0827 -0.1167 -0.1715 1.2252 ...
##  $ y2: num  -0.4161 -0.4672 0.0307 0.6556 0.4845 ...
##  $ y3: num  0.0471 0.4617 0.076 1.9888 1.1682 ...
##  $ y4: num  -0.4353 0.3373 0.3345 -0.4959 -0.0993 ...
##  $ y5: num  -0.245 0.979 1.506 -0.418 1.35 ...
##  $ y6: num  -0.27 0.919 1.394 -0.9 -0.184 ...
##  $ y7: num  -0.457 0.529 0.211 0.233 -0.661 ...
##  $ y8: num  -2.01 0.655 0.156 -1.677 0.194 ...
##  $ y9: num  -0.0866 -0.0947 -0.5407 -0.4327 1.3276 ...

model <- '
  # 測定方程式
    f1 =~ y1 + y2 + y3
    f2 =~ y4 + y5 + y6
    f3 =~ y7 + y8 + y9

  # 構造方程式
    f2 ~ f1
    f3 ~ f1 + f2
    
  # 共変関係
'

hist(data)

#dat = as.data.frame(scale(data))
fit <- sem(model, data=data)
#semPaths(object=fit, whatLabels="stand")
semPaths(fit, what = "stand")

e10 <- rnorm(1000,0,sqrt(0.59))
e11 <- rnorm(1000,0,sqrt(0.998))
f1 <- rnorm(1000,0,1)
f2 <- 0.64*f1+e10
f3 <- 0.014*f1+0.032*f2+e11
dat <- data.frame(f1=f1,f2=f2,f3=f3)

attach(dat)

## The following objects are masked _by_ .GlobalEnv:
## 
##     f1, f2, f3

model.RM <- '
  f2 ~ f1
  f3 ~ f1 + f2
'
fit <- sem(model.RM, data=dat)
semPaths(fit, what = "stand", style = "lisrel", layout = "tree", rotation = 2, nCharNodes = 0, nCharEdges = 0, fade = FALSE, edge.width = 0.2, label.scale = FALSE, label.cex = 1.2, theme = 'gray', asize = 6.0, node.width = 2.0, curve = 2.0)