Conjunto de dados

Os dados considerados para aplicação do método é denominado por longley, disponível no pacote datasets. Este conjunto de dados refere-se a um estudo de macroeconomia onde temos 7 variáveis avaliadas em 16 anos consecutivos, de 1947 a 1962.

No R o conjunto de dados tem estrutura conforme output abaixo:

str(longley)

## 'data.frame':    16 obs. of  7 variables:
##  $ GNP.deflator: num  83 88.5 88.2 89.5 96.2 ...
##  $ GNP         : num  234 259 258 285 329 ...
##  $ Unemployed  : num  236 232 368 335 210 ...
##  $ Armed.Forces: num  159 146 162 165 310 ...
##  $ Population  : num  108 109 110 111 112 ...
##  $ Year        : int  1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 ...
##  $ Employed    : num  60.3 61.1 60.2 61.2 63.2 ...

O modelo considerado aqui será o modelo gaussiano considerando Employed como variável de interesse e todas as demais como explicativas. Conforme Abaixo

\[ Employed \mid X \sim Normal(\hat{\mu}, \hat{\sigma}^2) \\ \hat{\mu} = X\underline{\beta} \]

onde a matriz X compreende todas as 6 variáveis explicativas e portanto tem dimensão 16x7 (6 \(\beta\)’s associados as covariáveis e intercepto \(\beta_0\)). Portanto para o modelo descrito serão estimados 8 parâmetros (7 de regressão e o \sigma^2).

Abaixo ajustamos os modelos:

• m0lm: Modelo gaussiano de regressão múltipla;
• m0ridge: Modelo gaussiano de regressão múltipla com penalização Ridge;
• m0lasso: Modelo gaussiano de regressão múltipla com penalização Lasso; e
• m0elnet: Modelo gaussiano de regressão múltipla com penalização dada por \(\alpha = 0.5\)

##----------------------------------------------------------------------
## Preditor adotado
preditor <- Employed ~ .
data <- longley

##-------------------------------------------
## Via LM (Linear Models)
m0lm <- lm(preditor, data = data)

##-------------------------------------------
## Via Reguralização
library(glmnet)

## É ncessário informar as matrizes X e y
X <- model.matrix(preditor, data = data)
y <- longley$Employed

## Regularização Ridge
m0ridge <- cv.glmnet(x = X, y = y, family = "gaussian", alpha = 0,
                     grouped = FALSE)

## Regularização Ridge
m0lasso <- cv.glmnet(x = X, y = y, family = "gaussian", alpha = 1,
                     grouped = FALSE)

## Regularização Elastic Net (com alpha = 0.5)
m0elnet <- cv.glmnet(x = X, y = y, family = "gaussian", alpha = 0.5,
                     grouped = FALSE)

Os objetos m0ridge, m0lasso e m0net tem classe cv.glmnet e para estas classes o pacote dispõe de vários métodos. Utilizaremos alguns deles para verificar o ajuste do modelos modelos. Note que utilzamos a função cv.glmnet o cv vem das iniciais Cross Validation, utilizando esta função o próprio algoritmo programada já sugere a melhor opção para \(\lambda\) a ser utilizada, segundo o critério implementado (Erro quadrático médio da validação cruzada).

Escolhendo a penalização adequada

Essas funções de estimação providas no pacote glmnet ajustam, por padrão, modelos com diferentes \(\lambda\)’s na penalização. É utilizado uma sequência de 100 valores escolhidos com base na escala da variável resposta e no método utilizado ($). Nos objetos de classe cv.glmnet temos a função plot que apresenta o critério utilizado para escolha do \(\lambda\).

par(mfrow = c(1, 3))
plot(m0ridge); title("Ridge", line = 2.5)
plot(m0lasso); title("Lasso", line = 2.5)
plot(m0elnet); title("alpha = 0.5", line = 2.5)

Com base nestes gráficos nota-se que o valor de \(\lambda\) sugerido, que minimiza o erro de validação cruzada é 0,3671, 0,0041, 0,0082 para os modelos com penalização Ridge, Lasso, e Elastic Net (\(\lambda = 0.5\)), ainda é apresentado (segunda linha vertical) o maior valor de \(\lambda\) que está no intervalo de um erro padrão do erro quadrático médio de validação que resultou nos valores 0,6414, 0,022, 0,0401.

Comparação de ajustes

Abaixo como comparação apresentamos os parâmetros estimados em ambos os métodos.

Primeiro apresentamos abaixo a comparação dos coeficientes estimados utilizando o \(\lambda\) que minimiza o erro de validação cruzada e o \(\lambda\) máximo pertencente ao intervalo de um erro padrão do erro de validação cruzada.

## cbind(coef(m0lm), coef(m0ridge), coef(m0lasso), coef(m0elnet))

coefRidge <- cbind("RIDGE-MIN" = coef(m0ridge, s = "lambda.min")[-2],
                   "RIDGE-1SE" = coef(m0ridge, s = "lambda.1se")[-2])
coefLasso <- cbind("LASSO-MIN" = coef(m0lasso, s = "lambda.min")[-2],
                   "LASSO-1SE" = coef(m0lasso, s = "lambda.1se")[-2])
coefElnet <- cbind("ELNET-MIN" = coef(m0elnet, s = "lambda.min")[-2],
                   "ELNET-1SE" = coef(m0elnet, s = "lambda.1se")[-2])

compareCoefs <- cbind(coefRidge, coefLasso, coefElnet)
rownames(compareCoefs) <- coef(m0elnet)@Dimnames[[1]][-2]

knitr::kable(compareCoefs, align = rep("c", 6), digits = 3)

Note na tabela que os valores, primeiramente olhando as colunas duas a duas (se foi pego o \(\lambda\) mínimo - -MIN ou o maior dentro do intervalo de 1 erro padrão - -1SE), dentre essas duas alternativas não vemos uma diferença muito discrepante das estimativas, porém perceba que sempre em -1SE temos uma estimativa mais próxima de zero do que a fornecida por -MIN, pois nestes casos a penalização é maior. Escolheremos o maior \(\lambda\), pois desejamos uma maior penalização. Agora comparando entre procedimentos de penalização RIDGE, LASSO e ELNET (com \(\lambda\) = 0,5) notamos maior diferença nas estimativas, note que no caso RIDGE as estimativas chegam próxima a zero, porém não são nulas, já nos casos LASSO e ELNET a nulidade já ocorre para alguns parâmetros.

Agora vamos a comparação com a tradicional metodologia já amplamente utilizada, a regressão gaussiana múltipla. Abaixo temos a comparação dos coeficientes estimatidas. Apresentamos somente as estimativas com penalização utilizando o maior \(\lambda\) sugerido.

coeflm <- cbind("LM" = coef(m0lm))
compareCoefs <- cbind(coeflm, compareCoefs[, seq(2, 6, 2)])
colnames(compareCoefs) <- c("LINEAR", "RIDGE", "LASSO", "ELNET")

knitr::kable(compareCoefs, align = rep("c", 4), digits = 3)

Os valores estimados são razoavelmente distintos, porém note que a regressão LASSO toma a nulidade do parâmetros relacionados as variáveis GNP.deflator, GNP e Population. Verificando o teste de significância para estes parâmetros no modelo linear gaussiano temos:

## summary(m0lm)

## Para usar kable
est = coef(summary(m0lm))[, "Estimate"]
sdt = coef(summary(m0lm))[, "Std. Error"]
tval = est/sdt
pval = 2*pt(abs(tval), df = m0lm$df.residual, lower = FALSE)

mysummary <- cbind(Estimate = est, Std.Error = sdt,
                   "t value" = tval, "P(>|t|)" = pval)

knitr::kable(mysummary, align = rep("c", 4), digits = 3)

Portanto, nota-se que a regressão por penalização LASSO fez, neste caso, o trabalho manual que seria de retirar estas as variáveis do modelo e reajustá-lo.