Regressão Múltipla

Um hospital da Alemanha deseja verificar se o nível de glicose no sangue de uma pessoa pode ser explicada pelas variáveis Idade (anos), Pressão Sanguínea (mm Hg) e IMC (kg/\(m^2\)). A partir de uma amostra de 768 indivíduo ajuste o modelo e avalie a significância das covariáveis.

Inicialmente, iremos fazer uma simples análise exploratória para compreender o comportamento dos dados. Observe, a partir do boxplot, a presenta de valores extremos em todas as variáveis estudadas. Pelo gráfico de dispersão abaixo, aparentemente, glicose e pressão sanguínea apresentam a menor correlação linear. Na verdade, de forma geral, não está muito evidente a relação inidividual de cada covariável com a glicose.

Para ajustar o modelo múltiplo, utilizamos a mesma função lm(). Agora, só precisamos adicionar as regressoras no preditor. Por exemplo, se temos três covariáveis \(X_1\), \(X_2\) e \(X_3\), definimos a função como: lm(Y ~ X1 + X2 + X3), onde o perador + é resposável pela inclusão das regressoras.

## 
## Call:
## lm(formula = Glucose ~ BloodPressure + Age + BMI, data = base.1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -131.784  -19.112   -2.026   18.401   84.459 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   68.18492    5.80137  11.753  < 2e-16 ***
## BloodPressure  0.06018    0.06033   0.998    0.319    
## Age            0.67280    0.09534   7.057 3.82e-12 ***
## BMI            0.81850    0.14390   5.688 1.83e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.13 on 764 degrees of freedom
## Multiple R-squared:  0.1154, Adjusted R-squared:  0.1119 
## F-statistic: 33.22 on 3 and 764 DF,  p-value: < 2.2e-16

Observe que a variável Pressão Sanguínea parece não possuir relação siginificativa com a Glicose. Podemos, então, ajustar um modelo sem esta variável.

## 
## Call:
## lm(formula = Glucose ~ Age + BMI, data = base.1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -132.310  -19.102   -1.977   18.310   83.766 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 70.29517    5.40191  13.013  < 2e-16 ***
## Age          0.69555    0.09257   7.514 1.61e-13 ***
## BMI          0.85891    0.13808   6.220 8.16e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.13 on 765 degrees of freedom
## Multiple R-squared:  0.1142, Adjusted R-squared:  0.1119 
## F-statistic: 49.33 on 2 and 765 DF,  p-value: < 2.2e-16

Note que as outras varíaveis continuam sendo relevantes para explicação da glocose. O coeficiente de determinação ajustado e a variância não se alteram, confirmando que a variável pressão sanguínea não era relevante.

Seleção de variáveis

Iremos comparar os modelos selecionados a partir dos métodos de seleção de variáveis. Para isso, dividiremos a amostra em duas partes: Treino e Teste.

Como exemplo, trabalharemos com os dados Hitters disponíveis na biblioteca ISLR.

O objetivo é prever o salário de um jogador de beisebol com base em várias estatísticas associadas ao desempenho no ano anterior. Primeiro, nota-se que a variável salário está ausente para alguns dos jogadores e precisamos desconsiderar esses \(\textit{missings}\) da base.

Regressão Stepwise

Existem diversas funções no \(\texttt{R}\) para realizar a regressão stepwise.

stepAIC(): disponível no pacote MASS, escolhe o melhor modelo a partir da medida AIC. Possui um argumento chamado direction, que pode assumir os seguintes valores: both (para seleção stepwise); backward (seleção backward); e forward (seleção forward). Retorna o melhor modelo.
regsubsets(): disponível no pacote leaps. Tem o argumento nvmax, que permite especificar o número máximo de preditores que devem ser incorporados no modelo. Retorna vários modelos com tamanhos diferentes até nvmax, isto é, precisamos comparar o desempenho dos diferentes modelos para escolher o melhor. Também tem o argumento method, que pode assumir as opções backward, forward e seqrep (combinação das seleções forward e backward).

Podemos, portanto, utilizar a função stepAIC() selecionar o subconjunto de covariáveis, identificando o melhor modelo a partir da medida AIC. Antes, podemos avaliar o modelo completo, contendo todas as regressoras no preditor. No modelo cheio é possível verificar as diversas variáveis que não parecem ser relevantes na explicação do salário. Vale lembrar que esse tipo de análise deve ser feita com cautela, pois pode existir colinearidade entre as regressoras e não avaliamos esse tipo de problema por enquanto.

## 
## Call:
## lm(formula = Salary ~ ., data = train.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -817.84 -170.49  -25.88  147.72 1867.60 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  177.89960   97.33814   1.828  0.06898 .  
## AtBat         -1.92129    0.68606  -2.800  0.00557 ** 
## Hits           6.74439    2.55568   2.639  0.00892 ** 
## HmRun          2.64519    6.56962   0.403  0.68761    
## Runs          -2.04517    3.15427  -0.648  0.51743    
## RBI            0.16848    2.77824   0.061  0.95170    
## Walks          5.72622    1.96034   2.921  0.00386 ** 
## Years         -6.36321   13.26437  -0.480  0.63191    
## CAtBat        -0.11994    0.14128  -0.849  0.39683    
## CHits          0.15081    0.71905   0.210  0.83408    
## CHmRun         0.24310    1.69601   0.143  0.88616    
## CRuns          1.33193    0.80452   1.656  0.09926 .  
## CRBI           0.50615    0.73307   0.690  0.49065    
## CWalks        -0.87123    0.34100  -2.555  0.01131 *  
## LeagueN       53.96746   88.57419   0.609  0.54297    
## DivisionW   -122.03711   43.66279  -2.795  0.00566 ** 
## PutOuts        0.36277    0.08426   4.305 2.53e-05 ***
## Assists        0.34642    0.23797   1.456  0.14691    
## Errors        -3.29725    4.62703  -0.713  0.47686    
## NewLeagueN   -17.42297   88.36679  -0.197  0.84388    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 320.8 on 216 degrees of freedom
## Multiple R-squared:  0.5344, Adjusted R-squared:  0.4935 
## F-statistic: 13.05 on 19 and 216 DF,  p-value: < 2.2e-16

Abaixo, realizamos a seleção passo a passo. O modelo selecionado tem 10 covariáveis, metade das regressoras observadas. É interessante avaliar a regressão stepwise com a função regsubsets(), onde é possível limitar a quantidade de covariáveis no modelo. Dessa forma, poderíamos comparar os modelos com 1, 2, 3, …, 10 regressoras. É possível que um modelo com menos variáveis seja tão bom quanto o modelo com 10 regressoras e é sempre preferível trabalharmos com modelos parcimoniosos (menor número de variáveis possível).

## 
## Call:
## lm(formula = Salary ~ AtBat + Hits + Walks + CAtBat + CRuns + 
##     CRBI + CWalks + Division + PutOuts + Assists, data = train.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -841.67 -178.01  -35.27  130.89 1909.04 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  162.14550   72.19186   2.246 0.025673 *  
## AtBat         -1.97866    0.58750  -3.368 0.000891 ***
## Hits           6.29063    1.80689   3.481 0.000599 ***
## Walks          5.44405    1.70903   3.185 0.001650 ** 
## CAtBat        -0.10379    0.05966  -1.740 0.083286 .  
## CRuns          1.35473    0.41191   3.289 0.001167 ** 
## CRBI           0.68379    0.22261   3.072 0.002390 ** 
## CWalks        -0.90070    0.27920  -3.226 0.001442 ** 
## DivisionW   -121.16513   42.18545  -2.872 0.004466 ** 
## PutOuts        0.37843    0.08027   4.714 4.25e-06 ***
## Assists        0.24937    0.16918   1.474 0.141880    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 316.1 on 225 degrees of freedom
## Multiple R-squared:  0.529,  Adjusted R-squared:  0.5081 
## F-statistic: 25.27 on 10 and 225 DF,  p-value: < 2.2e-16

Regressão LASSO

Usaremos o pacote glmnet para realizar a regressão lasso. A função principal neste pacote é glmnet(), que pode ser usada para ajustar modelos de regressão lasso, ridge e muito mais. Essa função tem uma sintaxe ligeiramente diferente de outras funções de ajuste de modelo que vimos até agora. Em particular, devemos passar uma matriz \(X\), bem como um vetor \(Y\), e não usamos a sintaxe Y~X. Usaremos novamente os dados Hitters com os valores ausentes removidos.

A função model.matrix() é particularmente útil para criar \(X\), pois não apenas produz uma matriz correspondente aos 19 preditores, mas também transforma automaticamente quaisquer variáveis qualitativas em variáveis . A última propriedade é importante porque glmnet() aceita apenas entradas numéricas e quantitativas.

A função glmnet() tem um argumento alfa que determina que tipo de modelo é ajustado. Se alfa = 1, então um modelo lasso é ajustado, e se alfa = 0, então, um modelo de regressão ridge é ajustado. Por padrão, a função glmnet() padroniza as variáveis para que fiquem na mesma escala.

Por padrão, a função glmnet() executa a regressão para um intervalo selecionado automaticamente de valores \(\lambda\). Podemos selecionar um desses modelos a partir da validação cruzada.

cv.glmnet() é a função principal para fazer validação cruzada, em conjunto com vários métodos de suporte, como plotagem e previsão. Definimos uma semente aleatória para que os resultados obtidos sejam reproduzíveis.

É possível selecionar \(\lambda\) e os coeficientes correspondentes do modelo. O resultado lambda.min é o valor de \(\lambda\) que fornece o erro quadrático médio mínimo.

## [1] 2.347514

Para visualizar os coeficientes do modelo estimado, basta utilizar a função coef().

## 20 x 1 sparse Matrix of class "dgCMatrix"
##                       s0
## (Intercept)  135.4834956
## AtBat         -1.4692550
## Hits           5.0091807
## HmRun          .        
## Runs           .        
## RBI            0.6256107
## Walks          4.3266247
## Years         -8.5990235
## CAtBat         .        
## CHits          .        
## CHmRun         0.5061689
## CRuns          0.8091727
## CRBI           0.2973508
## CWalks        -0.6459869
## LeagueN       33.0396107
## DivisionW   -122.6003028
## PutOuts        0.3602766
## Assists        0.1862368
## Errors        -2.2000969
## NewLeagueN     .

Comparação de modelos

Avaliaremos os modelos ajustados a partir de duas medidas verificam a acurácia de cada modelo, isto é, sua capacidade de acertar na previsão. Além disso, iremos calcular o AIC (medida que compara os modelos em termos de ajuste) dos modelos.

Critério de Informação de Akaike (AIC): medida de avalia a log-verossimilhança e uma penalização dada pela quantidade de parâmetros estimados. O modelo com menor AIC é preferível. \[ \textrm{AIC} = -2\log(\hat{L}) + 2q \] onde \(q\) é a quantidade de parâmetros e \(\hat{L}\) é a verossimilhança maximizada.

Para calcular as medidas de acurácia, as previsões foram feitas utilizando o dado de teste. Assim, o \(Y\) predito é comparado com o \(Y\) da base de teste. As medidas de comparação preditivas utilizadas são listadas abaixo:

Root Mean Squared Error (RMSE): mede a media das diferenças entre o valor verdadeiro e o valor predito. A raiz quadrada faz com que essa medida esteja na mesma escala da variável resposta. O modelo com menor RMSE é preferível. \[ \textrm{RMSE} = \sqrt{\frac{1}{n_{\textrm{teste}}} \sum |y_{\textrm{obs}}- y_{\textrm{pred}}}|^2 \]
Mean Absolute Percentage Error(MAPE): mede a porcentagem do erro em comparação ao valor verdadeiro. Fornece uma medida de erro padronizada. Por exemplo, se o erro fosse 10 e o valor verdadeiro fosse 100, então a porcentagem seria 10%. Por outro lado, se o erro fosse 100 e o valor verdadeiro de \(y\) fosse 1000, a medida ainda seria 10%. O modelo com menor MAPE é preferível. \[ \textrm{MAPE} = \frac{1}{n_{\textrm{teste}}} \sum \frac{|y_{\textrm{obs}}- y_{\textrm{pred}}|}{|y_{\textrm{obs}}|}*100 \]

Tabela - Comparação dos modelos ajustados

Modelo	RMSE	MAPE	AIC
FULL	287.6057	106.7649	3414.7
STEPWISE	288.7209	106.4917	3399.4
LASSO	301.2472	111.0805	–

Modelos Lineares I