Как я могу обучить модель gl mnet (семейство Пуассона) со смещенным членом, используя пакет кареток в R? - PullRequest
3 голосов
/ 08 апреля 2020

Я хочу смоделировать количество страховых претензий, используя Пуассона gl mnet. Данные, которые у меня есть, содержат количество заявок для каждой политики (которая является переменной ответа), некоторые характеристики политики (пол, регион и т. Д. c.), А также продолжительность политики (в годах). , Я хочу включить логарифмическую продолжительность как смещенный термин, как мы обычно делаем в актуарной науке. С функцией cv.glmnet пакета 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

Однако моя цель - обучить эту модель с использованием функции train пакета caret, потому что много преимуществ это дает. Действительно, проверка, предварительная обработка, а также выбор функций намного лучше с этим пакетом. Обучить базовый c gl mnet (без смещения) просто с помощью caret:

library(caret)

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

fit

Наивно, мы можем попытаться добавить аргумент offset в train function:

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 не заботится о предоставлении значения для аргумента newoffset в функции predict.glmnet.

В этой книге они показывают, как добавить смещенный член в модель GLM путем изменения исходного кода функции caret::train. Работает отлично. Однако функция predict.glm сильно отличается от функции predict.glmnet, поскольку не имеет аргумента newoffset. Я пытался изменить исходный код функции caret::train, но у меня возникли некоторые проблемы, потому что я недостаточно хорошо знаю, как эта функция работает.

Ответы [ 2 ]

2 голосов
/ 09 апреля 2020

Простой способ сделать это - передать столбец offset как часть x и в каждом вызове fit и predict передать как x столбцы x, которые не являются offset , В то время как offset / newoffset передайте столбец x, соответствующий offset.

. В следующем примере столбец offest для x также должен называться "offset". Это можно изменить относительно просто

Для создания функции мы просто будем использовать множество деталей из: https://github.com/topepo/caret/blob/master/models/files/glmnet.R

gl mnet является своеобразным, так как для этого требуется loop, остальное просто промыть и повторить из 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

Это не будет работать с предварительной обработкой в ​​каретке, так как мы передаем смещение как одну из функций. Однако он будет работать с рецептами , поскольку вы можете определить столбцы, для которых будут выполняться функции предварительной обработки, с помощью выборы . Более подробная информация: https://tidymodels.github.io/recipes/articles/Selecting_Variables.html

У меня не было времени, чтобы проверить мой код ошибки. Если возникнут какие-либо проблемы или если где-то будет ошибка, пожалуйста, прокомментируйте. Спасибо.

Вы также можете опубликовать проблему в caret github с просьбой добавить эту функцию (смещение / newoffset) в модель

1 голос
/ 08 апреля 2020

Я пытался изменить информацию о модели разными способами, но она с треском провалилась. Ниже я могу предложить одно решение, которое может быть не самым лучшим, но куда-то вас доставит, если ваши данные будут разумными. Вы можете прочитать больше здесь и здесь :

enter image description here

где tx - смещение. В gl mnet есть штрафной коэффициент, который вы можете ввести для каждого термина, и если вы позволите ему быть равным 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.gl mnet. Также следует отметить, что cv.gl mnet часто использует 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)

enter image description here

...