Logistic迴歸與R

一. 廣義線性模型

    R 中提供了擬合計算廣義線性模型的函數 glm(), 命令以下:

fitted.model <- glm(formula, family = family.generator, data = data.frame)

其中 formula 是擬合公式, , family 是分佈族, 即廣義線性模型的種類, 如 正態分佈,

Poisson分佈, 二項分佈 等, data 是數據框.

 

二. 二項分佈族

    在二項分佈族中, Logistic迴歸是最重要的模型. 在該回歸模型下, 響應變量 Y 是分類的,

定性變量, 常常是成功 或 失敗.

    對於響應變量 Y, 有 p 個自變量, 記爲 X1, X2, · · · , Xp. 在 p 個自變量的做用下出現成功的

條件機率記爲 P = P {Y = 1|X1, X2, · · · , Xp}, 那麼Logistic迴歸模型爲

    P = exp(Z) / 1 + exp(Z)

    Z = β0 + β1X1 + β2 + · · · + βpXp

    P = exp(β0 + β1X1 + β2 + · · · + βpXp) / 1 + exp(β0 + β1X1 + β2 + · · · + βpXp)

其中 β0 爲常數項或截距, 稱 β1, β2, · · ·, βp 爲Logistic迴歸係數

    從公式能夠看出, Logistic迴歸模型是一個非線性迴歸模型, 自變量 Xj( j = 1,2, · · · , p)

是連續變量, 也能夠是分類變量, 或啞變量 ( dummy variable ) 對自變量 Xj 任意取值,

Z 在 −∞ 到 +∞ 變化時, 公式的比值總在 0 到 1 之間變化, 這正是機率 P 的值域.

 

• 對公式做Logistic變換, Logistic迴歸模型能夠變成下列線性形式:

    logit(P) = ln( P / 1 - P ) = Z = β0 + β1X1 + β2 + · · · + βpXp

    從上式能夠看出, 咱們可以使用線性迴歸模型對參數進行估計, 這就是Logistic迴歸模型

屬於廣義線性模型的緣由.

    Logistic迴歸模型的公式爲:

fitted.model <- glm(formula, family = binomial, data = data.frame)

 

三. R. Norell 實驗

    解: 用數據框形式輸入數據, 再構造矩陣, 一列是成功( 響應 )的次數, 一列是失敗

( 不響應 )的次數,

而後再做Logistic迴歸.

    P = exp(Z) / 1 + exp(Z)

    Z = β0 + β1X

    logit(P) = ln( P / 1 - P ) = Z = β0 + β1X

代碼以下:

norell <- data.frame(
    x = 0 : 5, n = rep(70, 6), success = c(0, 9, 21, 47, 60, 63)
)

norell$Ymat = cbind(norell$success, norell$n - norell$success)

norell.glm <- glm(Ymat ~ x, family = binomial, data = norell)

summary(norell.glm)

 

Call:
glm(formula = Ymat ~ x, family = binomial, data = norell)

Deviance Residuals: 
      1        2        3        4        5        6  
-2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
norell$x      1.2459     0.1119   11.13   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 250.4866  on 5  degrees of freedom
Residual deviance:   9.3526  on 4  degrees of freedom
AIC: 34.093

Number of Fisher Scoring iterations: 4

即 β0 = −3.3010, β1 = 1.2459

Logistic迴歸模型爲:

    P = exp(−3.3010 + 1.2459X) / 1 + exp(−3.3010 + 1.2459X)

Ｘ : 電流強度(單位毫安)

 

在獲得Logistic迴歸模型後,

• 作預測

列如, 當電流爲 3.5 毫安時, 牛響應的機率爲多少?

pre <- predict(norell.glm, data.frame(x = 3.5))
exp(pre) / (1 + exp(pre))

0.742642 

即, 當電流爲 3.5 毫安時, 牛響應的機率爲74.26 %

 

• 作控制

列如, 當 50 % 的牛有響應, 電流強度 X 爲多少? 當 P = 0.5 時, ln( P / 1 - P ) = 0,

因此, X = - β0 / β1

X <- norell.glm$coefficients[[1]] / norell.glm$coefficients[[2]]

X = 2.649439 

即, 當電流爲 2.65 毫安時, 50 % 的牛有響應

 

• 響應比與Logistic迴歸曲線

# 電流強度, 橫座標的點, 0 : 5, 均勻分佈 100 個
d <- seq(0, 5, len = 100)
# 計算預測值
pre <- predict(norell.glm, data.frame(x = d))
# 響應比, 預測機率
p = exp(pre) / (1 + exp(pre))

norell$y = norell$success / norell$n
# 散點圖
plot(norell$x, norell$y)
# 預測曲線
lines(d, p)