Scikit-learn ou statsmodels? Avaliando meu modelo de regressão

O uso das duas ferramentas para avaliar premissas de um modelo de regressão linear

Se você já precisou treinar um modelo de regressão linear em python, é bem provável que conheça os benefícios de pelo menos uma dessas duas ferramentas. Antes de mais nada, o objetivo desse post não é justificar qual a melhor entres elas, mas tentar explorar seu uso diante de um determinado objetivo.

A fim de introduzir um pouco o tema, vou abreviar uma experiência profissional minha, na qual me peguei exatamente nesse dilema!

O primeiro modelo que eu precisei treinar em python “oficialmente” como cientista de dados foi um modelo de regressão linear múltipla. O modelo serviria para testar algumas hipóteses que haviam sido levantadas pelo time de negócios.

Em resumo, a regressão linear múltipla estima um valor esperado de uma variável dependente dados valores de outras variáveis independentes. O modelo visa ajustar uma equação linear aos dados e para isso o método utilizado mais comum é o Método de Mínimos Quadrados — ou Mínimos Quadrados Ordinários (MQO) — que calcula uma reta que se ajuste aos dados com o menor erro possível. A regressão linear múltipla é dada pela fórmula:

Onde Y é variável dependente que o modelo tentará prever; β0 é a constante que representa a interceptação da reta e o valor de Y quando todas as variáveis independentes forem nulas; β1 é a variação esperada em Y dado um incremento unitário em x1, mantendo-se constante as outras variáveis independentes; βk é a variação esperada em Y dado um incremento unitário em xk, mantendo-se constante
as outras variáveis independentes; e ε é o erro residual não explicado pelo modelo, dado pela diferença entre a variável dependente e a variável dependente prevista. Utilizando MQO, os valores de β são estimados de maneira que ε seja o menor possível.

Breve conceito apenas para contextualização, mas aproveito para compartilhar aqui um material do MIT bem completo e resumido sobre regressão.

Pois bem, na questão conceitual eu já tinha uma base em econometria e já havia rodado algumas regressões no R, e foi por isso que meu dilema começou. Como ainda não tinha muita proficiência em python, precisei pesquisar sobre as bibliotecas, e claro que me deparei com nossa queridinha e famosa scikit-learn.

O scikit-learn é uma biblioteca de aprendizado de máquina de software livre para a linguagem de programação python. Amplamente utilizada, a biblioteca possui vários algoritmos de aprendizado supervisionado e não supervisionado, como classificação, regressão e clusterização.

Com modelo treinado e orientada pelo meu chefe, precisei realizar e avaliar diversos testes estatísticos para validá-lo, e foi bem aí — com a cabeça no R ainda — que senti bastante falta do bom e velho summary(model). Nesse momento, um amigo de trabalho e que também estava envolvido no projeto, me apresentou o statsmodels.

Segundo documentação, o statsmodels é um módulo python que fornece várias funções para a estimativa de muitos modelos estatísticos diferentes, assim como testes estatísticos e a exploração de dados estatísticos.

Segundo levantamento feito no GitHub pelo blog do The Data Incubator, algumas tags atreladas as duas ferramentas evidenciam as diferenças entre elas: ao scikit-learn atrelam-se aprendizado de máquina e ciência de dados; e ao statsmodels atrelam-se econometria, modelos lineares generalizados, análise de séries temporais e modelos de regressão.

Em resumo, o statsmodels tem uma pegada mais econométrica, tornando possível realizar vários testes estatísticos necessários para avaliar a relevância dos parâmetros e do modelo em si, de maneira bem simples e com pouquíssimas linhas de código. Com o scikit-learn e com o auxílio de outras ferramentas, também é possível realizar esses testes, mas as vezes de maneira não tão intuitiva e mais “braçal”.

Como a minha demanda conversava muito mais com que o statsmodels poderia me oferecer, eu resolvi usá-lo para continuar a avaliação.

Então, vamos ver como elas funcionam? Let’s code!

Atenção! Como o objetivo desse post é realizar uma avaliação de um modelo de regressão linear utilizando essas duas ferramentas, não vou me estender com exploratória e pré-processamento, mas saibam que ignorar esse processo é um grande equívoco na vida real. Entender, tratar e manipular seus dados corretamente e um pré requisito para um bom desempenho de qualquer modelo.

Para o treino dessa regressão utilizaremos uma base de dados que contém informações sobre vinhos brancos, disponível no UCI Machine Learning Repository. Para esse caso, nossa variável dependente, ou seja, a variável de saída, será ’quality’.

Abaixo as importações necessárias:

Para o treino do modelo, todas as variáveis de entrada — que são numéricas — foram padronizadas, re-escalando os dados de maneira que o conjunto tenha média próxima a 0 e desvio próximo a 1 e amortizando valores discrepantes, para evitar que o algoritmo tenha viés para variáveis com maior ordem de grandeza. Para isso, a função StandardScaler() do pacote preprocessing da biblioteca scikit-learn foi usada, como mostrado abaixo:

Também antes do treino, plotamos uma matriz de correlação para identificar possível correlação significativa entre variáveis de entrada, mapeando assim a possibilidade de multicolinearidade, já que essas variáveis precisam ser linearmente independentes:

De acordo com a matriz, as variáveis independentes 'density' e 'residual_sugar' possuem considerável correlação positiva, mas seguiremos assim e mais à frente realizaremos esse diagnóstico.

Para que realmente exista um bom ajuste, devemos saber que um modelo de regressão linear deve obedecer algumas premissas. São elas:

1. Distribuição normal dos resíduos: Como alguns testes de ajuste assumem que os erros são normalmente distribuídos — como o teste F — , caso essa premissa não seja obedecida, os cálculos dos intervalos de confiança podem não ser fidedignos para as previsões do modelo.
2. Expectativa dos resíduos é 0: Intuitivamente, espera-se que a média dos erros seja 0, já que queremos o menor erro possível.
3. Independência dos resíduos: Não deve existir auto-correlação entre os resíduos, caso contrário, a distribuição do resíduo pode ser afetada, além de poder indicar subestimação/superestimação.
4. Homocedasticidade: Para qualquer valor de X, a variância do resíduo é a mesma, ou seja, os resíduos devem apresentar variação constante. Quando essa premissa não é obedecida, temos a presença de heterocedasticidade, o que torna mais difícil determinar o verdadeiro desvio padrão dos erros, colocando em dúvida o resultado dos intervalos de confiança.

Regressão Linear com scikit-learn

Como mencionado lá em cima, nessa minha experiência, comecei pelo scikit-learn. Aqui farei o mesmo. No código abaixo: definindo variáveis de entrada e saída; instanciando o modelo; ajustando o modelo aos dados; realizando a previsão e calculando o resíduo.

Visualizando os parâmetros (β1…βk) e o intercepto (β0) do modelo:

Calculando o coeficiente de determinação:

Coeficiente de determinação (R²) varia de 0 a 1, e indica o quanto da variabilidade dos dados foi explicada pelo modelo, que nesse caso foi 28%. No entanto, o R² sozinho não indica se um modelo de regressão é ou não adequado. É possível ter um valor baixo de R² para um bom modelo e vice e versa, e por isso é tão importante avaliar os resíduos também. Saiba mais sobre cálculo e avaliação do R² aqui e aqui.

Checando a média e desvio padrão do resíduo:

Plotando a distribuição dos resíduos:

Visualmente a distribuição se aproxima da distribuição normal, mas façamos alguns testes:

Entendendo as funções acima, que foram calculadas utilizando a biblioteca scipy:

1. Skewness: assimetria da distribuição. Quando mais próximo de 0 mais “perfeita” é a assimetria; para valores > 0 existe uma assimetria positiva, e negativa para valores < 0. Também é possível calcular a skewness com um método do pandas chamado skew().
2. Kurtosis: a curtose está associada ao achatamento da distribuição. A curtose de uma distribuição normal é 3. Para valores > 3 a distribuição é mais “alta” que a distribuição normal e para valores < 3 mais “achatada”. Também é possível calcular a curtose com um método do pandas chamado kurtosis(), porém o pandas calcula utilizando a definição de Fisher, que considera a curtose da normal igual a 0. Ou seja, esse resultado será expresso como um “excesso” de curtose, que é o valor do saldo após a subtração por 3. Acima calculamos das duas formas.
3. Teste Jarque-Bera e normaltest: ambos os testes tem como hipótese nula a normalidade. Sendo assim, para valores de p < 0,05 a normalidade é rejeitada.

Por mais que o gráfico de distribuição do resíduo se assemelhe a uma normal, os resultados acima não confirmam essa teoria. Com distribuição bem levemente e positivamente inclinada, talvez possamos dizer que temos uma distribuição próxima a normal, mas não exatamente.

Nesse momento decidi partir para o statsmodels. Além das outras premissas, avaliar multicolinearidade e a significância dos parâmetros da regressão com ele era mais simples.

Regressão Linear com statsmodels

Diferenças em relação ao scikit-learn: a) o statsmodels não adiciona automaticamente a constante (β0) ao modelo, então devemos fazer isso utilizando sm.add_constant(); b) y deve ser passado antes do x dentro do método sm.OLS() que estrutura o modelo.

Como mencionado lá em cima, se você vem do R, se sente em casa com a saída do summary():

Note que, exceto a média e o desvio padrão do resíduo, tudo o que foi explorado e calculado anteriormente (coeficientes, R², funções e testes de assimetria de distribuição do resíduo) se encontram no quadro. Além disso, é possível encontrar outras informações sobre o modelo, como R² ajustado e o teste F. Também é possível avaliar a significância dos coeficientes, observando o erro padrão, o p-valor e o intervalo de confiança para cada um deles.

Mais uma diferença legal de ser pontuada: ao contrário do scikit-learn, no statsmodels não precisamos calcular o resíduo “na mão”. Basta chamá-lo pela propriedade model.resid.

Seguindo com a avaliação, para verificar a independência dos resíduos podemos realizar o teste de Durbin-Watson, também presente no quadro.

Segundo documentação, a estatística do teste é aproximadamente igual a 2 * (1-r), onde r é a autocorrelação da amostra dos resíduos. Para r = 0, indicando que não há correlação, o resultado do teste é 2. O resultado sempre estará entre 0 e 4, no qual quanto mais próximo de 0, mais evidências de correlação positiva e quanto mais próximo de 4, mais evidências de correlação negativa.

Sendo assim, para o nosso caso, o resultado do teste aponta evidências de correlação positiva, o que sugere possibilidade do nosso modelo não ter sido configurado muito bem. Devemos olhar caso a caso, mas uma possibilidade seria avaliar a inserção de variáveis independentes importantes, ou até mesmo avaliar modelos alternativos.

Agora avaliaremos se os resíduos são homocedásticos. Para isso, podemos plotar um gráfico de resíduos versos valores previstos e para essa premissa ser verdadeira, o espalhamento dos resíduos não pode apresentar nenhum tipo de padrão nem tendência em relação a Y estimado. Além disso, podemos realizar o teste Goldfeld-Quandt, que tem como hipótese nula a homoscedasticidade, no qual para valores de p < 0,05 a homocedasticidade é rejeitada.

Com o gráfico observamos que distribuição dos resíduos é similar e não há um maior ou menor espalhamento para alguns valores de Y. Além disso, o resultado do teste confirma a hipótese de homocedasticidade.

Além das premissas enumeradas lá em cima, vamos checar também a presença de multicolinearidade no modelo. Quando estamos trabalhando com uma regressão linear múltipla, é muito importante checar a correlação entre as variáveis independentes, como vimos anteriormente. Caso as variáveis sejam muito correlacionadas entre elas, poderemos ter sérios efeitos nas estimativas de MQO dos coeficientes, no qual os cálculos sobre a relevância e o peso dos parâmetros são feitos com pouca confiança.

A detecção de multicolinearidade pode ocorrer através da VIF (Variance Inflation Factor), além da matriz de correlação. O fator de inflação basicamente mede o grau em que cada variável independente é explicada pelas demais variáveis independentes. Quanto maior for o fator de inflação, mais severa será a multicolinearidade. Veja mais sobre multicolinearidade e VIF aqui e aqui.

Vejamos o código e o resultado:

Uma regra prática e presente na literatura, sugere que um VIF > 10 indica multicolinearidade e pode estar influenciando as estimativas de mínimos quadrados. De acordo com nossa tabela, as variáveis 'density' e 'residual_sugar' encontram-se com VIF alto, exatamente as mesmas variáveis apontadas pela matriz de correlação. Para esse caso, excluiremos a variável 'density' do modelo.

Em seguida, assumindo que não existe mais multicolinearidade, vamos interpretar os resultados, olhando para as estatísticas dos coeficientes. Notamos que na tabela de saída do modelo, cada coeficiente possui seu teste t de significância individual. Esses testes são fundamentais para determinar se a variável independente é realmente importante para o modelo. O teste tem como hipótese nula a irrelevância do parâmetro (β = 0), no qual para valores de p < 0,05 a hipótese nula é rejeitada e podemos manter a variável no modelo. Lembrando que esses testes assumem que os erros são independentes e normalmente distribuídos.

No statsmodels é fácil identificar o p-valor dos coeficientes pois, além de presente na tabela, pode-se chamar a propriedade model.pvalue. Sendo assim, podemos rodar um código que treina o modelo em loops, retirando a cada loop as variáveis com p-valor > 0,05, até “sobrar” um modelo apenas com coeficientes significativos:

Agora temos um modelo sem a variável correlacionada (que já havia sido retirada) e com todos os coeficientes significativos. Lembrando que, como alteramos as variáveis do modelo inicial, deveríamos verificar novamente as análises e testes feitos até aqui. Visualizando o modelo final através do model.summary():

Para fins de exploração do conteúdo, dei continuidade nas avaliações mesmo sob o alarme da possibilidade de o modelo não ser o melhor possível e com brecha para melhorias. Na vida real, deveríamos considerar algumas sugestões, como: transformação de variáveis; interação entre variáveis; detecção de possíveis outliers; e como já apontado, inserção de mais variáveis importantes ou até mesmo avaliar outros modelos alternativos.

O código completo dessa análise encontra-se no meu GitHub.

Ah! Para quem ficou se perguntando se há como calcular o p-valor da significância dos coeficientes para um modelo do scikit-learn, compartilho aqui o link do código. Como disse, é possível, mas com muito mais trabalho.

Conclusões

Com o uso do scikit-learn e statsmodels realizamos algumas importantes avaliações quando o assunto é regressão linear. Para nossos testes, usamos principalmente o statsmodels, por oferecer funcionalidades bastante intuitivas para o que precisamos. De fato, o modelo deixou possibilidades de melhoria e uma mais profunda avaliação quanto a sua abordagem, mas como em todo processo de exploração de dados e modelos, devemos experimentar e testar muitas vezes, e para isso aproveitar o que cada ferramenta tem de melhor.

Obrigada pela leitura e qualquer sugestão é bem-vinda! (:

Referências

http://www.stat.yale.edu/Courses/1997-98/101/linreg.htm
https://blog.thedataincubator.com/2017/11/scikit-learn-vs-statsmodels/
http://labtrop.ib.usp.br/doku.php?id=cursos:planeco:roteiro:07-classrcmdr
https://towardsdatascience.com/verifying-the-assumptions-of-linear-regression-in-python-and-r-f4cd2907d4c0
https://towardsdatascience.com/testing-for-normality-using-skewness-and-kurtosis-afd61be860
https://encurtador.com.br/ehvB7
http://www.portalaction.com.br/analise-de-regressao/25-testes-individuais-e-intervalos-de-confianca-para-os-parametros