Por que o LASSO não encontra meu par perfeito de preditores em alta dimensionalidade?

20

Estou executando um pequeno experimento com regressão LASSO em R para testar se é capaz de encontrar um par preditor perfeito. O par é definido assim: f1 + f2 = resultado

O resultado aqui é um vetor predeterminado chamado 'idade'. F1 e f2 são criados pegando metade do vetor de idade e configurando o restante dos valores para 0, por exemplo: idade = [1,2,3,4,5,6], f1 = [1,2,3, 0,0,0] e f2 = [0,0,0,4,5,6]. Combino esse par preditivo com uma quantidade crescente de variáveis ​​criadas aleatoriamente por amostragem de uma distribuição normal N (1,1).

O que vejo é que, quando atingi 2 ^ 16 variáveis, o LASSO não encontra mais meu par. Veja os resultados abaixo.

Número de recursos por dobra por tamanho de dadosCoeficientes do par perfeito

Por que isso está acontecendo? Você pode reproduzir os resultados com o script abaixo. Percebi que quando escolho um vetor de idade diferente, por exemplo: [1: 193], o LASSO encontra o par em alta dimensionalidade (> 2 ^ 16).

O 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
fonte
comentário secundário: devido ao uso de 'createFolds', você também precisa do pacote 'caret' carregado.
IWS
2
Veja o teorema 2a de 'Wainwright: limiares agudos para recuperação de alta dispersão dimensional e barulhenta da esparsidade'. No regime em que você está, onde o suporte verdadeiro tem cardinalidade fixa 2 ep cresce com n fixo, parece provável que haja correlações muito altas se houver recursos suficientes, o que leva à baixa probabilidade de recuperação bem-sucedida do suporte que você percebe. (No entanto, desde os vetores não no verdadeiro suporte são muito pequenos (média 0 variância 1), parece que isso não pode ser a razão desde o recurso de idade verdade tem muito grandes entradas.)
user795305
1
@ Ben, acho que esta é a explicação correta e, dada a popularidade dessa pergunta, seria ótimo se você pudesse fornecer uma resposta que explique por que isso é verdade.
NRH 07/02
1
@ Maddenker ^sempre retorna um duplo para argumentos inteiros ou duplos em R. R também muda para duplos se ocorrer um estouro de número inteiro.
Roland
1
FYI: Adicionei um script atualizado na minha página do github. Neste script, uso menos amostras, o que já induz o problema em 2 ^ 5 variáveis. Isto permite tempos de corrida rápida e permite que você experimentar mais com os dados: github.com/sjorsvanheuveln/LASSO_pair_problem
Ansjovis86

Respostas:

4

Esse problema é bem conhecido por acadêmicos e pesquisadores. A resposta, no entanto, não é simples e pertence mais - na minha opinião - à otimização do que às estatísticas. As pessoas tentaram superar essas desvantagens incluindo uma penalidade adicional na cadeia, daí a regressão líquida elástica.p>n

O laço é uma ferramenta popular para regressão linear esparsa, especialmente para problemas em que o número de variáveis ​​excede o número de observação. Mas quando p> n, o critério do laço não é estritamente convexo e, portanto, pode não ter um minimizador exclusivo.

Como o @ben mencionou, quando você tem 2e16 covariáveis, não é diferente de algumas serem bastante semelhantes às verdadeiras covariáveis. Por isso, o ponto acima é relevante: o LASSO é indiferente a escolher um dos dois.

Talvez mais relevante e mais recentemente (2013), exista outro artigo de Candes sobre como, mesmo quando as condições estatísticas são ideais (preditores não correlacionados, apenas alguns grandes efeitos), o LASSO ainda produz falsos positivos, como o que você vê em seus dados:

Em cenários de regressão em que variáveis ​​explicativas têm correlações muito baixas e há relativamente poucos efeitos, cada um de grande magnitude, esperamos que o Lasso encontre as variáveis ​​importantes com poucos erros, se houver. Este artigo mostra que, em um regime de esparsidade linear - o que significa que a fração de variáveis ​​com efeito de não fuga - tende a uma constante, por menor que seja -, esse não pode realmente ser o caso, mesmo quando as variáveis ​​de design são estocásticas independentemente .

Mustafa S Eisa
fonte
Eu não sabia disso. Eu achava que o LASSO era uma ferramenta padrão e confiável para identificar modelos esparsos (ou pelo menos essa foi minha impressão lendo os dois livros de Hastie e Tibshirani e usando o método). Como você diz que o problema é conhecido, você sabe se há também soluções / e / ou abordagens alternativas?
DeltaIV 11/02/17
Se bem entendi, estes resultados parecem ser apenas sparsity linear, enquanto o problema na mão preocupações sub linear sparsity
user795305
@ Ben, claro, mas isso não torna o papel irrelevante. É uma das peças mais recentes da literatura que conheço que aborda esta questão. Acho que mostra algo simples: mesmo com condições estatísticas ideais, o LASSO não tem as melhores propriedades.
Mustafa S Eisa
@DeltaIV, LASSO é uma heurística de otimização convexa para fins de seleção de variáveis. No livro de Tibshirani, eles mostram que ele pode seguir um caminho semelhante ao da AIC ou dos métodos passo a passo, mas isso não é uma garantia. Na minha opinião, a maioria de seus problemas vem do fato de ser uma heurística e não a coisa real, mas você desiste disso para ganhar convexidade, que tem outras propriedades interessantes.
Mustafa S Eisa