線形回帰の95%CI:ポリゴンver.

線形モデルの95%CIをグレイのシェードで描きたいときは、polygon()を使います。ポアソンモデルを例に、グレイシェイドの95%CIを描きます(下図)。

#1. 統計モデルの構築
model <- glm(y~x, family=poisson)

#2. モデルの予測値を示すための説明変数のデータフレームを作る
new <- data.frame(x = seq(min(x),max(x),length=100))

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

y_pred <- exp(pred$fit) #予測値の平均値
up <- exp(pred$fit + 1.96*pred$se.fit) #Upper 95%CI:平均予測値に1.96*SEを加える
low <- exp(pred$fit - 1.96*pred$se.fit) #Lower 95%CI
# リンク関数がlogなので、predict()の返値をlogの逆関数exp()で変換する

#4. polygon()のx,yの値を設定する
X <- c(new$x, rev(new$x))
Y <- c(low, rev(up))
# rev()はベクターの並び順を逆転させる

#5. プロット
plot(y~x, type="n", ylim=c(min(Y,y),max(Y,y)),
     ann=F, axes=F)
box(bty="l")
polygon(X,Y,border = F, col=rgb(0,0,0,alpha=0.3))
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
lines(new$x, y_pred) # 平均予測値のラインを描く
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish abundance", side = 2, line = 3)

polygon()を後のほうに描いてしまうと、先に描いたポイントやラインも塗りつぶしてしまうので、先に描いておきます。また、プロットのy軸の範囲をylim=c(min(Y,y),max(Y,y))とすることで、データポイントおよび95%CIの最大・最小値をカバーできる描画領域を指定しています。

※polygon()補足
polygon()は、グラフ上のXY座標を指定し、その座標に囲まれた範囲を塗りつぶします。例えば、次のようなデータフレームがあるとします。

#ポリゴンの座標
X <- c(1,2,2,1); Y <- c(1,1,2,2)
data.frame(X,Y)
  X Y
1 1 1
2 2 1
3 2 2
4 1 2

これは、座標を(X, Y)と表記したとき、データ1が(1, 1)を示し、データ4が(1, 2)を示す、という意味なります。これをpolygon()に書かせてみます(右図)。polygon()では、第一引数でX座標を指定し、第二引数でY座標を指定します。

#ポリゴンを描画するプロット領域
plot(0, type="n", xlim=c(0,3),ylim=c(0,3),
     xlab="X", ylab="Y")
polygon(X,Y,col="black")
# colは塗りつぶしの色指定;borderは縁どりの有無を指定

SPONSOR LINK

Spread the love

コメントを残す

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