Curso Visualização de dados de interesse para a vigilância em saúde

Módulo 6 - Diagnóstico do modelo e avaliação da qualidade do ajuste

Vimos anteriormente que há inúmeras combinações de parâmetros ( $p, d, q$ , ou $P$ , $D$ e $Q$ ) que podemos assumir no processo de modelagem. Assim, surge a necessidade de comparar modelos e verificar qual ou quais se ajustaram bem aos dados.

Uma das formas mais conhecidas de comparação de modelos (e que estávamos fazendo anteriormente, mas sem explicá-la a fundo) é por meio do Critério de Informação de Akaike, conhecido como AIC. Formalmente, essa métrica é definida como:

A I C = - 2 l o g (L) + 2 (p + q + k + 1)

$AIC = - 2log(L) + 2(p+q+k+1)$

Em que $L$ é a verossimilhança do modelo (métrica que quantifica a bondade de ajuste do modelo), p é a ordem da parte autorregressiva, q é a ordem da parte de média móvel e k é definido como o número de parâmetros no modelo.

Ao comparar os modelos por meio do AIC, buscamos sempre o modelo que possui o menor valor de tal métrica. Ao olhar a fórmula, portanto, podemos notar que a métrica penaliza os modelos que possuem maior número de parâmetros (ou seja, que são mais complexos). Embora um modelo com mais parâmetros provavelmente se ajustará melhor aos dados, em alguns casos essa melhora não é relevante o suficiente quando comparada a um modelo mais simples. É esse o objetivo da penalização: optar pela simplicidade quando o ganho em ajuste não é relevante.

Quando ajustamos um modelo ARIMA no R vemos que ele nos fornece o AIC e, também, outras duas métricas. Voltemos a um exemplo da saída do modelo para a série temporal de SRAG:

arima_srag_auto

#> Series: srag_pr_ts 

#> ARIMA(1,0,0)(1,1,0)[12] 

#> 

#> Coefficients:

#>          ar1     sar1

#>       0.7863  -0.4943

#> s.e.  0.0852   0.1408

#> 

#> sigma^2 = 33762:  log likelihood = -319.47

#> AIC=644.95   AICc=645.49   BIC=650.56

O objeto arima_srag_auto armazena o modelo ARIMA ajustado automaticamente para os casos de SRAG no estado do Paraná.

As métricas de qualidade de ajuste aparecem na última linha da saída. Além do AIC, as outras métricas exibidas são o AICc (AIC corrigido) e o BIC (Critério de Informação Bayesiano). Todos seguem a mesma ideia: ao comparar modelos, buscamos diminuir esses valores para termos melhores ajustes. O AICc é similar ao AIC, porém com um termo de correção para casos em que a amostra de dados que temos é reduzida. Assim, penalizam-se modelos muito complexos em casos de pequenas amostras. O BIC também busca aplicar um peso maior a modelos com alta complexidade, tornando a penalização mais severa conforme a amostra aumenta.

Podemos, portanto, nos basear nessas métricas para comparar modelos ajustados à mesma série temporal. Vamos, por exemplo, comparar o ajuste do modelo para SRAG no Paraná com um outro modelo com outro conjunto de parâmetros, como um SARIMA(2,0,1)(1,0,1). Primeiro, vamos ajustá-lo:

arima_srag_2 <- Arima(y = srag_pr_ts, order = c(2,0,1), seasonal = c(1,0,1))

A função Arima() é utilizada para ajustar um modelo ARIMA sazonal à série temporal srag_pr_ts.
- order = c(2, 0, 1) define o modelo não sazonal com:
  . p = 2: dois termos autorregressivos.
  . d = 0: nenhuma diferenciação.
  . q = 1: um termo de médias móveis.
- seasonal = c(1, 0, 1) define o componente sazonal do modelo com:
  . P = 1: um termo autorregressivo sazonal.
  . D = 0: nenhuma diferenciação sazonal.
  . Q = 1: um termo de médias móveis sazonais.
O modelo ajustado é armazenado no objeto arima_srag_2.

Repare que agora usamos a função Arima(), com A maiúsculo, ao invés da arima() tradicional. A Arima() é a função para ajustes de modelos ARIMA do pacote {forecast}, e que possui mais opções, como a inclusão de termos sazonais.

Vamos chamar o modelo para exibir o seu resumo:

arima_srag_2

#> Series: srag_pr_ts 

#> ARIMA(2,0,1)(1,0,1)[12] with non-zero mean 

#> 

#> Coefficients:

#>          ar1      ar2      ma1    sar1     sma1      mean

#>       1.4941  -0.7247  -0.4804  0.3610  -0.1979  357.5365

#> s.e.  0.1504   0.1255   0.2008  0.7201   0.7071   56.0912

#> 

#> sigma^2 = 28749:  log likelihood = -390.97

#> AIC=795.94   AICc=798.1   BIC=810.6

Vemos que segundo os três critérios, o modelo com ajuste automático obtido anteriormente possui melhor ajuste, pois possui menor valor tanto para o AIC quanto para AICc ou BIC.

Os processos de ajuste automático que utilizamos anteriormente, portanto, seguem essa lógica: realiza-se um conjunto de ajustes de diferentes modelos, e estes são comparados de acordo com uma essas métricas (por padrão, o AIC).

Análise de resíduos

Assim como outros tipos de modelos estatísticos, ao escolher um modelo ARIMA para os dados é necessário verificar se o modelo está adequado. Para isso, deve-se atentar se o modelo atende aos pressupostos estabelecidos. Neste caso, o pressuposto que temos é de que os resíduos do modelo sejam independentes e distribuídos aproximadamente por uma distribuição normal, de média zero e uma variância constante.

Resíduos de um modelo são valores que representam a diferença entre os valores ajustados pelo modelo e os valores observados da série temporal. Aqui, chamamos os resíduos de $\alpha_t$ , onde:

α_{t} = z_{t} - {\hat{z}}_{t}

$\alpha_t = z_t - \hat{z}_t$

Ou seja, $\alpha_t$ representa a diferença de um valor observado na série no tempo $t$ ( $z_t$ ) e o que foi ajustado pelo modelo, para esse mesmo momento $t$ ( $\hat{z}_t$ ).

Assim, tem-se que o pressuposto em relação aos resíduos é:

α_{i} \sim N (0, σ_{α}^{2})

$\alpha_i \sim N(0, \sigma^2_{\alpha})$

Queremos que os resíduos ( $\alpha_i$ ) sigam uma distribuição normal centrada em zero, e com variância constante ( $\sigma^2_{\alpha}$ ). Ou seja, nós esperamos que o modelo capture toda a estrutura de autocorrelação dos dados. Ao olhar os resíduos, esperamos não haver nenhuma evidência de autocorrelação restante (chamada de autocorrelação residual).

Para analisar esse pressuposto, geralmente seguimos as seguintes etapas:

Fazemos um gráfico da série dos resíduos $\alpha_t$ , observando se não há uma tendência evidente (ou seja, se é estacionária), e se sua média está em torno de zero.
Se a série dos resíduos for de fato estacionária, calcula-se funções de autocorrelação (ACF) e autocorrelação parcial (PACF).
Se não há autocorrelação evidente, há indícios de que o modelo esteja adequadamente ajustado, pois não há nenhuma estrutura temporal “sobrando” nos dados. Caso haja, então é necessário rever e ajustar um outro modelo para a série trabalhada.

Vamos verificar se o modelo ajustado anteriormente para SRAG está adequado. Para isso, há uma função do pacote {forecast} que realiza esses passos diretamente, chamada checkresiduals(). Veja:

checkresiduals(arima_srag_auto)

Figura 40: Análise de resíduos de um modelo ARIMA para casos de SRAG no Paraná.

#> 

#>  Ljung-Box test

#> 

#> data:  Residuals from ARIMA(1,0,0)(1,1,0)[12]

#> Q* = 10.332, df = 10, p-value = 0.4119

#> 

#> Model df: 2.   Total lags used: 12

A função checkresiduals() é usada para diagnosticar os resíduos de um modelo ajustado, neste caso, o modelo armazenado em arima_srag_auto. Ela verifica se os resíduos do modelo ARIMA seguem as suposições de independência e normalidade.
Essa função gera três gráficos principais:
- A série temporal dos resíduos: Permite verificar se há padrões visíveis nos resíduos ao longo do tempo. O ideal é que eles estejam distribuídos de forma aleatória ao redor de zero e que não apresentem tendências ou variações sistemáticas.
- Um correlograma (ACF) dos resíduos: Mostra a autocorrelação dos resíduos em diferentes lags. O esperado é que os valores fiquem dentro das linhas pontilhadas. Ou seja, espera-se valores das autocorrelações de, no máximo, 0,25.
- Um histograma mostrando a distribuição dos resíduos: Permite avaliar se a distribuição dos resíduos se aproxima de uma distribuição normal.
Além desses gráficos, o teste de Ljung-Box também é aplicado para verificar se há autocorrelação significativa nos resíduos. Um p-valor não significativo (geralmente p>0,05) indica que não há autocorrelação residual relevante, reforçando a adequação do modelo.

O objetivo é verificar a adequação do modelo ajustado, garantindo que ele capturou toda a estrutura da série temporal sem deixar padrões significativos nos resíduos.

Caso ainda existam valores com autocorrelação superior a 0.25 ou -0.25, é indício de que o modelo não capturou toda a estrutura dos dados e não apresentou um bom ajuste. Dessa forma, um novo modelo deve ser ajustado. Isso pode envolver:

Testar outros parâmetros;
Aplicar transformações adicionais nos dados.

Agora vamos interpretar a Figura 40. Vê-se por meio da série dos resíduos (primeiro gráfico) e da distribuição dos resíduos em um histograma (terceiro gráfico) que sua maior concentração de valores está ao redor de zero. Sua variância também não parece aumentar nem diminuir de amplitude ao longo do tempo, o que indica uma variância estável.

Ao olhar para o segundo gráfico (ACF), que indica a autocorrelação dos resíduos, vê-se que não há autocorrelação relevante (acima ou abaixo da linha de referência) em nenhum lag. O p-valor para o teste de Ljung-Box também não foi significativo ao nível de 5%. Portanto, não há evidência de que tenha restado estrutura de autocorrelação temporal nos dados que não foi incorporada pelo modelo. Podemos dizer que o modelo está, portanto, adequado.

Contudo, se padrões diferentes do esperado fossem observados nos resíduos, ou se persistisse a presença de autocorrelação ao olhar o gráfico de ACF, um novo modelo que incorporasse essas estruturas seria necessário.

Com o modelo adequado, podemos seguir para realizar previsões com o modelo!

Modelos de previsão

Uma vez que ajustamos um modelo ARIMA e verificamos a qualidade de seu ajuste para uma série temporal, podemos utilizá-lo para realizar previsões para a série em um determinado número de “passos” (dias, semanas, meses) adiante. As previsões se dão com base nas estruturas modeladas (os termos de autocorrelação, de médias móveis, tendência e sazonalidade dos dados).

Vamos seguir com nosso exemplo de SRAG no Paraná. Anteriormente ajustamos um modelo ARIMA(1,0,0)(1,1,0) que pareceu adequado após a análise de resíduos, portanto vamos utilizá-lo para uma previsão de 12 meses à frente com a função predict():

pred_srag <- predict(arima_srag_auto, n.ahead = 12, se.fit = T) 

pred_srag <- data.frame(pred_srag)

pred_srag

#>         pred       se

#> 1   103.8709 183.7453

#> 2   149.4024 233.7406

#> 3   472.6541 259.8811

#> 4   777.6522 274.8005

#> 5   995.3029 283.6316

#> 6  1013.8055 288.9561

#> 7   588.6373 292.1993

#> 8   482.1176 294.1864

#> 9   345.9675 295.4081

#> 10  294.8000 296.1609

#> 11  260.4514 296.6254

#> 12  139.0994 296.9121

A função predict() é utilizada para gerar previsões baseadas no modelo ARIMA ajustado armazenado em arima_srag_auto. O argumento n.ahead = 12 indica que estão sendo previstas 12 observações à frente da série temporal original. O argumento se.fit = T também solicita os erros padrão das previsões.
O objeto pred_srag armazena o tibble resultante, que contém as previsões e seus respectivos erros padrão. Quando pred_srag é chamado, ele exibe as previsões geradas pelo modelo ARIMA para os próximos 12 períodos e outras informações associadas.
A função data.frame() transforma o objeto pred_srag em um data frame.

Obtemos um data.frame com duas variáveis: as estimativas pontuais para os 12 passos à frente, e o erro padrão associado a essas estimativas. A partir delas, podemos calcular o intervalo de confiança de 95% para a previsão:

pred_srag$inf <- pred_srag$pred - 1.96 * pred_srag$se

pred_srag$sup <- pred_srag$pred + 1.96 * pred_srag$se

Estas linhas de comando estão calculando os intervalos de confiança de 95% para as previsões armazenadas no objeto pred_srag.
- pred_srag$inf: Calcula o limite inferior do intervalo de confiança, subtraindo 1,96 vezes o erro padrão (se) da previsão (pred).
- pred_srag$sup: Calcula o limite superior do intervalo de confiança, somando 1,96 vezes o erro padrão (se) à previsão (pred).

Aqui, utilizamos os valores -1,96 e 1,96 pois são os valores dos percentis 2,5% e 97,5% de uma distribuição normal padrão. Ou seja, a partir deles delimitamos um intervalo de 95% de confiança para nossas estimativas.

Vejamos como a previsão se comporta no gráfico de série temporal na Figura 41.

Figura 41: Série temporal de SRAG no Paraná com previsão para um ano.

Vemos que o número esperado de casos para 2018 representa um padrão semelhante e aceitável em relação ao que foi visto nos anos anteriores, com um pico no inverno. É interessante perceber que o modelo estima um pico ainda maior que o ano anterior (2017). Contudo, os intervalos de confiança da previsão são amplos devida à alta incerteza associada.

<
1
2
3
4
5
6
7
8
Conclusão
>