Aula 02
RFN
24/02/2021
Obtenção de uma série da dados meteorológica
De forma geral, o primeiro passo para a análise de uma série temporal é a obtenção dos seus dados. Os valores, aqui utilizados, foram extraídos do do site do Instituto Nacional de Meteorogia - https://portal.inmet.gov.br/. No item Banco de dados meteorológico - https://bdmep.inmet.gov.br/.
Dados meteorológicos da cidade do Rio Grande - RS - Brasil
A série temporal escolhida, são dados de extraídos do Inmet para a a cidade de Rio Grande no período de 2001 à 2017. A tabela foi salva em .csv, desta forma para a leitura do arquivo pode-se usar a função read.csv()
dados<-read.csv('dadosTempRG2.csv',sep=';',dec='.')Após a leitura dos dados, verifica-se a estrutura dos dados com a função str().
str(dados)## 'data.frame': 28136 obs. of 8 variables:
## $ Data : Factor w/ 7034 levels "01/01/2002","01/01/2003",..: 3427 3427 3427 3427 3659 3659 3659 3659 3890 3890 ...
## $ Hora : int 0 600 1200 1800 0 600 1200 1800 0 600 ...
## $ TempArBulbo: num NA NA NA NA NA NA NA NA NA NA ...
## $ UmidadeAr : int NA NA NA NA NA NA NA NA NA NA ...
## $ VentoDir : int NA NA NA NA NA NA NA NA NA NA ...
## $ VentoRajada: num NA NA NA NA NA NA NA NA NA NA ...
## $ VentoVel : num NA NA NA NA NA NA NA NA NA NA ...
## $ Precpitacao: num NA NA NA NA NA NA NA NA NA NA ...O que observa-se? É possível verificar que o número de registros e quais são as variáveis. Especificamente neste conjunto de dados, percebe-se logo no início, a presença de diversos NA’s. Tópico que será abordado nesta aula.
A variável Data está como factor que é muitas vezes utilizado para expressar uma variável categórica presente em uma base de dados. No caso da variável Data o mais adequado é transformá-la em date.
dados$Data <- as.Date(dados$Data,"%d/%m/%Y")## Warning in strptime(x, format, tz = "GMT"): unable to identify current timezone 'U':
## please set environment variable 'TZ'str(dados)## 'data.frame': 28136 obs. of 8 variables:
## $ Data : Date, format: "2001-11-15" "2001-11-15" ...
## $ Hora : int 0 600 1200 1800 0 600 1200 1800 0 600 ...
## $ TempArBulbo: num NA NA NA NA NA NA NA NA NA NA ...
## $ UmidadeAr : int NA NA NA NA NA NA NA NA NA NA ...
## $ VentoDir : int NA NA NA NA NA NA NA NA NA NA ...
## $ VentoRajada: num NA NA NA NA NA NA NA NA NA NA ...
## $ VentoVel : num NA NA NA NA NA NA NA NA NA NA ...
## $ Precpitacao: num NA NA NA NA NA NA NA NA NA NA ...A primeira análise dos dados é feita com um summary(), para ter uma ideia geral dos dados.
summary(dados)## Data Hora TempArBulbo UmidadeAr
## Min. :2001-11-15 Min. : 0 Min. :-0.3 Min. : 19.00
## 1st Qu.:2006-09-08 1st Qu.: 450 1st Qu.:15.1 1st Qu.: 68.00
## Median :2011-07-02 Median : 900 Median :19.1 Median : 79.00
## Mean :2011-07-02 Mean : 900 Mean :18.7 Mean : 76.36
## 3rd Qu.:2016-04-25 3rd Qu.:1350 3rd Qu.:22.6 3rd Qu.: 87.00
## Max. :2021-02-16 Max. :1800 Max. :37.3 Max. :100.00
## NA's :3253 NA's :4507
## VentoDir VentoRajada VentoVel Precpitacao
## Min. : 1.0 Min. : 0.000 Min. : 0.000 Min. : 0.000
## 1st Qu.: 61.0 1st Qu.: 4.700 1st Qu.: 2.200 1st Qu.: 0.000
## Median :116.0 Median : 6.900 Median : 3.400 Median : 0.000
## Mean :143.4 Mean : 7.093 Mean : 3.607 Mean : 0.157
## 3rd Qu.:229.0 3rd Qu.: 9.200 3rd Qu.: 4.800 3rd Qu.: 0.000
## Max. :360.0 Max. :26.300 Max. :14.800 Max. :85.800
## NA's :5878 NA's :5423 NA's :5409 NA's :3652Pode-se buscar por dados faltantes usando a função missmap() do pacote Amelia.
library(Amelia)## Warning: package 'Amelia' was built under R version 3.6.3## Loading required package: Rcpp## ##
## ## Amelia II: Multiple Imputation
## ## (Version 1.7.6, built: 2019-11-24)
## ## Copyright (C) 2005-2021 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ##missmap(dados)
É possível verificar que somente as variáveis Data e Hora não contém NA’s. Desta forma vamos investigar os valores da variável TempArBulbo.
dados$TempArBulbo[3500:5500] # valor inicial tentando pegar o início do período dos dados faltantes## [1] NA NA 19.1 22.2 23.5 NA 20.3 22.1 20.8 21.2 20.6 22.1 21.4 20.5
## [15] 20.6 20.2 24.2 20.6 19.4 NA 26.1 21.9 21.4 23.4 29.1 24.2 22.8 23.6
## [29] 32.5 24.5 22.4 24.1 24.9 20.1 NA 21.8 25.3 NA NA 21.0 22.1 16.9
## [43] 12.6 18.7 22.6 20.1 NA 19.3 22.3 12.7 10.3 13.8 20.5 NA 15.9 20.4
## [57] 23.5 22.1 22.0 21.3 18.8 18.8 18.0 18.6 21.7 18.3 14.1 15.7 17.1 11.2
## [71] 9.9 13.7 NA NA 7.7 NA NA NA NA NA NA NA NA NA
## [85] NA NA NA NA NA NA 16.9 15.4 18.3 11.6 9.6 NA NA NA
## [99] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [113] NA NA NA NA NA NA 18.3 NA NA NA NA NA NA NA
## [127] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [141] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [155] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [169] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [183] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [197] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [211] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [225] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [239] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [253] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [267] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [281] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [295] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [309] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [323] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [337] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [351] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [365] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [379] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [393] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [407] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [421] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [435] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [449] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [463] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [477] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [491] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [505] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [519] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [533] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [547] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [561] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [575] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [589] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [603] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [617] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [631] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [645] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [659] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [673] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [687] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [701] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [715] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [729] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [743] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [757] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [771] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [785] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [799] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [813] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [827] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [841] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [855] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [869] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [883] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [897] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [911] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [925] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [939] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [953] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [967] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [981] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [995] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1009] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1023] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1037] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1051] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1065] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1079] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1093] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1107] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1121] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1135] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1149] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1163] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1177] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1191] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1205] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1219] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1233] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1247] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1261] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1275] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1289] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1303] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1317] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1331] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1345] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1359] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1373] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1387] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1401] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1415] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1429] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1443] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1457] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1471] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1485] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1499] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1513] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1527] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1541] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1555] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1569] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1583] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1597] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1611] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1625] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [1639] NA NA NA NA NA NA NA NA NA NA 18.5 11.9 10.1 10.8
## [1653] 13.4 10.0 5.9 8.7 16.4 11.6 10.7 13.4 17.5 12.7 12.9 14.0 20.2 16.8
## [1667] NA 17.1 22.3 20.2 17.9 16.2 16.4 15.7 15.4 17.6 21.8 16.7 17.1 16.2
## [1681] NA 18.4 18.1 18.6 25.3 20.4 19.0 21.2 24.4 20.7 19.4 20.1 22.6 18.6
## [1695] 20.0 20.2 23.0 19.8 17.3 18.5 22.6 18.2 18.9 19.7 24.0 19.3 18.2 20.0
## [1709] 24.6 16.7 15.1 16.7 18.4 17.9 17.6 17.7 22.1 19.5 18.2 17.9 22.1 18.2
## [1723] 17.5 16.7 17.2 17.3 16.7 17.6 19.5 17.3 14.9 15.8 NA NA NA 17.5
## [1737] 18.0 17.6 17.9 17.0 19.4 17.5 17.3 18.0 19.4 17.6 17.2 17.9 17.8 14.0
## [1751] 12.3 12.4 13.3 11.6 10.0 9.0 13.0 8.9 6.5 7.7 12.2 7.2 5.7 7.7
## [1765] 16.3 12.0 9.8 10.4 13.5 8.0 6.5 8.7 16.2 10.8 7.7 10.1 16.8 14.0
## [1779] 13.0 13.6 16.4 12.2 13.9 15.5 17.5 18.7 16.8 17.3 20.8 15.7 15.3 15.8
## [1793] 22.1 17.3 16.6 15.3 18.1 14.4 14.5 15.0 20.2 18.9 16.8 17.7 28.3 22.0
## [1807] 22.8 17.6 15.7 13.9 13.0 13.6 13.3 10.6 7.6 8.5 10.4 7.0 5.9 6.8
## [1821] 13.3 9.7 7.9 11.6 16.2 12.7 12.7 14.3 16.7 13.9 13.3 13.4 17.5 13.9
## [1835] 13.4 12.2 19.5 14.7 12.5 12.4 21.9 13.9 10.7 12.4 18.1 14.7 15.4 15.9
## [1849] 20.0 17.4 15.0 15.9 20.8 17.9 16.0 17.5 21.4 16.8 16.2 16.1 15.3 14.0
## [1863] 11.9 11.2 12.2 9.5 7.9 8.3 10.5 8.1 5.5 6.0 11.5 6.0 3.5 5.8
## [1877] 12.9 10.8 11.9 12.7 13.5 12.3 11.7 11.3 14.7 10.6 9.9 9.5 13.9 8.2
## [1891] NA 9.0 14.2 9.6 10.5 10.3 12.7 11.3 8.4 9.6 12.1 7.3 5.8 7.0
## [1905] 13.7 10.0 9.5 7.9 17.5 9.3 5.9 11.9 17.6 15.3 15.4 17.9 28.8 22.2
## [1919] 21.4 21.6 29.9 20.3 15.4 16.4 18.6 17.7 17.6 19.1 25.7 19.6 18.5 20.4
## [1933] 28.0 21.3 19.6 20.9 29.1 22.3 19.7 22.1 28.2 23.9 14.8 11.8 10.4 10.3
## [1947] 9.6 9.2 12.4 8.1 4.3 7.0 13.3 12.7 12.4 14.1 16.0 12.7 12.2 14.9
## [1961] 17.0 11.7 9.3 12.4 16.5 12.0 8.8 10.0 16.8 11.8 10.0 10.7 15.4 13.2
## [1975] 13.6 13.8 19.1 15.5 13.3 14.0 12.7 12.6 14.6 15.4 19.2 17.5 16.4 19.1
## [1989] 16.8 13.3 11.6 NA NA 13.0 7.7 13.5 18.6 15.9 15.0 17.0 18.8plot(dados$Data,dados$TempArBulbo)
Quando plotamos os dados de Temperatura x Data, encontramos uma faixa sem valores registados, neste momento temos duas alternativas:
descartar os primeiros anos da série (neste caso a série continua com vários anos de registros)
imputar os dados faltantes
A opção que será abordada é a imputação de dados multivariada
Seleção dos dados
Antes de realizar a imputação, por opção, serão selecionados para a análise os dados que se iniciam em 01Jan2002 e acabam em 31Dez2020. Para a manipulação dos dados serão utilizadas os pacotes dplyr, tidyr e tibble.
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(tibble)
dados2<-as_tibble(dados) # criando um tibble e preservando os dados originais
dadosSel <-dados2 %>% filter(Data > "2001-12-31")
dadosSel <-dadosSel %>% filter(Data < "2021-01-01")
# o objeto _dadosSel_ está recebendo os valores filtrados de _dados2_Agora a intenção é tentar visualizar o período dos dados faltantes e pode-se verificar que metade dos dados faltantes estão até meados de 2005.
dataTempNA<-dadosSel %>% filter(is.na(TempArBulbo)) %>% select(Data)
summary(dataTempNA)## Data
## Min. :2002-01-31
## 1st Qu.:2004-08-29
## Median :2005-03-18
## Mean :2007-12-26
## 3rd Qu.:2009-10-20
## Max. :2017-09-25Imputação de dados
Como primeiro passo será verificar o percentual de dados faltantes dos dados gerais e depois da variável temp.
#Percentual total de na's
sum(is.na(dadosSel))/(dim(dadosSel)[1]*7)*100 #7 é o número de variáveis de interesse, deixando de lado Data e Hora)## [1] 14.3732#Percentual de na's na Temperatura
sum(is.na(dadosSel$TempArBulbo))/(dim(dadosSel)[1])*100## [1] 11.60303Imputando dados usando o pacote mice
Imputa dados na temperatura
Aqui será usada uma variável completa, Hora, e uma variável com dados faltantes para ser imputada, TempArBulbo usando o pacote mice.
library(mice)## Warning: package 'mice' was built under R version 3.6.3##
## Attaching package: 'mice'## The following object is masked from 'package:stats':
##
## filter## The following objects are masked from 'package:base':
##
## cbind, rbind#cria variáveis auxiliares
imp<-dadosSel %>% select(Hora,TempArBulbo)
set.seed(279)
imputed = mice(imp, m=5)##
## iter imp variable
## 1 1 TempArBulbo
## 1 2 TempArBulbo
## 1 3 TempArBulbo
## 1 4 TempArBulbo
## 1 5 TempArBulbo
## 2 1 TempArBulbo
## 2 2 TempArBulbo
## 2 3 TempArBulbo
## 2 4 TempArBulbo
## 2 5 TempArBulbo
## 3 1 TempArBulbo
## 3 2 TempArBulbo
## 3 3 TempArBulbo
## 3 4 TempArBulbo
## 3 5 TempArBulbo
## 4 1 TempArBulbo
## 4 2 TempArBulbo
## 4 3 TempArBulbo
## 4 4 TempArBulbo
## 4 5 TempArBulbo
## 5 1 TempArBulbo
## 5 2 TempArBulbo
## 5 3 TempArBulbo
## 5 4 TempArBulbo
## 5 5 TempArBulboimp1 <- complete(imputed,1)
imp2 <- complete(imputed,2)
imp3 <- complete(imputed,3)
imp4 <- complete(imputed,4)
imp5 <- complete(imputed,5)
dadosImpMice<-dadosSel
dadosImpMice$TempArBulbo<-apply(cbind(imp1$TempArBulbo,imp2$TemArBulbo,imp3$TempArBulbo,imp4$TempArBulbo,imp5$TempArBulbo),1,mean)Comparando os dados anteriores a imputação e os dados posteriores de TempArBulbo tem-se:
#Dados Originais Selecionados
summary(dadosSel$TempArBulbo,na.rm = TRUE)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -0.30 15.00 19.00 18.64 22.60 37.30 3221#Dados Imputados com o pacote MICE
summary(dadosImpMice$TempArBulbo)#na.rm = TRUE nao precia ser utilizado pois os dados não contem mais NA's## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.30 15.20 18.80 18.57 22.20 37.30boxplot(dadosSel$TempArBulbo,dadosImpMice$TempArBulbo)
É possível verificar que onde a janela com NA’s era muito grande os valores da temperatura se aproximam mais da médias dos valores da temperatura.
plot.ts(dadosImpMice$TempArBulbo)
Imputando dados usando o pacote AMELIA
O pacote Amelia é uma outra alternativa para imputação múltipla. Ressalta-se que aqui só imputamos uma variável usando o pacote mice, o qual poderia ser usado para imputar todas as demais variáveis, para as demais variáveis faremos a imputação com o pacote Amelia tendo duas situações:
imputar a partir dos dadosSel
imputar a partir dos dadosImpMice
Imputar a partir de dadosSel usando o pacote Amelia
O objeto dadosSel contém todos os dados selecionados do período de 01/01/2002 até 31/12/2020. O primeiro cuidado é que todas as variáveis sejam numéricas.
glimpse(dadosSel)## Observations: 27,760
## Variables: 8
## $ Data <date> 2002-01-01, 2002-01-01, 2002-01-01, 2002-01-01, 2002-0...
## $ Hora <int> 0, 600, 1200, 1800, 0, 600, 1200, 1800, 0, 600, 1200, 1...
## $ TempArBulbo <dbl> 23.4, 22.5, 23.6, 29.4, 23.2, 20.7, 24.2, 26.2, 21.7, 2...
## $ UmidadeAr <int> 78, 89, 83, 50, 80, 92, 69, 64, 80, 92, 66, 62, 80, 93,...
## $ VentoDir <int> 62, 59, 317, 265, 283, 216, 220, 181, 183, 233, 211, 17...
## $ VentoRajada <dbl> 10.0, 4.2, 4.6, 12.7, 4.3, 4.9, 9.3, 9.9, 5.8, 5.0, 7.6...
## $ VentoVel <dbl> 5.5, 2.7, 3.0, 3.6, 2.1, 3.6, 6.3, 6.5, 2.2, 2.4, 4.4, ...
## $ Precpitacao <dbl> 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ...library(Amelia) #Carregar o pacote Amelia
dadosSel<-as.data.frame(dadosSel) #data.frame para o Amelia
#Determina a coluna e os limites inferiores e superiores para as imputações c(coluna,LimInf,LimSu)
lim.Umidade<-c(4,0,100)
lim.VentoDir<-c(5,1,360)
lim.VentoRajada<-c(6,0,27)
lim.VentoVel<-c(7,0,15)
lim.Precpitacao<-c(8,0,86)
limites<-rbind(lim.Umidade,lim.VentoDir,lim.VentoRajada,lim.VentoVel,lim.Precpitacao)
imputed <- amelia(dadosSel,bounds=limites) #imputando os dados## -- Imputation 1 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 2 --
##
## 1 2 3 4 5 6
##
## -- Imputation 3 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 4 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 5 --
##
## 1 2 3 4 5 6Entendendo o que aconteceu:
summary(imputed)##
## Amelia output with 5 imputed datasets.
## Return code: 1
## Message: Normal EM convergence.
##
## Chain Lengths:
## --------------
## Imputation 1: 7
## Imputation 2: 6
## Imputation 3: 7
## Imputation 4: 7
## Imputation 5: 6
##
## Rows after Listwise Deletion: 20371
## Rows after Imputation: 27760
## Patterns of missingness in the data: 13
##
## Fraction Missing for original variables:
## -----------------------------------------
##
## Fraction Missing
## Data 0.0000000
## Hora 0.0000000
## TempArBulbo 0.1160303
## UmidadeAr 0.1612032
## VentoDir 0.2105908
## VentoRajada 0.1942003
## VentoVel 0.1936960
## Precpitacao 0.1304035A maneira mais básica de fazer a imputação é escolher uma das cadeia e usá-las com os dados
dadosImpAmelia<-imputed$imputations$imp5
summary(dadosImpAmelia)## Data Hora TempArBulbo UmidadeAr
## Min. :2002-01-01 Min. : 0 Min. :-2.69 Min. : 19.00
## 1st Qu.:2006-10-01 1st Qu.: 450 1st Qu.:15.00 1st Qu.: 68.00
## Median :2011-07-02 Median : 900 Median :18.90 Median : 78.89
## Mean :2011-07-02 Mean : 900 Mean :18.62 Mean : 76.48
## 3rd Qu.:2016-04-01 3rd Qu.:1350 3rd Qu.:22.50 3rd Qu.: 87.00
## Max. :2020-12-31 Max. :1800 Max. :37.30 Max. :100.00
## VentoDir VentoRajada VentoVel Precpitacao
## Min. : 1.0 Min. : 0.000 Min. : 0.000 Min. : 0.0000
## 1st Qu.: 66.0 1st Qu.: 4.800 1st Qu.: 2.300 1st Qu.: 0.0000
## Median :130.0 Median : 7.000 Median : 3.500 Median : 0.0000
## Mean :147.1 Mean : 7.178 Mean : 3.701 Mean : 0.2581
## 3rd Qu.:228.0 3rd Qu.: 9.279 3rd Qu.: 4.923 3rd Qu.: 0.0000
## Max. :360.0 Max. :26.300 Max. :14.800 Max. :85.8000Imputar a partir de dadosImpMice usando o pacote Amelia
O objeto dadosImpMice será o novo ponto de partida, verificar se todas as variáveis sejam numéricas. Neste objeto a temperatura já foi imputada com o pacote mice
glimpse(dadosImpMice)## Observations: 27,760
## Variables: 8
## $ Data <date> 2002-01-01, 2002-01-01, 2002-01-01, 2002-01-01, 2002-0...
## $ Hora <int> 0, 600, 1200, 1800, 0, 600, 1200, 1800, 0, 600, 1200, 1...
## $ TempArBulbo <dbl> 23.4, 22.5, 23.6, 29.4, 23.2, 20.7, 24.2, 26.2, 21.7, 2...
## $ UmidadeAr <int> 78, 89, 83, 50, 80, 92, 69, 64, 80, 92, 66, 62, 80, 93,...
## $ VentoDir <int> 62, 59, 317, 265, 283, 216, 220, 181, 183, 233, 211, 17...
## $ VentoRajada <dbl> 10.0, 4.2, 4.6, 12.7, 4.3, 4.9, 9.3, 9.9, 5.8, 5.0, 7.6...
## $ VentoVel <dbl> 5.5, 2.7, 3.0, 3.6, 2.1, 3.6, 6.3, 6.5, 2.2, 2.4, 4.4, ...
## $ Precpitacao <dbl> 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ...dadosSel2<-as.data.frame(dadosImpMice) #data.frame para o Amelia
#Determina a coluna e os limites inferiores e superiores para as imputações c(coluna,LimInf,LimSu)
lim.Umidade<-c(4,0,100)
lim.VentoDir<-c(5,1,360)
lim.VentoRajada<-c(6,0,27)
lim.VentoVel<-c(7,0,15)
lim.Precpitacao<-c(8,0,86)
limites<-rbind(lim.Umidade,lim.VentoDir,lim.VentoRajada,lim.VentoVel,lim.Precpitacao)
imputed <- amelia(dadosSel2,bounds=limites) #imputando os dados## -- Imputation 1 --
##
## 1 2 3 4 5 6
##
## -- Imputation 2 --
##
## 1 2 3 4 5 6
##
## -- Imputation 3 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 4 --
##
## 1 2 3 4 5 6
##
## -- Imputation 5 --
##
## 1 2 3 4 5 6Entendendo o que aconteceu:
summary(imputed)##
## Amelia output with 5 imputed datasets.
## Return code: 1
## Message: Normal EM convergence.
##
## Chain Lengths:
## --------------
## Imputation 1: 6
## Imputation 2: 6
## Imputation 3: 7
## Imputation 4: 6
## Imputation 5: 6
##
## Rows after Listwise Deletion: 20371
## Rows after Imputation: 27760
## Patterns of missingness in the data: 13
##
## Fraction Missing for original variables:
## -----------------------------------------
##
## Fraction Missing
## Data 0.0000000
## Hora 0.0000000
## TempArBulbo 0.0000000
## UmidadeAr 0.1612032
## VentoDir 0.2105908
## VentoRajada 0.1942003
## VentoVel 0.1936960
## Precpitacao 0.1304035A maneira mais básica de fazer a imputação é escolher uma das cadeia e usá-las com os dados
dadosImpAmelia2<-imputed$imputations$imp5
summary(dadosImpAmelia2)## Data Hora TempArBulbo UmidadeAr
## Min. :2002-01-01 Min. : 0 Min. :-0.30 Min. : 19.00
## 1st Qu.:2006-10-01 1st Qu.: 450 1st Qu.:15.20 1st Qu.: 68.00
## Median :2011-07-02 Median : 900 Median :18.80 Median : 79.00
## Mean :2011-07-02 Mean : 900 Mean :18.57 Mean : 76.55
## 3rd Qu.:2016-04-01 3rd Qu.:1350 3rd Qu.:22.20 3rd Qu.: 87.00
## Max. :2020-12-31 Max. :1800 Max. :37.30 Max. :100.00
## VentoDir VentoRajada VentoVel Precpitacao
## Min. : 1.00 Min. : 0.000 Min. : 0.000 Min. : 0.0000
## 1st Qu.: 66.35 1st Qu.: 4.800 1st Qu.: 2.300 1st Qu.: 0.0000
## Median :130.54 Median : 7.000 Median : 3.508 Median : 0.0000
## Mean :147.16 Mean : 7.183 Mean : 3.713 Mean : 0.2521
## 3rd Qu.:227.36 3rd Qu.: 9.300 3rd Qu.: 4.968 3rd Qu.: 0.0000
## Max. :360.00 Max. :26.300 Max. :14.800 Max. :85.8000Separar os dados
A partir de agora, para simplificar a manipulação dos dados vamos salvar os dados completos em um novo objeto x.
x.Imp<-dadosImpAmelia2 #pode ser usado o dadosaImpAmeliaO conjunto de dados aqui utilizado é referente a Temperatura da cidade do Rio grande -RS e possui 4 registros diários que serão separados em 4 séries distintas.
dados00H<-x.Imp[x.Imp$Hora==0,]
dados06H<-x.Imp[x.Imp$Hora==600,]
dados12H<-x.Imp[x.Imp$Hora==1200,]
dados18H<-x.Imp[x.Imp$Hora==1800,]Transformando os objetos em séries temporais:
dados06HST<-ts(dados06H[,2:8],frequency=365,start=c(2002,1))
dados00HST<-ts(dados00H[,2:8],frequency=365,start=c(2002,1))
dados12HST<-ts(dados12H[,2:8],frequency=365,start=c(2002,1))
dados18HST<-ts(dados18H[,2:8],frequency=365,start=c(2002,1))Um forma de visualizar as séries dos diferentes horários do dia é fazer um gráfico com a função ts.plot().
ts.plot(dados00HST[,2], gpars=list(xlab="ano", ylab="temperatura", lty=c(1:4),main="Temp. na cidade do Rio Grande - RS - 00h - 2002 à 2020"))
Ajustando o gráfico, pode-se delimitá-lo com os valores máximo e mínimo da variável TempArBulbo.
minH<-min(x.Imp$TempArBulbo,na.rm = TRUE)
maxH<-max(x.Imp$TempArBulbo,na.rm = TRUE)E logo após plotar as séries com cores diferentes. Neste exemplo foi selecinado o primeiro ano de dados.
plot(dados18HST[1:365,2],type="l",ylim=c(minH,maxH), xlab="Ano 2002", ylab="Temperatura", main="Temp. na cidade do Rio Grande - RS - 2002")
points(dados00HST[1:365,2],type="l",lty=2,col='blue',)
points(dados06HST[1:365,2],type="l",lty=2,col='red')
points(dados12HST[1:365,2],type="l",lty=2,col='green')
#Adicionando legenda no gráfico
legend("bottomleft", legend=c('00H','06H','12H','18H'),
col=c("blue", "red", "green","black"), lty=2,cex=0.5)
Decomposição de uma série temporal
Uma das ideias iniciais é decompor a série temporal em seus componentes, para isso pode-se usar a função decompose() a qual realiza a decomposição sazonal clássica por médias móveis lidando com componentes sazonais aditivos ou multiplicativos.
x.decompClassic<-decompose(dados18HST[,2],type="multiplicative")
plot(x.decompClassic)
Neste gráfico pode-se verificar os dados observados, a tendência, a sazonalidade e a componente aleatória extraídos pela função decompose().
Outra função que realiza a decomposição de uma série temporal e está presente no pacót básico do R, stats é a fução stl(), a qual realiza a decomposição via Loess smothing, que é uma suavização de peso local. Esta função decompõe uma série temporal em componentes sazonais, de tendência e irregulares usando loess e traça os componentes separadamente, sendo que se o componente cíclico (se presente nos dados) está incluído no gráfico da tendência.
x.decompLoess<-decompose(dados18HST[,2],type="multiplicative")
plot(x.decompLoess)
Observa-se que não há quase diferenças entre s funções decompose() e stl().
Autocorrelação e autocorrelação parcial
Um outro resumo da série pode ser obtido com a função tsdisplay() do pacote forecast. (forecast::tsdisplay())
library(forecast)## Warning: package 'forecast' was built under R version 3.6.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zootsdisplay(dados18HST[,2]) 
# plot.type = c("partial", "histogram", "scatter", "spectrum"),
# ci.type = c("white", "ma") - tipo dos limites de confiança - assumir # entrada de ruído branco ou para lag k uma entrada MA (k-1)?No pacote ggplot esta função pode ser usada da seguinte forma:
library(ggplot2)
ggtsdisplay(dados18HST[,2])
#com outras opções de parâmetros
ggtsdisplay(dados18HST[,2],plot.type="scatter", theme=theme_bw())
Esta função plota uma série temporal junto com sua função de autocorrelação (acf ou FAC) e sua função de autocorrelação parcial (pacf ou facp) e um gráfico de dispersão defasado, ou ainda outras opções disponíveis na função (explore no help do R).
Mas o que são funções de autocorrelação e autocorrelação parcial?
A autocorrelação é a correlação de uma série temporal com atrasos em si mesma. Esta é uma métrica significativa porque, mostra se os estados anteriores (observações defasadas) da série temporal influenciam o estado atual. No gráfico de autocorrelação, se a autocorrelação cruzar a linha azul tracejada, significa que o atraso específico está significativamente correlacionado com a série atual.
acf18 <- acf(dados18HST[,2]) # autocorrelação
No gráfico da função acf() verifica-se que há autocorrelação significativa para todos os atrasos mostrados no eixo x.
A autocorrelação é muito usada para determinar se a série temporal é estacionária ou não. Uma série temporal estacionária fará com que a autocorrelação caia para zero rapidamente, mas para uma série não estacionária ela cairá gradualmente.
pacf18 <- pacf(dados18HST[,2]) # autorrelação parcial
A autocorrelação parcial é a correlação da série temporal com uma defasagem dela mesma, sendo removida a dependência linear de todas as defasagens entre elas.
Também é possível calcular a correlação cruzada entre duas séries temporais com a função ccf().
Extração da tendência e sazonalidade
As funções, vistas anteriormente, decompose() e forecast::stl() dividem a série temporal em suas componentes, a partir delas podemos extrair essas componentes da série temporal.
Caso se deseje estudar somente tendência de uma série pode-se trabalhar com a componente específica.
Para fazer o ajuste sazonal pode-se usar a função decompose():
x.decompClassic<-decompose(dados18HST[,2],type="multiplicative")
ajuste.sazonal<-dados18HST[,2]-x.decompClassic$seasonal
par(mfrow=c(2,1))
plot(dados18HST[1:365,2],type='l')
plot(ajuste.sazonal[5:370],type='l')
Usando a função stl() pode-se também remover a sazonalidade:
# Decompor a série
x.decompLoess <- stl(dados18HST[,2], s.window="periodic")
# pegar o componentes (sazonalidade)
ajuste.sazonal <- dados18HST[,2]-x.decompLoess$time.series[,"seasonal"]
par(mfrow=c(2,1))
plot(dados18HST[1:365,2],type='l')
plot(ajuste.sazonal[5:370],type='l')
Para remover a tendência pode-se utilizar uma regressão linear. Usando a regressão linear para modelar os dados da série temporal com índices lineares. Os resíduos do modelo resultante são uma representação da série temporal desprovida de tendência.
ajuste.tendencia <- lm(ajuste.sazonal~c(1:length(ajuste.sazonal)))
plot(resid(ajuste.tendencia), type="l") 
# resid(ajuste.tendencia) contém a série sem tendência
summary(ajuste.tendencia)##
## Call:
## lm(formula = ajuste.sazonal ~ c(1:length(ajuste.sazonal)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.162 -2.026 -0.074 1.891 14.477
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.048e+01 7.856e-02 260.7 <2e-16 ***
## c(1:length(ajuste.sazonal)) 2.079e-04 1.960e-05 10.6 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.272 on 6938 degrees of freedom
## Multiple R-squared: 0.01595, Adjusted R-squared: 0.01581
## F-statistic: 112.4 on 1 and 6938 DF, p-value: < 2.2e-16Vamos salvar os dados de trabalho em um arquivo chamado dados18HTS.csv.
write.csv2(dados18HST,'dados18HST.csv')..