線形回帰の95%CI:正規分布

lmやglmの予測値をプロットし、それらの信頼区間を合わせて表示したいケースは多いかと思います。正規分布のモデルであればpredict()の引数を調整することで簡単に描けますが、それ以外の確率分布となったときに案外描けなかったりします。以下のスクリプトでは正規分布を例に示していますが、この手順を覚えればほかの確率分布のモデルにも適用できます(多少の変更は必要)。プロットの基本的な調整でつくったサンプルデータを使ってスクリプトを書きます。以下のスクリプトから右図が描けます。

#1. 統計モデルの構築
model <- lm(y~x)

#2. モデルの予測値を示すための説明変数のデータフレームを作る
new <- data.frame(x = seq(min(x),max(x),length=100))
#予測用データフレーム内の説明変数の名称はモデルの変数名と一致する必要あり
#seq(min, max, length=サンプル数)となっており、最小値から最大値までの範囲で等間隔のベクターを生成

#3. predict()をつかって予測値を導く
pred <- predict(model,new,se.fit=T)

y_pred <- pred$fit #予測値の平均値
up <- pred$fit + 1.96*pred$se.fit #Upper 95%CI:平均予測値に1.96*SEを加える
low <- pred$fit - 1.96*pred$se.fit #Lower 95%CI

#4. プロット
plot(y~x, type="n", ann=F, axes=F)
box(bty="l")
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
lines(new$x, y_pred) # 平均予測値のラインを描く
lines(new$x, up, lty=2); lines(new$x, low, lty=2)# 95%CIを描く
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish density", side = 2, line = 3)

おまけ:matplot()を使った方法

Y <- cbind(up,y_pred,low)
matplot(new$x,Y, type="l", lty=c(2,1,2), col=rgb(0,0,0), ann=F, axes=F)
box(bty="l")
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish density", side = 2, line = 3)

なお、複数の説明変数がある重回帰モデルの場合、注目する説明変数以外は平均値(もしくは中央値)に固定したデータセットを作る必要があります(そうしない限り、ほかの変数の影響を受けた予測値になる)。以下にサンプルを示します。

model <- lm(y~x1+x2+x3)
new_x1 <- seq(min(x1),max(x1),length=100)
new <- data.frame(x1 = new_x1, x2 = median(x2), x3 = median(x3))

SPONSOR LINK

Spread the love

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です