欠測データの補完(調査観察データの統計科学 2.2節)

「調査観察データの統計科学」
http://www.amazon.co.jp/dp/4000069721/
の2.2節「欠測のメカニズム」にあった、例2.3「ランダムな欠測としての選抜効果」の例をRで試してみた。

本書P31の定数の決め方だと、P32のグラフのようにならないので、を求めるときの切片は無視してみた。

ソースは以下の通り。

# データの設定
N <- 1000
theta1 <- 50
theta2 <- 0.8
sigma <- 49
C <- 60

set.seed(1)
y1 <- theta1 + rnorm(N, 0, sqrt(100))
y2 <- theta2 * y1 + rnorm(N, 0, sqrt(sigma))
pcol <- ifelse(y1 >= C, "navy", "grey")
data <- data.frame("y1" = y1, "y2" = y2)


# データプロット
plot(y1,y2,
		type = "p",
		col = pcol,
		xlim = c(0,100),
		ylim = c(0, 100))
par(new = T)
abline(v = C, col = "black")


# 観測データの抽出
y11 <- subset(data, y2 >= C)$y1
y21 <- subset(data, y2 >= C)$y2
# 観測データの相関係数
cor1 <- cor(y11, y21, method = "pearson")


# 観測データの線形回帰によるパラメータ推定
ylm <- lm(y21~-1 + y11,data = data)
summary(ylm)
theta2_es <- coefficients(ylm)


# 欠測データ補完後の相関係数
cor2 <- (theta2_es * var(y1))/(sqrt(var(y1)) * sqrt(theta2_es^2 * var(y1) + sigma))


# データプロット
plot(y1,y2,
		type = "p",
		col = pcol,
		xlim = c(0,100),
		ylim = c(0,100))
par(new = T)
abline(v = C, col = "black")
par(new = T)
segments(0,0,100,100*cor1, col="green")
par(new = T)
segments(0,0,100,100*cor2, col="red")
legend(locator(1),c("欠測データのみ","欠測補完後"),lty = c(1,1), col=c("green","red"))

たしかに、欠測があるとして計算した方が、相関係数が高くなっている。