Proceso de construccion de rangos para el deficit hidrico climatico (def) en la primera temporada humeda (Abril a Septiembre -incluyendolos- para el 2000 y su revision.

Equipo DS DI

DataIntelligence
23-11-2021

1 Introducción

Para construir la parte física del IVCC utilizamos 5 variables extraídas de GEE: TerraClimate: Clima mensual y balance hídrico climático para superficies terrestres globales, Universidad de Idaho1.

  • def
  • pdsi
  • pr
  • tmmn
  • tmmx
Nombre Unidad Min Max Escala Description
def mm 0* 4548* 0.1 Déficit hídrico climático, derivado mediante un modelo unidimensional de balance hídrico del suelo

En éste documento analizaremos la integridad de def y seguiremos los siguientes pasos:

  • Corroboraremos que la tabla entregada sea la misma sobre la que trabajamos.
  • Verificaremos que los datos faltantes (0) imputados sean coherentes.
  • Revisaremos que los rangos sean los que deben ser implementando la lógica.

2 La tabla def original y sobre la que trabajamos

https://drive.google.com/drive/folders/1h47aFXMhgzE_0Dvuv5zprqVumXDSSfQM

def_original <- read.csv("C:/Users/chris/OneDrive/Documentos/GitHub/ds_datos_climaticos/comunas_y_subcuencas/COMUNA/MEDIA/TERRACLIMATE_MEDIA_def.csv")
# datatable(def_original, extensions = 'Buttons', escape = FALSE, rownames = TRUE,
#           options = list(dom = 'Bfrtip',
#           buttons = list('colvis', list(extend = 'collection',
#           buttons = list(
#           list(extend='copy'),
#           list(extend='excel',
#             filename = 'tabla'),
#           list(extend='pdf',
#             filename= 'tabla')),
#           text = 'Download')), scrollX = TRUE))

Son las mismas.

3 La imputación

3.1 La tabla def original

Observemos las estadísticas básicas de los primeras 10 columnas de la tabla original, junto con la última que nos servirá de control:

3.2 2000

def_original_11 <- def_original[, c(1:12)]
summary(def_original_11)
##   X200001_def      X200002_def       X200003_def      X200004_def       
##  Min.   :   0.0   Min.   :   0.00   Min.   :   0.0   Min.   :   0.0000  
##  1st Qu.: 908.4   1st Qu.:  34.92   1st Qu.: 347.7   1st Qu.:   0.0176  
##  Median :1526.2   Median : 842.16   Median :1007.2   Median : 475.5803  
##  Mean   :1314.6   Mean   : 731.22   Mean   : 792.0   Mean   : 375.5261  
##  3rd Qu.:1741.1   3rd Qu.:1314.63   3rd Qu.:1170.3   3rd Qu.: 584.5591  
##  Max.   :2004.2   Max.   :1573.06   Max.   :1566.8   Max.   :1344.3411  
##   X200005_def        X200006_def       X200007_def      X200008_def     
##  Min.   :   0.000   Min.   :   0.00   Min.   :  0.00   Min.   :   0.00  
##  1st Qu.:   0.000   1st Qu.:   0.00   1st Qu.:  0.00   1st Qu.:   0.00  
##  Median :   2.991   Median :   0.00   Median :  0.00   Median :  47.38  
##  Mean   : 142.911   Mean   :  47.94   Mean   : 68.67   Mean   : 157.27  
##  3rd Qu.: 192.315   3rd Qu.:   0.00   3rd Qu.: 37.96   3rd Qu.: 202.38  
##  Max.   :1161.226   Max.   :1026.38   Max.   :997.74   Max.   :1231.68  
##   X200009_def      X200010_def      X200011_def      X200012_def  
##  Min.   :   0.0   Min.   :   0.0   Min.   :   0.0   Min.   :   0  
##  1st Qu.:   0.0   1st Qu.: 149.2   1st Qu.: 271.2   1st Qu.: 470  
##  Median :   0.0   Median : 336.6   Median : 660.5   Median :1379  
##  Mean   : 104.8   Mean   : 424.8   Mean   : 723.1   Mean   :1091  
##  3rd Qu.:   0.0   3rd Qu.: 582.0   3rd Qu.:1270.2   3rd Qu.:1694  
##  Max.   :1482.3   Max.   :1753.0   Max.   :1828.6   Max.   :1983

3.3 La tabla def imputada

def_imp <- readxl::read_xlsx("C:/Users/chris/OneDrive/Documentos/GitHub/ds_datos_climaticos/imputaciones/comuna_medias_imputadas/excel/TERRACLIMATE_MEDIA_def_imp.xlsx")
# datatable(def_imp, extensions = 'Buttons', escape = FALSE, rownames = TRUE,
#           options = list(dom = 'Bfrtip',
#           buttons = list('colvis', list(extend = 'collection',
#           buttons = list(
#           list(extend='copy'),
#           list(extend='excel',
#             filename = 'tabla'),
#           list(extend='pdf',
#             filename= 'tabla')),
#           text = 'Download')), scrollX = TRUE))
def_imp_11 <- def_imp[, c(1:12)]
summary(def_imp_11)
##   X200001_def         X200002_def         X200003_def       
##  Min.   :   0.6213   Min.   :   0.0325   Min.   :   0.0076  
##  1st Qu.: 908.3645   1st Qu.: 137.9409   1st Qu.: 418.9063  
##  Median :1526.2338   Median : 864.4821   Median :1009.5254  
##  Mean   :1314.5973   Mean   : 760.3193   Mean   : 818.1768  
##  3rd Qu.:1741.1431   3rd Qu.:1314.6284   3rd Qu.:1170.2826  
##  Max.   :2004.1576   Max.   :1573.0612   Max.   :1566.7845  
##   X200004_def         X200005_def         X200006_def         X200007_def      
##  Min.   :   0.0176   Min.   :   0.0009   Min.   :   0.3924   Min.   :  0.0318  
##  1st Qu.: 171.2903   1st Qu.:   0.9005   1st Qu.: 880.5894   1st Qu.:  1.8554  
##  Median : 506.9764   Median : 139.1904   Median : 911.6222   Median : 64.8725  
##  Mean   : 449.3724   Mean   : 208.7234   Mean   : 803.3049   Mean   :225.2529  
##  3rd Qu.: 609.4491   3rd Qu.: 264.6537   3rd Qu.: 939.3866   3rd Qu.:383.3600  
##  Max.   :1344.3411   Max.   :1161.2257   Max.   :1026.3828   Max.   :997.7352  
##   X200008_def         X200009_def         X200010_def       
##  Min.   :   0.0002   Min.   :   0.0000   Min.   :   0.0127  
##  1st Qu.:   1.3305   1st Qu.:   0.0055   1st Qu.: 199.2226  
##  Median : 160.0805   Median :   0.8657   Median : 348.6122  
##  Mean   : 211.2668   Mean   : 434.2141   Mean   : 464.9004  
##  3rd Qu.: 262.9107   3rd Qu.:1099.9329   3rd Qu.: 613.8256  
##  Max.   :1231.6774   Max.   :1482.3231   Max.   :1752.9751  
##   X200011_def         X200012_def    
##  Min.   :   0.0053   Min.   :   0.0  
##  1st Qu.: 278.8127   1st Qu.: 489.6  
##  Median : 660.5409   Median :1378.9  
##  Mean   : 726.3515   Mean   :1102.0  
##  3rd Qu.:1270.1519   3rd Qu.:1694.0  
##  Max.   :1828.5873   Max.   :1983.1

Observamos que una vez aplicada la imputación, en algunas columnas obtenemos valores extremadamente distintos y otros en los cuales apenas cambia.

Uno de los campos que más difiere es: X200006_def; X202012_def apenas difiere. Especulamos que los registros ceros en el campo X200006_def son masivos:

Contemos los ceros en el dataframe original de 12 campos (meses) del 2000 y planteémoslo como porcentajes:

lapply(def_original_11, function(x){ length(which(x==0))/length(x)})
## $X200001_def
## [1] 0.005797101
## 
## $X200002_def
## [1] 0.1884058
## 
## $X200003_def
## [1] 0.08115942
## 
## $X200004_def
## [1] 0.2492754
## 
## $X200005_def
## [1] 0.4695652
## 
## $X200006_def
## [1] 0.9072464
## 
## $X200007_def
## [1] 0.573913
## 
## $X200008_def
## [1] 0.4115942
## 
## $X200009_def
## [1] 0.773913
## 
## $X200010_def
## [1] 0.1072464
## 
## $X200011_def
## [1] 0.07536232
## 
## $X200012_def
## [1] 0.04637681

Efectivamente, $X200006_def casi no tiene datos (0.9072464% son ceros) y X200001_def casi no tiene ceros (0.005797101%)

4 Boxplots

par(mfrow = c(1, 2))
par(cex.axis=0.5)
boxplot(def_original_11,  col=rainbow(length(unique(def_original_11))))
boxplot(def_imp_11,  col=rainbow(length(unique(def_imp_11))))

par(cex.axis=0.1)

Justificación para utilizar MICE:

4.1 MICE (Multivariate Imputation by Chained Equations)

MICE (Imputación multivariante mediante ecuaciones encadenadas) es uno de los paquetes más utilizados por los usuarios de R. La creación de múltiples imputaciones en comparación con una sola (como la media) se encarga de la incertidumbre en los valores perdidos.

MICE asume que los datos faltantes son Missing at Random (MAR). Cuando decimos que faltan datos al azar, queremos por ejemplo, que el NA de un registro se puede predecir a partir de otra información del mismo. Por ejemplo, si un niño no asiste a una evaluación educativa porque está (genuinamente) enfermo, esto podría ser predecible a partir de otros datos que tenemos sobre la salud del niño, pero no estaría relacionado con lo que hubiéramos medido si el niño no ha estado enfermo.

Por ejemplo, supongamos que tenemos \(X_1, X_2… .X_k\) variables. Si \(X_1\) tiene valores perdidos, entonces se regresará (en términos estadísticos) sobre las otras variables de \(X_2 \ a \ X_k\). Los valores perdidos en \(X_1\) serán reemplazados por los valores predichos obtenidos. De manera similar, si \(X_2\) tiene valores perdidos, las variables de \(X_1, X_3 \ a \ X_k\) se utilizarán en el modelo de predicción como variables independientes. Posteriormente, los valores perdidos se reemplazarán con valores predichos2.

(MAR) ocurre cuando la falta no es aleatoria, pero cuando se puede explicar completamente por variables donde hay información completa. Dado que MAR es un supuesto imposible de verificar estadísticamente, debemos basarnos en su razonabilidad. Un ejemplo es que los hombres tienen menos probabilidades de completar una encuesta de depresión, pero esto no tiene nada que ver con su nivel de depresión, después de tener en cuenta la masculinidad.

5 El código

Leemos nuestra base de datos ya imputada:

def_imp <- readxl::read_xlsx(“C:/Users/chris/OneDrive/Documentos/GitHub/ds_datos_climaticos/imputaciones/comuna_medias_imputadas/excel/TERRACLIMATE_MEDIA_def_imp.xlsx”)

Y renombramos el campo de comuna:

Leemos la base de datos de déficit hídrico completa:

datatable(def, extensions = 'Buttons', escape = FALSE, rownames = TRUE,
          options = list(dom = 'Bfrtip',
          buttons = list('colvis', list(extend = 'collection',
          buttons = list(
          list(extend='copy'),
          list(extend='excel',
            filename = 'tabla'),
          list(extend='pdf',
            filename= 'tabla')),
          text = 'Download')), scrollX = TRUE))

Para facilitar la lectura agreguemos el nombre de las comunas y regiones3:

codigos_comunales_actuales <- readxl::read_xls("CUT_2018_v04.xls")

debemos extraer los ceros iniciales del codigo comuna 2018

library("stringr")
codigos_comunales_actuales$`Código Comuna 2018` <- str_remove(codigos_comunales_actuales$`Código Comuna 2018`, "^0+")

y hacemos la unión:

def <- merge(codigos_comunales_actuales,def , all.x = TRUE)
def <- def[,-c(1)]
datatable(def, extensions = 'Buttons', escape = FALSE, rownames = TRUE,
          options = list(dom = 'Bfrtip',
          buttons = list('colvis', list(extend = 'collection',
          buttons = list(
          list(extend='copy'),
          list(extend='excel',
            filename = 'tabla'),
          list(extend='pdf',
            filename= 'tabla')),
          text = 'Download')), scrollX = TRUE))

Eliminemos las columnas no numéricas almacenándolas en el vector COM y dejemos sólo tres decimales en el dataframe numérico:

COM <- def[,c(1:6)]

Construyamos nuestro dataframe numerico

def <- def[,7:258]
options(scipen=999)

round_df <- function(df, digits) {
  nums <- vapply(df, is.numeric, FUN.VALUE = logical(1))

  df[,nums] <- round(df[,nums], digits = digits)

  (df)
}

def <- round_df(def, digits=3)

datatable(def, extensions = 'Buttons', escape = FALSE, rownames = TRUE,
          options = list(dom = 'Bfrtip',
          buttons = list('colvis', list(extend = 'collection',
          buttons = list(
          list(extend='copy'),
          list(extend='excel',
            filename = 'tabla'),
          list(extend='pdf',
            filename= 'tabla')),
          text = 'Download')), scrollX = TRUE))

5.1 Corrección

Corregimos por 0.1

def <- def*0.1

5.2 Unimos las columnas no numéricas

def <- cbind(def, COM)
datatable(def, extensions = 'Buttons', escape = FALSE, rownames = TRUE,
          options = list(dom = 'Bfrtip',
          buttons = list('colvis', list(extend = 'collection',
          buttons = list(
          list(extend='copy'),
          list(extend='excel',
            filename = 'tabla'),
          list(extend='pdf',
            filename= 'tabla')),
          text = 'Download')), scrollX = TRUE))

5.3 Seleccionemos las temporadas húmedas:

library(dplyr) 
library(tidyverse) 

df_new  <- def[, -grep("01_def", colnames(def))]
df_new  <- df_new[, -grep("02_def", colnames(df_new))]
df_new  <- df_new[, -grep("03_def", colnames(df_new))]
df_new  <- df_new[, -grep("10_def", colnames(df_new))]
df_new  <- df_new[, -grep("11_def", colnames(df_new))]
df_new  <- df_new[, -grep("12_def", colnames(df_new))]


datatable(df_new, extensions = 'Buttons', escape = FALSE, rownames = TRUE,
          options = list(dom = 'Bfrtip',
          buttons = list('colvis', list(extend = 'collection',
          buttons = list(
          list(extend='copy'),
          list(extend='excel',
            filename = 'tabla'),
          list(extend='pdf',
            filename= 'tabla')),
          text = 'Download')), scrollX = TRUE))

5.4 Añadimos una columna con las medias por comuna

df_new$mean <- rowMeans(df_new[,1:126])

df_new$mean <- round(df_new$mean, digits=3)

datatable(df_new, extensions = 'Buttons', escape = FALSE, rownames = TRUE,
          options = list(dom = 'Bfrtip',
          buttons = list('colvis', list(extend = 'collection',
          buttons = list(
          list(extend='copy'),
          list(extend='excel',
            filename = 'tabla'),
          list(extend='pdf',
            filename= 'tabla')),
          text = 'Download')), scrollX = TRUE))

5.5 Calculemos los máximos y los mínimos de mean:

Máximo:

max <- max(df_new$mean, na.rm = TRUE)
max
## [1] 121.53

Mínimo:

uno <- min(df_new$mean, na.rm = TRUE)
min <- min(df_new$mean, na.rm = TRUE)
uno
## [1] 10.548

5.6 Determinamos cinco rangos dividiendo la diferencia máximo mínimo de mean por 5 y establecemos el incremento:

incremento <- (max-min)/5
incremento 
## [1] 22.1964

5.7 Determinamos los límites de los rangos:

el mínimo y primero es:

uno
## [1] 10.548

el segundo es:

dos <- uno + incremento
dos
## [1] 32.7444

el tercero es:

tres <- dos  + incremento
tres
## [1] 54.9408

el cuarto es:

cuatro <- tres  + incremento
cuatro
## [1] 77.1372

el quinto es:

quinto <- cuatro  + incremento
quinto
## [1] 99.3336

y el sexto coincide (como debe ser) con el máximo:

sexto <- quinto  + incremento
sexto
## [1] 121.53

5.8 Determinamos los rangos

construímos el primero rango

rango <- c(1, 2, 3, 4, 5)
uno <- round(uno, digits = 3)
dos <- round(dos, digits = 3)
primer_rango <- paste(uno, "-", dos)
primer_rango 
## [1] "10.548 - 32.744"

construímos el segundo rango

dos <- round(dos, digits = 3)
tres <- round(tres, digits = 3)
segundo_rango <- paste(dos, "-", tres)
segundo_rango 
## [1] "32.744 - 54.941"

construímos el tercero rango

tres <- round(tres, digits = 3)
cuatro <- round(cuatro, digits = 3)
tercer_rango <- paste(tres, "-", cuatro)
tercer_rango 
## [1] "54.941 - 77.137"

construímos el cuarto rango

cuatro <- round(cuatro, digits = 3)
quinto <- round(quinto, digits = 3)
cuarto_rango <- paste(cuatro, "-", quinto)
cuarto_rango 
## [1] "77.137 - 99.334"

construímos el quinto rango

quinto <- round(quinto, digits = 3)
sexto <- round(sexto, digits = 3)
quinto_rango <- paste(quinto, "-", sexto)
quinto_rango 
## [1] "99.334 - 121.53"

lo desplegamos por completo:

rango <- c(1,2,3,4,5)
intervalos <- c(primer_rango, segundo_rango , tercer_rango , cuarto_rango , quinto_rango )
df <- data.frame(rango, intervalos)
print (df)
##   rango      intervalos
## 1     1 10.548 - 32.744
## 2     2 32.744 - 54.941
## 3     3 54.941 - 77.137
## 4     4 77.137 - 99.334
## 5     5 99.334 - 121.53

Añadimos una columna con los rangos:

Convención para determinar los rangos:

df_new$rango <- ifelse(df_new$mean > 10.547 & df_new$mean < 32.744 , 1, ifelse(df_new$mean > 32.745 & df_new$mean < 54.941 , 2, ifelse(df_new$mean > 54.942 & df_new$mean < 77.137 , 3, ifelse(df_new$mean > 77.138 & df_new$mean < 99.333 ,4 , ifelse(df_new$mean > 99.334  & df_new$mean <  121.54,5,"")))))

6 Tabla final

Recapitulemos:

Mientras mayor sea el valor del rango (1-5) mayor el déficit hídrico, el cual debe ser mayor en las comunas del norte de Chile y en particular en Valparaíso. Comunas de la misma región deben tener valores similares:

datatable(df_new, extensions = 'Buttons', escape = FALSE, rownames = TRUE,
          options = list(dom = 'Bfrtip',
          buttons = list('colvis', list(extend = 'collection',
          buttons = list(
          list(extend='copy'),
          list(extend='excel',
            filename = 'tabla.xlsx'),
          list(extend='pdf',
            filename= 'tabla.xlsx')),
          text = 'Download')), scrollX = TRUE))

7 La confirmación

se confirma