ANÁLISE DA SÉRIE TEMPORAL DO NÚMERO DE INTERNAÇÕES POR PROBLEMAS RESPIRATÓRIOS NO RIO DE JANEIRO
Relatório final do trabalho da disciplina de Séries Temporais (GET00127)
Resumo
Este trabalho teve como objetivo analisar a série temporal do número de internações hospitalares por doenças do aparelho respiratório no estado do Rio de Janeiro, entre janeiro de 2010 e dezembro de 2019. Com base nessa série, aplicaram-se métodos estatísticos clássicos de modelagem, como o modelo SARIMA e o modelo de Holt-Winters, com o intuito de capturar padrões estruturais e realizar previsões. A análise exploratória evidenciou forte sazonalidade e possível tendência decrescente ao longo do tempo. O modelo SARIMA \((1,1,1) \times (2,1,1)_{[12]}\) foi selecionado com base no menor valor do Critério de Informação Bayesiano (BIC) entre os modelos testados, apresentando bom desempenho preditivo e resíduos compatíveis com os pressupostos do modelo. A comparação com o modelo Holt-Winters mostrou que o SARIMA foi superior em termos de acurácia preditiva, com menor erro absoluto médio (MAE), menor erro percentual absoluto médio (MAPE) e menor autocorrelação dos resíduos. Os resultados obtidos reforçam a utilidade dos modelos de séries temporais como ferramenta de apoio à gestão em saúde pública, possibilitando projeções confiáveis para o planejamento hospitalar.
Introdução
Internações hospitalares por problemas respiratórios representam uma demanda significativa para o sistema público de saúde, exigindo planejamento adequado de recursos humanos, leitos, equipamentos e insumos. Antecipar essas demandas por meio de modelos preditivos pode fornecer subsídios valiosos para a gestão hospitalar e para políticas públicas mais eficazes, especialmente em contextos de alta sazonalidade e variações históricas consistentes.
A análise de séries temporais é uma ferramenta estatística fundamental nesse processo, pois permite identificar padrões estruturais nos dados como tendência, sazonalidade e ciclos, além de viabilizar a construção de modelos com capacidade preditiva. Métodos clássicos, como os modelos SARIMA, têm sido amplamente utilizados para esse tipo de tarefa buscando capturar padrões complexos e não lineares.
Este trabalho tem como objetivo realizar uma análise da série temporal do número de internações por doenças respiratórias no estado do Rio de Janeiro. Pretende-se investigar a estrutura da série, identificar suas componentes e comparar o desempenho de diferentes métodos de previsão, contribuindo para o aprimoramento de estratégias de alocação de recursos na saúde pública.
Materiais e Metodologias
Modelos SARIMA
Os modelos ARIMA (p,d,q) quando contém uma parcela sazonal (autorregressiva e de média móvel) se transformam em SARIMA (Seasonal Autorregressive Integrated Moving Average) (p,d,q)(P,D,Q)s, o qual considera a sazonalidade do tipo multiplicativa.
Ao analisar a notação dos modelos SARIMA, a parte “AR” representa a quantidade “p” de termos autorregressivos, ou seja, se for um AR (1), há um termo defasado em um período, e “p” é igual a 1. A parte “I” significa que a série precisa ser diferenciada “d” vezes para tornar-se estacionária. Já “MA” representa a quantidade de termos de média móvel “q”.
No caso da sazonalidade, pode ser necessária ainda uma diferenciação sazonal “D” e a inclusão de termos autorregressivos sazonais “P” e de médias móveis sazonais “Q”, com s=12, para o caso mensal, usualmente.
O modelo geral de sazonalidade multiplicativa SARIMA (p,d,q)(P,D,Q)s é dado pela equação a seguir, com base em Nelson, na qual se supõe que os resíduos \(a_t\) tenham média zero, variância constante, sejam normais e não autocorrelacionados serialmente:
Dizemos que \(\{Z_t\}\) é um modelo \(\text{ARIMA}(p,d,q) \times (P,D,Q)_s\) com período sazonal \(s\), sendo definido da seguinte forma:
\[ \Phi_P(B^s)\, \phi_p(B)\, (1 - B)^d\, (1 - B^s)^D\, Z_t = \theta_q(B)\, \Theta_Q(B^s)\, a_t \]
Onde:
- \(B\) é o operador defasagem (\(Z_{t-1} = BZ_{t}\));
- \(\phi_p(B)\) e \(\theta_q(B)\) são os polinômios autoregressivo e de médias móveis, respectivamente;
- \(\Phi_P(B^s)\) e \(\Theta_Q(B^s)\) são os polinômios sazonais de autoregressão e médias móveis;
- \(d\) e \(D\) são os graus de diferenciação não sazonal e sazonal;
- \(a_t\) é o ruído branco.
Modelo de Holt-Winters
O método de Holt-Winters (HW) é uma extensão dos modelos de suavização exponencial, adequado para séries temporais que apresentam não apenas tendência, mas também comportamento sazonal. Ele é amplamente utilizado por sua simplicidade e capacidade de fornecer previsões eficientes mesmo em contextos com padrões complexos.
Esse modelo considera três componentes principais da série: nível, tendência e sazonalidade, que são suavizados por meio de equações recursivas. Existem duas variações principais do modelo: com sazonalidade aditiva e com sazonalidade multiplicativa. A escolha entre elas depende da forma como o efeito sazonal se manifesta ao longo da série.
No caso da sazonalidade aditiva, assume-se que o padrão sazonal se mantém constante ao longo do tempo, independentemente do nível da série. Já na sazonalidade multiplicativa, o componente sazonal varia proporcionalmente ao nível da série, por exemplo, o efeito do aumento nas vendas durante o Natal tende a ser mais pronunciado quando o volume total de vendas também é maior.
A equação geral de previsão do modelo Holt-Winters com sazonalidade multiplicativa é dada por:
\[ \hat{Z}_{t+h} = (\bar{Z}_t + h \cdot \hat{T}_t) \cdot \hat{F}_{t+h-s}, \quad h = 1, 2, \ldots, s \]
Onde:
- \(\bar{Z}_t\) é a estimativa do nível da série no tempo \(t\);
- \(\hat{T}_t\) é a estimativa da tendência;
- \(\hat{F}_{t+h-s}\) é o fator sazonal estimado para o horizonte \(h\);
- \(s\) é o período da sazonalidade.
De forma geral, o modelo de Holt-Winters é eficaz na geração de previsões de curto prazo e permite atualizações contínuas à medida que novas observações se tornam disponíveis. Por esse motivo, é amplamente aplicado em contextos operacionais e de apoio à tomada de decisão, como o monitoramento de internações hospitalares ao longo do tempo.
Critério de Informação Bayesiano (BIC)
O Critério de Informação Bayesiano (BIC) é uma métrica amplamente utilizada para seleção de modelos estatísticos, especialmente em contextos onde diferentes estruturas de modelagem são comparadas. Ele combina a qualidade do ajuste, expressa pela verossimilhança, com uma penalização pela complexidade do modelo, de forma a evitar sobreajuste.
Sua formulação geral é dada por:
\[ \text{BIC} = -2 \cdot \ln(\hat{L}) + k \cdot \ln(N) \]
Onde:
- \(\hat{L}\) é o valor máximo da função de verossimilhança do modelo ajustado;
- \(k\) é o número total de parâmetros estimados;
- \(N\) é o número de observações da amostra.
O BIC favorece modelos que oferecem bom ajuste aos dados com o menor número possível de parâmetros, sendo preferido aquele que apresentar o menor valor desse critério. Por penalizar mais fortemente a complexidade o BIC tende a selecionar modelos mais parcimoniosos, o que é desejável em aplicações com risco de superparametrização.
Base de Dados
Os dados utilizados neste estudo referem-se ao número de internações hospitalares por doenças do aparelho respiratório no estado do Rio de Janeiro. A série contempla informações agregadas por mês, abrangendo o período de janeiro de 2010 a dezembro de 2019, totalizando 120 observações.
As informações foram obtidas por meio do sistema TABNET, disponível no portal do Departamento de Informática do Sistema Único de Saúde (DATASUS), que disponibiliza estatísticas oficiais de saúde pública no Brasil. A seleção da variável considerou todos os registros classificados sob o grupo de causas relacionadas ao aparelho respiratório, conforme os critérios da Classificação Internacional de Doenças (CID-10).
Optou-se por restringir a série até dezembro de 2019, de forma a evitar possíveis distorções associadas ao início da pandemia de COVID-19 em 2020. Esse evento impactou significativamente os padrões de morbidade e hospitalização, além de introduzir alterações operacionais no sistema de saúde, o que poderia comprometer a análise de regularidades históricas da série.
Dessa forma, a base adotada busca refletir um comportamento mais estável e representativo das internações em períodos anteriores à pandemia, permitindo a aplicação de métodos clássicos de modelagem e previsão de séries temporais sem interferências atípicas.
Resultados
Análise Descritiva
Abaixo pode ser visto o comportamento do número de internações de pessoas com doenças do aparelho respiratório ao decorrer do tempo.
A Figura 1 apresenta a série temporal do número mensal de internações por doenças do aparelho respiratório no estado do Rio de Janeiro, no período de janeiro de 2010 a dezembro de 2019. Observa-se um comportamento sazonal bastante evidente ao longo dos anos, com oscilações regulares que se repetem em períodos semelhantes. Além disso, nota-se uma possível tendência decrescente no nível da série, especialmente a partir de 2011, acompanhada de uma redução gradual na amplitude das variações, o que sugere que a série pode não apresentar variância constante ao longo do tempo. Embora tal comportamento tenha sido observado, o valor-p inferior a 0,01 obtido no Teste de Dickey-Fuller Aumentado sugere que a hipótese nula de não estacionaridade pode ser rejeitada ao nível de 5%, indicando que a série é estacionária.
A Figura 2 reforça a presença de sazonalidade, ao exibir de forma desagregada o perfil mensal das internações para cada ano. É possível verificar que os meses de abril a junho concentram, de forma recorrente, os maiores números de internações em comparação aos demais meses do ano. Esse padrão se mantém ao longo dos anos analisados, o que confirma a existência de um componente sazonal consistente na série. Tais características indicam que métodos de modelagem que considerem sazonalidade são apropriados para a análise dessa série temporal.
Modelagem SARIMA
Determinação das Ordens do SARIMA
A modelagem de um modelo SARIMA tem início com a identificação de suas ordens \((p,d,q) \times (P,D,Q)_{[s]}\). Como a série apresenta sazonalidade, a determinação dos parâmetros sazonais deve ser feita a partir da análise da série dessazonalizada. Já os componentes não sazonais \((p,d,q)\) podem ser definidos com base na série original, que, segundo o Teste de Dickey-Fuller Aumentado, pode ser considerada estacionária, dispensando, inicialmente, a necessidade de uma diferenciação não sazonal \((d = 0)\).
Pela Figura 3, observa-se que a aplicação de uma diferenciação sazonal com período de 12 meses foi eficaz na atenuação do componente sazonal da série. Assim, é razoável supor que apenas uma diferenciação sazonal seja necessária, ou seja, que a ordem sazonal de diferenciação seja dada por \(D = 1\) para fins de modelagem. O Teste de Dickey-Fuller resultou em um valor-p de 0,028, o que permite rejeitar a hipótese nula de não estacionaridade ao nível de significância de 5%. Dessa forma, pode-se considerar que a série diferenciada sazonalmente (D = 1) é estacionária em nível.
A seguir, foram analisadas as funções de autocorrelação (FAC) e autocorrelação parcial (FACP) em diferentes versões da série temporal original, diferenciada sazonalmente, com o objetivo de identificar as ordens apropriadas para um modelo SARIMA.
Na análise das autocorrelações da série original, a FAC apresenta autocorrelações significativas persistentes, inclusive em múltiplos de 12, o que indica a presença de sazonalidade anual. Já a FACP mostra cortes visíveis nas defasagens 1 e 2, sugerindo a presença de componentes autorregressivos de ordem, como AR(1) ou AR(2).
Após a aplicação de uma diferença sazonal com defasagem 12, observa-se uma redução das autocorrelações sazonais na FAC, o que indica que a sazonalidade foi atenuada. As autocorrelações remanescentes sugerem a possível presença de componentes de média móvel de ordem, como MA(1) ou MA(2).
A FACP da série diferenciada sazonalmente apresenta um corte acentuado no lag 1, o que indica a presença de um componente autorregressivo sazonal de ordem 1, ou seja, um termo SAR(1).
Na FAC da série diferenciada sazonalmente, observa-se que, a partir do lag 5, há um aumento considerável nas autocorrelações, sugerindo a presença de resquícios de sazonalidade ou de alguma estrutura não capturada pela diferenciação sazonal. Diante disso, optou-se por verificar a aplicação de uma diferenciação não sazonal, ou seja, considerar \(( d = 1 )\), a fim de melhorar verificar esse comportamento.
Assim, a série resultante após a aplicação da diferenciação não sazonal e sazonal pode ser visualizada a seguir:
Pelo gráfico acima, observa-se que ambas as séries apresentam um comportamento compatível com estacionaridade. Essa evidência é reforçada pelos resultados do Teste de Dickey-Fuller, que indicaram valores-p inferiores a 0,01 para ambas as séries, permitindo rejeitar a hipótese nula de não estacionaridade ao nível de 5%, podendo assim considerar ambas séries como estacionárias.
Na análise da FAC e da FACP após a aplicação de uma diferenciação não sazonal (\(d = 1\)), observa-se que as autocorrelações permanecem significativas em diversos lags, com destaque para um crescimento a partir do lag 5, o que pode indicar a presença de componentes sazonais ainda não totalmente eliminados. Além disso, a FACP apresenta cortes nos primeiros lags, sugerindo a presença de componentes autorregressivos de baixa ordem, como \(AR(1)\).
Ao aplicar tanto a diferenciação não sazonal (\(d = 1\)) quanto a sazonal (\(D = 1\)), observa-se uma redução nas autocorrelações da FAC, principalmente nas defasagens múltiplas de 12, indicando que a sazonalidade foi adequadamente atenuada. Além disso, as autocorrelações nos primeiros lags se mantêm dentro dos limites de significância, o que sugere ausência de componentes de média móvel de ou componentes de baixa ordem, como \(MA(1)\) ou \(MA(2)\)
A FACP da série com ambas as diferenciações apresenta baixa significância nos lags iniciais, reforçando a hipótese da não existência de um componente autorregressivo sazonal ou um componente de baixar ordem, como um termo \(SAR(1)\).
Ajuste do modelo SARIMA
Com base nas ordens identificadas para os componentes sazonais e não sazonais do modelo SARIMA, passa-se agora à etapa de ajuste e comparação entre diferentes especificações. Para isso, serão testadas combinações das ordens propostas, considerando a presença de tendência e sazonalidade identificadas na análise anterior. A seleção do melhor modelo será realizada com base no Critério de Informação Bayesiano (BIC), sendo escolhido aquele que apresentar o menor valor da métrica entre os candidatos testados. Esse critério favorece modelos mais parcimoniosos, penalizando a complexidade excessiva e contribuindo para uma escolha robusta do modelo final.
Além das combinações manuais de ordens previamente definidas, será também considerado como candidato o modelo selecionado automaticamente pela função auto.arima() do pacote forecast. Esse procedimento utiliza uma busca automatizada baseada em critérios de informação como o AIC, BIC ou AICc. A inclusão do auto.arima() como um dos modelos concorrentes visa oferecer uma referência adicional na comparação, podendo eventualmente indicar configurações que não foram exploradas manualmente.
Tabela 1 – BIC dos modelos SARIMA ajustados
| Modelo (d = 1 e D = 1) | BIC | Modelo (d = 0 e D = 1) | BIC |
|---|---|---|---|
| (1,1,1) × (2,1,1)[12] | 1421.92 | (1,0,0) × (2,1,0)[12] | 1442.2 |
| (1,1,1) × (0,1,0)[12] | 1634.244 | (1,0,2) × (2,1,1)[12] | 1612.60 |
| (1,1,1) × (1,1,1)[12] | 1596.58 | (1,0,1) × (2,1,1)[12] | 1615.21 |
| (1,1,1) × (1,1,0)[12] | 1596.58 | (1,0,0) × (2,1,1)[12] | 1611.89 |
| (1,1,1) × (2,1,0)[12] | 1602.702 | (1,0,1) × (2,1,0)[12] | 1622.99 |
| (1,1,1) × (2,1,1)[12] | 1593.184 | (2,0,2) × (2,1,1)[12] | 1624.58 |
Pela Tabela 1, observa-se que o modelo com o menor valor do Critério de Informação Bayesiano (BIC), entre os ajustados manualmente, foi o \(SARIMA(1,1,1) \times (2,1,1)_{[12]}\). Assim, esse modelo foi selecionado por apresentar o melhor desempenho de ajuste segundo o critério adotado.
Dessa forma, o modelo selecionado foi estimado como \(SARIMA(1,1,1) \times (2,1,1)_{[12]}\), com base na série. O modelo estimado segue a seguinte estrutura polinomial e apresenta os seguintes parâmetros:
\[ (1 - 0{,}2539 B^{12} - 0{,}3722 B^{24})(1 - 0{,}5525 B)(1 - B)(1 - B^{12}) Y_t = (1 - 0{,}8393 B)(1 - 0{,}7443 B^{12}) a_t \]
Tabela: Estimativas dos parâmetros e seus intervalos de confiança (95%)
| Parâmetro | Estimativa | Erro Padrão | IC 95% Inferior | IC 95% Superior |
|---|---|---|---|---|
| \(ar_1\) | 0,5525 | 0,1342 | 0,2894 | 0,8156 |
| \(ma_1\) | -0,8393 | 0,0857 | -1,0074 | -0,6712 |
| \(sar_1\) | -0,2539 | 0,1544 | -0,5566 | 0,0487 |
| \(sar_2\) | -0,3722 | 0,1208 | -0,6080 | -0,1364 |
| \(sma_1\) | -0,7443 | 0,2432 | -1,2200 | -0,2686 |
Análise de resíduos do modelo SARIMA ajustado
Com a seleção do modelo ajustado \(SARIMA(1,1,1) × (2,1,1)[12]\), resta agora analisar seus resíduos com o objetivo de verificar se os pressupostos do modelo são adequadamente atendidos. A validação dos resíduos é uma etapa fundamental para garantir a confiabilidade das inferências e previsões. Os principais pressupostos que devem ser verificados são:
A análise do gráfico de linhas dos resíduos mostra que eles oscilam de forma aproximadamente simétrica em torno de zero, sem apresentar indícios de estrutura temporal ou variações abruptas na variância, o que é compatível com o comportamento esperado de resíduos de um modelo bem ajustado.
Tabela: Resultados dos testes aplicados aos resíduos
| Teste | Hipótese nula | Valor-p (ou valor de referência) |
|---|---|---|
| Teste t para média | \(\mu = 0\) (média dos resíduos é zero) | 0,2348 |
| Shapiro-Wilk | Resíduos seguem distribuição normal | 0,2452 |
| Skewness (assimetria) | Assimetria ≈ 0 (simetria em torno de zero) | - 0.0192 |
| Ljung-Box | Ausência de autocorrelação | 0,5923 |
| Dickey-Fuller Aumentado (ADF) | Série não estacionária (possui raiz unitária) | 0,01 |
Os testes aplicados aos resíduos indicam que o modelo SARIMA ajustado apresenta um bom desempenho. O teste t não rejeita a hipótese de que a média dos resíduos seja zero (\(p = 0{,}2348\)), o que sugere que não há viés sistemático. O teste de Shapiro-Wilk (\(p = 0{,}2452\)) não rejeita a normalidade dos resíduos, o que é desejável para inferência estatística. A estatística de assimetria (skewness) próxima de zero (-0,019) reforça a ideia de simetria dos resíduos.
Além disso, o teste de Ljung-Box (\(p = 0{,}5923\)) não aponta evidências de autocorrelação nos resíduos, o que indica que o modelo conseguiu capturar adequadamente a estrutura temporal da série. Por fim, o teste de Dickey-Fuller Aumentado (\(p = 0{,}01\)) confirma que os resíduos são estacionários, como esperado em um modelo bem especificado.
Em conjunto, os seguintes gráficos apoiam a adequação do modelo ajustado.
Previsão do modelo SARIMA
Após a validação dos pressupostos do modelo ajustado e a constatação de que os resíduos apresentam comportamento adequado, isto é, ausência de autocorrelação, média próxima de zero e distribuição aproximadamente normal, procede-se à etapa de previsão. Com base no modelo SARIMA selecionado, realizou-se a projeção do número de internações para os 12 meses de 2019 que não foram incluidos para achar o modelo proposto , conforme apresentado a seguir.
Tabela: Métricas de desempenho do modelo SARIMA (treinamento)
| Métrica | Valor |
|---|---|
| ME | 35,71 |
| RMSE | 311,20 |
| MAE | 231,36 |
| MPE | 0,58% |
| MAPE | 4,96% |
| MASE | 0,40 |
| ACF1 | 0,051 |
As métricas indicam que o modelo SARIMA apresenta um bom desempenho preditivo no conjunto de treinamento. O erro absoluto médio (MAE) é moderado em relação à escala da série, e o erro percentual médio absoluto (MAPE) de aproximadamente 5% reforça a boa acurácia do modelo. Além disso, o valor baixo do MASE (< 1) indica que o modelo supera a previsão ingênua (média ou “naïve”), enquanto o ACF1 residual próximo de zero (0,051) sugere ausência relevante de autocorrelação nos resíduos.
Modelagem Holt-Winter
Previsão do modelo Holt-Winter
Tabela: Métricas de desempenho do modelo Holt-Winters (conjunto de treinamento)
| Métrica | Valor |
|---|---|
| ME | 51,13 |
| RMSE | 426,09 |
| MAE | 344,14 |
| MPE | 0,78% |
| MAPE | 7,39% |
| MASE | 0,60 |
| ACF1 | 0,180 |
Comparação de Modelos
As métricas obtidas para o modelo Holt-Winters no conjunto de treinamento indicam um desempenho razoável, mas inferior ao do modelo \(SARIMA(1,1,1)(2,1,1)_{[12]}\). O erro percentual absoluto médio (MAPE) de aproximadamente 7,4% é relativamente maior, e o valor do MASE acima de 0,5 sugere que o modelo não supera com ampla folga a previsão ingênua. Além disso, o valor da autocorrelação dos resíduos no primeiro lag (ACF1 = 0,18) indica possível presença de autocorrelação residual, o que pode sinalizar limitações na captura da estrutura temporal da série.
Conclusão
A análise realizada demonstrou que a série de internações por doenças respiratórias apresenta forte padrão sazonal e comportamento compatível com estacionaridade após a aplicação de diferenciações apropriadas. Entre os modelos avaliados, o SARIMA \((1,1,1) \times (2,1,1)_{[12]}\) apresentou o melhor desempenho segundo o Critério de Informação Bayesiano (BIC) e métricas de previsão como MAE e MAPE. O modelo de Holt-Winters, embora tenha fornecido previsões razoáveis, apresentou desempenho inferior, especialmente em relação à acurácia e à autocorrelação dos resíduos.
A validação dos pressupostos do modelo SARIMA confirmou que os resíduos se comportam como ruído branco, com média próxima de zero, distribuição aproximadamente normal e ausência de autocorrelação significativa. Esses resultados indicam que o modelo ajustado é estatisticamente adequado e confiável para fins preditivos.
Em síntese, o uso de modelos de séries temporais, especialmente o SARIMA, mostra-se uma ferramenta eficaz para o monitoramento e a previsão de internações hospitalares, podendo subsidiar ações estratégicas de planejamento no sistema de saúde pública. A capacidade de antecipar padrões sazonais e comportamentos recorrentes pode contribuir significativamente para uma alocação mais eficiente de recursos hospitalares.