# 羅吉斯回歸 (logistic regression)
###### tags: `R教學` `基礎統計`

**範例**：半致死率 (LD50) 的運算

 羅吉斯回歸是最早被發展出來的廣義線性回歸方程式，也是了解廣義線性回歸意義的最簡單案例。廣義線性回歸的概念就是把先用一個連接方程式把非線性的數據先轉換成線性之後，再以線性模型的方法做檢定。在R中的程式碼為glm()。
  羅吉斯回歸處理的目標是二項式分布的數據，就是平常我們最常遇到的機率數據。機率的上下界介於0 (0%) 與1 (100%) 之間，但是線性函數是沒有上下極限的，因此得到的資料會先藉由 link function 轉換成線性之後再以線性模型做分析。羅吉斯回歸的 link function 為：
$$ y = log_e (\frac {p}{1-p}) $$
藉此方程式將數據轉換成直線回歸的形式，因此轉換後的數值會符合 $y=\beta_0+ \beta_1 x$ 的動態。因此若要回頭把估計出來的直線轉換回原始的數值，就把估計出來的斜率 ($\beta_1$) 與截距 ($\beta_0$)上述公式的反函數來達成：
$$ p(x) = \frac {1}{1-e^{-(\beta_0 + \beta_1 x)}} $$
而這兩個公式也是本次的練習中會用到的轉換公式。

rm(list=ls())
library(MASS)

dat <- data.frame(
  dose = c(0, 0.1, 1, 10, 100, 1000, 10000, 100000),
  survival = c(100, 98, 92, 73, 61, 29, 8, 7),
  death =    c(0,    2,  8, 27, 39, 71, 92, 93)
)
dat$rate = dat$death / (dat$survival + dat$death)
dat_M <- dat
dat_M$dose <- dat_M$dose + 0.000001

plot(rate ~ log(dose), 
     data=dat_M)

G <- glm(as.matrix(dat_M[,3:2]) ~ log(dose), 
         ,family=binomial(link=logit)
         ,data=dat_M)
co <- G$coef
se.co <- summary(G)$coef[,2]
cov.co <- summary(G)$cov.scaled[1,2]
LD50 <- -co[1]/co[2]
LD90 <- (log(9)-co[1])/co[2]

se.ld50 <- abs(LD50) * sqrt( (se.co[1]/co[1])^2 +
                               (se.co[2]/co[2])^2 -
                               2*cov.co/(co[1]*co[2]) )

MASS::dose.p(G,p=c(0.5))
LD50 + c(-1,1) * se.ld50


jpeg(filename = "LD50.jpg",
     width = 510*2, height = 380*2, pointsize = 20, quality = 75)
plot(rate ~ log(dose), data=dat_M, 
     las=1, pch=19, xaxt="n", ylab = "Death rate", xlab = "dose")
axis(1, at=log(dat_M$dose), labels =dat$dose)
#axis(1, at=log(dat_M$dose), labels = formatC(dat_M$dose, digits = 1, format = "e"))
curve(1/(1+exp(-co[1]-co[2]*x)),add=TRUE,lwd=1)

points(LD50, 0.5, pch = 19, col = 2)
segments (-15, 0.5, LD50, 0.5, lty = 2)
arrows(LD50-se.ld50, 0.5, LD50+se.ld50, 0.5,
       code=3, length=0.05, angle=90, col=2)
dev.off()