class: center, middle, inverse, title-slide # Modelos de Regressão e Previsão (ACA228) ## Diagnóstico e avaliação do modelo. ### Prof. Carlos Trucíos
ctruciosm.github.io
carlos.trucios@facc.ufrj.br
### Faculdade de Administração e Ciências Contábeis, Universidade Federal de Rio de Janeiro --- layout: true <a class="footer-link" href="http://ctruciosm.github.io">ctruciosm.github.io — Carlos Trucíos (FACC/UFRJ)</a> ---
## Introdução Na aula anterior vimos a diferença entre **previsão**, **valores estimados**, **resíduos** e **resíduos de inovação**. -- - **Previsão** ( `\(\hat{y}_{T+h|T}\)`): usamos a informação disponível até o tempo `\(T\)` para prever as observações *fora da amostra* ( `\(T+1, T+3, \ldots, T+h\)`). - **Valores estimados** ( `\(\hat{y}_{t|t-1}\)`): usamos a informação disponível até o tempo `\(t-1\)` para "prever" (na verdade estimar) o que acontece no tempo `\(t\)`. .red[Contudo, todos os parâmetros estimados do modelo são estimados utilizando a informação até o tempo `\(T\)`. Por isso, o valor estimado não é uma previsão]. - **Resíduos:** é a diferença entre o valor verdadeiro e o valor estimado, `\(e_t = y_t - \hat{y}_{t|t-1}\)` - **Resíduos de inovação:** Se a série é transformada (por exemplo `\(y_t* = \log(y_t)\)`), os resídus de inovação são a diferença entre o valor verdadeiro (transformado) e o valor estimados da série transformada `\(\log(y_t) - \hat{y}_t*\)`. Se a série não for transformada, os resíduos e os resíduos de inovação são iguais. -- > Através dos resíduos de inovação saberemos se nosso modelo de previsão é **minimamente aceitável**. --- class: inverse, right, middle # Diagnóstico do modelo Quanto da dinâmica dos dados (dentro da amostra) foi capturada pelo modelo? --- ## Diagnóstico do modelo > Assim como o médico é capaz de fazer um diagnóstico do paciente baseando-se nos resultados dos exames médicos. Nós, somos capazes de fazer um diagnóstico do modelo através da análise dos resíduos de inovação. -- Em geral, queremos que: - **Os resíduos de inovação sejam não correlacionados.** Se os resídios de inovação forem correlacionados, significa que nosso modelo não capturou suficientemente a dinâmica dos dados e que existem informações nos dados que ainda precisam ser exploradas. * Para verificar isto, utilizamos a função de autocorrelação dos resíduos de inovação e testes "Portmanteau". -- - **Os resíduos de inovação tenham média zero**. Se a média não for zero, significa que nossas previsões serão viesadas (ou seja, tendemos a subestimar ou superestimar). * Para verificar isto, podemos fazer um histograma * Também podemos fazer um gráfico sequência dos resíduos de inovação. -- > **Todo modelo que não satisfaz essas duas caracteristicas básicas pode (e deve) ser melhorado!**. Isso não significa que, se os residuais de inovação satisfazem essas caracreristicas o modelo não possa ser melhorado, ele pode! mas essas duas caracteristica são o mínimo aceitavel. --- ## Diagnóstico do modelo Voltando ao nosso conjunto de dados sobre a temperatura de Lajeado/RS... .panelset[ .panel[.panel-name[Passo 0] ```r library(tsibble) library(fpp3) # Importando os dados uri <- "https://raw.githubusercontent.com/ctruciosm/ctruciosm.github.io/master/datasets/lajeado_rs.csv" lajeado_rs <- read.csv(uri, sep = ";") # Passando para um formato de séries temporais (Ver aula anterior) lajeado_rs <- lajeado_rs %>% mutate(ano_mes = yearmonth(ano_mes)) %>% select(ano_mes, temp_media) %>% as_tsibble(index = ano_mes) ``` A função `gg_tsresiduals()` do pacote `feasts` (que já é carregado quando utilizamos o pacote `fpp3`) nos ajudará a fazer o diagnóstico do modelo. A função fornece o **gráfico de sequência**, o **gráfico da função de autocorrelação** e o **histograma** dos resíduos de inovação. ] .panel[.panel-name[Modelo] Ajustando (ou treinando) o modelo ```r fit <- lajeado_rs %>% model(SNAIVE(temp_media~lag(12))) ``` Fazendo previsão h = 12 passos à frente e fazendo o gráfico ```r fit %>% forecast(h = 12) %>% autoplot(lajeado_rs) ``` ] .panel[.panel-name[Gráfico] ![](ACA228_ST_03_files/figure-html/unnamed-chunk-4-1.png)<!-- --> ] .panel[.panel-name[Diagnóstico] ```r fit %>% gg_tsresiduals() ``` ![](ACA228_ST_03_files/figure-html/unnamed-chunk-5-1.png)<!-- --> ] ] --- ## Diagnóstico do modelo: Testes Portmanteau Embora o gráfico da função de autocorrelação (também conhecido como **correlograma**) nos ajude a verificar se os resídios são ou não correlacionados, muitas vezes é dificil verificar isto visualmente. -- - Uma forma mais _formal_ de verificar isto é através de testes de hipóteses. -- - Testaremos um grupo de autocorrelações em lugar de apenas fazer o teste individualmente (daí o nome **portmanteau** que em francês significa _mala_). -- - Testamos se as primeiras `\(k\)` autocorrelações são diferentes do que nós esperamos (esperamos que os resíduais sejam um ruido branco, o que implica que as autocorrelações sejam todas zero). -- - Estudaremos dois testes bastante conhecidos e utilizados quando falamos em séries temporais: O teste **Box-Pierce** e o teste **Ljung-Box**. --- ## Diagnóstico do modelo: Testes Portmanteau .pull-left[ #### Teste Box-Pierce `$$Q = T \displaystyle \sum_{i=1}^k \hat{\rho}_i^2 \sim \chi^2_{k-q}$$` ] .pull-right[ #### Teste Ljung-Box `$$Q^{\ast} = T(T+2) \displaystyle \sum_{i=1}^k (T-i)^{-1} \hat{\rho}_i^2 \sim \chi^2_{k-q}$$` ] em que `\(\hat{\rho}_i = \widehat{\mathbb{Cov}}(\epsilon_t, \epsilon_{t-i})\)`, `\(k\)` é o número de desafagens utilizado (usualmente 10 para dados não sazonais ou `\(2 \times m\)` para dados com período sazonal `\(m\)`) e `\(q\)` é o número de parâmetros estimados no modelo. -- **Regra 1:** rejeitamos `\(H_0\)` se `\(Q > c\)` (ou `\(Q^{\ast} > c\)` dependendo do teste utilzado). Em que `\(c\)` é o quantil 1- `\(\alpha\)` da distribuição `\(\chi^2_{k-q}\)` -- **Regra 2:** rejeitamos `\(H_0\)` se o p-valor é menor do que o nível de significância `\(\alpha\)`. O p-valor é calculado através de `\(P(X > Q)\)` (ou `\(P(X > Q^{\ast})\)` dependo do teste utilizado), em que `\(X \sim \chi^2_{k-q}\)`. -- > Se a série for um ruido branco (que é o queremos dos resíduos), `\(\hat{\rho}_i \approx 0\)` e então `\(Q\)` (ou `\(Q^{\ast}\)`) serão pequenos. Valores muito grandes de `\(Q\)` (ou `\(Q^{\ast}\)`) levarão a rejeitar `\(H_0\)`. --- ## Diagnóstico do modelo: Testes Portmanteau ![](ACA228_ST_03_files/figure-html/unnamed-chunk-6-1.png)<!-- --> Se utilizarmos nossa **Regra 1** rejeitamos `\(H_0\)` se `\(Q\)` (ou `\(Q^{\ast}\)`) cair na área cinza. --- ## Diagnóstico do modelo: Testes Portmanteau No exemplo da temperatura de Lajeado/RS, temos sazonalidade mensal ( `\(m = 12\)`), assim `\(k = 2 \times 12\)` (`lag`). Como o modelo SNAIVE não estima parâmetros (apenas repete as últimas observações do período sazonal), `\(q = 0\)` (`dof`) .panelset[ .panel[.panel-name[Box-Pierce] ```r fit %>% augment() %>% features(.innov, box_pierce, lag = 24, dof = 0) ``` ``` ## # A tibble: 1 × 3 ## .model bp_stat bp_pvalue ## <chr> <dbl> <dbl> ## 1 SNAIVE(temp_media ~ lag(12)) 31.0 0.155 ``` ```r alpha = 0.05 qchisq(1 - alpha, 24) ``` ``` ## [1] 36.41503 ``` ] .panel[.panel-name[Ljung-Box] ```r fit %>% augment() %>% features(.innov, ljung_box, lag = 24, dof = 0) ``` ``` ## # A tibble: 1 × 3 ## .model lb_stat lb_pvalue ## <chr> <dbl> <dbl> ## 1 SNAIVE(temp_media ~ lag(12)) 38.2 0.0329 ``` ```r alpha = 0.05 qchisq(1 - alpha, 24) ``` ``` ## [1] 36.41503 ``` ] .panel[.panel-name[Quem é .innov?] **Lembre-se:** a função `augment()` fornece os valores estimados (.fitted), os resíduos (.resid) e os resíduos de inovação (.innov). E são esses resíduos de inovação que estão entrando como argumento na função `features()`. ```r fit %>% augment() ``` ``` ## # A tsibble: 79 x 6 [1M] ## # Key: .model [1] ## .model ano_mes temp_media .fitted .resid .innov ## <chr> <mth> <dbl> <dbl> <dbl> <dbl> ## 1 SNAIVE(temp_media ~ lag(12)) 2015 Jan 25.6 NA NA NA ## 2 SNAIVE(temp_media ~ lag(12)) 2015 Fev 24.8 NA NA NA ## 3 SNAIVE(temp_media ~ lag(12)) 2015 Mar 24.2 NA NA NA ## 4 SNAIVE(temp_media ~ lag(12)) 2015 Abr 21.2 NA NA NA ## 5 SNAIVE(temp_media ~ lag(12)) 2015 Mai 18.5 NA NA NA ## 6 SNAIVE(temp_media ~ lag(12)) 2015 Jun 16.2 NA NA NA ## 7 SNAIVE(temp_media ~ lag(12)) 2015 Jul 15.9 NA NA NA ## 8 SNAIVE(temp_media ~ lag(12)) 2015 Ago 21 NA NA NA ## 9 SNAIVE(temp_media ~ lag(12)) 2015 Set 17.8 NA NA NA ## 10 SNAIVE(temp_media ~ lag(12)) 2015 Out 19.4 NA NA NA ## # … with 69 more rows ``` ] .panel[.panel-name[Qual teste?] > Embora os dois testes sejam bastante utilizados, o teste de **Ljung-Box é preferido**. ] ] --- ## Diagnóstico do modelo: Portmanteau tests Quando utilizamos `features(.innov, ljung_box, lag = 24, dof = 0)` ou `features(.innov, box_pierce, lag = 24, dof = 0)`, precisamos definir o valor de `dof` (numero de parâmetros que foram estimados). -- **Como sei quantos parametros foram estimados no meu modelo?** -- A função `tidy()` fornece esta resposta. -- .panelset[ .panel[.panel-name[Modelos ajustados] ```r fit <- lajeado_rs %>% model(SNAIVE(temp_media~lag(12))) fit2 <- lajeado_rs %>% model(RW(temp_media ~ drift())) ``` ] .panel[.panel-name[Modelo 2] ```r tidy(fit2) ``` ``` ## # A tibble: 1 × 6 ## .model term estimate std.error statistic p.value ## <chr> <chr> <dbl> <dbl> <dbl> <dbl> ## 1 RW(temp_media ~ drift()) b -0.142 0.289 -0.492 0.624 ``` ] .panel[.panel-name[Modelo 1] ```r tidy(fit) ``` ``` ## # A tibble: 0 × 6 ## # … with 6 variables: .model <chr>, term <chr>, estimate <dbl>, ## # std.error <dbl>, statistic <dbl>, p.value <dbl> ``` ] .panel[.panel-name[Testes no modelo2] ```r fit2 %>% augment() %>% features(.innov, box_pierce, lag = 24, dof = 1) ``` ``` ## # A tibble: 1 × 3 ## .model bp_stat bp_pvalue ## <chr> <dbl> <dbl> ## 1 RW(temp_media ~ drift()) 173. 0 ``` ```r fit2 %>% augment() %>% features(.innov, ljung_box, lag = 24, dof = 1) ``` ``` ## # A tibble: 1 × 3 ## .model lb_stat lb_pvalue ## <chr> <dbl> <dbl> ## 1 RW(temp_media ~ drift()) 212. 0 ``` ] .panel[.panel-name[Modelo 2 Gráfico] ```r fit2 %>% forecast(h = 12) %>% autoplot(lajeado_rs) ``` ![](ACA228_ST_03_files/figure-html/unnamed-chunk-14-1.png)<!-- --> ] ] --- ## Diagnóstico do modelo: Resumo - Quando trabalhamos com séries temporais precisamos verificar se nosso modelo é **minimamente aceitável**. - Istó significa que os resíduos devem ser não correlaconados e ter média zero. - Para ver se essas duas características são satisfeitas ou não temos os gráficos (de sequência, de autocorrelações e histograma) que são obtidos através de ```r fit %>% gg_tsresiduals() ``` - De forma mais formal (e alternativa ao grárico das autocorrelações) são utiizados os testes Box-Pierce ou Ljung-Box, sendo este último preferivel. ```r fit %>% augment() %>% features(.innov, box_pierce, lag = k, dof = q) fit %>% augment() %>% features(.innov, ljung_box, lag = k, dof = q) ``` - Se a média dos resíduos de inovação for `\(\neq 0\)` (digamos `\(m\)`), as previsões serão viesadas. Uma forma simples de contornar este problema é substraindo `\(m\)` às previsões ( `\(\hat{y}_{T+h|T} - m\)`). - Já solucionar o problema de resíduos autocorrelacionados é mais complexo. --- class: inverse, right, middle # Avaliando as previsões --- ## Avaliando as previsões - O diagnóstico do modelo nos diz se o modelo capturou as informações contidas nos dados. **Contudo, não diz nada ao respeito das previões `\(\hat{y}_{T+h|T}\)`** -- - Na prática, podemos ter vários modelos considerados como **minimamente aceitáveis**. -- - Além do diagnostico do modelo, é importante saber qual o desempeneho do modelo _fora da amostra_. Isto só pode ser feito sabendo qual foi o desempenho do modelo com novos dados (dados que o modelo nunca viu). -- - Mas, como fazer isto se utilizamos todos os dados para ajustar (treinar) o modelo? -- - Dividimos os dados em dois conjuntos: dados de treinamento (ou dados _dentro da amostra_ ou _in-sample_) e dados de teste (ou dados _fora da amostra_ ou _out-of-sample_). -- ![](ACA228_ST_03_files/figure-html/unnamed-chunk-17-1.png)<!-- --> --- ## Avaliando as previsões .pull-left[ #### Dados de treinamento ou _in-sample_ - Usados para ajustar (treinar) o modelo ou seja, para estimar os parâmetros do modelo. - O modelo construido com os dados de treinamento é utilizado para fazer previsões (que posteriormente serão avaliadas com os dados de teste). ] -- .pull-right[ #### Dados de teste ou _out-of-sample_ - Simulam os novos dados, **nunca** são utilizados para estimar o modelo mas apenas para avaliar o desempenho do modelo. - São dados que o modelo nunca viu. - Comparamos a previsão feita com os dados de treinamento com as observações verdadeiras. ] -- .center[ .red[ > Na prática 20% ou 30% dos dados é utilizado como conjunto de teste. ]] -- - Dividimos o _dataset_ de `\(T\)` observações em `\(R\)` para dados de treinamento e `\(P\)` para dados de teste ( `\(T = R + P\)`). - Com os dados de treinamento calculamos as previsões `\(\hat{y}_{R+1|R}, \ldots, \hat{y}_{R+P|R}\)`. - Comparamos `\(\underbrace{\hat{y}_{R+1|R}, \ldots, \hat{y}_{R+P|R}}_{\text{Previsões}}\)` com `\(\underbrace{y_{R+1}, \ldots, y_{R+P}}_{\text{Dados observados}}\)` --- ## Avaliando as previsões ```r n <- nrow(lajeado_rs) # Obs. no dataset n_train <- round(n*0.7) # Obs. no conjunto de treinamento (70% train, 30% test) dados_train <- lajeado_rs[1:n_train,] # escolhemos as primeiras n_train observações. dados_teste <- lajeado_rs[-c(1:n_train), ] # escolhemos todas, menos as primeiras n_train observações. ``` -- Uma versão mais sofisticada (e mais usada) para avaliar o desempenho do modelo é atraves de _Time series cross-validation_ (validação cruzada para séries temporais). --- class: inverse, right, middle # Validação cruzada para Séries Temporais --- ## Validação cruzada para Séries Temporais - Não temos mais um único conjunto de dados de teste, mas vários. -- - O desempenho do modelo será calculado como a média do desempenho obtido em cada conjunto de teste. -- ![](ACA228_ST_03_files/figure-html/unnamed-chunk-19-1.png)<!-- --> --- ## Validação cruzada para Séries Temporais .center[ <div class="figure"> <img src="imagens/one_step_ahead.png" alt="Fonte: Livro 'Forecasting: Principles and Practice'" width="80%" /> <p class="caption">Fonte: Livro 'Forecasting: Principles and Practice'</p> </div> ] --- ## Validação cruzada para Séries Temporais .center[ <div class="figure"> <img src="imagens/h_step_ahead.png" alt="Fonte: Livro 'Forecasting: Principles and Practice'" width="80%" /> <p class="caption">Fonte: Livro 'Forecasting: Principles and Practice'</p> </div> ] --- ## Validação cruzada para Séries Temporais - Sejam `\(R^i\)` e `\(P^i\)` ( `\(T = R^i + P^i\)`) o número de observações nos conjuntos de treinamento e teste, respectivamente. ( `\(i = 1, \ldots, k\)`) -- - Suponha que, para cada conjunto de treinamento, fazemos a previsão `\(h = 1, \ldots, j\)` passos à frente. -- - Seja `\(\hat{y}_{R^i+h|R^i}\)` a previsão `\(h\)` passos à frente feita com o `\(i\)`-ésimo conjunto de treinamento. -- - Seja `\(e_{R^i+h} = \hat{y}_{R^i+h|R^i} - y_{R^i + h}\)` o erro de previsão (para a previsão `\(h\)` passos à frente) obtido com o modelo ajustado no `\(i\)`-ésimo conjunto de treinamento. -- .pull-left[ **RMSE** `$$\sqrt{mean(e_{R+h}^2)}$$` ] .pull-right[ **MAE** `$$mean(|e_{R+jh}|)$$` ] em que `\(e_{R+h} = (e_{R^1+h}, e_{R^2+h}, \ldots, e_{R^k+h})\)` é um vetor com todos os erros de previsão `\(h\)` passos a frente (obtidos com todos os conjuntos de treinamento). -- `\(e_{R+h}^2 = (e_{R^1+h}^2, e_{R^2+h}^2, \ldots, e_{R^k+h}^2)\)` e `\(|e_{R+h}| = (|e_{R^1+h}|, |e_{R^2+h}|, \ldots, |e_{R^k+h}|)\)`. -- > Versões análogas do MAE e RMSE podem também ser obtidas para os dados _dentro da amostra_ --- ## Validação cruzada para Séries Temporais .center[ <div class="figure"> <img src="imagens/h_step_ahead.png" alt="Fonte: Livro 'Forecasting: Principles and Practice'" width="80%" /> <p class="caption">Fonte: Livro 'Forecasting: Principles and Practice'</p> </div> ] --- ## Validação cruzada para Séries Temporais No **R**, utilizaremos a função `stretch_tsibble()` que nos ajudará a criar vários conjuntos de treinamento. O argumento `.init` define quantas observações utilizaremos no primeiros conjunto de treinamento e `.step` define o incremento no tamanho de amostra -- Para ilustrar o que acontece quando utilizamos a função `stretch_tsibble()` com `.init = 3` e `.step = 1`. ```r lajeado_rs_tr <- lajeado_rs %>% stretch_tsibble(.init = 3, .step = 1) glimpse(lajeado_rs_tr) ``` ``` ## Rows: 3,157 ## Columns: 3 ## Key: .id [77] ## $ ano_mes <mth> 2015 Jan, 2015 Fev, 2015 Mar, 2015 Jan, 2015 Fev, 2015 Mar,… ## $ temp_media <dbl> 25.6, 24.8, 24.2, 25.6, 24.8, 24.2, 21.2, 25.6, 24.8, 24.2,… ## $ .id <int> 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5,… ``` --- ## Validação cruzada para Séries Temporais Ajustar (treinar) um modelo com apenas 3 observações não parece ser algo muito sensato. Utilizaremos `.init = 50` para nosso exemplo. .panelset[ .panel[.panel-name[Dados de Treinamento] ```r lajeado_rs_tr <- lajeado_rs %>% stretch_tsibble(.init = 50, .step = 1) ``` ] .panel[.panel-name[Avaliando h = 1] ```r lajeado_rs_tr %>% model(RW(temp_media ~ drift())) %>% forecast(h = 1) %>% accuracy(lajeado_rs) ``` ``` ## # A tibble: 1 × 10 ## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1 ## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 RW(temp_media ~ drift(… Test -0.285 2.45 1.98 -2.29 10.5 1.44 1.38 0.356 ``` ] .panel[.panel-name[Avaliando h = 3] ```r lajeado_rs_tr %>% model(RW(temp_media ~ drift())) %>% forecast(h = 3) %>% group_by(.id) %>% mutate(h = row_number()) %>% ungroup() %>% accuracy(lajeado_rs, by = c("h", ".model")) ``` ``` ## # A tibble: 3 × 11 ## h .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1 ## <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 1 RW(temp_media ~ … Test -0.285 2.45 1.98 -2.29 10.5 1.44 1.38 0.356 ## 2 2 RW(temp_media ~ … Test -0.521 4.13 3.59 -4.98 18.8 2.62 2.32 0.756 ## 3 3 RW(temp_media ~ … Test -0.695 5.81 5.10 -8.11 26.9 3.73 3.27 0.793 ``` ] ] --- ## Exercícios Utilizando os dados de lajeado/RS (disponíveis [aqui](https://raw.githubusercontent.com/ctruciosm/ctruciosm.github.io/master/datasets/lajeado_rs.csv)) -- - Importe os dados (ver slide 5) - A variável `ano_mes` ainda não está como estrutura tempora. Transforme ela para formato temporal com a função `yearmonth` e convirta seu _dataset_ para um _dataset_ em formato de série temporal (ver slide 5) - Utilize a função `stretch_tsibble()` para criar vários conjuntos de treinamento. Defina `.init = 50` e `.step = 1`. (ver slide 22) - Ajuste um modelo `SNAIVE()` (`SNAIVE(temp_media~lag(12))`), faça a previsão `\(h = 3\)` passos à frente (para cada conjunto de treinamento) e finalmente avalie o desempenho do modelo com a função `accuracy()`. (ver slide 22) - Qual modelo teve um melhor desempenho _fora da amostra_, o modelo `SNAIVE` ou o modelo desenvolvido na aula? (compare o RMSE e o MAE). -- > **Entregar até antes de começar a próxima aula!** (carregar o código pelo Google Class) --- ### Referências - [Hyndman, R.J., & Athanasopoulos, G. (2021). Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia. OTexts.com/fpp3.](https://otexts.com/fpp3/). Chapters 5.4, 5.8, 5.10 - Shumway R. H., & Stoffer, D. S. (1999). Time Series Analysis and Its Applications. Springer. Chapter 1 -- Interessados em conhecer mais medidas para avaliar o desempenho do modelo? ver o artigo do Rob Hyndman. - Hyndman, R. J., & Koehler, A. B. (2006). Another look at measures of forecast accuracy. International Journal of Forecasting, 22(4), 679-688.[Link Artigo](https://www.sciencedirect.com/science/article/pii/S0169207006000239) -- [Link WP](https://robjhyndman.com/papers/mase.pdf)