如何使用带有 R 的 L-BFGS 执行逻辑回归？

Posted 2023-03-05

【中文标题】如何使用带有 R 的 L-BFGS 执行逻辑回归？

【英文标题】：How to perform logistic regression with L-BFGS with R?

【发布时间】：2015-07-21 16:30:23

【问题描述】：

我正在尝试使用带有 R 的 L-BFGS 执行逻辑回归。 这是我的数据集（390 obs。14 个变量，Y 是目标变量）

GEST    DILATE    EFFACE    CONSIS    CONTR    MEMBRAN    AGE    STRAT    GRAVID    PARIT    DIAB    TRANSF    GEMEL    Y
31           3       100         3        1         2     26         3         1        0       2         2       1     1
28           8         0         3        1         2     25         3         1        0       2         1       2     1
31           3       100         3        2         2     28         3         2        0       2         1       1     1
...

此数据集可在此处找到：http://tutoriels-data-mining.blogspot.fr/2008/04/rgression-logistique-binaire.html in "Données :prematures.xls"。 Y 是我在 Excel 中使用列“PREMATURE”创建的列，Y=IF(PREMATURE="positif";1;0)

我尝试过像https://stats.stackexchange.com/questions/17436/logistic-regression-with-lbfgs-solver 这样使用 optimx 包，代码如下：

install.packages("optimx")
  library(optimx)

vY = as.matrix(premature['Y'])
mX = as.matrix(premature[c('GEST','DILATE','EFFACE','CONSIS','CONTR','MEMBRAN','AGE','STRAT','GRAVID','PARIT','DIAB','TRANSF','GEMEL')])

#add an intercept to the predictor variables
mX = cbind(rep(1, nrow(mX)), mX)

#the number of variables and observations
iK = ncol(mX)
iN = nrow(mX)

#define the logistic transformation
logit = function(mX, vBeta) 
return(exp(mX %*% vBeta)/(1+ exp(mX %*% vBeta)) )

# stable parametrisation of the log-likelihood function
logLikelihoodLogitStable = function(vBeta, mX, vY) 
  return(-sum(
    vY*(mX %*% vBeta - log(1+exp(mX %*% vBeta)))
    + (1-vY)*(-log(1 + exp(mX %*% vBeta)))
  )  # sum
  )  # return 


# score function
likelihoodScore = function(vBeta, mX, vY) 
  return(t(mX) %*% (logit(mX, vBeta) - vY) )
    

# initial set of parameters (arbitrary starting parameters)
vBeta0 = c(10, -0.1, -0.3, 0.001, 0.01, 0.01, 0.001, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01)

optimLogitLBFGS = optimx(vBeta0, logLikelihoodLogitStable, method = 'L-BFGS-B',gr = likelihoodScore, mX = mX, vY = vY, hessian=TRUE)

这是错误：

Error in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower,  : Cannot evaluate function at initial parameters

【问题讨论】：

问题是什么？代码在哪里？有什么错误（如果有的话）？

为什么要使用 L-BFGS？您是否考虑过使用基本 glm 包？ stat.ethz.ch/R-manual/R-patched/library/stats/html/glm.html

我已经使用 glm 进行了逻辑回归，但现在我正在对 Spark（使用 L-BFGS）和 R 进行性能比较。

我的问题是：如何用 L-BFGS 和 R 进行逻辑回归？

你确定吗？我收到警告而不是错误：Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings. This can cause poor performance in optimization. It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.

【参考方案1】：

根据你的数据，我得到：

optimLogitLBFGS
#                p1         p2       p3         p4         p5         p6
# L-BFGS-B 9.720242 -0.1652943 0.525449 0.01681583 0.02781123 -0.3921004
#                 p7          p8         p9       p10        p11        p12
# L-BFGS-B -1.694412 -0.03461208 0.02759248 0.1993573 -0.6718275 0.02537887
#                 p13      p14   value fevals gevals niter convcode  kkt1  kkt2
# L-BFGS-B -0.8374338 0.625044 187.581    121    121    NA        1 FALSE FALSE
#          xtimes
# L-BFGS-B  0.044

使用的代码，稍作修改，不使用Excel直接创建vY：

library(readxl)  # used to read xls file
library(optimx)

premature <- read_excel("prematures.xls")

vY = as.matrix(premature['PREMATURE'])
# Recoding the response variable
vY = ifelse(vY == "positif", 1, 0)

mX = as.matrix(premature[c('GEST', 'DILATE', 'EFFACE', 'CONSIS', 'CONTR', 
                           'MEMBRAN', 'AGE', 'STRAT', 'GRAVID', 'PARIT', 
                           'DIAB', 'TRANSF', 'GEMEL')])

#add an intercept to the predictor variables
mX = cbind(rep(1, nrow(mX)), mX)

#the number of variables and observations
iK = ncol(mX)
iN = nrow(mX)

#define the logistic transformation
logit = function(mX, vBeta) 
  return(exp(mX %*% vBeta)/(1+ exp(mX %*% vBeta)) )


# stable parametrisation of the log-likelihood function
logLikelihoodLogitStable = function(vBeta, mX, vY) 
  return(-sum(
    vY*(mX %*% vBeta - log(1+exp(mX %*% vBeta)))
    + (1-vY)*(-log(1 + exp(mX %*% vBeta)))
  )  # sum
  )  # return 


# score function
likelihoodScore = function(vBeta, mX, vY) 
  return(t(mX) %*% (logit(mX, vBeta) - vY) )


# initial set of parameters (arbitrary starting parameters)
vBeta0 = c(10, -0.1, -0.3, 0.001, 0.01, 0.01, 0.001, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01)

optimLogitLBFGS = optimx(vBeta0, logLikelihoodLogitStable,
                         method = 'L-BFGS-B', gr = likelihoodScore, 
                         mX = mX, vY = vY, hessian=TRUE)

# Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower,  :
#  Parameters or bounds appear to have different scalings.
#  This can cause poor performance in optimization. 
#  It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.

【讨论】：

谢谢，我已经纠正了我的问题，这是因为 vY 被定义为我的数据集中的一个因素。您知道是否可以使用带有分类变量的 L-BFGS 执行逻辑回归？