1 南風原p226.

scratch-R様を参考にさせていただきました。
内容は「心理統計学の基礎」（南風原朝和先生）によります.

x1 <- c(12,12,7,17,14,9,10,13,15,12,12,15,11,14,17,17,16,15,15,10,12,9,12,12,19,11,14,15,15,15,16,15,12,10,11,12,15,13,15,12,12,12,13,17,13,11,14,16,12,12)
x2 <- c(2,2,2,3,2,2,3,3,3,1,3,3,2,2,4,2,4,3,4,2,2,1,2,2,4,2,3,2,3,3,2,3,2,2,3,1,2,3,2,2,2,3,3,3,2,3,2,4,2,2)
y <- c(6,11,11,13,13,10,10,15,11,11,16,14,10,13,12,15,16,14,14,8,13,12,12,11,16,9,12,13,13,14,12,15,8,12,11,6,12,15,9,13,9,11,14,12,13,9,11,14,16,8)
dat <- data.frame(x1, x2, y)
#dat <- data.frame(x1=scale(x1), x2=scale(x2), y=scale(y))
dat
x1 x2 y
12 2 6
12 2 11
7 2 11
17 3 13
14 2 13
9 2 10
10 3 10
13 3 15
15 3 11
12 1 11
12 3 16
15 3 14
11 2 10
14 2 13
17 4 12
17 2 15
16 4 16
15 3 14
15 4 14
10 2 8
12 2 13
9 1 12
12 2 12
12 2 11
19 4 16
11 2 9
14 3 12
15 2 13
15 3 13
15 3 14
16 2 12
15 3 15
12 2 8
10 2 12
11 3 11
12 1 6
15 2 12
13 3 15
15 2 9
12 2 13
12 2 9
12 3 11
13 3 14
17 3 12
13 2 13
11 3 9
14 2 11
16 4 14
12 2 16
12 2 8

1.1 データの3Dプロット.

res <- lm(y~x1+x2,data=dat)
options(rgl.printRglwidget = TRUE) 
library(rgl)
plot3d(x1, x2, y, type = "s", col = "blue", size = 1)

coefs <- coef(res)
a <- coefs["x1"]
b <- coefs["x2"]
c <- -1
d <- coefs["(Intercept)"]
planes3d(a, b, c, d, alpha = 0.5)
summary(res)
## 
## Call:
## lm(formula = y ~ x1 + x2, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0563 -1.4780 -0.0563  1.4751  4.9437 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   5.1716     1.6647   3.107   0.0032 **
## x1            0.3027     0.1467   2.064   0.0446 * 
## x2            1.1259     0.4713   2.389   0.0210 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.12 on 47 degrees of freedom
## Multiple R-squared:  0.3137, Adjusted R-squared:  0.2845 
## F-statistic: 10.74 on 2 and 47 DF,  p-value: 0.000144

1.2 (y,x2)の偏相関係数

x1で回帰した残差の相関係数. 即ち、x1で説明できない部分同士の相関である.

result_y <- lm(y~x1, data=dat)

result_x2 <- lm(x2~x1, data=dat)
cor(result_y$residuals,result_x2$residuals)
## [1] 0.3290749

残差の相関係数が偏相関係数.

1.3 相関・偏相関グラフ.

https://kusanagi.hatenablog.jp/entry/2014/06/01/181132

library(corpcor)
library(qgraph)
par(mfrow=c(1,2))
qgraph(cor(dat),edge.labels=T)
qgraph(cor2pcor(cor(dat)),edge.labels=T)

各相関は大きいので、他の変数の影響を取り除いた偏相関は少し小さくなっている。

1.4 psychパッケージ.

library(psych)
## 
## Attaching package: 'psych'
## The following object is masked from 'package:lavaan':
## 
##     cor2cov
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
describe(dat)
vars n mean sd median trimmed mad min max range skew kurtosis se
x1 1 50 13.20 2.4494897 13 13.200 2.9652 7 19 12 -0.0146969 -0.3366667 0.3464102
x2 2 50 2.48 0.7623808 2 2.425 0.7413 1 4 3 0.3366715 -0.4155983 0.1078169
y 3 50 11.96 2.5068071 12 12.075 2.9652 6 16 10 -0.3898264 -0.3968307 0.3545161 
cor(dat)
##           x1        x2         y
## x1 1.0000000 0.5376770 0.4799263
## x2 0.5376770 1.0000000 0.5014633
## y  0.4799263 0.5014633 1.0000000
#r12 <- cor(x1,x2)
#ry1 <- cor(x1,y)
#ry2 <- cor(x2,y)
#偏相関係数
library(ppcor)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
result <- pcor(dat[,1:3])
result
## $estimate
##           x1        x2         y
## x1 1.0000000 0.3913046 0.2882904
## x2 0.3913046 1.0000000 0.3290749
## y  0.2882904 0.3290749 1.0000000
## 
## $p.value
##             x1          x2          y
## x1 0.000000000 0.005432028 0.04455729
## x2 0.005432028 0.000000000 0.02095613
## y  0.044557288 0.020956131 0.00000000
## 
## $statistic
##          x1       x2        y
## x1 0.000000 2.915095 2.064053
## x2 2.915095 0.000000 2.389087
## y  2.064053 2.389087 0.000000
## 
## $n
## [1] 50
## 
## $gp
## [1] 1
## 
## $method
## [1] "pearson"

1.5 パス図.

パス係数は、標準化偏回帰係数である.

attach(dat)
## The following objects are masked _by_ .GlobalEnv:
## 
##     x1, x2, y
model.RM <- '
  y ~ x1 + x2

'
fit.RM <- sem(model.RM, data = dat, estimator = "ML")
par(family = "HiraKakuProN-W3") #日本語フォントの指定
semPaths(fit.RM, 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)

パス係数は偏回帰係数であり、偏相関係数とほぼ一致している。

1.5.1 理論値.

パス係数（標準偏回帰係数）は、偏相関係数から次のように計算される。

　　　\(\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]} }\).

今、\(r[y,x_1]=0.48\)、\(r[x_2,x_1]=0.54\)　より.

　　　\(\displaystyle \sqrt{\frac{1-r^2[y, x_1]} {1-r^2[x_2, x_1]} }=\sqrt{\frac{1-0.48^2}{1-0.54^2}}=\sqrt{\frac{0.7696}{0.7084}}=1.042\).

つまり、ほぼ等しいと言っても、この問題では 1.042 倍である。
実際、
　　　偏回帰係数は　\(0.3290\times 1.042=0.3429\).

と計算され、パス係数と一致している。

もし、\(r[y,x_1]=r[x_2,x_1]\)　x2,yとx1の相関係数が同程度なら. 標準偏回帰係数は偏相関係数に一致する。即ち、パス係数は偏相関係数に一致する。

1.6 単回帰は相関係数.

attach(dat)
## The following objects are masked _by_ .GlobalEnv:
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 3):
## 
##     x1, x2, y
model.RM <- '
  y ~ x1
  x2 ~ x1
'
fit.RM <- sem(model.RM, data = dat, estimator = "ML")
par(family = "HiraKakuProN-W3") #日本語フォントの指定
semPaths(fit.RM, 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)

  • x2,yをx1で単回帰すると、回帰係数は相関係数.

  • x2,yをx1で単回帰した回帰残差同士の相関は偏相関係数になる.

  • 重回帰における偏回帰係数は偏相関係数の\(\sqrt{\dfrac{1-r^2[y, x_1]} {1-r^2[x_2, x_1]} }\)倍である.

偏相関係数の数理的説明については、

  • 当サイト[相関分析からSEM・パス解析へ]でアップする予定である。
    また、
  • 「こころ」のための専門メディア 金子書房 様のサイトの中で、
    関西学院大学社会学部教授:吉田寿夫先生による.
    「偏回帰係数についてあらためて考える」に詳しい説明があある。.

2 偏相関係数型モデル.

アイスタット様「共分散構造分析（５/７）」で学ばせていただいたものです.

x1：身長、x2：学年、y：漢字個数（100個中の識別個数）.

x1 <- c(182,163,160,173,168,175,165,170,165,172,178,173,176,170,169,152,157,175,174,155)
x2 <- c(3,1,1,3,2,3,1,2,3,2,2,1,3,2,1,2,2,3,3,1)
y <- c(99,15,19,71,39,47,11,67,87,23,59,15,47,35,15,35,47,39,79,31)
dat <- data.frame(x1=x1, x2=x2, y)
dat
x1 x2 y
182 3 99
163 1 15
160 1 19
173 3 71
168 2 39
175 3 47
165 1 11
170 2 67
165 3 87
172 2 23
178 2 59
173 1 15
176 3 47
170 2 35
169 1 15
152 2 35
157 2 47
175 3 39
174 3 79
155 1 31 
#interaction.plot(x1,x2,y, type="b",  pch=c(1,2))

summaryで、偏回帰係数のP値をみると、身長は有意ではない。
学年による影響が大きいことがわかる。

2.1 データの3Dプロット.

res <- lm(y~x1+x2,data=dat)
options(rgl.printRglwidget = TRUE) 
library(rgl)
plot3d(x1, x2, y, type = "s", col = "blue", size = 1)

coefs <- coef(res)
a <- coefs["x1"]
b <- coefs["x2"]
c <- -1
d <- coefs["(Intercept)"]
planes3d(a, b, c, d, alpha = 0.5)
summary(res)
## 
## Call:
## lm(formula = y ~ x1 + x2, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.409  -7.709  -3.149  11.923  31.556 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.285533  90.518033  -0.080 0.936790    
## x1           0.004981   0.568863   0.009 0.993115    
## x2          24.607664   5.509660   4.466 0.000339 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.82 on 17 degrees of freedom
## Multiple R-squared:  0.6203, Adjusted R-squared:  0.5756 
## F-statistic: 13.89 on 2 and 17 DF,  p-value: 0.0002662
result_y <- lm(y~x1, data=dat)

result_x2 <- lm(x2~x1, data=dat)
cor(result_y$residuals,result_x2$residuals)
## [1] 0.7347715
library(psych)
describe(dat)
vars n mean sd median trimmed mad min max range skew kurtosis se
x1 1 20 168.60 7.9960517 170 169.0625 7.4130 152 182 30 -0.4906636 -0.7594455 1.7879715
x2 2 20 2.05 0.8255779 2 2.0625 1.4826 1 3 2 -0.0839707 -1.6013746 0.1846048
y 3 20 44.00 25.8212479 39 41.7500 29.6520 11 99 88 0.5406838 -0.8777964 5.7738066 
cor(dat)
##           x1        x2         y
## x1 1.0000000 0.5293962 0.4180589
## x2 0.5293962 1.0000000 0.7875929
## y  0.4180589 0.7875929 1.0000000
#r12 <- cor(x1,x2)
#ry1 <- cor(x1,y)
#ry2 <- cor(x2,y)
#偏相関係数
library(ppcor)
result <- pcor(dat[,1:3])
result
## $estimate
##             x1        x2           y
## x1 1.000000000 0.3575359 0.002123715
## x2 0.357535949 1.0000000 0.734771519
## y  0.002123715 0.7347715 1.000000000
## 
## $p.value
##           x1           x2            y
## x1 0.0000000 0.1328768030 0.9931154915
## x2 0.1328768 0.0000000000 0.0003393787
## y  0.9931155 0.0003393787 0.0000000000
## 
## $statistic
##             x1       x2           y
## x1 0.000000000 1.578498 0.008756323
## x2 1.578498197 0.000000 4.466275970
## y  0.008756323 4.466276 0.000000000
## 
## $n
## [1] 20
## 
## $gp
## [1] 1
## 
## $method
## [1] "pearson"
attach(dat)
## The following objects are masked _by_ .GlobalEnv:
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 3):
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 4):
## 
##     x1, x2, y
model.RM <- '
  y ~ x1 + x2

'
fit.RM <- sem(model.RM, data = dat, estimator = "ML")
par(family = "HiraKakuProN-W3") #日本語フォントの指定
semPaths(fit.RM, 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)

attach(dat)
## The following objects are masked _by_ .GlobalEnv:
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 3):
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 4):
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 5):
## 
##     x1, x2, y
model.RM <- '
  y ~ x1
  x2 ~ x1
'
fit.RM <- sem(model.RM, data = dat, estimator = "ML")
par(family = "HiraKakuProN-W3") #日本語フォントの指定
semPaths(fit.RM, 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)

2.2 相関・偏相関グラフ.

https://kusanagi.hatenablog.jp/entry/2014/06/01/181132

library(corpcor)
library(qgraph)
par(mfrow=c(1,2))
qgraph(cor(dat),edge.labels=T)
qgraph(cor2pcor(cor(dat)),edge.labels=T)

3 パス係数の標準解と非標準解.

岡田安功様によると、

パス係数.

  • パス(→)の上に書く数値.
  • 一方の変数(外生変数,内生変数)が他方の変数(内生変数)に与える影響の強さ
  • 重回帰分析の偏回帰係数のようなもの
  • 標準解は重回帰分析における標準化偏回帰係数のようなもの.
    (測定値の単位に影響されず比較 が可能、パス図に書くのは通常標準解).

4 まとめ.

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

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

4.1 B図.

attach(dat)
## The following objects are masked _by_ .GlobalEnv:
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 3):
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 4):
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 5):
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 6):
## 
##     x1, x2, y
model.RM <- '
  y ~ x2
  x1 ~  x2
'
fit.RM <- sem(model.RM, data = dat, estimator = "ML")
par(family = "HiraKakuProN-W3") #日本語フォントの指定
semPaths(fit.RM, 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)

4.2 C図.

attach(dat)
## The following objects are masked _by_ .GlobalEnv:
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 3):
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 4):
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 5):
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 6):
## 
##     x1, x2, y
## The following objects are masked from dat (pos = 7):
## 
##     x1, x2, y
model.RM <- '
  x2 ~ y
  x1 ~ y
'
fit.RM <- sem(model.RM, data = dat, estimator = "ML")
par(family = "HiraKakuProN-W3") #日本語フォントの指定
semPaths(fit.RM, 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)

参考サイト.

いちばんやさしい統計学.

偏相関係数と標準化偏回帰係数の違い.
荻原祐二様.

吉田寿夫様「因果関係の探究における標準偏回帰係数の意味」.