对数后的Logistic回归参数P值变化-R

如何解决对数后的Logistic回归参数P值变化-R

在R中计算逻辑回归时，我有一个问题，对我来说，这是没有意义的。我在模型中有一个参数，正数（分子量）。我有一个二进制响应变量，比方说A或B。我的数据表叫做df1。

str(df1)

data.frame':    1015 obs. of  2 variables:
 $ Protein_Class: chr  "A" "A" "A" "B" ...
 $ MW                  : num  47114 29586 26665 34284 104297 ...

我制作模型：

summary(glm(as.factor(df1[,1]) ~ df1[,2],family="binomial"))

结果是：

Call:
glm(formula = as.factor(df1[,family = "binomial")

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5556  -1.5516   0.8430   0.8439   0.8507  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  8.562e-01  1.251e-01   6.842  7.8e-12 ***
df1[,2]    -1.903e-07  3.044e-06  -0.063     0.95    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1239.2  on 1014  degrees of freedom
Residual deviance: 1239.2  on 1013  degrees of freedom
AIC: 1243.2

Number of Fisher Scoring iterations: 4

到目前为止，一切都很好。但是，当我取变量的对数时：

summary(glm(as.factor(df1[,1]) ~ log10(df1[,2]),family="binomial"))

Call:
glm(formula = as.factor(df1[,family = "binomial")

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8948  -1.4261   0.8007   0.8528   1.0469  

Coefficients:
                Estimate Std. Error z value Pr(>|z|)   
(Intercept)      -2.7235     1.1169  -2.438  0.01475 * 
log10(df1[,2])   0.8038     0.2514   3.197  0.00139 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1239.2  on 1014  degrees of freedom
Residual deviance: 1228.9  on 1013  degrees of freedom
AIC: 1232.9

Number of Fisher Scoring iterations: 4

p值已更改！怎么会这样？更重要的是，要使用哪个？我的理解是逻辑回归是基于等级的，而我所做的只是单调变换。注意，模型的AUROC曲线保持不变。转换过程中不会丢失零或负值。

我在这里错过了什么吗？有什么建议吗？

预先感谢

亚当

解决方法

有几件事情需要考虑。首先，您可以将搜索限制在1的一侧或另一侧。这就是降低x的幂-平方根，对数，逆等...-都具有相似的效果，但是在不同程度上。他们都吸收了大价值，散布了小价值。大于1的转换则相反，它们倾向于增大大值之间的价差，而减小小值之间的价差-所有这些通常都假设变量中没有非正值。那么，这确实是一个有关您想要什么样的转换以及之后的转换的问题-它必须要多么严重。

首先，您需要什么样的转换。我做了一些假数据来说明这一点：

library(dplyr)
library(tidyr)
library(ggplot2)
set.seed(1234)
x <- runif(1000,1,10000)
y.star <- -6 + log(x)
y <- rbinom(1000,plogis(y.star) )
df <- tibble(
  y=y,x=x,ystar=y.star)

接下来，由于这只是一个双变量关系，我们可以用一条黄土曲线将其绘制出来。不过，尤其是，我们想知道y相对于x的对数奇数。我们可以通过使用对数分位数函数qlogis()对黄土曲线的预测进行变换来做到这一点-这将概率取到对数形式。然后，我们可以进行绘制了。

lo <- loess(y ~ x,span=.75)
df <- df %>% mutate(fit = predict(lo),fit = case_when(
               fit < .01 ~ .01,fit > .99 ~ .99,TRUE ~ fit))
ggplot(df) + 
  geom_line(aes(x=x,y=qlogis(fit)))

这看起来像是类日志关系。然后，我们可以实现一些不同的转换并将其绘制-平方根，对数和负逆。

lo1 <- loess(y ~ sqrt(x),span=.5)
lo2 <- loess(y ~ log(x),span=.5)
lo3 <- loess(y ~ I(-(1/x)),span=.5)

df <- df %>% mutate(fit1 = predict(lo1),fit1 = case_when(
                      fit1 < .01 ~ .01,fit1 > .99 ~ .99,TRUE ~ fit1))
df <- df %>% mutate(fit2 = predict(lo2),fit2 = case_when(
                      fit2 < .01 ~ .01,fit2 > .99 ~ .99,TRUE ~ fit2))
df <- df %>% mutate(fit3 = predict(lo3),fit3 = case_when(
                      fit3 < .01 ~ .01,fit3 > .99 ~ .99,TRUE ~ fit3))

接下来，我们需要转换数据，以便绘图看起来正确：

plot.df <- df %>% 
  tidyr::pivot_longer(cols=starts_with("fit"),names_to="var",values_to="vals") %>% 
  mutate(x2 = case_when(
    var == "fit" ~ x,var == "fit1" ~ sqrt(x),var == "fit2" ~ log(x),var == "fit3" ~ -(1/x),TRUE ~ x),var = factor(var,labels=c("Original","Square Root","Log","Inverse")))

然后，我们可以进行绘图：

ggplot(plot.df,aes(x=x2,y=vals)) + 
  geom_line() + 
  facet_wrap(~var,scales="free_x")

在这里，日志似乎是最线性的-不足为奇，因为我们使用y.star来创建变量log(x)。如果我们想测试这些不同的可能性，罗切斯特的政治学家凯文·克拉克（Kevin Clarke）提出了配对符号测试，以评估非嵌套模型之间的差异。 here上有一篇论文。我编写了一个名为clarkeTest的程序包，该程序包在R中实现了此功能。因此，我们可以使用它来测试各种不同的选择：

m0 <- glm(y ~ x,data=df,family=binomial)
m1 <- glm(y ~ sqrt(x),family=binomial)
m2 <- glm(y ~ log(x),family=binomial)
m3 <- glm(y ~ I(-(1/x)),family=binomial)

针对平方根测试原稿：

library(clarkeTest)
> clarke_test(m0,m1)
# 
# Clarke test for non-nested models
# 
# Model 1 log-likelihood: -309
# Model 2 log-likelihood: -296
# Observations: 1000
# Test statistic: 400 (40%)
# 
# Model 2 is preferred (p = 2.7e-10)

这表明平方根比原始的未转换变量好。

clarke_test(m0,m2)
# 
# Clarke test for non-nested models
# 
# Model 1 log-likelihood: -309
# Model 2 log-likelihood: -284
# Observations: 1000
# Test statistic: 462 (46%)
# 
# Model 2 is preferred (p = 0.018)

以上显示日志比未转换的变量要好。

> clarke_test(m0,m3)
# 
# Clarke test for non-nested models
# 
# Model 1 log-likelihood: -309
# Model 2 log-likelihood: -292
# Observations: 1000
# Test statistic: 550 (55%)
# 
# Model 1 is preferred (p = 0.0017)

以上说明未转换的变量比负的逆变量更可取。然后，我们可以测试两种模型与原始模型之间的差异。

> clarke_test(m1,m2)
# 
# Clarke test for non-nested models
# 
# Model 1 log-likelihood: -296
# Model 2 log-likelihood: -284
# Observations: 1000
# Test statistic: 536 (54%)
# 
# Model 1 is preferred (p = 0.025)

这表明，就单个对数似然而言，平方根比对数转换好。

另一个选择是对可能的转换进行网格搜索，并每次查看AIC。首先，我们需要创建一个函数来处理变换功率= 0的情况，在这里我们应该用对数代替。然后，我们可以为每个不同的转换运行一个模型，并获得AIC。

grid <- seq(-1,by=.1)

trans <- function(x,power){
  if(power == 0){
    tx <- log(x)
  }else{
    tx <- x^power
  }
  tx
}
mods <- lapply(grid,function(p)glm(y ~ trans(x,p),family=binomial))
aic.df <- tibble(
  power = grid,aic = sapply(mods,AIC))

接下来，我们可以将AIC绘制为功率的函数。

ggplot(aic.df,aes(x=power,y=aic)) + 
  geom_line()

这告诉我们大约-.25是适当的转换参数。请注意，Clarke测试结果与AIC之间存在差异，因为AIC基于 overall 对数似然，而Clarke测试基于个体对数中的差异-可能性。

我们会发现，这种新提议的转换也比平方根差。

m4 <- glm(y ~ I(x^-.25),family=binomial)
clarke_test(m1,m4)
# 
# Clarke test for non-nested models
# 
# Model 1 log-likelihood: -296
# Model 2 log-likelihood: -283
# Observations: 1000
# Test statistic: 559 (56%)
# 
# Model 1 is preferred (p = 0.00021)

因此，如果您有几个不同的候选人，并且您喜欢Clarke测试背后的想法，则可以使用它来找到合适的转换。如果您不考虑候选人，那么总是可以进行网格搜索。

对数后的Logistic回归参数P值变化-R

如何解决对数后的Logistic回归参数P值变化-R

解决方法

相关推荐