線形回帰の95%CI:ポアソン分布

ポアソン分布は、非負の離散値(カウントデータなど)に対して適用されることの多い確率分布です。ポアソン分布を仮定した線形モデルでは、母数である平均ラムダが負の値をとれないため、リンク関数にlogを使います。そのため、predict()の返値は対数スケールの値になっており、予測値曲線をプロットする際にはもとのスケールに戻す必要があります。記事末尾のサンプルデータを使い、ポアソンモデルの予測曲線をプロットします(右図)。

#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()で変換する
# SEの演算は指数変換を施す前に行う
# ダメな例:exp(pred$fit) + exp(1.96*pred$se.fit)

#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 abundance", side = 2, line = 3)

基本的なスクリプトは正規分布の場合と変わりませんが、辺値を変換するステップが必須となります。また、予測値の95%CIを出す際に、exp()で変換後に演算すると負の値が現れたりしますので注意が必要です。

※注 predict(model, new, se.fit=T, type=”response”)とすると、実測値スケールで予測値を返してくれますが、このSEをもとに四則演算を行うと上記の「ダメな例」と同様のことをすることになります。

N <- 50; a <- 1.4; b <- 0.7
x <- rnorm(N,0,1)

lambda <- y <- NULL
for(i in 1:N){
  lambda[i] <- exp(a + b*x[i])
  y[i] <- rpois(1,lambda[i])
}

SPONSOR LINK

Spread the love

コメントを残す

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