Econometria: dica muito útil de visualizar resultados com pacote performance

Abstract

The example is a simplified exercise using a linear multiple regression with the package performance.

Licença

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

License: CC BY-SA 4.0

Citação

Sugestão de citação: FIGUEIREDO, Adriano Marcos Rodrigues. Econometria: dica muito útil de visualizar resultados com pacote performance. Campo Grande-MS,Brasil: RStudio/Rpubs, 2021. Disponível em http://rpubs.com/amrofi/dica_performance.

1 Introdução

Usarei o pacote mpg para exemplo do pacote performance de Lüdecke, Ben-Shachar, et al. (2021). Instruções interessantes estão em <https://youtu.be/EPIxQ5i5oxs>, e <https://easystats.github.io/performance/>.

2 Dados e pacote

library(tidyverse)
library(performance)
mpg

3 Modelo com data.frame

Seja um modelo usual, linear de \(hwy=f(disp,class,cyl)\) para os dados do dataset mpg. Posso acessar o \(R^2\) com uso da funçãor2(), ou a performance geral do modelo pela função model_performance(). Esta última fornece mais indicadores do que apenas o \(R^2\) colocado antes.

model_lm <- lm(hwy ~ displ + class + cyl, data = mpg)

model_lm

Call:
lm(formula = hwy ~ displ + class + cyl, data = mpg)

Coefficients:
    (Intercept)            displ     classcompact     classmidsize  
        39.0048          -0.5779          -3.2478          -2.9486  
   classminivan      classpickup  classsubcompact         classsuv  
        -6.9398         -10.2181          -2.6336          -9.0290  
            cyl  
        -1.3306  

r2(model_lm)

# R2 for Linear Regression
       R2: 0.809
  adj. R2: 0.803

model_performance(model_lm)

A função check_model(model_lm) permitirá olhar graficamente:

check_model(model_lm)  # requer o pacote see e ggrepel e qqplotr

Posso comparar com um segundo modelo:

model2 <- glm(am ~ wt + cyl, data = mtcars, family = binomial)
r2(model2)

# R2 for Logistic Regression
  Tjur's R2: 0.705

model_performance(model2)

Podemos comparar os modelos, fazendo (neste caso, fiz um novo modelo lm para comparar modelos de mesma classe):

model_lm2 <- lm(hwy ~ displ + class, data = mpg)
2

[1] 2

model_lm2

Call:
lm(formula = hwy ~ displ + class, data = mpg)

Coefficients:
    (Intercept)            displ     classcompact     classmidsize  
         38.953           -2.298           -5.312           -4.947  
   classminivan      classpickup  classsubcompact         classsuv  
         -8.799          -11.923           -4.699          -10.585  

compare_performance(model_lm, model_lm2)

Podemos também fazer essa comparação fazendo a combinação de funções e o plot. Em todos os casos (indicadores de performance), o ajuste foi melhor para o model_lm.

compare_performance(model_lm, model_lm2, rank = T)

plot(compare_performance(model_lm, model_lm2))

Vou comparar fazendo um modelo de fixed effects, e com datasets iguais. Observe que ele avisará caso os modelos venham de datasets diferentes.

model_lmer <- lme4::lmer(hwy ~ displ + class + cyl + (1 | cty), data = mpg)

model_lmer

Linear mixed model fit by REML ['lmerMod']
Formula: hwy ~ displ + class + cyl + (1 | cty)
   Data: mpg
REML criterion at convergence: 846.322
Random effects:
 Groups   Name        Std.Dev.
 cty      (Intercept) 7.656   
 Residual             1.157   
Number of obs: 234, groups:  cty, 21
Fixed Effects:
    (Intercept)            displ     classcompact     classmidsize  
        31.7544           0.1468          -1.5995          -1.0515  
   classminivan      classpickup  classsubcompact         classsuv  
        -3.3315          -5.7680          -2.0567          -5.1225  
            cyl  
        -0.3089  

r2(model_lmer)

# R2 for Mixed Models

  Conditional R2: 0.979
     Marginal R2: 0.068

model_performance(model_lmer)

plot(compare_performance(model_lm, model_lmer))

Portanto, neste caso, o modelo lmer é melhorar para previsões, conforme valores menores de AIC (868<1130) e BIC (906<1164) para o modelo lmer. Também olhando o R2 ajustado, para o model_lm = 0.803 enquanto o R2 (cond) = 0.979.

compare_performance(model_lm, model_lmer)

check_model(model_lmer)

4 Modelo com tidy

library(tidymodels)

# * Linear Regression ----
model_lm_tidy <- linear_reg() %>%
    set_engine("lm") %>%
    fit(hwy ~ displ + class + cyl, data = mpg)

check_model(model_lm_tidy)

model_lm2_tidy <- linear_reg() %>%
    set_engine("lm") %>%
    fit(hwy ~ displ + class, data = mpg)
compare_performance(model_lm_tidy, model_lm2_tidy, rank = T)

plot(compare_performance(model_lm_tidy, model_lm2_tidy))

check_model(model_lm2_tidy)

5 Testes individuais

É possível realizar alguns testes de pressupostos de modo individualizado, mas ele não mostra os resultados de modo evidente, mas apenas indica se tem problemas. Vejamos.

check_autocorrelation(model_lm)

Warning: Autocorrelated residuals detected (p < .001).

check_collinearity(model_lm)  # VIF

check_heteroscedasticity(model_lm)  # BP test

Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).

plot(check_heteroscedasticity(model_lm))

Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).

check_normality(model_lm)  # shapiro test

Warning: Non-normality of residuals detected (p < .001).

plot(check_normality(model_lm))

Warning: Non-normality of residuals detected (p < .001).

check_outliers(model_lm)

OK: No outliers detected.

plot(check_outliers(model_lm))

check_singularity(model_lm)

[1] FALSE

Vejamos como obter os resultados dos testes. O valor resultante é o p-valor do teste realizado. Normalmente os artigos científicos exigirão também que se relate o valor da estatística de teste.

x <- check_heteroscedasticity(model_lm)

Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).

print(x)

[1] 1.749816e-09
attr(,"object_name")
[1] "model_lm"
attr(,"class")
[1] "check_heteroscedasticity"     "see_check_heteroscedasticity"
[3] "numeric"                     

6 Índice geral de performance

É possível também gerar um indice geral de performance, com modelos de classes diferentes. Vejamos. Sejam quatro modelos distintos: linear, binomial, mixed, poisson.

m1 <- lm(mpg ~ wt + cyl, data = mtcars)  #linear
model_performance(m1)

m2 <- glm(vs ~ wt + mpg, data = mtcars, family = "binomial")  #binominal
model_performance(m2)

library(lme4)
m3 <- lme4::lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)  #mixed
model_performance(m3)

counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome <- gl(3, 1, 9)
treatment <- gl(3, 3)
m4 <- glm(counts ~ outcome + treatment, family = poisson())  #poisson
model_performance(m4)

Posso olhar cada um dos resultados anteriores, ou fazer um índice geral de performance como no chunk a seguir. Veja o aviso de que os datasets de origem são diferentes. Atentar para a instrução do pacote acerca deste indicador: calculation is based on normalizing all indices (i.e. rescaling them to a range from 0 to 1), and taking the mean value of all indices for each model.

compare_performance(m1, m2, m3, m4, rank = TRUE)

A coluna do Performance-Score é a última desta comparação (o leitor precisará passar as colunas para o lado para visualizá-la).

Como ao gerar o html de output ele está formatando toda a tabela de saída do compare_performance, sugiro ao leitor que olhe a saída do print da comparação.

print(compare_performance(m1, m2, m3, m4, rank = TRUE))

# Comparison of Model Performance Indices

Name |   Model |      AIC |      BIC |   RMSE |  Sigma | Performance-Score
--------------------------------------------------------------------------
m2   |     glm |   31.298 |   35.695 |  0.359 |  0.934 |           100.00%
m4   |     glm |   56.761 |   57.747 |  3.043 |  1.132 |            96.21%
m1   |      lm |  156.010 |  161.873 |  2.444 |  2.568 |            92.46%
m3   | lmerMod | 1755.628 | 1774.786 | 23.438 | 25.592 |             0.00%

A interpretação é a seguinte: os indicadores utilizados para comparação, considerando que são modelos distintos, serão neste caso: AIC, BIC, RMSE e Sigma.

O modelo m2 teve menores valores de cada um deles comparativamente aos demais, em todos os casos (Performance-Score =100%). O modelo m4 foi o segundo melhor conform eo Performance-Score = 96.21%.

O plot de teia de aranha pode auxiliar na visualização entre os componentes do indice.

plot(compare_performance(m1, m2, m3, m4, rank = TRUE))

Se comparássemos apenas m2 e m4, outros indicadores entrariam no computo do indicador geral.

print(compare_performance(m2, m4, rank = TRUE))

# Comparison of Model Performance Indices

Name | Model |    AIC |    BIC |  RMSE | Sigma | Score_log | Score_spherical | Performance-Score
------------------------------------------------------------------------------------------------
m2   |   glm | 31.298 | 35.695 | 0.359 | 0.934 |   -14.903 |           0.095 |            66.67%
m4   |   glm | 56.761 | 57.747 | 3.043 | 1.132 |    -2.598 |           0.324 |            33.33%

plot(compare_performance(m2, m4, rank = TRUE))

Referências

Referências de pacotes utilizados - lista todos os pacotes que foram abertos em algum procedimento (mesmo que dentro de rotinas de outros pacotes). Ou seja, chamamos apenas o tidyverse, performance e tidymodels, mas estes usam muitos outros!

Dancho, Matt (2021). easystats: Quickly investigate model performance. <https://www.business-science.io/r/2021/07/13/easystats-performance-check-model.html?utm_source=Business+Science+-+Combined+List&utm_campaign=54dc4a32b4-R_TIPS_041_EASYSTATS&utm_medium=email&utm_term=0_a4e5b7c52f-54dc4a32b4-308012619&mc_cid=54dc4a32b4&mc_eid=b83cd76e52>. Accessed on July, 15th, 2021.

Allaire, JJ, Yihui Xie, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, Hadley Wickham, Joe Cheng, Winston Chang, and Richard Iannone. 2021. Rmarkdown: Dynamic Documents for r. https://CRAN.R-project.org/package=rmarkdown.

Barnier, Julien. 2021. Rmdformats: HTML Output Formats and Templates for Rmarkdown Documents. https://github.com/juba/rmdformats.

Bates, Douglas, and Martin Maechler. 2021. Matrix: Sparse and Dense Matrix Classes and Methods. http://Matrix.R-forge.R-project.org/.

Bates, Douglas, Martin Maechler, Ben Bolker, and Steven Walker. 2021. Lme4: Linear Mixed-Effects Models Using Eigen and S4. https://github.com/lme4/lme4/.

Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.

Bray, Andrew, Chester Ismay, Evgeni Chasnovski, Ben Baumer, and Mine Cetinkaya-Rundel. 2021. Infer: Tidy Statistical Inference. https://CRAN.R-project.org/package=infer.

Henry, Lionel, and Hadley Wickham. 2020. Purrr: Functional Programming Tools. https://CRAN.R-project.org/package=purrr.

Kuhn, Max. 2020a. Dials: Tools for Creating Tuning Parameter Values. https://CRAN.R-project.org/package=dials.

———. 2020b. Modeldata: Data Sets Used Useful for Modeling Packages. https://CRAN.R-project.org/package=modeldata.

———. 2021a. Tune: Tidy Tuning Tools. https://CRAN.R-project.org/package=tune.

———. 2021b. Workflowsets: Create a Collection of Tidymodels Workflows. https://CRAN.R-project.org/package=workflowsets.

Kuhn, Max, and Davis Vaughan. 2021a. Parsnip: A Common API to Modeling and Analysis Functions. https://CRAN.R-project.org/package=parsnip.

———. 2021b. Yardstick: Tidy Characterizations of Model Performance. https://CRAN.R-project.org/package=yardstick.

Kuhn, Max, and Hadley Wickham. 2020. Tidymodels: A Collection of Packages for Modeling and Machine Learning Using Tidyverse Principles. https://www.tidymodels.org.

———. 2021a. Recipes: Preprocessing Tools to Create Design Matrices. https://CRAN.R-project.org/package=recipes.

———. 2021b. Tidymodels: Easily Install and Load the Tidymodels Packages. https://CRAN.R-project.org/package=tidymodels.

Lüdecke, Daniel, Mattan S. Ben-Shachar, Indrajeet Patil, Philip Waggoner, and Dominique Makowski. 2021. “performance: An R Package for Assessment, Comparison and Testing of Statistical Models.” Journal of Open Source Software 6 (60): 3139. https://doi.org/10.21105/joss.03139.

Lüdecke, Daniel, Dominique Makowski, Mattan S. Ben-Shachar, Indrajeet Patil, Philip Waggoner, and Brenton M. Wiernik. 2021. Performance: Assessment of Regression Models Performance. https://easystats.github.io/performance/.

Müller, Kirill, and Hadley Wickham. 2021. Tibble: Simple Data Frames. https://CRAN.R-project.org/package=tibble.

R Core Team. 2021. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Robinson, David, Alex Hayes, and Simon Couch. 2021. Broom: Convert Statistical Objects into Tidy Tibbles. https://CRAN.R-project.org/package=broom.

Silge, Julia, Fanny Chow, Max Kuhn, and Hadley Wickham. 2021. Rsample: General Resampling Infrastructure. https://CRAN.R-project.org/package=rsample.

Vaughan, Davis. 2021. Workflows: Modeling Workflows. https://CRAN.R-project.org/package=workflows.

Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.

———. 2021a. Forcats: Tools for Working with Categorical Variables (Factors). https://CRAN.R-project.org/package=forcats.

———. 2021b. Stringr: Simple, Consistent Wrappers for Common String Operations.

———. 2021c. Tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.

———. 2021d. Tidyverse: Easily Install and Load the Tidyverse. https://CRAN.R-project.org/package=tidyverse.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, Hiroaki Yutani, and Dewey Dunnington. 2021. Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://CRAN.R-project.org/package=ggplot2.

Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2021. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.

Wickham, Hadley, and Jim Hester. 2020. Readr: Read Rectangular Text Data. https://CRAN.R-project.org/package=readr.

Wickham, Hadley, and Dana Seidel. 2020. Scales: Scale Functions for Visualization. https://CRAN.R-project.org/package=scales.

Xie, Yihui. 2014. “Knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC. http://www.crcpress.com/product/isbn/9781466561595.

———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

———. 2016. Bookdown: Authoring Books and Technical Documents with R Markdown. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/bookdown.

———. 2021a. Bookdown: Authoring Books and Technical Documents with r Markdown. https://CRAN.R-project.org/package=bookdown.

———. 2021b. Knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.

Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.

Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.