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)