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
<- read.csv("C:/Users/chris/OneDrive/Documentos/GitHub/ds_datos_climaticos/comunas_y_subcuencas/COMUNA/MEDIA/TERRACLIMATE_MEDIA_def.csv") def_original
# 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[, c(1:12)]
def_original_11 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
<- readxl::read_xlsx("C:/Users/chris/OneDrive/Documentos/GitHub/ds_datos_climaticos/imputaciones/comuna_medias_imputadas/excel/TERRACLIMATE_MEDIA_def_imp.xlsx") def_imp
# 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[, c(1:12)]
def_imp_11 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:
<- readxl::read_xls("CUT_2018_v04.xls") codigos_comunales_actuales
debemos extraer los ceros iniciales del codigo comuna 2018
library("stringr")
$`Código Comuna 2018` <- str_remove(codigos_comunales_actuales$`Código Comuna 2018`, "^0+") codigos_comunales_actuales
y hacemos la unión:
<- merge(codigos_comunales_actuales,def , all.x = TRUE)
def <- def[,-c(1)] def
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:
<- def[,c(1:6)] COM
Construyamos nuestro dataframe numerico
<- def[,7:258] def
options(scipen=999)
<- function(df, digits) {
round_df <- vapply(df, is.numeric, FUN.VALUE = logical(1))
nums
<- round(df[,nums], digits = digits)
df[,nums]
(df)
}
<- round_df(def, digits=3)
def
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*0.1 def
5.2 Unimos las columnas no numéricas
<- cbind(def, COM) def
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)
<- 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))]
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
$mean <- rowMeans(df_new[,1:126])
df_new
$mean <- round(df_new$mean, digits=3)
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.5 Calculemos los máximos y los mínimos de mean:
Máximo:
<- max(df_new$mean, na.rm = TRUE)
max max
## [1] 121.53
Mínimo:
<- min(df_new$mean, na.rm = TRUE)
uno <- min(df_new$mean, na.rm = TRUE)
min uno
## [1] 10.548
5.6 Determinamos cinco rangos dividiendo la diferencia máximo mínimo de mean por 5 y establecemos el incremento:
<- (max-min)/5
incremento 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:
<- uno + incremento
dos dos
## [1] 32.7444
el tercero es:
<- dos + incremento
tres tres
## [1] 54.9408
el cuarto es:
<- tres + incremento
cuatro cuatro
## [1] 77.1372
el quinto es:
<- cuatro + incremento
quinto quinto
## [1] 99.3336
y el sexto coincide (como debe ser) con el máximo:
<- quinto + incremento
sexto sexto
## [1] 121.53
5.8 Determinamos los rangos
construímos el primero rango
<- c(1, 2, 3, 4, 5)
rango <- round(uno, digits = 3)
uno <- round(dos, digits = 3)
dos <- paste(uno, "-", dos)
primer_rango primer_rango
## [1] "10.548 - 32.744"
construímos el segundo rango
<- round(dos, digits = 3)
dos <- round(tres, digits = 3)
tres <- paste(dos, "-", tres)
segundo_rango segundo_rango
## [1] "32.744 - 54.941"
construímos el tercero rango
<- round(tres, digits = 3)
tres <- round(cuatro, digits = 3)
cuatro <- paste(tres, "-", cuatro)
tercer_rango tercer_rango
## [1] "54.941 - 77.137"
construímos el cuarto rango
<- round(cuatro, digits = 3)
cuatro <- round(quinto, digits = 3)
quinto <- paste(cuatro, "-", quinto)
cuarto_rango cuarto_rango
## [1] "77.137 - 99.334"
construímos el quinto rango
<- round(quinto, digits = 3)
quinto <- round(sexto, digits = 3)
sexto <- paste(quinto, "-", sexto)
quinto_rango quinto_rango
## [1] "99.334 - 121.53"
lo desplegamos por completo:
<- c(1,2,3,4,5)
rango <- c(primer_rango, segundo_rango , tercer_rango , cuarto_rango , quinto_rango )
intervalos <- data.frame(rango, intervalos)
df 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:
$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,""))))) df_new
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
https://developers.google.com/earth-engine/datasets/catalog/IDAHO_EPSCOR_TERRACLIMATE?hl=en#bands↩︎
https://www.analyticsvidhya.com/blog/2016/03/tutorial-powerful-packages-imputing-missing-values/↩︎
http://www.subdere.gov.cl/documentacion/c%C3%B3digos-%C3%BAnicos-territoriales-actualizados-al-06-de-septiembre-2018↩︎