检查逻辑回归中的线性关系。

3

为了检查逻辑回归中的线性关系 ->

independent1independent2变量是否与depdendent的对数几率呈线性相关?

我想优化这个(工作中的)计算:

这是代码:

# Check Linearity  ---------------------------------------------------------

# quartiles of independent1
quantile(df$independent1, probs=c(0, 0.25, 0.5, 0.75, 1))


table(df$dependent[df$independent1<52])
table(df$dependent[df$independent1>=52 & df$independent1 < 60])
table(df$dependent[df$independent1>=60 & df$independent1 < 73])
table(df$dependent[df$independent1>=73 & df$independent1 < 91])

p1 <- mean(df$dependent[df$independent1<52])
p2 <- mean(df$dependent[df$independent1>=52 & df$independent1 < 60])
p3 <- mean(df$dependent[df$independent1>=60 & df$independent1 < 73])
p4 <- mean(df$dependent[df$independent1>=73 & df$independent1 < 91])
probs <- c(p1, p2, p3, p4)

# calculate the log-odds
logits <- log(probs/(1-probs))

# quartiles of independent1
q <- quantile(df$independent1, probs=seq(0,1,0.25))

# calculate median independent1 for each of the 4 groups
meds <- c( median(df$independent1[ df$independent1<q[2]]),
           median(df$independent1[ df$independent1>=q[2] & df$independent1<q[3]]),
           median(df$independent1[ df$independent1>=q[3] & df$independent1<q[4]]),
           median(df$independent1[ df$independent1>=q[4]])
)

plot(meds, logits, main="xxx",
     xlab = "independent1",
     ylab = "log-odds(dependent|independent1)", las=1)

对于一个变量来说,这可能还可以。但是我有更多的自变量。那么我该如何针对每个自变量(在这个例子中是independent1independent2)优化这段代码(检查和绘图)?

我的数据框:

df <- structure(list(dependent = c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), independent1 = c(84, 
49, 54, 75, 49, 70, 75, 42, 60, 72, 80, 73, 51, 61, 59, 78, 45, 
38, 78, 65, 91, 60, 39, 31, 42, 72, 41, 77, 73, 74, 39, 86, 71, 
55, 43, 75, 80, 75, 67, 74, 46, 70, 57, 66, 57, 72, 46, 52, 53, 
76, 57, 86, 67, 71, 57, 50, 76, 61, 41, 57, 62, 41, 64, 82, 53, 
75, 59, 38, 54, 56, 68, 63, 73, 26, 75, 76, 81, 46, 77, 53, 59, 
66, 51, 72, 80, 70, 39, 57, 62, 85, 84, 57, 73, 55, 70, 78, 66, 
69, 60, 51, 72, 68, 60, 62, 64, 44, 50, 59, 45, 81, 54, 68, 75, 
66, 54, 45, 52, 87, 44, 77, 49, 84, 68, 76, 82, 44, 58, 55, 69, 
33, 48, 62, 60, 76, 56, 73, 55, 58, 53, 53, 60, 52, 60, 41, 39, 
36, 38, 59, 54, 64), independent2 = c(23, 25, 34, 25, 31, 25, 
32, 19, 25, 28, 22, 18, 30, 26, 25, 25, 25, 19, 24, 27, 23, 28, 
39, 27, 30, 28, 22, 28, 25, 23, 18, 27, 27, 19, 25, 27, 26, 26, 
21, 26, 23, 28, 37, 32, 24, 32, 26, 23, 24, 27, 28, 25, 24, 22, 
34, 23, 35, 20, 29, 29, 21, 29, 25, 26, 23, 33, 25, 26, 29, 27, 
26, 28, 19, 22, 29, 22, 26, 35, 32, 29, 26, 23, 31, 30, 27, 28, 
23, 27, 34, 22, 24, 28, 21, 25, 18, 32, 21, 24, 31, 31, 24, 30, 
27, 23, 16, 26, 26, 19, 38, 21, 32, 34, 28, 19, 30, 24, 26, 24, 
40, 26, 15, 26, 28, 22, 25, 26, 31, 24, 26, 42, 26, 30, 28, 21, 
21, 19, 22, 20, 26, 31, 22, 25, 21, 20, 27, 27, 26, 29, 22, 24
)), row.names = c(NA, -150L), class = c("tbl_df", "tbl", "data.frame"
))

3
这个方法的统计功效不会很大。你可能想要加入一个二次多项式项。 - IRTFM
@IRTFM 感谢您的快速回复。我的目的是逐个变量检查逻辑回归变量的所有假设。对于连续变量,如果假设不成立,我将对变量进行二值化处理。 - TarJae
4
哎呀,统计学有效性的威胁变得更加严重了。将呈现非线性关联证据的变量二分化会隐藏(而不是调查)可能有用的信息。如果你采纳了我的建议,你可以保留满足你选择的“非线性”规则的二次项。 - IRTFM
@IRTFM,你能给我一个答案吗?我会非常感激的。 - TarJae
与 IRTFM 所述的统计功率无关。但我正在寻找一种编码解决方案(更直接的方式)来检查逻辑回归中的线性关系。 - TarJae
1个回答

1

我将演示一种有些不同且明显更有效的方法来拆分一个将用于逻辑回归模型的变量:

df$q41 <- with(df, cut(independent1, quantile(independent1), include = TRUE))
 # creates 4 level factor into roughly equally sized groups
 table(df$q41)
#--------------------
#[26,52] (52,60] (60,73] (73,91] 
#     39      37      39      35 

#Examine for "eyeball" trends in the log-odds of dependent 
fit1.q41 <- glm(dependent~q41+0, data=df, fam="binomial")
fit1.q41
#---------------------------
Call:  glm(formula = dependent ~ q41 + 0, family = "binomial", data = df)

Coefficients:
q41[26,52]  q41(52,60]  q41(60,73]  q41(73,91]  
    -3.638      -2.862      -2.918      -2.048  

Degrees of Freedom: 150 Total (i.e. Null);  146 Residual
Null Deviance:      207.9 
Residual Deviance: 65.52    AIC: 73.52

我选择移除截距项,因为它的存在阻碍了将最低组的系数与上三个组在同一刻度下观察。这些系数只是我创建的分组的逻辑回归值。对比如下:

> logits
[1] -3.555348 -2.740840 -2.970414 -2.169054
> coef(fit1.q41)
q41[26,52] q41(52,60] q41(60,73] q41(73,91] 
 -3.637586  -2.862201  -2.917771  -2.047693 

我随后尝试自动化这个过程,但遇到了一个问题,因为一个四分位组中的事件数量很少。在 independent2 中最低四分位数的极低系数是由于该类别中缺乏任何事件或 "1" 导致的。(估计对数几率为 -19.566069 确实指向了比例为 0。)

 lapply( df[-1], function(x){cat(str(x)); IVq <- cut(x, quantile(x), include = TRUE); logits<-coef( summary(glm(df$dependent~IVq+0, fam="binomial"))); logits})
 num [1:150] 84 49 54 75 49 70 75 42 60 72 ...
 num [1:150] 23 25 34 25 31 25 32 19 25 28 ...
$independent1
            Estimate Std. Error   z value     Pr(>|z|)
IVq[26,52] -3.637586  1.0130639 -3.590678 3.298191e-04
IVq(52,60] -2.862201  0.7270292 -3.936845 8.256004e-05
IVq(60,73] -2.917771  0.7259663 -4.019155 5.840732e-05
IVq(73,91] -2.047693  0.5312796 -3.854266 1.160776e-04

$independent2
             Estimate   Std. Error     z value     Pr(>|z|)
IVq[15,23] -19.566069 1639.9716035 -0.01193074 9.904809e-01
IVq(23,26]  -3.091042    0.7229988 -4.27530783 1.908734e-05
IVq(26,28]  -2.397895    0.7385489 -3.24676555 1.167245e-03
IVq(28,42]  -1.856298    0.4808846 -3.86017349 1.133066e-04


> lapply( df[-1], function(x){ IVq <- cut(x, quantile(x), include = TRUE); table(IVq, df$dependent) })
$independent1
         
IVq        0  1
  [26,52] 38  1
  (52,60] 35  2
  (60,73] 37  2
  (73,91] 31  4

$independent2
         
IVq        0  1
  [15,23] 43  0
  (23,26] 44  2
  (26,28] 22  2
  (28,42] 32  5

无论如何,我认为我已经展示了一种更具R风格的方法来计算四分位数内的logits。这也为您提供了一种模型比较方法来检查非线性偏差,并展示可能存在的陷阱。如果您有更多事件,您可能会考虑观察四分位因子在简单线性模型之上对空模型偏离的偏差变化......甚至使用poly更有力地创建您的比较模型。
在过去,当处理具有足够事件数量的数据集时,我选择基于从event==1子集中计算出的分位数进行拆分,而不是让拆分基于整个数据集。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接