２次元正規分布のデータから２通りの方法で回帰直線を引く

タイトルの通りです。

• lmで単回帰直線
• 平均と分散共分散からE[Y|X=x] *1を計算

の２通りで回帰直線を引きます。

コード

library(MASS)
library(scales)
options(digits = 3)

# 平均
mx <- 10 
my <- 22
# 分散
Sx <- sqrt(9)
Sy <- sqrt(16)
Sxy <- sqrt(9.6)
# サンプルサイズ
n <- 10000

# 平均ベクトル
M <- c(mx, my)

# 分散共分散行列
Sig <- c(Sx^2,Sxy,Sxy,Sy^2) %>% matrix(2,2)
z <- MASS::mvrnorm(n, mu = M, Sigma = Sig )  

# dev.off()
# 散布図プロット
plot(z[,2] ~ z[,1], xlab = "x", ylab = "y",col = c(alpha('orange', 0.2)) )

# 単回帰直線(青)
lm.z <- lm(z[,2] ~ z[,1])
coef(lm.z)
abline(lm.z, col ="blue", lwd = 2, lty = 2)

# E[Y|X = x]で直線を引く(赤)
list("Intercept" = my - mx * (Sxy/(Sx*Sy))*(Sy/Sx), "Slope" = (Sxy/(Sx*Sy))*(Sy/Sx) ) %>% unlist()
abline(my - mx * (Sxy/(Sx*Sy))*(Sy/Sx), (Sxy/(Sx*Sy))*(Sy/Sx), col="red", lwd = 3, lty = 2)

結果

だいたい一致しました。青が単回帰直線、赤がE[Y|X = x]です。

• lmでの係数

> coef(lm.z)
(Intercept)      z[, 1] 
     18.552       0.344

• E[Y|X = x]での係数

> list("Intercept" = my - mx * (Sxy/(Sx*Sy))*(Sy/Sx), "Slope" = (Sxy/(Sx*Sy))*(Sy/Sx) ) %>% unlist()
Intercept     Slope 
   18.557     0.344 

参考

改訂版 日本統計学会公式認定 統計検定2級対応 統計学基礎 のp.86 (2.9.18)

*1:母数から計算しているのでおそらく母回帰直線、と思われる。