独習統計学24講 - 驚異のアニヲタ社会復帰の予備

MikuHatsune2017-08-25

読んだ。

作者: 鶴田陽和
出版社/メーカー: 朝倉書店
発売日: 2013/05/29
メディア: 単行本
この商品を含むブログ (3件) を見る

すべての医療系学生・研究者に贈る独習統計学応用編24講 ─分割表・回帰分析・ロジスティック回帰─

作者: 鶴田陽和
出版社/メーカー: 朝倉書店
発売日: 2016/11/02
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

COI：ラボにあった。
　
分布や極簡単な線形代数的計算を使いながら、文系や初学者でも統計が勉強できる、ということで実例もいくつか具体的にやりながらやっていた。
　
応用編で、いままでよくわかってなかったけどなんかよくわかったような気になったのが、回帰分析を行ったあとの信頼区間についてである。ここで、単純な単回帰 $Y=\alpha + \beta X$ を考えるとき、次の4つの場合について推定を考えている。
1：回帰係数 $\alpha$ , $\beta$ の信頼区間と仮説検定
2：ある $X_0$ に対する「回帰直線の $Y$ 座標 $\eta_0=\alpha + \beta X_0$ の信頼区間
3：目的変数 $Y$ の信頼区間（「予測区間」と呼んでいる）
4 回帰式 $Y=\alpha + \beta X$ の信頼区間（=回帰直線が通る範囲）
　
適当にデータを作成する。

x <- 1:10
b0 <- c(0, 0.5)
eps <- rnorm(length(x), mean=0, sd=b0[2]) # 誤差
Y <- b[1] + b[2] * x # 真の回帰
y <- Y + eps         # 実際に得られるデータ

lm0 <- lm(y ~ x)     # 回帰推定

plot(x, Y, type="l", lty=3)
points(x, y, pch=15)
abline(lm0)

summary(lm0)

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.4143 -0.2917 -0.2228  0.2794  0.7716 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.10317    0.29861   0.345    0.739    
x            0.46034    0.04813   9.565 1.18e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4371 on 8 degrees of freedom
Multiple R-squared:  0.9196,	Adjusted R-squared:  0.9095 
F-statistic:  91.5 on 1 and 8 DF,  p-value: 1.182e-05

　
1 の、回帰モデルのなかでなんの説明変数が有意に効いているか効いていないかは、各係数の推定値と信頼区間をみればよい。ここでは、傾きは0.5 で設定していて、有意にはなっている。一方で切片は0 にしたので、有意に0 ではないかというとそうではない結果になっている（かと言って、回帰の結果を見て本当に切片が0 かは本当はわからない）。
　
2 は、 $X$ のなかのある特定の値 $X_0$ に対する回帰直線の $Y$ 座標、ということで求めるらしい。
　
3 は、母集団のパラメータではなく変量 $Y$ の値の推測をしたいという状況で、推定（estimation, 母数に対して）ではなく予測（prediction, 変量に対して）という。これは、2 に対して更に誤差が入ったものなので、2 より広い。
　
4 は、回帰直線全体としてみたときに線が通る範囲である。2 や3 はある $X$ だけを考えていたが、4 では $X$ 全体という制約で動かないといけないので、範囲としては狭いはず。回帰直線上の同時信頼区間になる、
ほとんどの統計ソフトは、1-3 の計算は実装されているが、4 を実装しているものは少ない、ということで、R で確かめてみる。
ちなみに1-4 の計算方法は以下の通りである。ここで $p$ は有意水準。

推定	パラメータ	信頼区間
回帰係数	$\beta$	$\hat{\beta}\pm t_{n-2}(p/2)\frac{\hat{\sigma}}{\sqrt{S_{xx}}}$
回帰直線の $Y$ 座標	$\eta_0$	$(\hat{\alpha}+\hat{\beta}X_0)\pm t_{n-2}(p/2)\sqrt{\frac{1}{n}+\frac{(X_0-\bar{X})^2}{S_{xx}}}\sigma$
目的変数	$Y_0$	$(\hat{\alpha}+\hat{\beta}X_0)\pm t_{n-2}(p/2)\sqrt{1+\frac{1}{n}+\frac{(X_0-\bar{X})^2}{S_{xx}}}\sigma$
回帰直線	$Y=\alpha+\beta X$	$(\hat{\alpha}+\hat{\beta}X_0)\pm \sqrt{2F_{2, n-2}(p)}\sqrt{\frac{1}{n}+\frac{(X_0-\bar{X})^2}{S_{xx}}}\sigma$

ここで、 $\hat{\sigma}$ は誤差の不偏分散（ $n-2$ になる）である。
　

sig <- 0.05                         # 有意水準
Sxx <- sum((x - mean(x))^2)         # 平方和
Sxy <- sum((x-mean(x))*(y-mean(y))) # xy の積和
b <- Sxy/Sxx                        # 回帰直線の傾きは最小二乗法からこれで得られる
a <- mean(y) - mean(x) * b          # 切片
n <- length(x)                      # 標本数
est <- c(cbind(1, x) %*% c(a, b))   # 推定されたY
resid <- y - est                    # 残差
delta <- sqrt(sum((resid-mean(resid))^2)/(n-2)) # 本当はsigma
tdist <- qt(1-sig/2, n-2)           # t 分布

# 係数の信頼区間
conf_interval <- tdist * delta / sqrt(Sxx)

# ita にあたる信頼区間
Y0_interval <- tdist * sqrt(1 + 1/n + ((x-mean(x))^2)/Sxx) * delta

# 目的変数の信頼区間
Y_interval <- tdist * sqrt(1/n + ((x-mean(x))^2)/Sxx) * delta

# 回帰直線の信頼区間
reg_interval <- sqrt(2*qf(sig, 2, n-2)) * sqrt(1/n + ((x-mean(x))^2)/Sxx) * delta

intervals <- list(Y_interval, Y0_interval, reg_interval)
yy <- mapply(function(z) outer(z, c(1, -1), "*") + est, intervals, SIMPLIFY=FALSE)

px <- c(x, rev(x))
cols <- c(grey(0.8), "skyblue", "hotpink")
par(mar=c(4.5, 5, 2, 2), cex.lab=2, cex.axis=2, las=1)
plot(x, Y, type="n", ylim=range(yy))
for(i in c(2,1,3)){
  polygon(px, c(yy[[i]][,1], rev(yy[[i]][,2])), border=NA, col=cols[i])
}
lines(x, Y, col="black", lty=3, lwd=3)
points(x, y, pch=15, col="red")
legend("bottomright", 
       legend=c("真の線", "観測したデータ", "回帰直線のY 座標", "目的変数 Y0", "回帰直線"),
       col=c("black", "red", cols), pch=c(NA, rep(15, 4)), pt.cex=2.5, cex=1.5, lty=c(3, rep(NA, 4)), lwd=5)

　
ここで、R でこの回帰分析に対してpredict を区間出力で実行すると

predict(lm0, as.data.frame(x), interval="confidence")

         fit         lwr      upr
1  0.5635107 -0.02895294 1.155974
2  1.0238546  0.52137751 1.526332
3  1.4841984  1.06160441 1.906792
4  1.9445422  1.58493039 2.304154
5  2.4048860  2.08133054 2.728442
6  2.8652298  2.54167436 3.188785
7  3.3255737  2.96596185 3.685185
8  3.7859175  3.36332351 4.208511
9  4.2462613  3.74378425 4.748738
10 4.7066051  4.11414144 5.299069

となる。これは3 の目的変数 $Y$ の予測区間になるので、各 $X$ を独立に見て $Y$ を見たときの予測区間になる。そのため、回帰直線全体として通る区間をみたいのであれば、F分布を使って4 を地道に計算する必要がある。
　
てこ比（leverage）
外れ値のような極端な値や、回帰範囲の端の値は回帰の性能に影響しやすい。これをてこ比 $h_{ij}$ という、測定値 $Y_i$ が変化すると $\hat{Y}_i$ の推定にどれくらい影響するか定量する指標を考える。
$h_{ij}=\frac{1}{n}+\frac{(X_i-\bar{X})(X_j-\bar{X})}{S_{xx}}$
を正規化して2SD 超えるとかなら注意した方がいいらしい。