A way of thinking

筆者個人の思考過程です。意見には個人差があります。

二項分布をfamilyとするGLMで確率が0%や100%はモデル予測に使われていないのか?という素朴な疑問(例えば,応答変数が生存率だった場合,10匹中0匹生き残ったとか,10匹とも生き残ったという状況)。logit変換の式を見ると,確率pが1または0の時は分母が0になるので,オッズが計算できない*1。結果として(回帰結果が変わるので),パラメータ推定に使っている模様。まぁよくよく考えると,下のprobaの式の状態で探せば,問題ないということか。しかしまぁ,まだこういうシンプルな話ですら理解してないっことが身にしみます。

# 二項分布のGLMで0や100%は推定に使われていないのか?という疑問
# 簡単なコードと疑似データで調べてみる。
set.seed(99)
# 化学物質の濃度のようなxを想定
x <- seq(1, 5, length=20)
a <- 5
b <- -1.5
# 応答変数は生存率とする
(proba <- 1/(1 + exp(-(a + b*x))))

N.sur<- rbinom (length(proba), 10, proba)
N.dead <- 10-N.sur

(d0 <- data.frame(x, N.sur, N.dead))

m.all <- glm(cbind(N.sur, N.dead) ~ x, d0, family=binomial)
summary(m.all)

# 生存数0のデータを除いてみる
d1 <- d0[-which(d0$N.sur == 0),]
m1 <- glm(cbind(N.sur, N.dead) ~ x, d1, family=binomial)
summary(m1)


# 全部生存のデータを除いてみる
d2 <- d0[-which(d0$N.sur == 10),]
m2 <- glm(cbind(N.sur, N.dead) ~ x, d2, family=binomial)
summary(m2)

*1:これらが上の疑問と関係あるかよくわからないのですが,まぁやってみたという感じです