如何使用R中的插入符号包训练一个带偏移项的glmnet模型(Poisson族)？

7eumitmz 于 2022-12-20 发布在其他

关注(0)|答案(2)|浏览(150)

我想使用Poisson glmnet对保险索赔计数进行建模。我手头的数据包含每份保单的索赔数量（这是响应变量），关于策略的一些特性（性别、地区等）以及政策的持续时间（以年为单位）。我想把对数持续时间作为一个补偿项，就像我们在精算学中通常做的那样。使用glmnet包中的cv.glmnet函数，它是直截了当：

library(tidyverse)
library(glmnet)

n <- 100

dat <- tibble(
 nb_claims = rpois(n, lambda = 0.5),
 duration = runif(n),
 x1 = runif(n),
 x2 = runif(n),
 x3 = runif(n)
)

fit <- cv.glmnet(
  x = dat %>% dplyr::select(x1, x2, x3) %>% as.matrix(),
  y = dat %>% pull(nb_claims),
  family = "poisson",
  offset = dat %>% pull(duration) %>% log()
)

fit

然而，我的目标是使用caret包中的train函数来训练这个模型，因为它提供了许多优点，实际上，使用这个包，验证、预处理以及特征选择都要好得多，使用caret训练一个基本的glmnet（没有偏移项）非常简单：

library(caret)

fit <- caret::train(
  x = dat %>% dplyr::select(x1, x2, x3) %>% as.matrix(),
  y = dat %>% pull(nb_claims),
  method = "glmnet",
  family = "poisson"
)

fit

我们可以天真地尝试在train函数中添加offset参数：

fit <- caret::train(
  x = dat %>% dplyr::select(x1, x2, x3) %>% as.matrix(),
  y = dat %>% pull(nb_claims),
  method = "glmnet",
  family = "poisson",
  offset = dat %>% pull(duration) %>% log()
)

fit

不幸的是，这段代码抛出了错误Error : No newoffset provided for prediction, yet offset used in fit of glmnet，这是因为caret::train函数没有注意为predict.glmnet函数中的newoffset参数给予。
在这本书中，他们展示了如何通过修改caret::train函数的源代码来向GLM模型添加偏移项。它工作得很完美。然而，predict.glm函数与predict.glmnet函数有很大不同，因为它没有newoffset参数。我尝试修改caret::train函数的源代码。但是我遇到了一些麻烦，因为我不太清楚这个函数是如何工作的。

来源：https://stackoverflow.com/questions/61104205/how-can-i-train-a-glmnet-model-poisson-family-with-an-offset-term-using-the-ca

2条答案

按热度按时间

0dxa2lsx1#

执行此操作的简单方法是将offset列作为x的一部分传递，并在每个fit和predict调用中作为x的x列（不是offset）传递，而作为offset/newoffset传递与offset对应的x列。
在下面的示例中，x的offest列也需要命名为“offset”，这可以相对容易地更改
要创建函数，我们只需使用以下代码中的许多部分：https://github.com/topepo/caret/blob/master/models/files/glmnet.R
glmnet是特殊的，因为它需要一个loop，其余的只是冲洗和reapeat从https://topepo.github.io/caret/using-your-own-model-in-train.html#illustrative-example-1-svms-with-laplacian-kernels
将在全文中指定family = "poisson"，以将此采用代码从https://github.com/topepo/caret/blob/master/models/files/glmnet.R更改为

glmnet_offset <- list(type = "Regression",
                      library = c("glmnet", "Matrix"),
                      loop = function(grid) {
                        alph <- unique(grid$alpha)
                        loop <- data.frame(alpha = alph)
                        loop$lambda <- NA
                        submodels <- vector(mode = "list", length = length(alph))
                        for(i in seq(along = alph)) {
                          np <- grid[grid$alpha == alph[i],"lambda"]
                          loop$lambda[loop$alpha == alph[i]] <- np[which.max(np)]
                          submodels[[i]] <- data.frame(lambda = np[-which.max(np)])
                        }
                        list(loop = loop, submodels = submodels)
                      }) 

glmnet_offset$parameters <- data.frame(parameter = c('alpha', 'lambda'),
                                       class = c("numeric", "numeric"),
                                       label = c('Mixing Percentage', 'Regularization Parameter'))

glmnet_offset$grid <- function(x, y, len = NULL, search = "grid") {
  if(search == "grid") {
    init <- glmnet::glmnet(Matrix::as.matrix(x[,colnames(x) != "offset"]), y,
                           family = "poisson",
                           nlambda = len+2,
                           alpha = .5,
                           offset = x[,colnames(x) == "offset"])
    lambda <- unique(init$lambda)
    lambda <- lambda[-c(1, length(lambda))]
    lambda <- lambda[1:min(length(lambda), len)]
    out <- expand.grid(alpha = seq(0.1, 1, length = len),
                       lambda = lambda)
  } else {
    out <- data.frame(alpha = runif(len, min = 0, 1),
                      lambda = 2^runif(len, min = -10, 3))
  }
  out
}

所以x[,colnames(x) != "offset"]是x，而offset是x[,colnames(x) == "offset"]

glmnet_offset$fit <- function(x, y, wts, param, last, ...) {

  theDots <- list(...)

  ## pass in any model weights
  if(!is.null(wts)) theDots$weights <- wts

  if(!(class(x)[1] %in% c("matrix", "sparseMatrix")))
    x <- Matrix::as.matrix(x)
  modelArgs <- c(list(x = x[,colnames(x) != "offset"],
                      y = y,
                      alpha = param$alpha,
                      family = "poisson",
                      offset = x[,colnames(x) == "offset"]),
                 theDots)

  out <- do.call(glmnet::glmnet, modelArgs)
  if(!is.na(param$lambda[1])) out$lambdaOpt <- param$lambda[1]
  out
}

glmnet_offset$predict <- function(modelFit, newdata, submodels = NULL) {
  if(!is.matrix(newdata)) newdata <- Matrix::as.matrix(newdata)
    out <- predict(modelFit,
                   newdata[,colnames(newdata) != "offset"],
                   s = modelFit$lambdaOpt,
                   newoffset = newdata[,colnames(newdata) == "offset"],
                   type = "response") #important for measures to be appropriate

  if(is.matrix(out)) out <- out[,1]
  out
  if(!is.null(submodels)) {
      tmp <- as.list(as.data.frame(predict(modelFit,
                                          newdata[,colnames(newdata) != "offset"],
                                          s = submodels$lambda,
                                          newoffset = newdata[,colnames(newdata) == "offset"],
                                          type = "response"),
                                   stringsAsFactors = TRUE))
    out <- c(list(out), tmp)

  }
  out

}

出于某种原因，我不明白，但它不工作，没有prob插槽

glmnet_offset$prob <- glmnet_offset$predict

glmnet_offset$tags = c("Generalized Linear Model", "Implicit Feature Selection",
                       "L1 Regularization", "L2 Regularization", "Linear Classifier",
                       "Linear Regression")

glmnet_offset$sort = function(x) x[order(-x$lambda, x$alpha),]
glmnet_offset$trim = function(x) {
  x$call <- NULL
  x$df <- NULL
  x$dev.ratio <- NULL
  x
}

library(tidyverse)
library(caret)
library(glmnet)

n <- 100
set.seed(123)
dat <- tibble(
  nb_claims = rpois(n, lambda = 0.5),
  duration = runif(n),
  x1 = runif(n),
  x2 = runif(n),
  x3 = runif(n)
)

x = dat %>%
  dplyr::select(-nb_claims) %>%
  mutate(offset = log(duration)) %>%
  dplyr::select(-duration) %>%
  as.matrix

fit <- caret::train(
  x = x,
  y = dat %>% pull(nb_claims),
  method = glmnet_offset,
)

fit
100 samples
  4 predictor

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 100, 100, 100, 100, 100, 100, ... 
Resampling results across tuning parameters:

  alpha  lambda        RMSE       Rsquared    MAE      
  0.10   0.0001640335  0.7152018  0.01805762  0.5814200
  0.10   0.0016403346  0.7152013  0.01805684  0.5814193
  0.10   0.0164033456  0.7130390  0.01798125  0.5803747
  0.55   0.0001640335  0.7151988  0.01804917  0.5814020
  0.55   0.0016403346  0.7150312  0.01802689  0.5812936
  0.55   0.0164033456  0.7095996  0.01764947  0.5783706
  1.00   0.0001640335  0.7152033  0.01804795  0.5813997
  1.00   0.0016403346  0.7146528  0.01798979  0.5810811
  1.00   0.0164033456  0.7063482  0.01732168  0.5763653

RMSE was used to select the optimal model using the smallest value.
The final values used for the model were alpha = 1 and lambda = 0.01640335.

predict(fit$finalModel,  x[,1:3], newoffset = x[,4]) #works

这不适用于插入符号中的预处理，因为我们将偏移作为特性之一传递。但是它适用于recipes，因为您可以通过selections定义要执行预处理功能的列。https://tidymodels.github.io/recipes/articles/Selecting_Variables.html
我还没有时间检查我的代码错误。如果有任何问题发生或如果有一个错误的地方，请评论。谢谢。
您还可以在插入符号github中发布一个问题，要求将此特性（offset/newoffset）添加到模型中

赞(0）回复(0）举报 2022-12-20

nhjlsmyf2#

我尝试了很多方法来改变模型信息，但都失败得很惨。下面我可以提出一个解决方案，可能不是最好的，但如果你的数据是明智的，它会让你到达某个地方。
在Poisson / negative binom..回归中，因子中的偏移被引入回归，您可以阅读更多here和here：

其中tx是偏移量。在glmnet中，可以为每个项引入一个惩罚因子，如果某项的惩罚因子为0，基本上不会惩罚它，它总是包含在内。我们可以将其用于偏移量，只有在使用有意义的数据集时才能看到这种效果（注意，在示例数据集中，偏移量是没有意义的数字）。
下面我使用来自MASS的保险索赔数据集：

library(tidyverse)
library(glmnet)
library(MASS)

dat <- Insurance
X =  model.matrix(Claims ~ District + Group + Age,data=dat)
Y = dat$Claims
OFF = log(dat$Holders)

fit_cv <- cv.glmnet(
  x = X,
  y = Y,
  family = "poisson",
  offset = OFF
)

现在使用插入符号，我将在没有任何训练的情况下拟合它，并使用从cv. glmnet中的拟合获得的相同lambda。您还应该注意的一点是，cv.glmnet经常使用lambda.1se而不是lambda.min：

fit_c <- caret::train(
  x = cbind(X,OFF),
  y = Y,
  method = "glmnet",
  family = "poisson",
  tuneGrid=data.frame(lambda=fit_cv$lambda.1se,alpha=1),
  penalty=c(rep(1,ncol(X)),0),
  trControl = trainControl(method="none")
)

我们可以看到预测的差异有多大：

p1 = predict(fit_cv,newx=X,newoffset=OFF)
p2 = predict(fit_c,newx=cbind(X,OFF))

plot(p1,p2)

赞(0）回复(0）举报 2022-12-20

我来回答

如何使用R中的插入符号包训练一个带偏移项的glmnet模型(Poisson族)？

2条答案

相关问题

热门标签

最新问答