Na área de planejamento e gerenciamento de recursos hídricos, séries de dados completas são necessárias para muitas variáveis, como chuva, vazão, evapotranspiração e temperatura (ELSHORBAGY, PANU e SIMONOVIC, 2000). Específicamente dados de chuva são fundamentais na hidrologia, pois trata-se da principal entrada de água no ciclo hidrológico de bacias hidrográficas. Entretanto, comumente ocorrem falhas em séries históricas.
A existência de falhas em dados hidrológicos pode ser atribuída a inúmeros fatores, como interrupção de medições por falhas em equipamentos (por ação de fenômenos naturais como furacões, descargas elétricas, etc.), perdas acidentais de dados em sistemas computacionais, dentre outros fatores de causa humana ou natural (ELSHORBAGY; PANU; SIMONOVIC, 2000).
Dados faltantes em séries históricas de precipitação (bem como outras variáveis apresentadas) se apresentam como um obstáculo para análise dos dados, bem como em estudos de modelagem hidrológica (SUN et al., 2017), por isso algumas técnicas de preenchimento (estimativa) de falhas hidrológicas foram desenvolvidas.
Dentre as técnicas de estimativa de falhas hidrológicas estão: ponderação regional, regressão linear, ponderação regional baseada em regressão linear e, mais recentemente, por meio de redes neurais artificiais (NKUNA; ODIYO, 2011), bem como por um método Bayesiano (SUN et al., 2017).
O objetivo deste tutorial é desenvolver uma rotina em linguagem R para preencher automaticamente falhas hidrológicas pelo método de regressão linear. De modo que, a técnica de estimativa de falhas será explicada ao longo da rotina, bem como a entrada e saída dos dados.
O método de preenchimento de falha por regressão linear é considerado um método simplificado onde se utiliza regressões simples ou múltiplas para gerar a informação no período com falha. Na regressão linear simples, as precipitações do posto com falhas (Y) e de um posto vizinho (X) são correlacionadas. As estimativas dos dois parâmetros da equação podem ser obtidas graficamente ou através do critério de mínimos quadrados. Para o ajuste da regressão linear simples, correlaciona-se o posto com falhas (Y) com outro vizinho (X). A correlação produz uma equação, cujos parâmetros podem ser estimados por métodos como o de mínimos quadrados, ou graficamente através da plotagem cartesiana dos pares de valores (X, Y), traçando-se a reta que melhor representa os pares de pontos. Uma vez definida a equação semelhante à apresentada abaixo, as falhas podem ser preenchidas:
\({f(x)}\)=\(a\)+\(b\).\(x\)
Onde: a se refere ao coeficiente linear e o b se refere ao coeficiente angular. Existem alguns fatores que podem limitar a utilização do método como o valor do r² ser menor que 0,7 e r inferior a 0,84.
#install.packages("readxl")
#install.packages("dplyr")
#install.packages("glue")
#install.packages("writexl")
library(readxl)
library(dplyr)
library(glue)
library(writexl)A primeira etapa é a importação da base de dados. Nesse caso a série histórica (de dados mensais ou anuais) com dados hidrológicos de postos de coleta de apoio (aqueles que servirão para estimar o dado faltante), bem como os dados da estação de coleta que apresenta falhas na série histórica. Para isso, a biblioteca readxl foi usada.
A biblioteca readxl dispõe da função read_excel que faz a leitura de planilhas em arquivo de extensão .xlsx. Dessa forma a função foi usada com objetivo de facilitar uma preparação anterior dos dados em planilhas desse formato. Além disso, como default, a função read_excel reconhece as falhas de dados (missing data ou Not Available) na base de dados, e as preenche como NA.
Assim o código chunk de entrada dos dados é escrito da seguinte forma:
precipitacao <- read_excel("preenchimento_falhas.xlsx", col_types = "numeric")
precipitacaoObserve que, o usuário manipula o valor do arquivo (file = “exemplo.xlsx”). No arquivo “preenchimento_falhas.xlsx” há somente uma planilha (sheet), mas, se houver mais de uma, é necessário informar como argumento da função read_excel qual planilha será usada (usando o argumento sheet = ). É válido notar que os nomes das colunas do data.frame usado neste tutorial são adotadas didaticamente da seguinte forma:
colnames(precipitacao)## [1] "Ano" "Falha" "apoio_1" "apoio_2" "apoio_3"Isso porquê nesse código, que automatiza o preenchimento de falhas (onde o usuário manipula apenas a entrada e saída dos dados), foi necessário padronizar as colunas do data.frame de entrada, sendo necessariamente as colunas 1 e 2, a data da informação (anual ou mensal) e a série da estação que apresenta dados faltantes a serem estimados.
Após importação, a rotina deve realizar as regressões (quantas forem possíveis de acordo com o número de estações de apoio). A análise de regressão consiste no ajuste de uma equação que demonstra o comportamento da variação de duas variáveis, uma em função de outra. Nesse caso trata-se de uma regressão linear, matematicamente expressa da seguinte forma:
\(\hat{Y}\)=\(\hat{\beta_{0}}\)+\(\hat{\beta_{1}}\).\(X\)
Ou ainda:
\({f(x)}\)=\(a\)+\(b\).\(x\)
Portanto o ajuste do modelo de regressão exige a obtenção (ou estimativa) dos parâmetros \(\hat{\beta_{0}}\) e \(\hat{\beta_{1}}\). A função usada para ajustar a regressão linear é a lm, da biblioteca nativa em R, stats. Entretanto, como dito anteriormente, o número de ajustes de regressões depende do número de séries históricas de estações de apoio existentes. Portanto, como o objetivo é automatizar o preenchimento de falhas, foi usado um loop descrito doravante.
A saída do loop, usado para ajustar i vezes o modelo de regressão linear, é um data.frame composto de 4 colunas: a primeira para identificar a estação que é correlacionada com a estação com falhas, a segunda com o coeficiente de determinação \(r^2\), as colunas 3 e 4 para \(\hat{\beta_{0}}\) e \(\hat{\beta_{1}}\), respectivamente. Arbitrariamente o data.frame com os parâmetros de ajuste da regressão, foi denominado parameters, como mostra o código chunk abaixo:
parameters <- data.frame(estacao = colnames(precipitacao[,-1:-2]),
r_sqrd = NA,
a = NA,
b = NA)
parameters## estacao r_sqrd a b
## 1 apoio_1 NA NA NA
## 2 apoio_2 NA NA NA
## 3 apoio_3 NA NA NACriado o data.frame que receberá os parâmetros, pode-se então realizar o ajuste dos modelos de regressão usando a função lm (stats), usando a função for. Assim, a variável que faz mudar a posição no vetor é a variável i, e o vetor é descrito como: 3:ncol(precipitacao); indicando que as colunas relacionadas serão da 3ª (lembre-se que, as colunas 1 e 2 são padronizadas para data da informação e série histórica com falhas, respectivamente) até quantas colunas houver no data.frame da base de dados.
O modelo de regressão linear é calculado argumentando na função lm o data.frame precipitacao como matriz (usando a função as.matrix): padronizando o Y (variável no eixo das abcissas) como a coluna 2 do data.frame, e no eixo X a cada coluna de estação de apoio (i). Para recuperar o valor de \(r^2\) foi usada a função summary.lm (que sumariza a saída da função como uma lista), especificando o \(r^2\) pelo argumento r.squared (que vai retorna um vetor com este valor). Os coeficientes \(\hat{\beta_{0}}\) e \(\hat{\beta_{1}}\), foram extraídos (na forma de vetor), usando a função coefficients. Em seguida, os vetores com \(r^2\), \(\hat{\beta_{0}}\) e \(\hat{\beta_{1}}\) são concatenados no data.frame (criado logo acima) parameters:
for (i in 3:ncol(precipitacao)) {
reglin <- lm(as.matrix(precipitacao[,2]) ~ as.matrix(precipitacao[,i]),
data = precipitacao)
r_sqrd <- summary.lm(reglin)$r.squared
a_b <- coefficients(reglin)
parameters[i-2,-1] <- c(r_sqrd, a_b)
}Observe que a saída é um data.frame parameters preenchido com os valores necessários doravante:
parameters## estacao r_sqrd a b
## 1 apoio_1 0.7230519 7.459739 1.0234817
## 2 apoio_2 0.9250707 -24.492741 1.1794832
## 3 apoio_3 0.8700884 -6.081400 0.9892785str(parameters)## 'data.frame': 3 obs. of 4 variables:
## $ estacao: Factor w/ 3 levels "apoio_1","apoio_2",..: 1 2 3
## $ r_sqrd : num 0.723 0.925 0.87
## $ a : num 7.46 -24.49 -6.08
## $ b : num 1.023 1.179 0.989Nesta etapa, serão utilizadas prioritariamente funções do pacote dplyr, desenvolvido para manipulação de dados em data frame.
Obtidos n modelos de regressão, só poderão ser usados para estimativa de falhas, aquelas séries de dados que se correlacionam com um coeficiente de correlação: \(r \ge 0.84\); ou com um coeficiente de determinação: \(r^2 \ge 0.7\). Para atender esse critério, é aplicado um filtro ao data frame parameters, gerando um novo data frame chamado filter1, apenas com linhas de estações com \(r^2 \ge 0.7\):
filter1 <- parameters %>% filter(r_sqrd>=0.7)Veja que, no caso da base de dados deste tutorial, todas as estações apresentam coeficiente de determinação \(r^2 \ge 0.7\) (comparar com o data.frame preenchido parameters):
filter1## estacao r_sqrd a b
## 1 apoio_1 0.7230519 7.459739 1.0234817
## 2 apoio_2 0.9250707 -24.492741 1.1794832
## 3 apoio_3 0.8700884 -6.081400 0.9892785str(filter1)## 'data.frame': 3 obs. of 4 variables:
## $ estacao: Factor w/ 3 levels "apoio_1","apoio_2",..: 1 2 3
## $ r_sqrd : num 0.723 0.925 0.87
## $ a : num 7.46 -24.49 -6.08
## $ b : num 1.023 1.179 0.989Agora que temos apenas aqueles modelos ajustados com \(r^2 \ge 0.7\), o ideal é usar o modelo que melhor se ajusta, ou seja, o que apresenta maior \(r^2\). Paralelo a isso, pode ocorrer que não exista nenhum modelo que atende ao critério \(r^2 \ge 0.7\). Como então arranjar as linhas em filter1, verificando a existência destas? A maneira empregada foi usar a função de condição ifelse.
A função ifelse exige nos argumentos: uma condição, uma ação para condição verdadeira, e uma ação para uma condição falsa; nesse caso, os argumentos foram colocados da seguinte forma análoga: Se o número de linhas do data frame filter1 for maior ou igual a 1 (ou seja, se houver ao menos 1 modelo que atende ao critério \(r^2 \ge 0.7\)) então use a função arrange (do pacote dplyr) e organize os modelos em ordem decrescente de \(r^2\), se não, emita uma mensagem de erro (usando a função stop) afirmando: Não há nenhum coeficiente de determinação maior do que 0,7:
ifelse(nrow(filter1)>=1, arrange1 <- arrange(filter1, desc(r_sqrd)),
stop("There isn't anyone determination coefficient larger than 0,7"))## [[1]]
## [1] apoio_2 apoio_3 apoio_1
## Levels: apoio_1 apoio_2 apoio_3Agora, se a condição for verdadeira, teremos um novo data frame denominado arrange1 no qual, os modelos ajustados que atendem ao critério \(r^2 \ge 0.7\) estão arranjados em ordem decrescente de \(r^2\), facilitando o próximo passo. Observe o arrange1 abaixo e compare com filter1 (logo acima).
arrange1str(arrange1)## 'data.frame': 3 obs. of 4 variables:
## $ estacao: Factor w/ 3 levels "apoio_1","apoio_2",..: 2 3 1
## $ r_sqrd : num 0.925 0.87 0.723
## $ a : num -24.49 -6.08 7.46
## $ b : num 1.179 0.989 1.023O procedimento de arrumar as linhas em ordem decrescente de \(r^2\) visa facilitar esta etapa: coletar os parâmetros do modelo que apresentou melhor ajuste. Esta etapa foi feita gerando um novo data frame com as mesmas variáveis de parameters, filter1 e arrange1, entretanto com apenas uma linha. Esse procedimento foi realizado usando a função slice do pacote dplyr, que seleciona uma linha pela posição:
best_fit <- arrange1 %>% slice(1)
best_fit## estacao r_sqrd a b
## 1 apoio_2 0.9250707 -24.49274 1.179483str(best_fit)## 'data.frame': 1 obs. of 4 variables:
## $ estacao: Factor w/ 3 levels "apoio_1","apoio_2",..: 2
## $ r_sqrd : num 0.925
## $ a : num -24.5
## $ b : num 1.18Feito isso, didaticamente foram inseridas duas funções com o objetivo de gerar uma mensagem indicando os parâmetros do modelo da série histórica que melhor se ajusta. Essa etapa foi realizada usando a função glue, da biblioteca análoga: glue.
A função glue é um sinônimo da função paste (nativa do R), entretanto apresenta uma sintaxe mais fluídica: variáveis são inseridos dentro do texto através de chaves {}. De modo que, a função preenche a variável info, logo descrita:
info <- glue("A estação de apoio que melhor se relaciona com a estação com falhas é {best_fit[[1,1]]}.
Esta estação apresentou um coeficiente de determinação r² = {signif(best_fit[[1,2]], digits = 4)}.
De modo que, os coeficientes linear e angular (interceptor e slope) foram, respectivamente, de: a = {signif(best_fit[[1,3]], digits = 6)} e b = {signif(best_fit[[1,4]], digits = 6)}.")
info## A estação de apoio que melhor se relaciona com a estação com falhas é apoio_2.
## Esta estação apresentou um coeficiente de determinação r² = 0.9251.
## De modo que, os coeficientes linear e angular (interceptor e slope) foram, respectivamente, de: a = -24.4927 e b = 1.17948.Ainda assim, a mensagem é exportada em um arquivo de texto (em extensão .txt), para consultas posteriores por parte do analista:
write(info, file = "info_regressao.txt")Agora que a rotina possui os parâmetros de melhor ajuste do modelo de regressão linear, é necessário gerar uma função que reconheça esses valores e que faça a estimativa (cálculo) dos valores faltantes na estação com falhas. A função foi denominada extrapolacao e foi criada usando a função function, definindo as variáveis linear (\(\hat{\beta_{0}}\)), angular (\(\hat{\beta_{1}}\)) e chuva_mm (\(X\)), seguindo a estrutura da equação:
\(\hat{Y}\)=\(\hat{\beta_{0}}\)+\(\hat{\beta_{1}}\).\(X\)
extrapolacao <- function(linear,angular,chuva_mm) {
linear+(angular*chuva_mm)
}Nesta etapa o papel da rotina é: identificar linhas com falhas na coluna 2 do data frame precipitacao. Para isso foi é necessário usar a função filter (pacote dplyr) e operadores lógicos do R: a função is.na, como mostra o código abaixo:
precit_falha <- precipitacao %>% filter(is.na(precipitacao[,2]))
precit_falha## # A tibble: 2 x 5
## Ano Falha apoio_1 apoio_2 apoio_3
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2001 NA 89.8 80 92.7
## 2 2003 NA 88.6 150. 175.De modo que possuindo um data frame apenas com as linhas com falhas na coluna 2, o procedimento é estimar o valor usando a função gerada na etapa anterior. Como temos um data frame com os parâmetros de melhor ajuste (lembre-se do best_fit), e a função extrapolacao lê esses parâmetros como variáveis, define-se que os valores da coluna 2 do data frame precit_falha, deverão ser calculados usando a fução.
A variável chuva_mm ( da função) é obtida com o valor correspondente ao NA da coluna 2 na coluna que apresentar o melhor ajuste. Por isso, é criado o vetor n, que obtém o primeiro valor do data frame best_fit (o nome da estação que melhor se ajustou) e gera um texto. Esse texto entra na função entre os operadores [[]] para indicar a coluna onde serão obtidos os valores de chuva_mm:
n <- paste(first(best_fit))
precit_falha[,2] <- extrapolacao(best_fit[1,3],best_fit[1,4],precit_falha[[n]])
precit_falhaAgora que calculamos as falhas, vamos preenchê-las no data frame original para então exportá-lo. Para isso, novamente aplicaremos o comando filter no data frame precipitacao, de modo a selecionar as linhas nas quais não há NA na coluna 2, usando o operador Booleano “!”:
precipitacao <- precipitacao %>% filter(!is.na(precipitacao[,2]))
precipitacaoAgora, temos o data frame precipitacao sem as falhas da coluna 2, vamos unir com o data frame precit_falha com as falhas da coluna 2 estimadas, usando o comando bind_rows. Feito isso, as linhas são reordenadas em ordem crescente da data:
precipitacao <- bind_rows(precipitacao, precit_falha)
precip_preench <- arrange(precipitacao, Ano)
precip_preenchA última etapa é a exportação do data frame completo, ou seja, com as falhas da coluna 2 preenchidas. Foi escolhido uma exportação também em arquivo excel (extensão .xlsx), de modo que a função capaz de executar isso foi a função write_xlsx do pacote writexl, conforme segue código abaixo:
write_xlsx(precip_preench,"falhas_preenchidas.xlsx")ELSHORBAGY, Amin A.; PANU, U. S.; SIMONOVIC, S. P.. Group-based estimation of missing hydrological data: I. Approach and general methodology. Hydrological Sciences Journal, [s.l.], v. 45, n. 6, p.849-866, dez. 2000. Informa UK Limited. http://dx.doi.org/10.1080/02626660009492388
SUN, Siao et al. A Bayesian method for missing rainfall estimation using a conceptual rainfall–runoff model. Hydrological Sciences Journal, [s.l.], v. 62, n. 15, p.2456-2468, 31 out. 2017. Informa UK Limited. http://dx.doi.org/10.1080/02626667.2017.1390317
NKUNA, T.R.; ODIYO, J.O.. Filling of missing rainfall data in Luvuvhu River Catchment using artificial neural networks. Physics And Chemistry Of The Earth, Parts A/b/c, [s.l.], v. 36, n. 14-15, p.830-835, jan. 2011. Elsevier BV. http://dx.doi.org/10.1016/j.pce.2011.07.041
citation("readxl")##
## To cite package 'readxl' in publications use:
##
## Hadley Wickham and Jennifer Bryan (2019). readxl: Read Excel
## Files. R package version 1.3.1.
## https://CRAN.R-project.org/package=readxl
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {readxl: Read Excel Files},
## author = {Hadley Wickham and Jennifer Bryan},
## year = {2019},
## note = {R package version 1.3.1},
## url = {https://CRAN.R-project.org/package=readxl},
## }citation("dplyr")##
## To cite package 'dplyr' in publications use:
##
## Hadley Wickham, Romain François, Lionel Henry and Kirill Müller
## (2019). dplyr: A Grammar of Data Manipulation. R package version
## 0.8.0.1. https://CRAN.R-project.org/package=dplyr
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {dplyr: A Grammar of Data Manipulation},
## author = {Hadley Wickham and Romain François and Lionel Henry and Kirill Müller},
## year = {2019},
## note = {R package version 0.8.0.1},
## url = {https://CRAN.R-project.org/package=dplyr},
## }citation("glue")##
## To cite package 'glue' in publications use:
##
## Jim Hester (2018). glue: Interpreted String Literals. R package
## version 1.3.0. https://CRAN.R-project.org/package=glue
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {glue: Interpreted String Literals},
## author = {Jim Hester},
## year = {2018},
## note = {R package version 1.3.0},
## url = {https://CRAN.R-project.org/package=glue},
## }citation("writexl")##
## To cite package 'writexl' in publications use:
##
## Jeroen Ooms (2018). writexl: Export Data Frames to Excel 'xlsx'
## Format. R package version 1.1.
## https://CRAN.R-project.org/package=writexl
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {writexl: Export Data Frames to Excel 'xlsx' Format},
## author = {Jeroen Ooms},
## year = {2018},
## note = {R package version 1.1},
## url = {https://CRAN.R-project.org/package=writexl},
## }