Imputação em séries temporais

Compartilhando o código em R

Tutorial transformação digital nº 6

Autor

Afiliação

Pataca, Luiz Carlos Moutinho (luiz.pataca@fjp.mg.gov.br)

Fundação João Pinheiro

Data de Publicação

10 de novembro de 2023

Resumo

Neste tutorial são abordados alguns conceitos sobre séries temporais e procedimentos para imputação de dados em séries temporais utilizando a linguagem e ambiente de computação estatística R. São apresentados alguns exemplos práticos com os respectivos códigos na linguagem R.

Contato: transformacao.digital@fjp.mg.gov.br

1 Introdução

Segundo o NIST/SEMATECH série temporal é uma sequência ordenada de valores de uma variável em intervalos de tempo igualmente espaçados (NIST/SEMATECH, s.d.). As séries temporais são usadas em estatística, processamento de sinais, reconhecimento de padrões, econometria, finanças, previsão do tempo, previsão de terremotos, eletroencefalografia, engenharia de controle, astronomia, engenharia de comunicações e em qualquer domínio da ciência aplicada e engenharia que envolva medições temporais (Shabou, s.d.).

O problema de imputação de dados é bem conhecido e coberto por vários pacotes disponíveis para a linguagem R. Entretanto, para o caso das séries temporais univariadas, são poucos os pacotes disponíveis capazes de preencher valores ausentes. A razão para isso reside no fato de que a maioria dos algoritmos de imputação depende de correlações entre atributos, enquanto a imputação de séries temporais univariadas precisa empregar dependências de tempo. Este fato torna a imputação em séries temporais um caso especial de imputação onde os pacotes existentes para imputação utilizados para outros dados, que não sejam séries temporais, levam a resultados insatisfatórios quando utilizados em séries temporais (Moritz, s.d.).

2 O pacote imputeTS

O pacote imputeTS foi desenvolvido especificamente para imputação em séries temporais univariadas e pode ser encontrado no CRAN. É um pacote de fácil utilização e que oferece vários algorítmos para imputação em séries temporais numéricas univariadas com dados equidistantes no tempo. Os algoritmos de imputação disponíveis incluem: ‘Média’, ‘LOCF’, ‘Interpolação’, ‘Média móvel’, ‘Decomposição sazonal’, ‘Suavização de Kalman em modelos de série temporal estrutural’, ‘Suavização de Kalman em modelos ARIMA’. Para cada um dos algoritmos existe uma função específica. Além das funções de imputação o pacote imputeTS oferece algumas funcionalidades de criação de gráficos e análise da distribuição dos valores ausentes (Moritz e Bartz-Beielstein 2017).

# Define os pacotes necessários
pacotes <-
  c("imputeTS", "gt")

# Verifica se os pacotes estão instalados e instala se necessário
install.packages(setdiff(x = pacotes,
                         y = rownames(installed.packages())))

# Carrega os pacotes
lapply(X = pacotes,
       FUN = library,
       character.only = TRUE)

3 Exemplos

3.1 Importação dos dados

Utilizaremos em nossos exemplos os dados de número de policiais militares em exercício efetivo dos municípios de Contagem, João Monlevade e Ouro Fino, do Estado de Minas Gerais, no período de 2000 a 2021. Os dados estão disponiveis no Google Drive no link https://drive.google.com/file/d/1pZl6cOr8NQ3kkIXrkpDARoYH4hFeqfzb/view?usp=drive_link. A tabela listando as funções e opções dos algoritmos de imputação do pacote imputeTS está disponível no link https://drive.google.com/file/d/1-QP49yN87htqusNyPkKwpVevZWKoig-c/view?usp=drive_link.

# Importar dados de arquivo texto com tabulação como separador de campo
dados <- read.delim2(file = "dadosSegurancaMG.tsv")

# Visualizar os dados de policiais por município
gt(head(x = dados[c("IBGE7", "ANO", "P_PM", "P_PM1")], 5))

IBGE7	ANO	P_PM	P_PM1
3100104	2000	1074,33	6
3100203	2000	721,29	31
3100302	2000	667,40	20
3100401	2000	1296,33	3
3100500	2000	574,45	20

# Importar tabela com os tipos de algoritmos de imputação
algoritmos <- read.delim2(file = "algoritmos.tsv")

# Visualizar os algoritmos e opções de imputação
gt(algoritmos)

Funcao	Opcao	Descrição
na_interpolation	linear	Imputação por interpolação linear
na_interpolation	spline	Imputação por interpolação spline
na_interpolation	stine	Imputação por interpolação de Stineman
na_kalman	StructTS	Imp. modelo estrutural e suavização de Kalman
na_kalman	auto.arima	Imp. por ARIMA e suavização de Kalman
na_locf	locf	Imputação pela última observação realizada
na_locf	nocb	Imputação pela próxima observação retrocedida
na_ma	simple	Imputação por média móvel simples
na_ma	linear	Imputação por média móvel ponderada linear
na_ma	exponential	Imp. por média móvel ponderada exponencial
na_mean	mean	Imputação pelo valor da média
na_mean	median	Imputação pelo valor da mediana
na_mean	mode	Imputação pelo valor da moda
na_random	NA	Imputação por amostra aleatória
na_seadec	NA	Imputação pela decomposição sazonal
na_seasplit	NA	Imputação pela decomposição sazonal dividida

Além dos algoritmos apresentados existe a imputação por um valor definido na_replace(), e a função na_remove() que não pode ser considerado um método de imputação e sim a remoção dos valores NA.

O pacote imputeTS só executa imputação sobre valores ausentes, valores NA, não executando operações sobre valores suspeitos. No caso de existência de valores suspeitos, antes, estes devem ser substituídos por NA na base de dados.

3.2 Criando a série temporal

Para transformar os dados univariados em uma série temporal univariada com intervalo de tempo de um ano entre as medidas utilizaremos a função ts() do pacote stats que é um dos pacotes básicos que são distribuidos com a versão R do CRAN.

# Define os códigos IBGE7 para as cidades de Contagem, João Monlevade e Ouro Fino.
ibge7 <- c(Contagem = 3118601L, JoaoMonlevade = 3136207L, OuroFino = 3146008L)
cidades <- c("Contagem", "João Monlevade", "Ouro Fino")
indice <- dados$IBGE7 %in% ibge7
dados <- dados[indice, c("IBGE7", "ANO", "P_PM1", "P_PM")]

# Selecionando os dados da cidade de Contagem
j <- 1
dadosCidade <- dados[which(dados$IBGE7 == ibge7[j]), ]

# Cria a série temporal para a quantidade de policiais na cidade de Contagem
serieCidade <- ts(data = dadosCidade$P_PM1, start = 2000)

3.3 Vizualizando a série temporal

plot(serieCidade,
    xlab = "Ano",
    ylab = "Número de policiais militares",
    main = paste0("Policiais militares na cidade de ", cidades[j]),
    type = "l"
)

Vamos remover algumas amostras para criar uma série temporal com dados ausentes.

# Definindo as amostras que serão removidas
anomalas <- NULL
anomalas$Index <- c(6, 15, 22)

# Criando a série temporal com os dados ausentes
serieCidadeNA <- serieCidade
serieCidadeNA[anomalas$Index] <- NA

O pacote imputeTS possui a função statsNA() que apresenta um sumário da distribuição dos dados ausentes.

# Apresenta a distribuição dos dados ausentes
statsNA(serieCidadeNA)

[1] "Length of time series:"
[1] 22
[1] "-------------------------"
[1] "Number of Missing Values:"
[1] 3
[1] "-------------------------"
[1] "Percentage of Missing Values:"
[1] "13,6%"
[1] "-------------------------"
[1] "Number of Gaps:"
[1] 3
[1] "-------------------------"
[1] "Average Gap Size:"
[1] 1
[1] "-------------------------"
[1] "Stats for Bins"
[1] "  Bin 1 (6 values from 1 to 6) :      1 NAs (16,7%)"
[1] "  Bin 2 (6 values from 7 to 12) :      0 NAs (0%)"
[1] "  Bin 3 (6 values from 13 to 18) :      1 NAs (16,7%)"
[1] "  Bin 4 (4 values from 19 to 22) :      1 NAs (25%)"
[1] "-------------------------"
[1] "Longest NA gap (series of consecutive NAs)"
[1] "1 in a row"
[1] "-------------------------"
[1] "Most frequent gap size (series of consecutive NA series)"
[1] "1 NA in a row (occurring 3 times)"
[1] "-------------------------"
[1] "Gap size accounting for most NAs"
[1] "1 NA in a row (occurring 3 times, making up for overall 3 NAs)"
[1] "-------------------------"
[1] "Overview NA series"
[1] "  1 NA in a row: 3 times"

Além da função statsNA() existem duas funções de construção de gráficos para análise dos dados ausentes: ggplot_na_distribution() que constroi um gráfico de linhas da distribuição dos valores ausentes e ggplot_na_distribution2() que constroi um gráfico de barras por intervalo dos valores ausentes.

3.4 Imputando os valores ausentes na série temporal

Utilizaremos no nosso exemplo o algoritmo de imputação por Suavização de Kalman em modelos de série temporal estrutural. Para tanto executaremos a função na_kalman() com a opção StructTS.

# Executa a imputação utilizando a função na_kalman com a opção StructTS
i <- 4
serieCidadeImp <- lapply(serieCidadeNA, FUN = algoritmos[i, 1], algoritmos[i, 2])

# Cria a sequência de anos para o eixo X
anos <-
    seq(from = as.Date(paste(tsp(serieCidade)[1], "-1-1", sep = "")),
        to = as.Date(paste(tsp(serieCidade)[2], "-1-1", sep = "")),
        by = "years")
# Cria a tabela dos dados imputados dos valores imputados
valoresImputados <- data.frame(NULL)
vetorOriginal <- round(serieCidade[anomalas$Index], digits = 0)
vetorImputado <- round(serieCidadeImp[[1]][anomalas$Index],
                       digits = 0)
vetorValores  <- c(rbind(vetorOriginal, vetorImputado))
names(vetorValores) <- 
    paste0("Y", rep(x = format(anos[anomalas$Index], "%Y"), each = 2), c("Verdadeiro", "Imp"))
valoresImputados <-
    rbind(valoresImputados,
          c(Municipio = cidades[j],as.vector(algoritmos[i, c(1, 2)]),
            vetorValores))

# Visualiza os vaores imputados
gt(valoresImputados)

Municipio	Funcao	Opcao	Y2005Verdadeiro	Y2005Imp	Y2014Verdadeiro	Y2014Imp	Y2021Verdadeiro	Y2021Imp
Contagem	na_kalman	StructTS	1709	1708	5	1496	1494	1282

3.5 Visualizando o gráfico da série temporal com os dados imputados

Os gráficos da série temporal são construídos com a função ggplot_na_imputations() . Esta função utiliza o pacote ggplot2 para gerar os gráficos. Os principais argumentos da função são o x_with_na que corresponde a série temporal com os valores ausentes e o argumento x_with_imputations que corresponde a série temporal com os valores imputados. Também pode ser apresentado como argumento à função, usando o argumento x_with_true, a série temporal com os dados verdadeiros, que no nosso caso está disponível e corresponde aos dados originais.

# Cria a sequência de anos para o eixo X
g <- ggplot_na_imputations(
            x_with_na = serieCidadeNA,
            x_with_imputations = serieCidadeImp[[1]],
            x_with_truth = serieCidade,
            title = "Valores imputados",
            subtitle = paste(
                'Função de imputação "',
                algoritmos[i, 1],
                '", opção "',
                algoritmos[i, 2],
                sep = ''
            ),
            xlab = "Ano",
            ylab = "Número de policiais militares",
            label_imputations = "Valores imputados",
            label_known = "Valores conhecidos",
            label_truth = "Valores reais",
            x_axis_labels = anos
        )
print(g)

Podemos verificar que o valor imputado para o ano de 2014 ficou longe do valor original. Mas é pertinente observar aqui que o valor original desvia muito dos outros valores da série de dados.

3.6 Executando para todos os três municípios

O código abaixo executa a imputação para os três municípios (Contagem, João Monlevade e Ouro Fino) utilizando todos os algoritmos e opções disponíveis no pacote imputeTS

valoresImputados <- data.frame(NULL)
anomalas <- NULL
anomalas$Index <- c(6, 15, 22)
for (j in 1:3) {
    dadosCidade <- dados[which(dados$IBGE7 == ibge7[j]),]
    serieCidade <- ts(data = dadosCidade$P_PM1, start = 2000)
    serieCidadeNA <- serieCidade
    serieCidadeNA[anomalas$Index] <- NA
    for (i in 1:nrow(algoritmos)) {
        if (!is.na(algoritmos[i, 2])) {
            if (algoritmos[i, 1] == "na_ma") {
                serieCidadeImp <-
                    lapply(serieCidadeNA,
                           FUN = algoritmos[i, 1],
                           weighting = algoritmos[i, 2])
            } else {
                serieCidadeImp <-
                    lapply(serieCidadeNA, FUN = algoritmos[i, 1], algoritmos[i, 2])
            }
        } else {
            serieCidadeImp <- lapply(serieCidadeNA, FUN = algoritmos[i, 1])
        }
        vetorOriginal <- round(serieCidade[anomalas$Index], digits = 0)
        vetorImputado <- round(serieCidadeImp[[1]][anomalas$Index],
                               digits = 0)
        vetorValores  <- c(rbind(vetorOriginal, vetorImputado))
        names(vetorValores) <-
            paste0("Y", rep(x = format(anos[anomalas$Index], "%Y"), each = 2), c("Verdadeiro", "Imp"))
        valoresImputados <-
            rbind(valoresImputados,
                  c(Municipio = cidades[j], as.vector(algoritmos[i, c(1, 2)]),
                    vetorValores))
    }
}
gt(valoresImputados)

Municipio	Funcao	Opcao	Y2005Verdadeiro	Y2005Imp	Y2014Verdadeiro	Y2014Imp	Y2021Verdadeiro	Y2021Imp
Contagem	na_interpolation	linear	1709	1709	5	1465	1494	1492
Contagem	na_interpolation	spline	1709	1710	5	1481	1494	1492
Contagem	na_interpolation	stine	1709	1709	5	1465	1494	1492
Contagem	na_kalman	StructTS	1709	1708	5	1496	1494	1282
Contagem	na_kalman	auto.arima	1709	1716	5	1501	1494	1241
Contagem	na_locf	locf	1709	1709	5	1502	1494	1492
Contagem	na_locf	nocb	1709	1709	5	1428	1494	1492
Contagem	na_ma	simple	1709	1710	5	1495	1494	1358
Contagem	na_ma	linear	1709	1710	5	1487	1494	1365
Contagem	na_ma	exponential	1709	1710	5	1478	1494	1372
Contagem	na_mean	mean	1709	1575	5	1575	1494	1575
Contagem	na_mean	median	1709	1622	5	1622	1494	1622
Contagem	na_mean	mode	1709	1709	5	1709	1494	1709
Contagem	na_random	NA	1709	1721	5	1088	1494	1273
Contagem	na_seadec	NA	1709	1709	5	1465	1494	1492
Contagem	na_seasplit	NA	1709	1709	5	1465	1494	1492
João Monlevade	na_interpolation	linear	132	132	122	116	120	119
João Monlevade	na_interpolation	spline	132	131	122	116	120	119
João Monlevade	na_interpolation	stine	132	132	122	116	120	119
João Monlevade	na_kalman	StructTS	132	132	122	116	120	118
João Monlevade	na_kalman	auto.arima	132	132	122	116	120	119
João Monlevade	na_locf	locf	132	132	122	118	120	119
João Monlevade	na_locf	nocb	132	132	122	115	120	119
João Monlevade	na_ma	simple	132	128	122	122	120	116
João Monlevade	na_ma	linear	132	129	122	120	120	117
João Monlevade	na_ma	exponential	132	130	122	119	120	118
João Monlevade	na_mean	mean	132	125	122	125	120	125
João Monlevade	na_mean	median	132	130	122	130	120	130
João Monlevade	na_mean	mode	132	132	122	132	120	132
João Monlevade	na_random	NA	132	117	122	118	120	133
João Monlevade	na_seadec	NA	132	132	122	116	120	119
João Monlevade	na_seasplit	NA	132	132	122	116	120	119
Ouro Fino	na_interpolation	linear	38	38	38	36	46	46
Ouro Fino	na_interpolation	spline	38	38	38	37	46	46
Ouro Fino	na_interpolation	stine	38	38	38	36	46	46
Ouro Fino	na_kalman	StructTS	38	38	38	38	46	37
Ouro Fino	na_kalman	auto.arima	38	-194	38	-204	46	-93
Ouro Fino	na_locf	locf	38	38	38	38	46	46
Ouro Fino	na_locf	nocb	38	38	38	33	46	46
Ouro Fino	na_ma	simple	38	37	38	36	46	37
Ouro Fino	na_ma	linear	38	38	38	37	46	38
Ouro Fino	na_ma	exponential	38	38	38	37	46	41
Ouro Fino	na_mean	mean	38	38	38	38	46	38
Ouro Fino	na_mean	median	38	38	38	38	46	38
Ouro Fino	na_mean	mode	38	38	38	38	46	38
Ouro Fino	na_random	NA	38	3	38	20	46	50
Ouro Fino	na_seadec	NA	38	38	38	36	46	46
Ouro Fino	na_seasplit	NA	38	38	38	36	46	46

Referências

Moritz, Steffen. s.d. «Comparison of Different Methods for Univariate Time Series Imputation in R».

Moritz, Steffen, e Thomas Bartz-Beielstein. 2017. «imputeTS: Time Series Missing Value Imputation in R». The R Journal 9 (1): 207. https://doi.org/10.32614/RJ-2017-009.

NIST/SEMATECH. s.d. «6.4.1. Definitions, Applications and Techniques». https://www.itl.nist.gov/div898/handbook/pmc/section4/pmc41.htm.

Shabou, Saif. s.d. Time Series with R. https://s-ai-f.github.io/Time-Series/.

Citação

BibTeX

@online{luiz carlos moutinho (luiz.pataca@fjp.mg.gov.br)2023,
  author = {Luiz Carlos Moutinho (luiz.pataca@fjp.mg.gov.br) , Pataca},
  title = {Tutorial transformação digital n\textsuperscript{o} 6},
  date = {2023-11-10},
  langid = {pt},
  abstract = {Neste tutorial são abordados alguns conceitos sobre séries
    temporais e procedimentos para imputação de dados em séries
    temporais utilizando a linguagem e ambiente de computação
    estatística R. São apresentados alguns exemplos práticos com os
    respectivos códigos na linguagem R.}
}

Por favor, cite este trabalho como:

Luiz Carlos Moutinho (luiz.pataca@fjp.mg.gov.br), Pataca. 2023. “Tutorial transformação digital n^o 6.” November 10, 2023.