Sabatina,

assunto: Regularização

Abril de 2016


Banco de dados


data("longley")

Um dataset de 7 variáveis econômicas, observadas de 1947 até 1962 (n = 16)

  • GNP.deflator
    • deflator de preços implícito do GNP
  • GNP
    • Produto Nacional Bruto
  • Unemployed
    • Número de desempregados (1954 = 100)
  • Armed.Forces
    • número de pessoas nas forças armadas
  • Population
    • população ‘não institucionalizada’ \(\geq\) 14 anos de idade
  • Year
    • a ano (tempo)
  • Employed
    • número de pessoas empregadas

summary(longley)
  GNP.deflator         GNP          Unemployed     Armed.Forces     Population         Year         Employed    
 Min.   : 83.00   Min.   :234.3   Min.   :187.0   Min.   :145.6   Min.   :107.6   Min.   :1947   Min.   :60.17  
 1st Qu.: 94.53   1st Qu.:317.9   1st Qu.:234.8   1st Qu.:229.8   1st Qu.:111.8   1st Qu.:1951   1st Qu.:62.71  
 Median :100.60   Median :381.4   Median :314.4   Median :271.8   Median :116.8   Median :1954   Median :65.50  
 Mean   :101.68   Mean   :387.7   Mean   :319.3   Mean   :260.7   Mean   :117.4   Mean   :1954   Mean   :65.32  
 3rd Qu.:111.25   3rd Qu.:454.1   3rd Qu.:384.2   3rd Qu.:306.1   3rd Qu.:122.3   3rd Qu.:1958   3rd Qu.:68.29  
 Max.   :116.90   Max.   :554.9   Max.   :480.6   Max.   :359.4   Max.   :130.1   Max.   :1962   Max.   :70.55  
lattice:: splom(longley, pscales = 0, col = 1, type = c("p", "g", "smooth"), xlab = NULL
                , main = "Todos os possíveis gráficos de dispersão 2 x 2"
                , sub = "Curvas de tendência estimadas por suavização loess")

No gráfico acima podemos observar que várias variáveis são altamente correlacionadas, como por exemplo: Employed x Year, Year x Population, Employed x GNP, Year x GNP, Population x GNP, Year x GNP.deflator, Population x GNP.deflator, GNP x GNP.deflator


Regularização


library(glmnet)

Por default a função glmnet utiliza um modelo linear Gaussiano


matriz de covariáveis

(x <- model.matrix(Employed ~ ., longley)[, -1])
     GNP.deflator     GNP Unemployed Armed.Forces Population Year
1947         83.0 234.289      235.6        159.0    107.608 1947
1948         88.5 259.426      232.5        145.6    108.632 1948
1949         88.2 258.054      368.2        161.6    109.773 1949
1950         89.5 284.599      335.1        165.0    110.929 1950
1951         96.2 328.975      209.9        309.9    112.075 1951
1952         98.1 346.999      193.2        359.4    113.270 1952
1953         99.0 365.385      187.0        354.7    115.094 1953
1954        100.0 363.112      357.8        335.0    116.219 1954
1955        101.2 397.469      290.4        304.8    117.388 1955
1956        104.6 419.180      282.2        285.7    118.734 1956
1957        108.4 442.769      293.6        279.8    120.445 1957
1958        110.8 444.546      468.1        263.7    121.950 1958
1959        112.6 482.704      381.3        255.2    123.366 1959
1960        114.2 502.601      393.1        251.4    125.368 1960
1961        115.7 518.173      480.6        257.2    127.852 1961
1962        116.9 554.894      400.7        282.7    130.081 1962

vetor resposta

y <- longley$Employed

Na função glmnet temos o argumento alpha

  • alpha = 0, regressão ridge

  • alpha = 1 (default), regressão lasso

Pra cada tipo de regressão (lasso e ridge) são ajustados três modelos, cada um com um \(\lambda\) diferente, e assim verificamos como isso impacta a estimação

Por default é usada uma sequência de 100 valores pra \(\lambda\), onde o maior valor é o menor valor para o qual todos os coeficientes são zero, e o menor é 0.0001, caso o número de observações seja maior que o número de variáveis, como acontece aqui


lasso


fit <- glmnet(x, y)
fit2 <- glmnet(x, y, lambda = seq(3, .001, length.out = 100))
fit3 <- glmnet(x, y, lambda = seq(3, .01, length.out = 100))

par(mfrow = c(1, 3))
plot(fit, xvar = "lambda", ylab = "Coeficientes", las = 1, sub = "Lambda default")
plot(fit2, xvar = "lambda", ylab = "Coeficientes", las = 1, sub = "Lambda: seq(3, 0.001, length.out = 100)")
plot(fit3, xvar = "lambda", ylab = "Coeficientes", las = 1, sub = "Lambda: seq(3, 0.01, length.out = 100)")

No gráfico da esquerda é usado o \(\lambda\) default

Na abscissa temos o log de \(\lambda\), na ordenada temos os valores do coeficientes

Cada curva é uma variável e no topo dos gráficos temos o número de coeficientes diferentes de zero

O maior log de \(\lambda\) é próximo de 1, logo, o maior \(\lambda\) é próximo de 3

3 é o maior \(\lambda\) considerado nos outros modelos

O número de \(\lambda\)s é fixo, 100

O que foi aumentando (leia-se, ficando mais próximo de zero) é o \(\lambda\) mínimo

No gráfico do meio ele é de 0.001 e no gráfico da direira é de 0.01, lembrando que o default é de 0.0001, ou seja, fomos aumentando o \(\lambda\) mínimo por uma razão de 10


Observações:

No \(\lambda\) default e no \(\lambda\) mínimo de 0.001 acabamos com todos os coeficientes diferentes de zero, quando aumentamos o \(\lambda\) mínimo pra 0.01 terminamos com 3 coeficientes de zero


ridge


fit4 <- glmnet(x, y, alpha = 0)
fit5 <- glmnet(x, y, alpha = 0, lambda = seq(3000, .001, length.out = 100))
fit6 <- glmnet(x, y, alpha = 0, lambda = seq(3000, .01, length.out = 100))

par(mfrow = c(1, 3))
plot(fit4, xvar = "lambda", ylab = "Coeficientes", las = 1, sub = "Lambda default")
plot(fit5, xvar = "lambda", ylab = "Coeficientes", las = 1, sub = "Lambda: seq(3000, 0.001, length.out = 100)")
plot(fit6, xvar = "lambda", ylab = "Coeficientes", las = 1, sub = "Lambda: seq(3000, 0.01, length.out = 100)")

O maior log de \(\lambda\) é próximo de 8, logo, o maior \(\lambda\) é próximo de 3000

3000 é o maior \(\lambda\) considerado nos outros modelos


Observações:

Para todos os \(\lambda\)s terminamos com todos os coeficientes diferentes de zero, na verdade os coeficientes se mostram diferentes de zero em todo o range da sequência de \(\lambda\)


E olhando para as estimativas dos coeficientes? (pra ambas regressões)


round(cbind(coefficients(fit, s = 0.0001), coefficients(fit4, s = 0.0001)), 5)
7 x 2 sparse Matrix of class "dgCMatrix"
                       1          1
(Intercept)  -2776.77868 -370.26268
GNP.deflator    -0.02405    0.08373
GNP             -0.00743    0.01072
Unemployed      -0.01615   -0.00687
Armed.Forces    -0.00932   -0.00166
Population      -0.19185    0.12030
Year             1.47226    0.21049
round(cbind(coefficients(fit2, s = 0.001), coefficients(fit5, s = 0.001)), 5)
7 x 2 sparse Matrix of class "dgCMatrix"
                       1           1
(Intercept)  -2546.62372 -2531.12692
GNP.deflator    -0.00981    -0.01257
GNP             -0.00638    -0.00550
Unemployed      -0.01581    -0.01570
Armed.Forces    -0.00899    -0.00897
Population      -0.15003    -0.15417
Year             1.35094     1.34321
round(cbind(coefficients(fit3, s = 0.01), coefficients(fit6, s = 0.01)), 5)
7 x 2 sparse Matrix of class "dgCMatrix"
                       1           1
(Intercept)  -1758.48727 -1345.39001
GNP.deflator     .           0.03914
GNP              .           0.00875
Unemployed      -0.01377    -0.01274
Armed.Forces    -0.00678    -0.00730
Population       .          -0.04652
Year             0.93628     0.72385

Pegando o \(\lambda\) de 0.0001

Quando mudamos o tipo de regressão e deixamos o \(\lambda\) default as estimativas diferem consideravelmente! Maaaas alguns padrões ainda podem ser traçados, como por exemplo: O maior coeficiente ainda é o da variável Year, seguido de Population (com um detalhe, na regressão lasso o coeficiente é negativo, na regressão ridge ele é positivo!)

Com base nas estimativas podemos dizer que existe uma relação positiva com Year e uma relação negativa com Unemployed e Armed.Forces, i.e, conforme o ano passa o número de pessoas empregadas aumenta (Employed) e conforme o desemprego e o número de pessoas nas forças armadas aumenta, o número de pessoas desempregados diminui. Contudo, esses dois últimos coeficientes são pequenos, logo, essa relação não é muito grande

Pegando o \(\lambda\) de 0.001

Aqui as estimativas são muito mais próximas. O maior coeficiente ainda pertence ao Year e a menor a Popupation, i.e., conforme a população aumenta o número de pessoas empregadas diminui. Com exceção da variável Year, conforme as variáveis aumentam o valor da resposta diminui

Pegando o \(\lambda\) de 0.01

Na regressão lasso apenas três coeficientes foram diferentes de zero, Unemployed, Armed.Forces e Year. Na regressão ridge todos foram diferentes de zero. Nessas três variáveis existe uma certa conformidade de valores


Considerações:

A regressão ridge é conhecida por falhar na parcimônia do modelo, incluindo um número maior de variáveis no modelo. Isso pode justificar o comportamento observado com \(\lambda\) = 0.01

Em geral, quando não deixamos o \(\lambda\) no default os resultados obtidos com as duas regressões são similares, apontando para as mesmas direções

Lembrando que esse conjunto de dados é famoso por sua alta colinearidade

Com \(\lambda\)s maiores as regressão divergem em número de variáveis diferentes de zero, com \(\lambda\)s menores algumas estimativas se mostram diferentes e com \(\lambda\)s intermediários (\(\lambda\) de 0.001, nesse caso), as estimativas se mostram muito mais similares