Introdução

A previsão de eventos climáticos, como a ocorrência de chuva, é uma tarefa importante em diversas áreas, como agricultura, planejamento urbano e gestão de recursos hídricos. Este relatório apresenta uma análise exploratória e preditiva da precipitação diária na cidade de Niquelândia (GO), utilizando uma série histórica de dados meteorológicos entre 2001 e 2018.

O foco inicial da análise é identificar padrões na ocorrência de chuva, tanto ao longo do tempo quanto em função de variáveis climáticas, como temperatura, umidade, pressão atmosférica e vento. Para isso, técnicas de modelagem baseadas em cadeias de Markov são aplicadas, possibilitando a construção de matrizes de transição diárias e mensais para estimar a probabilidade de continuidade ou interrupção de períodos chuvosos.

Além disso, são geradas novas variáveis derivadas (como diferenças e razões entre temperaturas e umidades), com o objetivo de enriquecer a base preditiva. Por fim, avaliamos a acurácia e o poder discriminativo dessas variáveis na previsão de chuva utilizando curvas ROC, AUC e outras métricas clássicas de classificação (sensibilidade, especificidade, acurácia, kappa, youden), destacando quais indicadores apresentam maior capacidade de previsão.

Leitura dos Dados

Leitura dos pacotes:

## Pacotes Utilizados
library(data.table)
library(dplyr)
library(markovchain)
library(pROC)
library(ggplot2)
library(cutpointr)
library(knitr)
library(kableExtra)
library(caret)
library(purrr)
# Leitura dos Dados
data <- fread("dados_A004_D_2001-05-30_2018-06-20.csv", 
               sep = ";", skip = 10, header = TRUE, dec = ',', 
               na.string = "null", drop = "V12") %>% na.omit()

# Renomear variáveis
names(data) <- c("data", "prec.tot", "pressao.med", "temp.orv", 
                  "temp.max", "temp.med", "temp.min", "ur.med", "ur.min", 
                  "vento.max", "vento.med")
# Checar os dados
str(data)
## Classes 'data.table' and 'data.frame':   4628 obs. of  11 variables:
##  $ data       : IDate, format: "2001-06-02" "2001-06-03" ...
##  $ prec.tot   : num  1.6 0 0 0 0 0 0 0 0 0 ...
##  $ pressao.med: num  940 940 941 940 941 ...
##  $ temp.orv   : num  14 15.2 13.8 14.5 13 11 12.3 13.1 13.1 11.6 ...
##  $ temp.max   : num  31.1 29.7 30.6 31.5 29.8 30.8 29.8 28.6 28.9 30 ...
##  $ temp.med   : num  24.5 24.4 23.8 24.8 24.1 23.3 23.1 22.9 22.5 22.8 ...
##  $ temp.min   : num  19.5 19.9 16.9 19.7 20 16.2 17 17.6 17.1 17.3 ...
##  $ ur.med     : num  53.5 57.1 55.3 55.1 51.6 48.5 52.1 55.1 56.8 51.9 ...
##  $ ur.min     : int  30 40 32 31 32 22 31 36 35 23 ...
##  $ vento.max  : num  11.8 13 8.4 8.1 12.5 8.5 10.4 12.6 11.2 10 ...
##  $ vento.med  : num  2.6 2.8 2.4 2 3.3 2.1 2.6 3.9 2.7 2.7 ...
##  - attr(*, ".internal.selfref")=<externalptr>

Com todas as variáveis originais no formato correto, podemos criar e categorizar algumas variáveis:

# Função para categorização da variável precipitação. cutpoint = > 2.5 ml
dados <-    data %>% mutate( dia        =   lubridate::day(data), 
                           mes          =   lubridate::month(data),
                           ano          =   lubridate::year(data),
                           tri          =   quarter(data),
                           ndias        =   lubridate::days_in_month(data), 
                           prec.H    =      ifelse(prec.tot > 2.5, 1, 0), 
                           temp.raz  =   temp.min / temp.max,
                           temp.dif  =   temp.max - temp.min,
                           ur.raz    =   ur.min / ur.med,
                           ur.dif    =   ur.med - ur.min,
                           idseq     =   cumsum(c(1, diff(data) != 1))) %>%  
  group_by(idseq) %>%
  mutate(prec.A = as.factor(lead(prec.H)), ene = n()) %>% 
  ungroup() %>%  
  dplyr::select(-data, -idseq, -dia, -ano, -tri, -ndias, -ene) %>% 
  dplyr::select(prec.H, prec.A, everything())
dados <- dados[!apply(dados, 1, function(x) any(is.na(x))), ]

str(dados)
## tibble [4,504 × 17] (S3: tbl_df/tbl/data.frame)
##  $ prec.H     : num [1:4504] 0 0 0 0 0 0 0 0 0 0 ...
##  $ prec.A     : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ prec.tot   : num [1:4504] 1.6 0 0 0 0 0 0 0 0 0 ...
##  $ pressao.med: num [1:4504] 940 941 941 942 941 ...
##  $ temp.orv   : num [1:4504] 14 13.8 13 11 12.3 13.1 13.1 11.6 11.1 11.4 ...
##  $ temp.max   : num [1:4504] 31.1 30.6 29.8 30.8 29.8 28.6 28.9 30 30.3 30.5 ...
##  $ temp.med   : num [1:4504] 24.5 23.8 24.1 23.3 23.1 22.9 22.5 22.8 22.9 23.4 ...
##  $ temp.min   : num [1:4504] 19.5 16.9 20 16.2 17 17.6 17.1 17.3 15.7 17.2 ...
##  $ ur.med     : num [1:4504] 53.5 55.3 51.6 48.5 52.1 55.1 56.8 51.9 49.8 48.5 ...
##  $ ur.min     : int [1:4504] 30 32 32 22 31 36 35 23 23 26 ...
##  $ vento.max  : num [1:4504] 11.8 8.4 12.5 8.5 10.4 12.6 11.2 10 10.1 9.5 ...
##  $ vento.med  : num [1:4504] 2.6 2.4 3.3 2.1 2.6 3.9 2.7 2.7 2.2 1.7 ...
##  $ mes        : num [1:4504] 6 6 6 6 6 6 6 6 6 6 ...
##  $ temp.raz   : num [1:4504] 0.627 0.552 0.671 0.526 0.57 ...
##  $ temp.dif   : num [1:4504] 11.6 13.7 9.8 14.6 12.8 11 11.8 12.7 14.6 13.3 ...
##  $ ur.raz     : num [1:4504] 0.561 0.579 0.62 0.454 0.595 ...
##  $ ur.dif     : num [1:4504] 23.5 23.3 19.6 26.5 21.1 19.1 21.8 28.9 26.8 22.5 ...

Aplicando a função acima, temos as seguintes variáveis na base de dados:

  • Prec.H: Se choveu hoje (>2.5 mm de chuva), então é = 1, 0 c.c..
  • Prec.A: Se choveu no dia seguinte (>2.5mm de chuva), então = 1, 0 c.c.
  • Prec.tot: Total em mm de precipitação (chuva) no dia.
  • Pressao.med: Pressão média medida no dia.
  • Temp.orv : Temperatura (ºC) do ponto de orvalho no dia.
  • Temp.max: Temperatura (ºC) máxima registrada no dia.
  • Temp.min: Temperatura (ºC) mínima registrada no dia.
  • Temp.med: Temperatura (ºC) média do dia.
  • Ur.med : Umidade relativa % do ar média no dia.
  • Ur.min: Umidade relativa % mínima medida no dia.
  • Vento.max: Velocidade máxima de vento (km/h) registrada no dia.
  • Vento.med: Velocidade média do vento (km/h) no dia.
  • Mês: Número do mês que foi feito a medição (1 = Janeiro, 2 = Fevereiro, …)
  • Temp.raz: Razão de temperatura entre a máxima e mínima.
  • Temp.dif: Range de temperatura no dia.
  • Ur.raz: Razão da umidade relativa calculada no dia.
  • Ur.dif: Range de umidade relativa no dia.
dados <- dados[!apply(dados, 1, function(x) any(is.na(x))), ]

Removemos todas as linhas que possuem NA, devido a diferença de dia.

Usando proportion table, temos uma matriz com as probabilidades de transição:

  • 0 -> 0 : Não choveu hoje, não chove amanhã
  • 0 -> 1: Não choveu hoje, chove amanhã
  • 1 -> 0: Choveu hoje, não chove amanhã
  • 1 -> 1 : Choveu hoje, chove amanhã.
table(dados$prec.H, dados$prec.A)
##    
##        0    1
##   0 3365  390
##   1  390  359
round(prop.table(table(dados$prec.H, dados$prec.A), margin = 1),3)
##    
##         0     1
##   0 0.896 0.104
##   1 0.521 0.479

Temos então que se NÃO está chovendo hoje, a probabilidade de ainda continuar sem chuva amanhã é de 89.6%. E se ESTÁ chovendo hoje, a probabilidade de continuar chovendo amanhã é de 47.9%.

Isso da um indício do clima da região, que indica certa aderência ao estado de seca/chuva. Em dias secos, é pouco provavel que tenha chuva (10,4%) e em dias chuvosos, praticamente a mesma chance existe para continuar chovendo ou não.

Para calcular a matriz de transição em n passos, usaremos o pacote MarkovChain. Iremos simular um n = 30, onde a matriz de transição está estabilizada.

seqano <- c(dados$prec.H)
mcseq <- markovchainFit(seqano)
mcseqvalores <- mcseq$estimate ; n <- 30
matriz_transicao_n <- mcseqvalores^n
matriz_transicao_n
## MLE Fit^30 
##  A  2 - dimensional discrete Markov Chain defined by the following states: 
##  0, 1 
##  The transition matrix  (by rows)  is defined as follows: 
##           0         1
## 0 0.8336664 0.1663336
## 1 0.8336664 0.1663336

Independente do dia, temos que a probabilidade de CHUVA é de 16.6% e de NÃO CHUVA é de 83.4%. Portanto, a cidade de Niquelândia em Goiás é uma cidade predominantemente não chuvosa.

Um breve pesquisa sobre a cidade nos mostra que o inverno é seco e quente, e no verão é umido e chuvoso. Iremos analisar esse padrão mês a mês:

# Gráfico de Precipitação por Mês
dados %>%
  group_by(mes) %>%
  summarise(prec_total = sum(prec.tot, na.rm = TRUE)) %>%
  ggplot(aes(x = factor(mes), y = prec_total)) +
  geom_col(fill = "darkgreen") +
  labs(
    title = "Precipitação Total por Mês",
    x = "Mês",
    y = "Precipitação Total (mm)"
  ) +
  theme_minimal()
Gráfico de Precipitação Total em mm por mês - Niquelândia / GO

Gráfico de Precipitação Total em mm por mês - Niquelândia / GO

Pelo gráfico, é possível observar que os dias mais chuvosos concentram-se especialmente no verão, enquanto os períodos mais secos ocorrem no inverno.

Uma possível extensão dessa análise é a construção de matrizes de transição mensais, permitindo estimar a probabilidade de ocorrência ou ausência de chuva em função do mês específico.

# Função p/ calcular as matrizes de transição mês a mês
matrizes_mensais <- lapply(1:12, function(i) {
  dados_mes <- dados %>%
    filter(mes == i & !is.na(prec.H) & !is.na(prec.A))
  if (nrow(dados_mes) > 0) {
    trans <- prop.table(table(dados_mes$prec.H, dados_mes$prec.A), margin = 1)
    return(trans)
  } else {
    return(NULL)
  }
})
names(matrizes_mensais) <- month.name # Adiciona os nomes dos meses

## Função para calcular matriz transição n mês a mês
n_passos <- 30
matrizes_npassos_mensais <- lapply(1:12, function(m) {
  dados_mes <- dados %>% filter(mes == m)
  dados_mes <- dados_mes[!is.na(dados_mes$prec.H), ]
  if (nrow(dados_mes) > 1) {
    seq_mes <- as.character(dados_mes$prec.H) 
    fit <- markovchainFit(data = seq_mes)
    matriz_1passo <- fit$estimate
    matriz_npassos <- matriz_1passo^n_passos
    return(matriz_npassos)
  } else {
    return(NULL)
  }
})
names(matrizes_npassos_mensais) <- month.name

Montando as saídas em uma tabela:

MÊS Matriz p^1 Matriz p^n … n=30 Probabilidade no Mês
Janeiro (Verão) 0.83 0.16
0.40 0.60		 0.71 0.29
0.71 0.29		 Prob. de Chover: 29%
Prob. de Não Chover: 71%
Fevereiro (Verão) 0.80 0.20
0.46 0.54		 0.72 0.28
0.72 0.28		 Prob. de Chover: 28%
Prob. de Não Chover: 72%
Março (Verão/Outono) 0.75 0.25
0.51 0.49		 0.68 0.32
0.68 0.32		 Prob. de Chover: 32%
Prob. de Não Chover: 68%
Abril (Outono) 0.88 0.12
0.72 0.28		 0.85 0.15
0.85 0.15		 Prob. de Chover: 15%
Prob. de Não Chover: 85%
Maio (Outono) 0.97 0.03
0.85 0.15		 0.96 0.04
0.96 0.04		 Prob. de Chover: 04%
Prob. de Não Chover: 96%
Junho (Outono/Inverno) 0.99 0.01
1.00 0.00		 0.98 0.02
0.98 0.02		 Prob. de Chover: 02%
Prob. de Não Chover: 98%
Julho (Inverno) 0.99 0.01
NA NA		 1.00 0.00
1.00 0.00		 Prob. de Chover: 0%
Prob. de Não Chover: 100%
Agosto (Inverno) 0.99 0.01
0.66 0.34		 0.99 0.01
0.99 0.01		 Prob. de Chover: 01%
Prob. de Não Chover: 99%
Setembro (Inverno/Primavera) 0.94 0.06
0.76 0.24		 0.95 0.05
0.95 0.05		 Prob. de Chover: 05%
Prob. de Não Chover: 95%
Outubro (Primavera) 0.81 0.19
0.60 0.40		 0.79 0.21
0.79 0.21		 Prob. de Chover: 21%
Prob. de Não Chover: 79%
Novembro (Primavera) 0.68 0.32
0.46 0.54		 0.58 0.42
0.58 0.42		 Prob. de Chover: 42%
Prob. de Não Chover: 58%
Dezembro (Primavera / Verão) 0.73 0.27
0.47 0.53		 0.63 0.37
0.63 0.37		 Prob. de Chover: 37%
Prob. de Não Chover: 63%

Pela tabela, observamos que, nos meses finais da primavera e durante o verão, há maior probabilidade de chuva na cidade de Niquelândia. Por outro lado, nos meses de inverno, predomina a seca e a chance de ocorrência de chuva é mínima.

Nota-se, por exemplo, que no mês de julho houve apenas um único dia com precipitação, seguido por um longo período de estiagem.

Mesmo com o uso apenas de estatísticas descritivas, é possível identificar o padrão das precipitações no município. Essas informações são úteis para a agricultura, especialmente no planejamento do plantio de culturas ou na escolha do período mais adequado para a semeadura.

Agora, partiremos para uma análise mais avançada, identificando as melhores variáveis que podem prever a ocorrência de chuva no dia seguinte.

Análise utilizando outras métricas

Apesar de termos, em números, a probabilidade de chuva para qualquer mês do ano, ainda não sabemos exatamente quão precisa essa informação é — ou seja, o quanto ela acerta ou erra.

Para isso, precisamos analisar outras métricas, como sensibilidade (prever chuva e realmente chover), especificidade (prever ausência de chuva e realmente não chover), taxas de falso positivo e falso negativo, além de índices como Youden e Kappa e curvas ROC.

Até agora, estávamos avaliando apenas a probabilidade de chover no dia seguinte com base no dia anterior. Agora, vamos estender essa análise para outras variáveis preditoras do nosso banco de dados e escolher um ponto de corte que maximize as métricas definidas previamente.

Treinamento / Teste

Para evitar vieses, faremos a divisão do nosso banco de dados em conjunto de treinamento e teste, com proporção de 70% para 30%. Dessa forma, utilizamos 70% dos dados para treinar o modelo e calcular as métricas, e testaremos o desempenho da predição nos 30% restantes, garantindo uma avaliação mais realista e robusta.

set.seed(107732) #seed do meu RA
indices <- createDataPartition(dados$prec.A, p = 0.7, list = FALSE)
treinamento <- dados[indices, ] 
teste <- dados[-indices, ]

Comparações via Área de Curva ROC

Com os datasets devidamente separados, iremos analisar todas as áreas de curva ROC para cada variável no banco de treinamento.

rocs <- pROC::roc(prec.A ~., data = treinamento, quiet = T)
sort(sapply(rocs, pROC::auc))
##         mes   vento.max      ur.dif   vento.med    temp.med    temp.max 
##   0.5379659   0.5845123   0.6111576   0.6175298   0.6594360   0.6622072 
##    temp.min      prec.H      ur.raz pressao.med    prec.tot    temp.dif 
##   0.6758036   0.6802699   0.6878281   0.7502925   0.7541165   0.7750345 
##    temp.raz      ur.min    temp.orv      ur.med 
##   0.7919908   0.8070952   0.8258067   0.8352656

Observamos que as variáveis com maior área de curva ROC são:

  1. ur.med : umidade relativa média do dia (83.52%)
  2. temp.orv: temperatura de orvalho. (82.58%)
  3. ur.min: umidade relativa minima observada. (80.70%)
  4. temp.raz: razão de temperatura (máx/min). (79.19%)

Poderíamos, a partir deste ponto, selecionar a variável preditora que apresentou a maior área sob a curva ROC (AUC) e analisar suas métricas em detalhes. No entanto, é importante destacar que a maior AUC nem sempre implica a melhor escolha prática. Uma variável com excelente desempenho global pode falhar justamente nos tipos de erro que mais nos preocupam.

Por isso, antes de escolher qualquer variável ou ponto de corte, é fundamental definir claramente nosso objetivo de estudo e, principalmente, quais tipos de erro são aceitáveis e quais devem ser evitados a todo custo. Essa definição orienta diretamente se devemos priorizar sensibilidade, especificidade ou um equilíbrio entre ambas — e é especialmente crítica em aplicações reais como previsão de chuva em contextos de alto impacto.

Objetivo

Nosso objetivo com este estudo é prever se choverá ou não em um determinado dia, com base nas variáveis observadas no dia anterior. A ideia central é construir um modelo que aprenda padrões relevantes nos dados históricos e seja capaz de antecipar eventos de chuva com boa precisão.

Além de simplesmente prever a ocorrência de chuva, buscamos entender a qualidade dessa previsão por meio de métricas como sensibilidade e especificidade. No contexto do clima de Goiás — caracterizado por longos períodos de estiagem alternados com chuvas intensas e concentradas — é especialmente importante priorizar a sensibilidade. Isso significa maximizar a capacidade do modelo de detectar corretamente os dias em que vai chover, mesmo que isso implique em alguns falsos alarmes.

A lógica por trás dessa escolha é prática: prever chuva e ela não ocorrer (falso positivo) pode gerar apenas pequenos incômodos, enquanto não prever chuva quando ela ocorre (falso negativo) pode causar impactos mais graves, como prejuízos na agricultura, alagamentos urbanos e falhas na preparação logística. Assim, o modelo será otimizado não apenas para acurácia geral, mas com foco em minimizar os riscos reais associados aos erros de previsão.

Análise via pacote CutpointR

var_alvo <- "prec.A"
predictor_vars <- setdiff(names(treinamento[,-1]), var_alvo)
results <- map_dfr(predictor_vars, function(var) {
  cp <- cutpointr(
    data = treinamento,
    x = !!sym(var),
    class = !!sym(var_alvo),
    pos_class = 1,  # Chuva
    method = maximize_metric,
    metric = youden)
  
  # Extraindo as métricas:
  tibble(
    variable = var,
    optimal_cutoff = cp$optimal_cutpoint,
    sensitivity = cp$sensitivity,
    specificity = cp$specificity,
    youden = cp$youden,
    kappa = cp$kappa
  )
})

Temos:

kable(results)
variable optimal_cutoff sensitivity specificity youden
prec.tot 0.2000000 0.6800000 0.7957398 0.4757398
pressao.med 938.1000000 0.7942857 0.5831114 0.3773972
temp.orv 17.5000000 0.8990476 0.6565234 0.5555710
temp.max 28.8000000 0.3714286 0.8904526 0.2618812
temp.med 24.0000000 0.6266667 0.6101179 0.2367846
temp.min 19.2000000 0.9314286 0.3963484 0.3277770
ur.med 67.6000000 0.8266667 0.7014074 0.5280740
ur.min 37.0000000 0.8285714 0.6413085 0.4698799
vento.max 10.3000000 0.4438095 0.6854317 0.1292412
vento.med 1.9000000 0.7504762 0.4351464 0.1856226
mes 10.0000000 0.4857143 0.7911754 0.2768896
temp.raz 0.6447368 0.7295238 0.7048307 0.4343545
temp.dif 10.6000000 0.6342857 0.7755801 0.4098658
ur.raz 0.6577181 0.4819048 0.8394827 0.3213875
ur.dif 27.2000000 0.4571429 0.7078737 0.1650166

Para critério de escolha, iremos escolher os três maiores índices de Youden (que balanceia sensibilidade com especificidade. Iremos analisar um a um.

Escolhas

temp.orv

cp.temp.orv <- cutpointr(
  data = treinamento,
  x = temp.orv,
  class = prec.A,
  pos_class = 1,
  method = maximize_metric,
  metric = youden)
summary(cp.temp.orv)
## Method: maximize_metric 
## Predictor: temp.orv 
## Outcome: prec.A 
## Direction: >= 
## 
##     AUC    n n_pos n_neg
##  0.8258 3154   525  2629
## 
##  optimal_cutpoint youden    acc sensitivity specificity  tp fn  fp   tn
##              17.5 0.5556 0.6969       0.899      0.6565 472 53 903 1726
## 
## Predictor summary: 
##     Data Min.    5% 1st Qu. Median     Mean 3rd Qu.  95% Max.       SD NAs
##  Overall  0.8  6.90    11.5   16.5 15.14759    19.1 20.3 21.3 4.491920   0
##        0  0.8  6.60    10.7   15.1 14.38429    18.6 20.0 21.3 4.488488   0
##        1  7.0 16.02    18.6   19.4 18.96990    19.9 20.5 21.2 1.676103   0
plot_roc(cp.temp.orv)

plot_sensitivity_specificity(cp.temp.orv)

Temperatura de Orvalho

  • AUC: 0.826
  • Sensibilidade: 0.899
  • Especificidade: 0.657
  • FN: 53 (melhor)

Excelente para prever corretamente os dias com chuva. Apesar de gerar mais falsos positivos, é a variável mais sensível, ou seja, a que mais “captura” os eventos reais de chuva.

ur.med

cp.ur.med <- cutpointr(
  data = treinamento,
  x = ur.med,
  class = prec.A,
  pos_class = 1,
  method = maximize_metric,
  metric = youden)
summary(cp.ur.med)
## Method: maximize_metric 
## Predictor: ur.med 
## Outcome: prec.A 
## Direction: >= 
## 
##     AUC    n n_pos n_neg
##  0.8353 3154   525  2629
## 
##  optimal_cutpoint youden    acc sensitivity specificity  tp fn  fp   tn
##              67.6 0.5281 0.7223      0.8267      0.7014 434 91 785 1844
## 
## Predictor summary: 
##     Data Min.     5% 1st Qu. Median     Mean 3rd Qu.  95% Max.       SD NAs
##  Overall 18.6 29.665    46.0   61.0 59.98554    76.0 85.5 94.5 17.77696   0
##        0 18.6 28.600    43.1   57.1 56.66105    71.1 82.3 90.5 16.99772   0
##        1 33.8 56.020    72.0   79.1 76.63333    84.4 89.4 94.5 10.92695   0
plot_roc(cp.ur.med)

plot_sensitivity_specificity(cp.ur.med)

Umidade relativa média:

  • AUC: 0.835 (melhor)
  • Sensibilidade: 0.827
  • Especificidade: 0.701
  • FN: 91

Apresenta o melhor desempenho global (maior AUC), com bom equilíbrio entre sensibilidade e especificidade. Embora não seja tão sensível quanto temp.orv, é mais equilibrada, resultando em menos falsos positivos. É uma boa escolha para quando temos que equilibrar sensibilidade + especificidade.

prec.total

cp.prec.tot <- cutpointr(
  data = treinamento,
  x = prec.tot,
  class = prec.A,
  pos_class = 1,
  method = maximize_metric,
  metric = youden)
summary(cp.prec.tot)
## Method: maximize_metric 
## Predictor: prec.tot 
## Outcome: prec.A 
## Direction: >= 
## 
##     AUC    n n_pos n_neg
##  0.7541 3154   525  2629
## 
##  optimal_cutpoint youden    acc sensitivity specificity  tp  fn  fp   tn
##               0.2 0.4757 0.7765        0.68      0.7957 357 168 537 2092
## 
## Predictor summary: 
##     Data Min. 5% 1st Qu. Median     Mean 3rd Qu.   95% Max.        SD NAs
##  Overall    0  0       0      0 2.760875     0.4 18.47 89.8  8.483208   0
##        0    0  0       0      0 1.621833     0.0 10.40 89.8  6.234897   0
##        1    0  0       0      2 8.464762    11.4 41.52 88.6 14.105175   0
plot_roc(cp.prec.tot)

plot_sensitivity_specificity(cp.prec.tot)

Precipitação total no dia anterior

  • AUC: 0.754 (pior)
  • Sensibilidade: 0.680 (pior)
  • Especificidade: 0.796 (melhor)
  • FN: 168 (pior)

Apesar de ser eficaz para prever dias sem chuva (alta especificidade), falha consideravelmente em prever os dias em que chove. Essa característica pode ser perigosa, pois muitos eventos de chuva passariam despercebidos. Poderia ser usada como filtro auxiliar, mas não como preditor principal.

Conclusão

No contexto climático da cidade de Niquelândia, Goiás — onde prever corretamente a ocorrência de chuva é mais importante do que evitar falsos alarmes (prever chuva e não chover) — a variável temp.orv se destaca pela alta sensibilidade e menor número de falsos negativos. A ur.med é uma alternativa mais equilibrada e robusta, podendo ser combinada com temp.orv para aumentar a performance do modelo.

Já prec.tot, embora intuitiva, apresenta desempenho inferior, especialmente na detecção de eventos de chuva, e não deve ser usada isoladamente para previsão.

Aplicando na base Teste

Como escolhemos a variável temperatura de orvalho como nossa principal preditora para chuva/não chuva, iremos aplica-la na base teste (30% dos dados) e ver como ela se comporta:

preds <- predict(cp.temp.orv, teste)
confusionMatrix(table(preds, teste$prec.A), positive = "1")
## Confusion Matrix and Statistics
## 
##      
## preds   0   1
##     0 769  25
##     1 357 199
##                                           
##                Accuracy : 0.717           
##                  95% CI : (0.6922, 0.7409)
##     No Information Rate : 0.8341          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.3585          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.8884          
##             Specificity : 0.6829          
##          Pos Pred Value : 0.3579          
##          Neg Pred Value : 0.9685          
##              Prevalence : 0.1659          
##          Detection Rate : 0.1474          
##    Detection Prevalence : 0.4119          
##       Balanced Accuracy : 0.7857          
##                                           
##        'Positive' Class : 1               
##
  • Sensibilidade: 0.888 → Excelente → o modelo acerta ~88% dos dias com chuva
  • Especificidade: 0.683 → OK → o modelo acerta ~68% dos dias sem chuva
  • Acurácia: 0.717 → moderada
  • Kappa: 0.36 → fraco
  • PPV (valor preditivo positivo): 0.36 → ruim
  • NPV (valor preditivo negativo): 0.97 → excelente
  • Prevalência amostral de chuva: 16.6% dos dias

O modelo está muito bom para detectar dias com chuva (alta sensibilidade).

O modelo é ótimo para detectar dias com chuva pela sua alta sensibilidade. O alto valor de predição negativa mostra que quando o modelo diz que não vai chover, ele quase sempre acerta. Ele está fazendo o que se espera dele: capturando bem os dias com chuva (sensível), mas com muitos alarmes falsos (PPV baixo).

Isso é justificado pelo contexto regional da cidade, onde o erro mais custoso é não prever chuva quando ela acontece.

Para finalizar, tentarei usar duas variáveis regressoras usando modelo logístico simples. Irei escolher temperatura de orvalho e umidade relativa média.

Regressão Logística Simples

mod <- glm(prec.A ~ temp.orv + ur.med, data = treinamento, family = binomial)
probs <- predict(mod, newdata = teste, type = "response")
preds_comb <- ifelse(teste$temp.orv >= 17.5 & teste$ur.med >= 67.6, 1, 0)
confusionMatrix(factor(preds_comb), factor(teste$prec.A), positive = "1")
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 829  36
##          1 297 188
##                                           
##                Accuracy : 0.7533          
##                  95% CI : (0.7294, 0.7761)
##     No Information Rate : 0.8341          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.3924          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.8393          
##             Specificity : 0.7362          
##          Pos Pred Value : 0.3876          
##          Neg Pred Value : 0.9584          
##              Prevalence : 0.1659          
##          Detection Rate : 0.1393          
##    Detection Prevalence : 0.3593          
##       Balanced Accuracy : 0.7878          
##                                           
##        'Positive' Class : 1               
##

Observamos uma diminuição pequena na sensibilidade e uma melhora na especificidade. Apesar disso, acredito que o modelo apenas com uma variável seria ideal, pois esse diminuiu nossa sensibilidade, o que não é ideal.

Poderiamos nos aprofundar em melhorar o modelo, combinando variáveis, analisando residuos, etc. porém não é o assunto desse trabalho.

Conclusão e Resultados

A análise revelou que Niquelândia apresenta um padrão climático bem definido, com longos períodos de estiagem no inverno e alta concentração de chuvas durante o verão. As matrizes de transição mensais reforçam essa sazonalidade, mostrando probabilidades de ocorrência de chuva próximas a 0% em julho e superiores a 30% em meses como janeiro e março.

As cadeias de Markov se mostraram úteis para modelar a persistência dos estados climáticos (chuva ou não chuva) ao longo do tempo. A matriz de transição em múltiplos passos indicou que, no longo prazo, há uma tendência de estabilização em torno de 16% de dias chuvosos, independentemente do estado inicial.

Complementando a análise probabilística, métricas de avaliação como curva ROC e AUC foram utilizadas para mensurar a capacidade de diferentes variáveis meteorológicas em prever a ocorrência de chuva no dia seguinte. Variáveis como temperatura de orvalho, umidade relativa e precipitação total no dia apresentaram AUCs superiores a 0.75, indicando um bom poder discriminativo. Além disso, a análise das métricas de sensibilidade e especificidade permitiu avaliar o desempenho dos modelos de forma mais robusta, considerando o desequilíbrio entre dias chuvosos e secos.

Em resumo, mesmo com abordagens estatísticas relativamente simples, foi possível extrair informações valiosas sobre o comportamento climático da região.

Referências:

(Tettey et al. 2017)

(Rotondi 2010)

(Thiele and Hirschfeld 2021)

Rotondi, Michael A. 2010. “To Ski or Not to Ski: Estimating Transition Matrices to Predict Tomorrow’s Snowfall Using Real Data.” Journal of Statistics Education 18 (03).
Tettey, Meshach, Francis T. Oduro, David Adedia, and Daniel A. Abaye. 2017. “Markov Chain Analysis of the Rainfall Patterns of Five Geographical Locations in the South Eastern Coast of Ghana.” Earth Perspectives 4 (1). https://doi.org/10.1186/s40322-017-0042-6.
Thiele, Christian, and Gerrit Hirschfeld. 2021. “Cutpointr: Improved Estimation and Validation of Optimal Cutpoints in r” 98. https://doi.org/10.18637/jss.v098.i11.