Использование glmnet для поиска оптимизированной модели для заданного числа предикторов - PullRequest
Новая
       6

Использование glmnet для поиска оптимизированной модели для заданного числа предикторов

0 голосов
/ 09 января 2019

Я пытаюсь использовать LASSO для функции, немного отличающейся от первоначально разработанной. У меня есть 22 различных задания на тесте, которые при усреднении дают окончательный результат. Я хочу посмотреть, какая комбинация ограниченного числа задач лучше всего предскажет общий балл в надежде создать краткую форму этого теста.

Я использую glmnet для запуска лассо, и он работает как положено. Затем я могу легко найти модель по заданной лямде с 

coef(cvfit, s = s)

Однако мне интересно, можно ли было бы указать n предикторов, имеющих ненулевые коэффициенты, а не параметр штрафа?

Я создал очень неэффективный способ сделать это, как показано ниже, извлекая модели из сетки тестовых лямбд, но мне было интересно, есть ли более эффективный способ сделать это 

nvar <- list()
coeffs <- list()

for(j in 1:20000) {

  s <- j / 20000

  coeffs[j] <- coef(cvfit, s = s) ##Get coefficient list at given lamda

  nvar[j] <- sum(as.vector(coef(cvfit, s = s)) != 0) - 1 ##Count number of variables with non-zero coeff and subtract one because intercept is always non-zero

}

nvar <- unlist(nvar)

getlamda <- function(numvar = 4) {

  min.lambda <- min(lambdas[nvar == numvar]) / 20000 ##Find the smallest lambda which resulted in the given number of non-zero coefficients

  coeffs[min.lambda]

}

Ответы [ 2 ]

0 голосов
/ 10 января 2019

После игры с решением Blended выше я понял, что есть еще более простой способ сделать это.

Использование набора данных Бостона, использованного в примере: 

library(dplyr)
library(glmnet)

(boston <- MASS::Boston %>% tbl_df())

tr_x <- model.matrix(medv ~ ., data = boston)[,-1]
tr_y <- boston$medv
cvfit <- glmnet(tr_x, tr_y)

Объект cvfit уже имеет все компоненты, которые нам нужны, чтобы найти ответ для заданного числа переменных. df - это число степеней свободы, и это число переменных параметров, в которых мы заинтересованы. lambda - лямбда для каждой модели.

Таким образом, мы можем создать простую функцию, которая возвращает наилучшую модель для заданного числа переменных. 

get_nparam <- function(mod, numvar) {

  coef(mod, s = with(cvfit, min(lambda[df == numvar])))

}

get_nparam(cvfit, 4)

#14 x 1 sparse Matrix of class "dgCMatrix"
#                       1
#(Intercept) 15.468034114
#crim         .          
#zn           .          
#indus        .          
#chas         .          
#nox          .          
#rm           3.816165372
#age          .          
#dis          .          
#rad          .          
#tax          .          
#ptratio     -0.606026131
#black        0.001518042
#lstat       -0.495954410
#

Еще раз спасибо Blender за другое решение и за то, что он помог мне в этом.

0 голосов
/ 10 января 2019

Вы можете использовать rowSums(). 

(boston <- MASS::Boston %>% tbl_df())
#> # A tibble: 506 x 14
#>       crim    zn indus  chas   nox    rm   age   dis   rad   tax ptratio
#>  *   <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <dbl>   <dbl>
#>  1 0.00632  18    2.31     0 0.538  6.58  65.2  4.09     1   296    15.3
#>  2 0.0273    0    7.07     0 0.469  6.42  78.9  4.97     2   242    17.8
#>  3 0.0273    0    7.07     0 0.469  7.18  61.1  4.97     2   242    17.8
#>  4 0.0324    0    2.18     0 0.458  7.00  45.8  6.06     3   222    18.7
#>  5 0.0690    0    2.18     0 0.458  7.15  54.2  6.06     3   222    18.7
#>  6 0.0298    0    2.18     0 0.458  6.43  58.7  6.06     3   222    18.7
#>  7 0.0883   12.5  7.87     0 0.524  6.01  66.6  5.56     5   311    15.2
#>  8 0.145    12.5  7.87     0 0.524  6.17  96.1  5.95     5   311    15.2
#>  9 0.211    12.5  7.87     0 0.524  5.63 100    6.08     5   311    15.2
#> 10 0.170    12.5  7.87     0 0.524  6.00  85.9  6.59     5   311    15.2
#> # ... with 496 more rows, and 3 more variables: black <dbl>, lstat <dbl>,
#> #   medv <dbl>

Для вышеуказанного набора данных ( Бостонский корпус ) рассмотрим medv ~ .. 

library(glmnet)
tr_x <- model.matrix(medv ~ ., data = boston)[,-1]
tr_y <- boston$medv
cvfit <- glmnet(tr_x, tr_y)

head(t(coef(cvfit)))
    #> 6 x 14 sparse Matrix of class "dgCMatrix"
    #>    [[ suppressing 14 column names '(Intercept)', 'crim', 'zn' ... ]]
    #>                                                        
    #> s0 22.53281 . . . . . .         . . . . . .  .         
    #> s1 23.60072 . . . . . .         . . . . . . -0.08439977
    #> s2 23.67264 . . . . . 0.1278413 . . . . . . -0.15358093
    #> s3 21.44649 . . . . . 0.5694424 . . . . . . -0.19698136
    #> s4 19.42057 . . . . . 0.9714620 . . . . . . -0.23654740
    #> s5 17.57464 . . . . . 1.3377669 . . . . . . -0.27259852

Полагаю, вы проделали эту процедуру.

Примечания

  1. Может быть удобно транспонировать матрицу коэффициентов, чтобы каждая переменная становилась каждым столбцом.
  2. Для t(coef(cvfit)), rowSums(t(coef(cvfit)) != 0) находит номер ненулевого элемента для каждой переменной.
  3. Далее мы просто сопоставляем numvar с этим rowSums и находим значение коэффициента.

Обозначим, что от s0 до s5 лямбда s0 больше s5 - больше штрафуется. 

head(cvfit$lambda)
#> [1] 6.777654 6.175546 5.626927 5.127046 4.671574 4.256564

Подмножество коф с numvar

На основании этих фактов 

get_nparam <- function(mod, numvar) {
  beta <- coef(mod)
  non_zero <- rowSums(t(beta)[,-1] != 0) # ignore intercept
  min_lam <- which(non_zero == numvar) # numvar non-zero coef
  t(beta)[dplyr::last(min_lam),] # last index = smallest lambda
}

С помощью этой функции вы можете получить 

get_nparam(cvfit, 4)
#>  (Intercept)         crim           zn        indus         chas 
#> 15.468034114  0.000000000  0.000000000  0.000000000  0.000000000 
#>          nox           rm          age          dis          rad 
#>  0.000000000  3.816165372  0.000000000  0.000000000  0.000000000 
#>          tax      ptratio        black        lstat 
#>  0.000000000 -0.606026131  0.001518042 -0.495954410

rm, ptratio, black и lstat отличны от нуля, а остальные равны нулю.

...