Análisis exploratorio de datos de temperatura

1. Contextualización

Este informe presenta un Análisis Exploratorio de Datos (EDA) para una serie de tiempo de temperatura en Brasil, específicamente en la región sureste en la estación “Ecologia agricola”. El objetivo es comprender la estructura de los datos y preparar el conjunto para modelos de pronóstico.

Los datos fueron extraídos del Climate Weather Surface of Brazil - Hourly, el cual contiene registros horarios de 623 estaciones meteorológicas en Brasil entre los años 2000 y 2021. Estos datos fueron obtenidos del Instituto Nacional de Meteorología de Brasil (INMET). Estos fueron descargados desde Kaggle, Climatología de la superficie de Brasil con frecuencia horaria.

El conjunto de datos descargado corresponde al archivo southeast.csv, que tiene un tamaño aproximado de 2.5 GB. Contiene información de múltiples estaciones meteorológicas, con variables como temperatura, humedad, velocidad del viento y radiación solar, entre otras.

1.1 Descripción del conjunto de datos

Las principales variables de interés incluyen:

Variable Formato Descripción
Data (YYYY-MM-DD) Fecha de registro
Hora (HH:00) Hora de observación
TEMPERATURA DO AR (°C) double Temperatura del aire en grados Celcius

Para reducir el volumen de información y enfocarnos en una única serie de tiempo, se seleccionaron los datos de la estación ECOLOGIA AGRICOLA, identificada con el código A601. Esta estación ha estado operativa desde el 07/05/2000, y se encuentra en:

  • Latitud: -22.8
  • Longitud: -43.6833
  • Altitud: 33 metros sobre el nivel del mar

Para realizar esta selección, se empleó localamente una consulta a SQL:

SELECT “index”, “Data”, “Hora”, “TEMPERATURA DO AR - BULBO SECO, HORARIA (°C)” FROM “southeast” WHERE “station_code” = ‘A601’ ORDER BY “Data” ASC, “Hora” ASC;

Finalmente, se exportó esta consulta como “soudeste_a601.csv”.

2. Análisis preliminar

A continuación, se hace un carague de los datos para una revisión general, con el propósito de realizar un análisis perliminar. Para esto se cargan las librerías que utilizarán durante el desarrollo del EDA.

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(skimr)
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(fpp3)
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.1 ──
## ✔ tsibble     1.1.6     ✔ feasts      0.4.1
## ✔ tsibbledata 0.4.1     ✔ fable       0.4.1
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()    masks base::date()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval()  masks lubridate::interval()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ tsibble::setdiff()   masks base::setdiff()
## ✖ tsibble::union()     masks base::union()

Para empezar, se cargan los datos dentro del entorno de trabajo.

temp = read_csv("soudeste_a601.csv") %>% 
  rename(Temperatura = `TEMPERATURA DO AR - BULBO SECO, HORARIA (°C)`) %>% 
  select(-index)
## Rows: 183936 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl  (2): index, TEMPERATURA DO AR - BULBO SECO, HORARIA (°C)
## date (1): Data
## time (1): Hora
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Posteriormente, los datos se procesan para que tengan un formato POSIXct para trabajar con el tiempo de modo que tanto la fecha como la hora sean una misma variable llamada “Date”.

tsibble_temp = temp %>% 
  mutate(Date = as.POSIXct(paste(Data, Hora), format = "%Y-%m-%d %H:%M:%S")) %>% 
  select(Date, Temperatura) %>% 
  as_tsibble(index = Date)

head(tsibble_temp)
## # A tsibble: 6 x 2 [1h] <?>
##   Date                Temperatura
##   <dttm>                    <dbl>
## 1 2000-05-07 00:00:00       -9999
## 2 2000-05-07 01:00:00       -9999
## 3 2000-05-07 02:00:00       -9999
## 4 2000-05-07 03:00:00       -9999
## 5 2000-05-07 04:00:00       -9999
## 6 2000-05-07 05:00:00       -9999

A simple vista, esta primera aproximación exhibe el hecho de que existen datos de temperatura erróneos con valores de -9.999. Para esto se emplea la función skim() que muesta un resumen estadísitco de los datos crudos.

skim(tsibble_temp$Date) %>% 
  select(-n_missing, -complete_rate)
Data summary
Name tsibble_temp$Date
Number of rows 183936
Number of columns 1
_______________________
Column type frequency:
POSIXct 1
________________________
Group variables None

Variable type: POSIXct

skim_variable min max median n_unique
data 2000-05-07 2021-04-30 23:00:00 2010-11-02 23:30:00 183936

En cuanto al tiempo, se observa que la serie de tiempo comenzó el 07 de mayo del 2000 a las 00:00:00 y su finalización fue el 30 de abril del 2021 a las 23:00:00. Dado que las observaciones tienen una frecuencia horaria la cantidad total de datos tomados fue de 183936.

skim(tsibble_temp$Temperatura)
Data summary
Name tsibble_temp$Temperatura
Number of rows 183936
Number of columns 1
_______________________
Column type frequency:
numeric 1
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
data 0 1 -624.99 2466.09 -9999 19.9 23 26.3 41.3 ▁▁▁▁▇

Por otro lado, en cuanto a los valores de temperatura, anque no hay NA’s sí existen valores erróneos. Esto ha hecho que la distribución de los datos sea significativamente diferente y además ha distorsionado medidas como la media que son susecpible a datos anómalos.

2.1 Tratamiento de outliners y NA’s

Para un tratamiento automático de estos valores erróneos y los faltantes se emplea la librería forecast. Sin embargo, es preciso convertir la serie de tiempo a un objeto ts, simultaneamente, se emplea la función tsclean para limpiar estos errores. Finalmente, se reescribe la variable tsibble_temp con los valores imputados.

Por medio de la función tsoutliers() se puede aprecira que el total de datos erróneos es de 11.906.

length(tsoutliers(ts(tsibble_temp$Temperatura))$index)
## [1] 11906

A partir de esto se limipian los datos con la función tsclean().

temp_clean = tsclean(ts(tsibble_temp$Temperatura))

tsibble_temp = 
  tibble(
    Date = tsibble_temp$Date, 
    Temperatura = temp_clean
  ) %>% 
    as_tsibble(index = Date)

autoplot(tsibble_temp, Temperatura)+
  labs(title = "Temperatura del aire",
       subtitle = "Sudeste brasileño - estación A601",
       y = "Grados celcius [ºC]",
       x = "Fecha [1h]")

Como se puede observar, en la gráfica de la serie de tiempo, ya no se hayan estos datos erróneos de -9.999 ºC.

2.2 Estadísticas descriptiva

A continuación, se presenta un resumen estadístico para la temperatura de la zona de estudio.

skim(tsibble_temp$Temperatura) %>% 
  select(-n_missing, -complete_rate)
Data summary
Name tsibble_temp$Temperatura
Number of rows 183936
Number of columns 1
_______________________
Column type frequency:
numeric 1
________________________
Group variables None

Variable type: numeric

skim_variable mean sd p0 p25 p50 p75 p100 hist
data 23.9 4.54 8.6 20.8 23.5 26.8 41.3 ▁▅▇▂▁

3. Visualización temporal

En el caso de los datos de series temporales, el gráfico obvio con el que empezar es un gráfico de tiempo. Este se representó anteriormente donde se expresó la temperatura con frecuencia horaria. Sin embargo, a continuación, se realizará un tratamiento para gráficar las temperaturas en diferentes periodos, para así pasar de horas a días, semanas, meses y finalmente años.

3.1 Agrupación diaría

tsibble_temp %>% 
  mutate(Diario =  as_date(Date)) %>% 
  index_by(Diario) %>%  
  summarise(Temp_promedio = mean(Temperatura)) %>% 
  autoplot(Temp_promedio)

3.2 Agrupación semanal

tsibble_temp %>% 
  mutate(Semanal =  yearweek(Date)) %>% 
  index_by(Semanal) %>%  
  summarise(Temp_promedio = mean(Temperatura)) %>% 
  autoplot(Temp_promedio)

3.3 Agrupación mensual

tsibble_temp %>% 
  mutate(Mensual =  yearmonth(Date)) %>% 
  index_by(Mensual) %>%  
  summarise(Temp_promedio = mean(Temperatura)) %>% 
  autoplot(Temp_promedio)

3.4 Agrupación anual

tsibble_temp %>% 
  mutate(Año =  year(Date)) %>% 
  index_by(Año) %>%  
  summarise(Temp_promedio = mean(Temperatura)) %>% 
  autoplot(Temp_promedio)+
  geom_smooth(method = "lm", se= FALSE)
## `geom_smooth()` using formula = 'y ~ x'

4. Análisis de patrones recurrentes

El gráfico a continuación es similar a un gráfico temporal, salvo que los datos se representan en función de las “estaciones” individuales en las que se observaron.

4.1 Análisis horario de la temperatura

tsibble_temp %>% 
  gg_season(Temperatura, period = "day")

En este gráfico cada linea representa un día y el eje x (Date) exhibe el comportamiento del variable en una frecuencia horaria. De modo que se puede observar la variación de la temperatura en el aire a lo lagro de un día. Una de las primeras característica notables de la serie, es que hay una especie de ruido de lineas diarías aglutinadas y con un comportamiento totalmente anómalo al rededor de los 35 y 30 grados celcius. Esto se debe probablemente a un error en la imputación de datos faltantes, ya que como se pudo observar en las gráficas de visualización temporal con sus agrupaciones durante el año 2005 la estación dejo de tomar datos, generando NA’s que fueron reemplazados automáticamente con la función tsclean().

Para obsrevar esta tendencia horaria más fácilmente también se optó por la realización de una gráfica de cajas y bigotes agrupadas por hora. De esta manera se podrá obsevar el cambio de la temperatura durante el ciclo diarío.

tsibble_temp %>% 
  mutate(hora = hour(Date)) %>% 
  index_by(hora) %>% 
  ggplot(aes(x = factor(hora), y = Temperatura)) +
  geom_boxplot()

Siguiendo el mismo patron que el gráfico anterior, se puede apreciar que la temperatura del aire alcanza su pico máixmo a las 16 horas del día. Esto se debe a la dinámica del balance de energía en la superficie terrestre y la atmósfera. Es correcto considerar que la radiación solar alcanza su máximo al mediodía. Sin embargo, la temperatura del aire no alcanza su pico inmediatamente, ya que la roca continental y la atmósfera siguen acumulando calor después de este instante. La superficie terrestre absorbe la radiación y la reemite como calor (radiación infrarroja), calentando el aire por conducción, a través del contacto directo entre el aire y la roca; o convección, que implica el movimiento de calor a través de un medio fluido (el aire). Este proceso continúa hasta que la energía saliente iguala a la energía entrante, lo que según nos muestra los datos, suele ocurrir entre 15:00 y 17:00 horas, aproximadamente.

4.1 Análisis semanal de la temperatura

tsibble_temp %>% 
  gg_season(Temperatura, period = "week")

A nivel semanal la estacionalidad es evidente, y repite el patron diarío que se observo anteriormente. Esto viene determinado por el ciclo solar.

4.2 análisis mensual de la temperatura

tsibble_temp %>% 
  gg_season(Temperatura, period = "year")

También se realiza un gráfico de cajas y bigotes para apreciar la variación estacional a nivel de cada mes.

tsibble_temp %>% 
  mutate(mes = month(Date, label = TRUE)) %>% 
  index_by(mes) %>% 
  ggplot(aes(x = mes, y = Temperatura)) +
  geom_boxplot()

En este caso, el patrón que se observa en el sudeste de Brasil, las temperaturas aumentan en la época seca de diciembre, enero y febrero, eventualmente, disminuyen de forma progresiva.

En primer lugar, el verano austral que comprende los meses de diciembre, enero y febrero, el Sol está casi perpendicular sobre el hemisferio sur, lo que implica mayor radiación solar y días más largos. Esta mayor insolación calienta la superficie terrestre y el aire, elevando las temperaturas. La atmósfera también es más seca en muchas áreas, lo que permite una mayor absorción de calor por parte del suelo y reduce el enfriamiento por nubosidad o precipitaciones.

Por otro lado, la Zona de Convergencia Intertropical (ZCIT) juega un papel clave en la distribución de temperaturas en Brasil, especialmente en el norte y noreste del país. Sin embargo, en el sudeste de Brasil, su influencia es indirecta, afectando la circulación de los vientos, la humedad y, en consecuencia, las temperaturas a lo largo del año.

La ZCIT se ubica entre los 14°N y 2°S a lo largo del año y se desplaza hacia el sur en ciertas épocas. En el verano austral (diciembre, enero y febrero), la ZCIT aún no ha alcanzado su punto más al sur y se mantiene más alejada del sudeste. Como resultado, el sudeste de Brasil queda bajo la influencia de sistemas de alta presión subtropicales, como el Anticiclón del Atlántico Sur, que inhibe la formación de nubes y favorece un clima más seco y caluroso. La menor humedad relativa reduce la nubosidad y las lluvias, permitiendo que la superficie terrestre absorba más calor y generando temperaturas más elevadas.

En el otoño e invierno (mayo a septiembre - MJJAS), la ZCIT ya ha migrado hacia el Hemisferio Norte, alejándose aún más del sudeste de Brasil. Esto provoca un cambio en la circulación atmosférica que facilita la llegada de masas de aire frío de origen polar, lo que enfría la región.

El aumento de temperatura en el sudeste de Brasil durante la época seca de diciembre a febrero se debe a varios factores. En primer lugar, la radiación solar es más intensa en esta época, ya que el Sol está casi perpendicular sobre el hemisferio sur, aumentando la cantidad de energía recibida. Además, con la ZCIT aún alejada, hay menos nubosidad, lo que permite una mayor absorción de calor en la superficie. La presencia del Anticiclón del Atlántico Sur también contribuye a mantener un tiempo seco y cálido, ya que dificulta la formación de lluvias.

tsibble_temp %>% 
  mutate(Mensual =  yearmonth(Date)) %>% 
  index_by(Mensual) %>%  
  summarise(Temp_promedio = mean(Temperatura)) %>% 
  gg_subseries(Temp_promedio)

A nivel anual también se repite este patron estacional de los meses, sin embargo, nuevamente se observa que durante el año 2005 la imputación de datos no resultó eficiente al captura el patron estacional mesual ya que durante los meses de agosto, septiembre y especialmente octubre se logra apreciar un pico desaforado y totalmente anómalo.

tsibble_temp %>% 
  mutate(Mensual =  yearmonth(Date)) %>% 
  index_by(Mensual) %>%  
  summarise(Temp_promedio = mean(Temperatura)) %>% 
  gg_season(Temp_promedio)

Este error de imputación de datos se puede volver a observar en este gráfico, el cual suaviza la serie al tomar los promedios mensuales para obsrevar el comportamiento estacional a nivel mensual. Aquí se puede apreciar también como existe un cambio importante en la pendiente del cambio de la temperatura durante la época de transición entre la temporada seca y de lluvia comprendida entre marzo abril y mayo. Dando paso así al invierno austral, es decir, junio, julilo, agosto.

5. Descomposición de la serie

Sea \(Y_t\) una serie de tiempo. Su descomposición aditiva se expresa como: \(Y_t = T_t + S_t + E_t\)

donde:

  • \(T_t\) es la componente de tendencia.
  • \(S_t\) es la componente estacional.
  • \(E_t\) es la componente de error o irregular

STL es un método versátil y robusto para descomponer series temporales. STL es un acrónimo de “Descomposición estacional y de tendencias utilizando Loess”, mientras que loess es un método para estimar relaciones no lineales. El método STL fue desarrollado por RB Cleveland et al. (1990) , y luego ampliado para manejar múltiples patrones estacionales por Bandara et al. (2022).

STL destaca por su flexibilidad y robustez, especialmente en el análisis de series con estacionalidad compleja o cambiante. Su capacidad para adaptarse a diversos tipos de estacionalidad, permitir la variación del componente estacional y ajustar la suavidad de la tendencia.

5.1 Descomposición STL

En el caso de la serie de tiempo original, al ser de frecuencia horaria, las componentes estacionales pueden dividirse en frecuencias diarías, semanales y anuales.

dcmp = tsibble_temp %>% 
  model(stl = STL(Temperatura))


components(dcmp) %>% 
  autoplot()

Como se puede observar la serie de tiempo se ha descompuesto en cuatro componentes de estacionalidad principales, las cuales también fueron analizadas anteriormente, la estacionalidad semanal, diaria y anual.

5.2 Descomposición clasica

Para simplificar la serie, como se hizo en la subsección 3.4, se hará un suavizado de la misma para quedarnos sólo con el componente estacional anual.

tsibble_temp %>% 
  mutate(Mensual =  yearmonth(Date)) %>% 
  index_by(Mensual) %>%  
  summarise(Temp_promedio = mean(Temperatura)) %>%
  select(Mensual, Temp_promedio) %>% 
  model(classical_decomposition(Temp_promedio, type = "additive")) %>% 
  components() %>% 
  autoplot() 
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).

6. Gráfica de autocorrelación

Así como la correlación mide el grado de una relación lineal entre dos variables, la autocorrelación mide la relación lineal entre valores rezagados de una serie temporal.

Los coeficientes de autocorrelación para los datos de temperatura se pueden calcular utilizando la función ACF()

6.1 Autocorrelación de frecuencia horaria

tsibble_temp %>% 
  ACF(Temperatura)  %>% 
  autoplot()

En este gráfico se puede notar el claro patrón estacional de los datos,los picos tienden a estar separados por una frecuencia de 24 horas (un día) al igual que los valles.

Las líneas azules discontinuas indican si las correlaciones son significativamente diferentes de cero. En este claro, evidentemente, no es así. La serie es fuertemente estacional y además, los datos de temperatura tienen una tendencia creciente, lo cual implica que las las autocorrelaciones para los pequeños rezagos tienden a ser grandes y positivas porque las observaciones cercanas en el tiempo también están cercanas en valor. Por lo tanto, la ACF de una serie temporal con tendencia tiende a tener valores positivos que disminuyen lentamente a medida que aumentan los rezagos.

6.2 Autocorrelación de frecuencia mensual

tsibble_temp %>% 
  mutate(Mensual =  yearmonth(Date)) %>% 
  index_by(Mensual) %>%  
  summarise(Temp_promedio = mean(Temperatura)) %>%
  select(Mensual, Temp_promedio) %>% 
  ACF(Temp_promedio)  %>% 
  autoplot()

Para este caso también se obsreva un comportamiento similar, es decir, la presencia de una tendencia creciente que acorta los valores en retargo alejados y la notable componente estacional que separa los valles y crestas con una frencuencia de 12 meses (un año).

7. Prueba de hipotesis

7.1 Dickey-Fuller Aumentado (ADF)

#install.packages("urca")
library(urca)

adf_test = ur.df(tsibble_temp$Temperatura, type = "drift", selectlags = "AIC")

summary(adf_test)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.2656  -0.4598  -0.0394   0.4740   8.3142 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.283383   0.012917   99.36   <2e-16 ***
## z.lag.1     -0.053690   0.000531 -101.10   <2e-16 ***
## z.diff.lag   0.516984   0.001996  259.02   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.025 on 183931 degrees of freedom
## Multiple R-squared:  0.2802, Adjusted R-squared:  0.2802 
## F-statistic: 3.581e+04 on 2 and 183931 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -101.1039 5111.002 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78

7.2 Prueba KPSS (Kwiatkowski-Phillips-Schmidt-Shin)

#install.packages("tseries")
library(tseries)

# Aplicar la prueba KPSS
kpss_test <- kpss.test(tsibble_temp$Temperatura)
## Warning in kpss.test(tsibble_temp$Temperatura): p-value smaller than printed
## p-value
# Ver resultados
print(kpss_test)
## 
##  KPSS Test for Level Stationarity
## 
## data:  tsibble_temp$Temperatura
## KPSS Level = 1.573, Truncation lag parameter = 26, p-value = 0.01

7.3 Gráfica de densidad de los datos

tsibble_temp %>% 
  ggplot(aes(x = Temperatura, y = after_stat(density))) +
  geom_histogram(fill = "lightblue", colour = "grey60", size = .2, bins = 25) +
  geom_density()+
  labs(
    title = "Distribución temporal de la temperatura del aire",
    subtitle = "Sudeste brasileño - estación A601",
    x = "Temperatura (ºC)", 
    y = "Densidad"
  )
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

7.4 QQplot

tsibble_temp %>% 
  ggplot(aes(sample = Temperatura)) +
  geom_qq(alpha = 0.1) +
  geom_qq_line()+
  labs(
    x= "Cuantiles teóricos", 
    y= "Cuantiles de la muestra"
  )

ks.test(tsibble_temp$Temperatura, "pnorm", mean(tsibble_temp$Temperatura, na.rm = TRUE), sd(tsibble_temp$Temperatura, na.rm = TRUE))
## Warning in ks.test.default(tsibble_temp$Temperatura, "pnorm",
## mean(tsibble_temp$Temperatura, : ties should not be present for the one-sample
## Kolmogorov-Smirnov test
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  tsibble_temp$Temperatura
## D = 0.040367, p-value < 2.2e-16
## alternative hypothesis: two-sided

Dado que el valor de p es tan pequeño se rechaza la hipotesis nula y se concluye que los datos de temperatura no siguen una distribución normal.