myscratchpad

##　心理統計学の基礎p226 より。感謝

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)
library(psych)
describe(dat)[2:4]
cor(dat)

r12 <- cor(x1,x2)
ry1 <- cor(x1,y)
ry2 <- cor(x2,y)

# 部分相関 (semipartial correlation, part correlation)
## x2からx1の影響を除き、そのx2とyとの相関を求める
ry21 <- (ry2 -ry1*r12)/sqrt(1-r12^2)
ry21
## x1からx2の影響を除き、そのx1とyとの相関を求める
ry12 <- (ry1 -ry2*r12)/sqrt(1-r12^2)
ry12
## 関数にしておこう
spcor <- function(x,y,el) {
  rxe <- cor(x,el)
  ryx <- cor(x,y)
  rye <- cor(el,y)
  cat("xからelの影響を除き、そのxとyの相関\n")
  return((ryx -rye*rxe)/sqrt(1-rxe^2))
}
spcor(x2, y, x1)
ry21

# 偏相関 (partial correlation)
## x2からx1の影響を除き、yからもx1の影響を除き、そのうえでx2とyとの相関をとる
pry21 <- (ry2-ry1*r12)/((sqrt(1-ry1^2))*(sqrt(1-r12^2)))
pry21
## 部分相関から偏相関を求める
ry21
pry21*(sqrt(1-ry1^2))
## 偏相関から部分相関を求める
pry21
ry21/sqrt(1-ry1^2)
## 関数にしておこう
pcor <- function(x,y,el) {
  rxe <- cor(x,el)
  ryx <- cor(x,y)
  rye <- cor(el,y)
  cat("xからelの影響を除き、yからelの影響を除き、そのxとyの相関\n")
  return((ryx-rye*rxe)/((sqrt(1-rye^2))*(sqrt(1-rxe^2))))
}
pcor(x2, y, x1)
pry21