2016-07-14 13 views
8

Mi aspetto LASSO senza penalizzazione ($ \ lambda = 0 $) per ottenere le stesse stime del coefficiente (o molto simili) di un adattamento OLS. Tuttavia, ottengo diverse stime dei coefficienti a R mettendo gli stessi dati (x, y) in

  • glmnet(x, y , alpha=1, lambda=0) per LASSO misura senza penalizzazione e
  • lm(y ~ x) per OLS misura.

Perché è quello?

+3

Invece di concentrarsi su una particolare funzione in R, sarebbe meglio spiegare perché si pensa che i due attacchi dovrebbero essere molto simili. Per esempio. Dì che LASSO senza penalizzazioni non dovrebbe dare altro se non un adattamento OLS, se è questo che intendi. Puoi anche approfondire il motivo per cui la pensi così usando le formule. –

+0

Ho pensato che fosse abbastanza ovvio LASSO senza penalizzazione e OLS dovrebbe dare gli stessi risultati. Mi stavo chiedendo perché due algoritmi mi danno stime diverse. –

+0

Ciò che è ovvio per te potrebbe non essere ovvio per tutti gli altri, quindi nel caso sia meglio essere il più espliciti e il più precisi possibile. –

risposta

0

Da help glmnet: Si noti inoltre che per "gaussian", glmnet standardizza y per avere varianza unità prima di calcolare la sua sequenza lambda (e quindi non standardizza i coefficienti risultanti); se si desidera effettuare una repro- /confrontare i risultati con altri software, è meglio fornire un y standardizzato.

+0

La differenza tra i coefficienti lm e glmnet si riduce, poiché il valore assoluto dei coefficienti si riduce. Ho ancora le stesse differenze quando non standardizzo i coefficienti. –

+1

C'è un altro avvertimento nel file di aiuto, in particolare la descrizione del parametro lambda, dicendo che l'algoritmo potrebbe avere problemi se si fornisce uno scalare e non un vettore. Non sono sicuro che ciò causi solo un problema di velocità o che in realtà possa influenzare le stime. – tomka

1

Ho eseguito con la "prostata" dataset di esempio del libro di Hastie il codice successivo:

out.lin1 = lm(lpsa ~ . , data=yy) 
out.lin1$coeff    
out.lin2 = glmnet(as.matrix(yy[ , -9]), yy$lpsa, family="gaussian", lambda=0, standardize=T ) 
coefficients(out.lin2) 

e il risultato dei coefficienti sono simili. Quando usiamo l'opzione standardize i coefficienti restituiti da glmnet() sono nelle unità originali delle variabili di input. Si prega di verificare che si stia utilizzando la famiglia "gaussiana"

+0

L'aggiunta di family = "gaussian" non ha modificato il risultato –

+0

È possibile includere il codice R e la manipolazione dei dati? –

3

Si sta utilizzando la funzione errata. x dovrebbe essere la matrice del modello. Non il valore del predittore grezzo. Quando lo fai, ottieni esattamente gli stessi risultati:

x <- rnorm(500) 
y <- rnorm(500) 
mod1 <- lm(y ~ x) 

xmm <- model.matrix(mod1) 
mod2 <- glmnet(xmm, y, alpha=1, lambda=0) 

coef(mod1) 
coef(mod2) 
0

Ho lo stesso problema. Penso che glmnet non possa gestire serie non stazionarie, cioè quando la serie sembra integrata o una passeggiata casuale. Quando simulo i dati stazionari, i risultati di glmnet e OLS sono abbastanza vicini. Ma, in teoria, glmnet con lambda = 0 dovrebbe dare gli stessi risultati di OLS, indipendentemente dal fatto che le serie siano integrate o meno.

Il codice seguente utilizza il tasso di disoccupazione per contea in California dal 1990 al 1999 dallo Bureau of Labor Statistics Local Area database diminuendo il periodo di tempo. Ho inserito una copia CSV di tali dati here per comodità. Il codice regredisce il valore di una contea sui valori passati di tutte le altre contee. La traiettoria della disoccupazione nelle contee n. 30, 34 e 36 (Contea di Orange, Contea di Sacramento e Contea di San Bernardino) sembra integrata. OLS fornisce un coefficiente autoregressivo, il coefficiente di regressione corrispondente ai valori passati della contea sul lato sinistro, per essere inferiore a 1. Ma glmnet restituisce un valore maggiore di 1. Un coefficiente autoregressivo maggiore di 1 potrebbe dare un percorso esplosivo. Il codice enuncia i valori predefiniti per famiglia, standardizzazione, pesi e intercetta. Stabilisce inoltre criteri di convergenza molto più rigidi rispetto ai valori predefiniti.

L'esecuzione di questo codice per la contea # 30 (Contea di Orange) fornisce un coefficiente OLS di 0,9 e un coefficiente LASSO di 1,2. Il suggerimento di not_bonferroni dà esattamente gli stessi risultati (e non si applicherebbe in un problema di Big K con più regressori di osservazioni).

county_wide <- read.csv(file = "county_wide.csv") 

# Problematic counties: #30, #34, #36 
# All three look like their path is integrated rather than stationary 
selected_county <- 30 

# Get dimensions 
num_entities <- dim(county_wide)[2] 
num_observations <- dim(county_wide)[1] 

# Dependent variable: most recent observations of selected county 
Y <- as.matrix(county_wide[1:(num_observations - 1), selected_county]) 

# Independent variables: lagged observations of all counties 
X <- as.matrix(county_wide[2:num_observations, ]) 

# Plot the county to show that it is integrated 
plot(county_wide[, selected_county]) 

# Run OLS, which adds an intercept by default 
ols <- lm(Y ~ X) 
ols_coef <- coef(ols) 

# Control glmnet settings 
glmnet.control(factory = T) 
glmnet.control(fdev = 1e-20) 
glmnet.control(devmax = 0.99999999999999999) 

# run glmnet with lambda = 0 and spelling out the 
# default values for arguments, e.g. intercept 
library("glmnet") 
lasso0 <- glmnet(y = Y, 
       x = X, 
       intercept = T, 
       lambda = 0, 
       weights = rep(1, times = num_observations - 1), 
       alpha = 1, 
       standardize = T, 
       family = "gaussian") 
lasso_coef <- coef(lasso0) 

# compare OLS and LASSO 
comparison <- data.frame(ols = ols_coef, 
         lasso = lasso_coef[1:length(lasso_coef)] 
) 
comparison$difference <- comparison$ols - comparison$lasso 

# Show average difference 
mean(comparison$difference) 

# Show the two values for the autoregressive parameter 
comparison[1 + selected_county, ] 

# Note how different these are: glmnet returns a coefficient above 1, ols returns 
# a coefficient below 1!! 


# not_bonferroni's suggested solution returns exactly the same 
# results with these data 
mod1 <- lm(Y ~ X) 

xmm <- model.matrix(mod1) 
mod2 <- glmnet(xmm, Y, alpha = 1, lambda = 0, intercept = T) 

coef(mod1)[selected_county + 1] # Index +1 for the intercept 
coef(mod2)[selected_county + 2] # Index +2 for the intercepts of OLS and LASSO 
Problemi correlati