如何在glm中指定一种特殊的优化方法

Question

在glm的文档中，有一个method选项，在文档中它说

用户提供的拟合函数既可以作为函数提供，也可以作为命名函数的字符串提供，函数采用与glm.fit相同的参数。如果指定为字符串，则从stats命名空间中查找。

我想使用二项式非负最小二乘优化器，使系数为非负，并将其加到1。这个优化器的一个例子是在带有选项method = "method.NNLS"的SuperLearner包中。下面是一个可重复的例子：

library(SuperLearner)
# binary outcome
set.seed(1)
N <- 200
X <- matrix(rnorm(N*10), N, 10)
X <- as.data.frame(X)
Y <- rbinom(N, 1, plogis(.2*X[, 1] + .1*X[, 2] - .2*X[, 3] +
                           .1*X[, 3]*X[, 4] - .2*abs(X[, 4])))

SL.library <- c("SL.glmnet", "SL.glm", "SL.knn", "SL.mean")

# least squares loss function
test.NNLS <- SuperLearner(Y = Y, X = X, SL.library = SL.library,
                          verbose = TRUE, method = "method.NNLS", family = binomial())
> test.NNLS

Call:  
SuperLearner(Y = Y, X = X, family = binomial(), SL.library = SL.library, method = "method.NNLS",  
    verbose = TRUE) 


                   Risk      Coef
SL.glmnet_All 0.2460486 0.0000000
SL.glm_All    0.2507033 0.2423697
SL.knn_All    0.2508500 0.3493301
SL.mean_All   0.2475494 0.4083002

请注意，Coefs是非负的，和为1。有办法在glm中使用相同的优化器吗？还是使用相同约束的不同优化器？我试过：

glm(Y ~ as.matrix(X), family = "binomial", method = "method.NNLS")

但这不起作用。

Answer 1

解决这个问题的一种方法(不是唯一的方法，可能不是最有效或最紧凑的方法)是对变换后的空间中的参数进行一般的极大似然估计；这种变换可能最常见的是加性对数比变换(例如在compositions::alr()中)。这个解决方案使用bbmle包，它是围绕optim()的包装；在本例中，它实际上没有提供比直接使用optim()更多的优势

library(bbmle)
library(compositions)
linkfun <- function(mu) as.numeric(compositions::alr(mu))
invlink <- function(eta) as.numeric(compositions::alrInv(eta))

## negative log-likelihood function: first convert the parameter vector
## from the constrained (alr) space to the unconstrained space, then
## convert the linear predictor from an unconstrained space to the
## constrained (probability, (0,1)) space
lfun <- function(beta) {
    -sum(dbinom(Y, size=1, 
                prob = plogis(as.matrix(X) %*% invlink(beta)),
                log=TRUE))
}
npar <- ncol(X)

使用mle2

## extra setup for using mle2 with a score function that takes a vector
## parameter
pnm <- paste0("beta", seq(npar))
parnames(lfun) <- pnm
    
m1 <- mle2(lfun, start = setNames(rep(0, npar-1), pnm[-1]), vecpar = TRUE)
invlink(coef(m1))
##  [1] 4.897617e-01 2.221793e-01 1.274830e-05 2.198596e-05 7.541545e-02
##  [6] 1.175966e-01 1.590407e-05 4.662563e-02 4.825376e-02 1.170077e-04

或

opt <- optim(par = rep(0, npar-1), fn = lfun, method = "BFGS")
invlink(opt$par)

附注1

我的模型fit与您用来生成数据的方法不匹配(这将以~ 0 + x1 + x2 + x3 + x3:x4 + abs(x4)的形式编码)，但是您的SuperLearner示例也不匹配，所以我不担心它。

附注2

我的猜测(没有仔细考虑)是，上面的方法将not给出与使用SuperLearner相同的结果；查看代码，SuperLearner所做的是找到(加权的)最小二乘问题的非负解，然后通过将系数向量除以其和将系数转换到组合空间。这与求解参数为非负且和为1的最小偏差/最大似然问题不同，如果您想要SuperLearner解，可以通过向glm传递适当的method=参数来实现，但是我不知道如何用glm来实现约束的最大似然解。