Módulo 3 - Prática em R

A partir desse momento vamos aprender a realizar análises de séries temporais de interesse para a vigilância utilizando a linguagem R. Caso você ainda não tenha conhecimento intermediário de R, pode fazer nosso curso introdutório. Ele está disponível online e gratuitamente em <https://sites.google.com/view/cursos-de-vigilancia/curso-inicial>. Recomendamos que tenha instalada a versão mais recentes do R em seu computador.

Seguindo no curso, vamos utilizar como exemplo a série temporal mensal de casos de SRAG no Paraná. Os dados foram obtidos do INFOGRIPE, e consolidados em um arquivo csv. Podemos ler o arquivo diretamente da web utilizando a função read.csv do R:

# Carregando o arquivo csv

srag_pr_mes <- read.csv(file = "https://raw.githubusercontent.com/ 
joaohmorais/dados-mod-temp/main/csv/srag_PR_mes.csv")



# Vendo os dados

head(srag_pr_mes)
#>    ano mes    n

#> 1 2013   1   54

#> 2 2013   2   66

#> 3 2013   3  108

#> 4 2013   4  291

#> 5 2013   5  718

#> 6 2013   6 1037

  • A função read.csv() é utilizada para carregar dentro do R um arquivo CSV que, nesse caso, está disponível no repositório github do curso. O arquivo é lido diretamente da internet e é atribuído ao objeto srag_pr_mes, que agora armazena os dados do arquivo. O arquivo contém dados sobre Síndrome Respiratória Aguda Grave (SRAG) por mês no Paraná de 2013 a 2017.

  • Já a função head() exibe as primeiras seis linhas do conjunto de dados armazenado em srag_pr_mes. Ele é útil para uma visualização inicial dos dados, permitindo verificar rapidamente a estrutura e as primeiras observações sem exibir o conjunto completo.

Os dados estão num formato tradicional de data.frame:

# Verificando a classe do objeto

class(srag_pr_mes)
#> [1] "data.frame"
  • O comando class() retorna o tipo do objeto, ou seja, ele exibe a classe ou as classes do objeto, no caso srag_pr_mes.
  • data.frame é uma estrutura de dados amplamente usada em R que organiza os dados em formato tabular, similar a uma planilha, com linhas e colunas.

No R, o formato mais comum para dados de séries temporais é o formato time series - ts. Portanto, vamos converter os dados. Utilizaremos a função ts(), que transforma os objetos em série temporal. Para esta função, devemos especificar os seguintes argumentos:

  • data: Os dados (nesse caso, o número de casos);
  • start: O início da série;
  • frequency: Sua frequência (diária = 365, mensal = 12, etc.).
# Convertendo os dados para um objeto de classe ts

srag_pr_ts <- ts(

  data = srag_pr_mes$n, # variavel que indica o número de casos

  start = c(2013, 1), # início da série - Janeiro de 2013

  frequency = 12 # frequência da série: mensal

)



class(srag_pr_ts)
#> [1] "ts"

  • ts(): Esta função transforma um objeto para a classe ts (time series, série temporal). Neste caso, o objeto ts criado com base nos dados srag_pr_mes está sendo armazenado no objeto srag_pr_ts.

    • srag_pr_mes$n: Esta é a coluna que contém o número de casos de SRAG em cada mês, que estamos usando como os dados principais da série.

    • start = c(2013, 1): Define o ponto inicial da série temporal, neste caso, janeiro de 2013, ou seja, primeiro mês de 2013. O formato c(ano, mês) indica o início da série.

    • frequency = 12: Define a frequência da série como 12 observações por ano, ou seja, uma série mensal.

Pronto! Agora temos um objeto de série temporal adequado para nossas primeiras análises no R. Antes de começar, vamos ver como é a estrutura do objeto ts?

# Vendo os dados

srag_pr_ts
#>       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec

#> 2013   54   66  108  291  718 1037 1378  664  362  248  120   67

#> 2014   80   83  128  296  477  514  455  263  178  180   86   94

#> 2015   77   77  180  255  319  373  428  274  231  101  115   99

#> 2016   69  123  590 1104 1189 1362  708  525  350  238  219  126

#> 2017   83  132  324  432  785  657  459  430  334  344  296  148

Percebam que o objeto ts srag_pr_ts tem os números de casos de SRAG por ano nas linhas e por meses nas colunas.

Agora, vamos fazer um gráfico da série temporal srag_pr_ts utilizando a função plot(). Este gráfico mostra a evolução do número de casos de SRAG ao longo do tempo, facilitando a identificação de tendências, sazonalidade ou outros padrões nos dados.

# Fazendo o gráfico da série temporal

plot(srag_pr_ts, 

     ylab = 'Casos de SRAG', # nomeando o eixo y

     xlab = 'Tempo') # nomeando o eixo x

Partindo para as análises, vamos primeiro testar a autocorrelação da série por meio da função de autocorrelação - acf():

acf(srag_pr_ts)

  • O comando acf() calcula e plota a função de autocorrelação (ACF) para um objeto ts, neste caso, a série temporal srag_pr_ts. Esse gráfico é utilizado para verificar a correlação entre os valores da série ao longo de diferentes defasagens (lags), indicando o quanto os valores da série em diferentes períodos estão correlacionados com valores anteriores.

Outra abordagem para identificação de autocorrelação é por meio do teste estatístico de Ljung-Box. Nele, tem-se que:

H0:ρh=0
H1:ρh≠0

sendo ρh a correlação da série com ela mesma (autocorrelação) defasada em h lags. No R, podemos realizar o teste de Ljung-Box com a função Box.test().

Box.test(srag_pr_ts, lag = 20, type = "Ljung-Box")
#> 

#>  Box-Ljung test

#> 

#> data:  srag_pr_ts

#> X-squared = 146.35, df = 20, p-value < 2.2e-16

O comando Box.test() realiza o teste de Ljung-Box para avaliar a independência de uma série temporal armazenada como objeto ts, no caso, de srag_pr_ts.

  • lag = 20: especifica que o teste deve incluir até 20 defasagens. Em outras palavras, ele analisa a presença de autocorrelação nas primeiras 20 defasagens da série.

  • type = "Ljung-Box": define o tipo do teste como “Ljung-Box”, teste comum para verificar autocorrelação em séries temporais.

Vê-se, tanto por meio do gráfico (onde vemos altas correlações nos primeiros lags) quanto por meio do teste (em que rejeitamos a hipótese nula de independência), que há evidências de dependência temporal na série de SRAG mensal no Paraná.


Breve análise descritiva da série

O comando summary() fornece medidas resumo sobre a série temporal. Com ele obtemos os valores mínimo, máximo, mediana e intervalos interquartílicos.

summary(srag_pr_ts)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 

#>    54.0   118.8   268.5   358.4   456.0  1378.0

O comando summary() exibe um resumo estatístico básico, neste caso da série temporal srag_pr_ts.

Vemos que a mediana de casos de SRAG por mês é 268.5, variando de 54 até 1378 casos por mês no período. Para entender como o número mensal de casos se distribui podemos fazer um histograma:

hist(srag_pr_ts, breaks = 10)

O comando hist() gera um histograma, neste caso, da série temporal srag_pr_ts.

  • breaks = 10: Especifica que o histograma terá 10 intervalos. O número de intervalos afeta a granularidade do gráfico, onde um número maior de intervalos fornece mais detalhes, enquanto um número menor de intervalos dá uma visão mais generalizada da distribuição.

Vemos que a maior concentração se dá até 200 casos por mês, mas alguns meses têm mais de 1000 casos.


Analisando seus componentes

Vimos que as séries são compostas de tendência, sazonalidade, ciclo e ruído aleatório. Vamos tentar identificar graficamente cada um desses componentes.


Tendência

Um modelo de regressão linear simples pode ser um indicador para uma possível tendência linear na série. Vamos verificar se é o caso:

plot(srag_pr_ts, main = "SRAG no Paraná - 2013-2017")

abline(reg = lm(formula = srag_pr_ts ~ time(srag_pr_ts)),

       col = "red")

  1. O comando plot() está sendo utilizado para plotar a série temporal srag_pr_ts.
  • O argumento main define o título do gráfico como “SRAG no Paraná - 2013-2017”.
  1. O comando abline() adiciona uma linha ao gráfico criado por plot() no comando anterior. Neste caso, abline() está adicionando uma linha de regressão linear obtida pela função lm():

  • lm(formula = srag_pr_ts ~ time(srag_pr_ts)) gera um modelo de regressão linear que usa o tempo time(srag_pr_ts) para prever os valores da série srag_pr_ts.

  • col = "red" define a cor da linha da regressão linear como vermelha, destacando-a no gráfico.

Vê-se que a linha tem uma leve e quase imperceptível inclinação, indicando que a série não possui uma tendência muito forte. Em caso de tendências não lineares ao longo do tempo, a função lowess() também pode ser uma boa indicadora. Vamos verificar.

# Vamos utilizar um pacote chamado Kendall para calcular a tendência

library(Kendall)



# Visualizando a tendência

plot(srag_pr_ts, main = "SRAG no Paraná - 2013-2017")

lines(lowess(time(srag_pr_ts), srag_pr_ts),

      lwd = 3,

      col = "red")

O pacote Kendall é usado para exibir a tendência na série temporal srag_pr_ts. Caso precise, instale-o através do comando install.packages("Kendall").

  • plot(srag_pr_ts, main = "SRAG no Paraná - 2013-2017"): Plota a série temporal com o título “SRAG no Paraná - 2013-2017”.

  • lines(lowess(time(srag_pr_ts), srag_pr_ts), lwd = 3, col = 2): Adiciona uma linha de suavização (lowess) sobre os dados para visualizar a tendência, com a linha de espessura 3 (lwd = 3) e na cor vermelha (col = "red").

A linha lowess facilita a observação da tendência subjacente, ajudando a destacar variações de longo prazo na série temporal.

Vemos que agora a linha de tendência se adapta conforme o tempo, indicando maior inclinação nos anos recentes da série (2016 e 2017). Contudo, ainda vemos uma tendência fraca.


Sazonalidade

Uma maneira simples de verificar a existência de sazonalidade é visualizar como a distribuição de casos se dá em cada mês. Dessa forma, vemos meses em que o número de casos tende a ser maior ou menor. Vamos, então, criar um boxplot para verificar a sazonalidade.

# Criando um boxplot para verificar a sazonalidade

boxplot(formula = srag_pr_ts ~ cycle(srag_pr_ts),

        xlab = "Mês",

        ylab = "Casos de SRAG")

A função boxplot() neste caso gera boxplots para verificar a sazonalidade na série temporal srag_pr_ts, mostrando a variação dos casos de SRAG ao longo dos meses.

  • srag_pr_ts ~ cycle(srag_pr_ts): Agrupa os dados da série temporal srag_pr_ts por ciclo mensal. A função cycle(srag_pr_ts) extrai o mês de cada observação, possibilitando uma comparação mês a mês.

  • xlab = "Mês" e ylab = "Casos de SRAG": Colocar o rótulo nos eixos X e Y como “Mês” e “Casos de SRAG”, respectivamente.

Por meio do boxplot obtido, é nítida a maior distribuição de casos no final do outono e inverno, como o período de maior a agosto. Também podemos obter uma visualização com o monthplot():

# Utilizando a função monthplot() para visualizar a sazonalidade

monthplot(srag_pr_ts)

  • monthplot(srag_pr_ts) gera um gráfico que exibe a média dos valores de srag_pr_ts para cada mês, destacando padrões sazonais ao longo dos anos.

Logo, percebemos que a sazonalidade se manifesta mais fortemente na série de SRAG do que a tendência. Podemos verificar isto por meio da função de decomposição da série.


Decomposição da série

Há um método clássico no R para decomposição da série chamado de decompose(), que realiza a decomposição da série em seus componentes utilizando médias móveis. Veja:

# Decomposição da série usando a função decompose()

plot(decompose(srag_pr_ts))

  • A função decompose() decompõe uma série temporal, no caso srag_pr_ts, em suas componentes (tendência, sazonalidade e resíduos).

  • O comando plot() plota cada um desses componentes separadamente para facilitar a análise visual.

Uma das características da aplicação do método de médias móveis é justamente a perda dos valores iniciais e finais da série (veremos mais a frente o motivo). Há outra abordagem de decomposição, chamada de STL - Seasonal and Trend Decomposition Using Loess. Como o nome diz, a decomposição é feita por meio do método Loess (que também veremos mais adiante), e realiza uma decomposição com mais robustez, mais sensível a outros tipos de sazonalidade e que lida melhor com outliers.

# Decomposição da série usando a função stl()

plot(stl(srag_pr_ts, s.window = 12))

  • A função stl() faz uma decomposição STL (Seasonal and Trend decomposition using Loess) de uma série temporal, no caso srag_pr_ts, separando suas componentes (tendência, sazonalidade e resíduos). s.window = 12 define uma janela de suavização sazonal de 12 períodos.

  • plot(): Plota o resultado da decomposição, exibindo gráficos separados para cada componente.

Como visto anteriormente, vê-se a série original (data, no topo), seu componente sazonal (seasonal, segundo gráfico), tendência (trend, terceiro gráfico) e resíduos (remainder, último gráfico). Ao lado direito, a barra de escala indica a importância para composição da série de cada termo (quanto menor, mais relevante). Vê-se que a sazonalidade e o componente aleatório são, nesse caso, os componentes mais presentes na série. Ao olhar o componente aleatório, observa-se os maiores picos em 2013 e 2016, onde o número de casos atingiu patamares maiores do que o comum.

Podemos salvar o objeto decomposto:

# Salvando o objeto decomposto e visualizando seus itens

srag_pr_decomp <- stl(srag_pr_ts, s.window = 12)



head(srag_pr_decomp$time.series)
#>       seasonal    trend  remainder

#> [1,] -285.0417 407.0146  -67.97293

#> [2,] -261.9554 409.8916  -81.93621

#> [3,] -100.8078 412.7687 -203.96086

#> [4,]  109.6366 415.6457 -234.28225

#> [5,]  328.4953 413.3284  -23.82369

#> [6,]  433.6265 411.0110  192.36244

Esse código utiliza a decomposição STL para separar a série temporal srag_pr_ts em suas componentes.

  • srag_pr_decomp <- stl(srag_pr_ts, s.window = 12): Aplica a decomposição STL (Seasonal-Trend decomposition using Loess) à série srag_pr_ts, com uma janela sazonal com 12 períodos, e armazena no objeto srag_pr_decomp.

  • head(srag_pr_decomp$time.series): Exibe as primeiras linhas da decomposição, com três componentes: seasonal (sazonalidade), trend (tendência) e remainder (resíduos), para rápida visualização dos dados de cada parte.

Vemos que a sazonalidade é a primeira coluna, tendência a segunda e o componente aleatório a terceira. Vamos extraí-las para objetos próprios:

srag_pr_tendencia <- srag_pr_decomp$time.series[,2]

srag_pr_sazonal <- srag_pr_decomp$time.series[,1]

srag_pr_aleatorio <- srag_pr_decomp$time.series[,3]

Esses comandos extraem as componentes de tendência, sazonalidade e aleatoriedade de uma decomposição de série temporal armazenada em srag_pr_decomp.

  • srag_pr_tendencia <- srag_pr_decomp$time.series[,2]: Armazena a componente de tendência da decomposição no objeto srag_pr_tendencia.

  • srag_pr_sazonal <- srag_pr_decomp$time.series[,1]: Armazena a componente sazonal (variação periódica) no objeto srag_pr_sazonal.

  • srag_pr_aleatorio <- srag_pr_decomp$time.series[,3]: Armazena a componente aleatória (ruído) no objeto srag_pr_aleatorio.

Essas variáveis permitem analisar cada componente separadamente para entender como a tendência, sazonalidade e ruído contribuem para o comportamento da série.

par(mfrow = c(3, 1))

plot(srag_pr_tendencia, main = "Componente de Tendência")

plot(srag_pr_sazonal, main = "Componente Sazonal")

plot(srag_pr_aleatorio, main = "Componente Aleatório")

Esse comando cria uma visualização conjunta dos três componentes de uma decomposição da série temporal, cada um em seu próprio gráfico:

  • par(mfrow = c(3, 1)): Divide a área de plotagem em 3 linhas e 1 coluna, exibindo três gráficos verticalmente.

  • plot(srag_pr_tendencia, main = "Componente de Tendência"): Plota o componente de tendência da série.

  • plot(srag_pr_sazonal, main = "Componente Sazonal"): Plota o componente sazonal.

  • plot(srag_pr_aleatorio, main = "Componente Aleatório"): Plota o componente aleatório (residual).

Esse layout facilita a análise visual de cada componente individual da série temporal.

E como a série é formada pela soma dos componentes, podemos obter a série original adicionando-os:

plot(

  srag_pr_tendencia + srag_pr_sazonal + srag_pr_aleatorio,

  main  = "Série recomposta",

  ylab = "Casos de SRAG"

)

Esse comando cria um gráfico da série recomposta somando os objetos com as componentes: srag_pr_tendencia, srag_pr_sazonal e srag_pr_aleatorio de srag_pr_ts.

  • srag_pr_tendencia + srag_pr_sazonal + srag_pr_aleatorio: Soma as componentes de tendência, sazonalidade e aleatoriedade da série para reconstruir a série original.

  • main = "Série recomposta": Define o título do gráfico como “Série recomposta”.

  • ylab = "Casos de SRAG": Define o rótulo do eixo Y como “Casos de SRAG”.

Chegamos ao final da nossa primeira prática com R. Agora, vamos seguir com o curso, dedicando nosso tempo a entender as técnicas de transformações e suavizações aplicadas às séries temporais, especialmente quando os dados estão muito “ruidosos”. Preparados? Vamos lá!