線形回帰の95%CI:二項分布

二項分布は、「複数回の試行のうち、何回成功したか」のようなデータ(上限値が分かっている非負の離散値)に対して適用されることが多い確率分布です(ベルヌーイ分布の一般化)。生態学のデータでは、種の在/不在データに適用されることが多いと思います。二項分布のパラメーターは、ある事象Xが生じる確率pであり、pは0-1の範囲に収まる必要があります。そのため、リンク関数にはlogit関数(logit(p) = log(p/(1-p)))が用いられます。二項分布モデルの場合、predict()の返値をlogitスケールになっていますので、予測値をプロットするためには逆ロジット関数(inverse logit: ilogit(x) = exp(x)/(1+exp(x))で変換する必要があります。その他、注意点などについてはポアソン分布の場合と同様です。サンプルデータは記事末尾にあります。

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

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

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

#4. inverse logit functionを作る
ilogit <- function(x){ exp(x)/(1+exp(x)) }
# Rのデフォルト関数にilogit関数はない

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

#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("Probability of fish present", side = 2, line = 3)
ilogit <- function(x){ exp(x)/(1+exp(x)) }
N <- 50
a <- 0.4; b <- 1.5
x <- rnorm(N,0,1)

p <- y <- NULL
for(i in 1:N){
  p[i] <- ilogit(a + b*x[i])
  y[i] <- rbinom(1,1,p[i])
}

SPONSOR LINK

Spread the love

コメントを残す

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