Pourquoi LASSO ne trouve-t-il pas ma paire de prédicteurs parfaite à haute dimensionnalité?

20

J'exécute une petite expérience avec la régression LASSO dans R pour tester s'il est capable de trouver une paire de prédicteurs parfaite. La paire est définie comme ceci: f1 + f2 = résultat

Le résultat ici est un vecteur prédéterminé appelé «âge». F1 et f2 sont créés en prenant la moitié du vecteur d'âge et en définissant le reste des valeurs à 0, par exemple: âge = [1,2,3,4,5,6], f1 = [1,2,3, 0,0,0] et f2 = [0,0,0,4,5,6]. Je combine cette paire de prédicteurs avec une quantité croissante de variables créées aléatoirement en échantillonnant à partir d'une distribution normale N (1,1).

Ce que je vois, c'est quand je frappe 2 ^ 16 variables, LASSO ne trouve plus ma paire. Voir les résultats ci-dessous.

Nombre de fonctionnalités par pli par taille de donnéesCoefficients de la paire parfaite

Pourquoi cela arrive-t-il? Vous pouvez reproduire les résultats avec le script ci-dessous. J'ai remarqué que lorsque je choisis un vecteur d'âge différent, par exemple: [1: 193], alors LASSO trouve la paire à haute dimensionnalité (> 2 ^ 16).

Le script:

## Setup ##
library(glmnet)
library(doParallel)
library(caret)

mae <- function(errors){MAE <- mean(abs(errors));return(MAE)}
seed = 1
n_start <- 2 #start at 2^n features
n_end <- 16 #finish with 2^n features
cl <- makeCluster(3)
registerDoParallel(cores=cl)

#storage of data
features <- list()
coefs <- list()
L <- list() 
P <- list() 
C <- list() 
RSS <- list() 

## MAIN ##
for (j in n_start:n_end){
  set.seed(seed)
  age <- c(55,31,49,47,68,69,53,42,58,67,60,58,32,52,63,31,51,53,37,48,31,58,36,42,61,49,51,45,61,57,52,60,62,41,28,45,39,47,70,33,37,38,32,24,66,54,59,63,53,42,25,56,70,67,44,33,50,55,60,50,29,51,49,69,70,36,53,56,32,43,39,43,20,62,46,65,62,65,43,40,64,61,54,68,55,37,59,54,54,26,68,51,45,34,52,57,51,66,22,64,47,45,31,47,38,31,37,58,66,66,54,56,27,40,59,63,64,27,57,32,63,32,67,38,45,53,38,50,46,59,29,41,33,40,33,69,42,55,36,44,33,61,43,46,67,47,69,65,56,34,68,20,64,41,20,65,52,60,39,50,67,49,65,52,56,48,57,38,48,48,62,48,70,55,66,58,42,62,60,69,37,50,44,61,28,64,36,68,57,59,63,46,36)
  beta2 <- as.data.frame(cbind(age,replicate(2^(j),rnorm(length(age),1,1))));colnames(beta2)[1] <-'age'

  f1 <- c(age[1:96],rep(0,97)) 
  f2 <- c(rep(0,96),age[97:193])
  beta2 <- as.data.frame(cbind(beta2,f1,f2))

  #storage variables
  L[[j]] <- vector()
  P[[j]] <- vector()
  C[[j]] <- list()
  RSS[[j]] <- vector()

  #### DCV LASSO ####
  set.seed(seed) #make folds same over 10 iterations
  for (i in 1:10){

    print(paste(j,i))
    index <- createFolds(age,k=10)
    t.train <- beta2[-index[[i]],];row.names(t.train) <- NULL
    t.test <- beta2[index[[i]],];row.names(t.test) <- NULL

    L[[j]][i] <- cv.glmnet(x=as.matrix(t.train[,-1]),y=as.matrix(t.train[,1]),parallel = T,alpha=1)$lambda.min #,lambda=seq(0,10,0.1)
    model <- glmnet(x=as.matrix(t.train[,-1]),y=as.matrix(t.train[,1]),lambda=L[[j]][i],alpha=1)
    C[[j]][[i]] <- coef(model)[,1][coef(model)[,1] != 0]
    pred <- predict(model,as.matrix(t.test[,-1]))
    RSS[[j]][i] <- sum((pred - t.test$age)^2)
    P[[j]][i] <- mae(t.test$age - pred)
    gc()
  }
}

##############
## PLOTTING ##
##############

#calculate plots features
beta_sum = unlist(lapply(unlist(C,recursive = F),function(x){sum(abs(x[-1]))}))
penalty = unlist(L) * beta_sum
RSS = unlist(RSS)
pair_coefs <- unlist(lapply(unlist(C,recursive = F),function(x){
  if('f1' %in% names(x)){f1 = x['f1']}else{f1=0;names(f1)='f1'}
  if('f2' %in% names(x)){f2 = x['f2']}else{f2=0;names(f2)='f2'}
  return(c(f1,f2))}));pair_coefs <- split(pair_coefs,c('f1','f2'))
inout <- lapply(unlist(C,recursive = F),function(x){c('f1','f2') %in% names(x)})
colors <- unlist(lapply(inout,function(x){if (x[1]*x[2]){'green'}else{'red'}}))
featlength <- unlist(lapply(unlist(C,recursive = F),function(x){length(x)-1}))

#diagnostics
plot(rep(n_start:n_end,each=10),pair_coefs$f1,col='red',xaxt = "n",xlab='n/o randomly generated features (log2)',main='Pair Coefficients',ylim=c(0,1),ylab='pair coefficients');axis(1, at=n_start:n_end);points(rep(n_start:n_end,each=10),pair_coefs$f2,col='blue');axis(1, at=n_start:n_end, labels=(n_start:n_end));legend('bottomleft',fill=c('red','blue'),legend = c('f1','f2'),inset=.02)
plot(rep(n_start:n_end,each=10),RSS+penalty,col=colors,xaxt = "n",xlab='n/o randomly generated features (log2)',main='RSS+penalty');axis(1, at=n_start:n_end, labels=(n_start:n_end));legend('topleft',fill=c('green','red'),legend = c('Pair Selected','Pair not Selected'),inset=.02)
plot(rep(n_start:n_end,each=10),penalty,col=colors,xaxt = "n",xlab='n/o randomly generated features (log2)',main='Penalty');axis(1, at=n_start:n_end, labels=(n_start:n_end));legend('topleft',fill=c('green','red'),legend = c('Pair Selected','Pair not Selected'),inset=.02)
plot(rep(n_start:n_end,each=10),RSS,col=colors,xaxt = "n",xlab='n/o randomly generated features (log2)',main='RSS');axis(1, at=n_start:n_end, labels=(n_start:n_end));legend('topleft',fill=c('green','red'),legend = c('Pair Selected','Pair not Selected'),inset=.02)
plot(rep(n_start:n_end,each=10),unlist(L),col=colors,xaxt = "n",xlab='n/o randomly generated features (log2)',main='Lambdas',ylab=expression(paste(lambda)));axis(1, at=n_start:n_end, labels=(n_start:n_end));legend('topleft',fill=c('green','red'),legend = c('Pair Selected','Pair not Selected'),inset=.02)
plot(rep(n_start:n_end,each=10),featlength,ylab='n/o features per fold',col=colors,xaxt = "n",xlab='n/o randomly generated features (log2)',main='Features per Fold');axis(1, at=n_start:n_end, labels=(n_start:n_end));legend('topleft',fill=c('green','red'),legend = c('Pair Selected','Pair not Selected'),inset=.02)
plot(penalty,RSS,col=colors,main='Penalty vs. RSS')
Ansjovis86
la source
commentaire mineur: en raison de l'utilisation de 'createFolds', vous devez également charger le package 'caret'.
IWS
2
Voir le théorème 2a de «Wainwright: seuils nets pour une récupération de densité élevée et bruyante». Dans le régime dans lequel vous êtes, où le véritable support a une cardinalité fixe 2 et p croît avec n fixe, il semble probable qu'il puisse y avoir des corrélations très élevées s'il y a suffisamment de fonctionnalités, ce qui conduit à la faible probabilité de récupération du support réussie que vous remarquez. (Cependant, comme les vecteurs qui ne sont pas dans le vrai support sont assez petits (moyenne 0 variance 1), il semble que ce ne soit pas la raison puisque la fonction de l'âge réel a de très grandes entrées.)
user795305
1
@Ben, je pense que c'est l'explication correcte, et étant donné la popularité de cette question, ce serait formidable, si vous pouviez fournir une réponse qui explique pourquoi il en est ainsi.
NRH
1
@Maddenker ^renvoie toujours un double pour les arguments entiers ou doubles dans R. R passe également en double en cas de dépassement d'entier.
Roland
1
FYI: J'ai ajouté un script mis à jour sur ma page github. Dans ce script, j'utilise moins d'échantillons, ce qui induit déjà le problème avec 2 ^ 5 variables. Cela permet des temps d'exécution rapides et vous permet d'expérimenter davantage avec les données: github.com/sjorsvanheuveln/LASSO_pair_problem
Ansjovis86

Réponses:

4

Ce problème est bien connu des universitaires et des chercheurs. Cependant, la réponse n'est pas simple et se rapporte davantage - à mon avis - à l'optimisation qu'à la statistique. Les gens ont tenté de surmonter ces inconvénients en incluant une pénalité de crête supplémentaire, d'où la régression nette élastique. Cet article de Tibshirani concerne le problème (c'est-à-dire le nombre de covariables plus grand que le nombre d'observations):p>n

Le lasso est un outil populaire pour la régression linéaire clairsemée, en particulier pour les problèmes dans lesquels le nombre de variables dépasse le nombre d'observations. Mais lorsque p> n, le critère du lasso n'est pas strictement convexe, et donc il peut ne pas avoir un minimiseur unique.

Comme @ben l'a mentionné, lorsque vous avez des covariables 2e16, ce n'est pas sans rappeler que certaines sont assez similaires aux vraies covariables. D'où la raison pour laquelle le point ci-dessus est pertinent: LASSO est indifférent à choisir l'un ou l'autre.

Peut-être plus pertinent et plus récemment (2013), il y a un autre article de Candes sur la façon dont, même lorsque les conditions statistiques sont idéales (prédicteurs non corrélés, seulement quelques grands effets), le LASSO produit toujours des faux positifs, comme ce que vous voyez dans vos données:

Dans les paramètres de régression où les variables explicatives ont des corrélations très faibles et où il y a relativement peu d'effets, chacun de grande ampleur, nous nous attendons à ce que le Lasso trouve les variables importantes avec peu d'erreurs, le cas échéant. Cet article montre que dans un régime de rareté linéaire --- ce qui signifie que la fraction de variables avec un effet non-disparaissant tend vers une constante, même petite --- cela ne peut pas vraiment être le cas, même lorsque les variables de conception sont stochastiquement indépendantes .

Mustafa S Eisa
la source
Je ne le savais pas. Je pensais que LASSO était un outil standard et fiable pour identifier un modèle clairsemé (ou du moins, c'était mon impression en lisant les deux livres Hastie et Tibshirani, et en utilisant la méthode moi-même). Puisque vous dites que le problème est bien connu, savez-vous s'il existe également des solutions et / ou des approches alternatives?
DeltaIV
Si je comprends bien, ces résultats semblent ne concerner que la rareté linéaire, alors que le problème à résoudre concerne la sous-linéarité
user795305
@Ben, bien sûr, mais cela ne rend pas le papier hors de propos. C'est l'un des documents les plus récents que je connaisse qui touche à cette question. Je pense que cela montre quelque chose de simple: même avec des conditions statistiques idéales, LASSO n'a pas les meilleures propriétés.
Mustafa S Eisa
@DeltaIV, LASSO est une heuristique d'optimisation convexe à des fins de sélection de variables. Dans le livre de Tibshirani, ils montrent qu'il peut suivre un chemin similaire à l'AIC ou aux méthodes par étapes, mais ce n'est pas une garantie. À mon avis, la plupart de ses problèmes viennent du fait que c'est une heuristique et non la vraie chose, mais vous abandonnez pour gagner en convexité qui a d'autres belles propriétés.
Mustafa S Eisa