Introducción

En este taller aprenderemos comandos tanto básicos como avanzados para:

  1. Visualizar datos de suelos
  2. Procesar datos
  3. Desarrollar modelos predictivos
  4. Graficar un mapa

Antes de comenzar, debemos descargar varios paquetes (librerías o bibliotecas). Son varias ya que veremos varias herramientas que nos ofrece R para procesar datos de suelo en específico.

install.packages(c("Hmisc","dplyr","soilDB","soiltexture","ggplot2","sf",
                   "caret","rnaturalearth","skimr","cowsay","randomForest","kernlab",
                   "MASS"),
                 repos = "http://cran.us.r-project.org")

Descarga de datos y visualización

Cargaremos datos de suelos chilenos muestreados por la National Cooperative Soil Survey, USDA. El comando bbox de la función `fetchKSSL corresponde al bounding box para Chile, correspondiente a (Longitud mínima, Latitud mínima, Longitud máxima, Latitud máxima). Para mas información de los datos descargados: https://ncsslabdatamart.sc.egov.usda.gov

library(soilDB)
## Warning: package 'soilDB' was built under R version 3.5.2
## Warning: package 'aqp' was built under R version 3.5.2
suelos.chile<-fetchKSSL(bbox=c(-75.64, -55.61, -66.96, -17.58))

Se cargaron 25 pedones muestreados a nivel nacional. El objeto suelos.chile es de tipo SoilProfileCollection (SPC). Con la función plot() los podemos visualizar de forma rápida

plot(suelos.chile)

Los objetos SPC son útiles para visualización. Además es posible transformar datos propios en este tipo de objetos. Para más información de las opciones de los paquetes aqp y soilDB: http://ncss-tech.github.io/AQP/. Ahora graficaremos la arcilla en cada horizonte.

par(mar=c(0,0,4,0)) #ajuste margen
plot(suelos.chile, color='clay', col.label="") #colorearemos según % arcilla
title("Arcilla (%)", line = 3) #título del gráfico

Procesamiento de datos

Primero, transformamos el objeto SPC en un data.frame, lo cual nos permitirá trabajar con el de forma más directa

data.chile<-as(suelos.chile, 'data.frame')

Podemos explorar los datos de forma directa en la consola o en RStudio.

summary(data.chile)
View(data.chile)

Nos enfocaremos en las siguientes variables cuantitativas:

Nombre base de datos Variables
sand, silt, clay Arena, Limo y Arcilla
oc Carbono orgánico
p_nz Fosforo New Zealand
db_od Densidad Aparente
base_sum Suma de bases
ph_h2o pH en agua
w3cld Capacidad de Campo
w15l2 Punto de Marchitez Permanente

Ahora limpiamos la base de datos y vemos correlacion entre variables. Es esta parte usaremos pipes ( %>%, más info: https://r4ds.had.co.nz/pipes.html). Removemos datos repetidos con la función distinct(), luego seleccionamos variables numéricas de interés (select()) y finalmente vemos la correlación estadística entre ellas.

library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.2
library(Hmisc)
## Warning: package 'Hmisc' was built under R version 3.5.2
data.chile %>% 
  distinct(pedon_key,labsampnum,sand,silt,clay, .keep_all= TRUE) %>% 
  select(sand,silt,clay,oc,p_nz,db_od,base_sum,ph_h2o,w3cld,w15l2) %>% #seleccionamos variables
  as.matrix() %>% # los ordenamos como matriz
  rcorr()   # vemos correlacion y significancia, ademas hay muchos datos incompletos
##           sand  silt  clay    oc  p_nz db_od base_sum ph_h2o w3cld w15l2
## sand      1.00 -0.67 -0.69 -0.19 -0.04 -0.15    -0.32   0.02  0.23 -0.27
## silt     -0.67  1.00 -0.08  0.44  0.32 -0.37    -0.03  -0.07  0.06  0.38
## clay     -0.69 -0.08  1.00 -0.15 -0.28  0.54     0.41   0.04 -0.35 -0.01
## oc       -0.19  0.44 -0.15  1.00  0.31 -0.57     0.42  -0.54  0.40  0.72
## p_nz     -0.04  0.32 -0.28  0.31  1.00 -0.80    -0.33  -0.33  0.67  0.75
## db_od    -0.15 -0.37  0.54 -0.57 -0.80  1.00     0.45   0.35 -0.73 -0.60
## base_sum -0.32 -0.03  0.41  0.42 -0.33  0.45     1.00   0.66 -0.21 -0.02
## ph_h2o    0.02 -0.07  0.04 -0.54 -0.33  0.35     0.66   1.00 -0.25 -0.24
## w3cld     0.23  0.06 -0.35  0.40  0.67 -0.73    -0.21  -0.25  1.00  0.58
## w15l2    -0.27  0.38 -0.01  0.72  0.75 -0.60    -0.02  -0.24  0.58  1.00
## 
## n
##          sand silt clay oc p_nz db_od base_sum ph_h2o w3cld w15l2
## sand      180  180  180 75  167   152      133    176   147   174
## silt      180  180  180 75  167   152      133    176   147   174
## clay      180  180  180 75  167   152      133    176   147   174
## oc         75   75   75 75   74    61       74     75    56    75
## p_nz      167  167  167 74  167   139      120    163   134   161
## db_od     152  152  152 61  139   152      111    152   147   148
## base_sum  133  133  133 74  120   111      133    133   106   133
## ph_h2o    176  176  176 75  163   152      133    176   147   170
## w3cld     147  147  147 56  134   147      106    147   147   143
## w15l2     174  174  174 75  161   148      133    170   143   174
## 
## P
##          sand   silt   clay   oc     p_nz   db_od  base_sum ph_h2o w3cld 
## sand            0.0000 0.0000 0.0947 0.6304 0.0719 0.0002   0.8217 0.0056
## silt     0.0000        0.2848 0.0000 0.0000 0.0000 0.7208   0.3788 0.4351
## clay     0.0000 0.2848        0.1901 0.0003 0.0000 0.0000   0.5810 0.0000
## oc       0.0947 0.0000 0.1901        0.0079 0.0000 0.0002   0.0000 0.0020
## p_nz     0.6304 0.0000 0.0003 0.0079        0.0000 0.0002   0.0000 0.0000
## db_od    0.0719 0.0000 0.0000 0.0000 0.0000        0.0000   0.0000 0.0000
## base_sum 0.0002 0.7208 0.0000 0.0002 0.0002 0.0000          0.0000 0.0329
## ph_h2o   0.8217 0.3788 0.5810 0.0000 0.0000 0.0000 0.0000          0.0022
## w3cld    0.0056 0.4351 0.0000 0.0020 0.0000 0.0000 0.0329   0.0022       
## w15l2    0.0003 0.0000 0.8446 0.0000 0.0000 0.0000 0.7876   0.0017 0.0000
##          w15l2 
## sand     0.0003
## silt     0.0000
## clay     0.8446
## oc       0.0000
## p_nz     0.0000
## db_od    0.0000
## base_sum 0.7876
## ph_h2o   0.0017
## w3cld    0.0000
## w15l2

Hay varias opciones para manejar datos ausentes, rellenar con promedios, usar funciones de pedotransferencia, etc. Para este taller, filtraremos el dataset solo con datos completos (sin variables ausentes)

dataset<-data.chile %>% 
  distinct(pedon_key,labsampnum,sand,silt,clay, .keep_all= TRUE) %>% 
  select(sand,silt,clay,oc,p_nz,db_od,base_sum,ph_h2o,w3cld,w15l2) %>%
  filter(complete.cases(.)) # solamente nos quedamos son filas con todos los datos en las columnas

Un paquete que nos ayuda a visualizar datos es skimRr. Revisemos el dataset filtrado en la consola.

library(skimr)
## Warning: package 'skimr' was built under R version 3.5.2
skim_to_wide(dataset)
## # A tibble: 10 x 13
##    type  variable missing complete n     mean  sd    p0    p25   p50  
##    <chr> <chr>    <chr>   <chr>    <chr> <chr> <chr> <chr> <chr> <chr>
##  1 inte… p_nz     0       54       54    91.59 17.01 10    92.5  98   
##  2 nume… base_sum 0       54       54    " 3.… " 3.… " 0 … " 0.… " 2.…
##  3 nume… clay     0       54       54    10.68 21.91 " 0 … " 0.… " 1.…
##  4 nume… db_od    0       54       54    " 0.… " 0.… " 0.… " 0.… " 0.…
##  5 nume… oc       0       54       54    " 4.… " 4.… " 0.… " 0.… " 3.…
##  6 nume… ph_h2o   0       54       54    " 5.… " 0.… " 5 … " 5.… " 5.…
##  7 nume… sand     0       54       54    58.37 25.73 " 3.… "47.… "67.…
##  8 nume… silt     0       54       54    30.95 15.82 " 1.… "19.… "30.…
##  9 nume… w15l2    0       54       54    "27.… 10.88 " 2 … 21.92 28.05
## 10 nume… w3cld    0       54       54    84.97 44.81 "16.… "46.… 75.65
## # … with 3 more variables: p75 <chr>, p100 <chr>, hist <chr>

Siguiendo con la visualización, podemos hacer un triángulo de textura con los datos filtrados

library(soiltexture)
TT.plot(class.sys = "USDA.TT",
        tri.data = dataset,
        css.names = c("clay","silt","sand"),
        cex.axis=1,cex.lab=1)

Modelos predictivos

Ahora avanzaremos con el código para modelos predictivos. Por tiempo, veremos dos ejemplos: a) Un modelo de regresión lineal (LM), con un proceso de selección de variables (step-wise), b) Máquina de vectores de soporte (Support Vector Machine) y c) Árbol aleatorio (Random Forest). Sin embargo, la sintáxis general puede ser replicada para modelos de regresión utilizando redes neuronales, árboles de decisión, etc. Para mas detalles del paquete carety sus funcionalidades ver: http://topepo.github.io/caret/. Esta tabla resume a grandes rasgos características de cada modelo. Sin embargo, el paquete caret permite hacer ajustes y preprocesamiento para optimizar el uso de estos modelos.

Características LM SVM RF
Parámetros de ajuste NO SI (C y sigma) SI (mtry)
Preprocesamiento NO SI NO
Seleccion variables SI (Step-wise) NO NO

Paso previo: Partición de la base de datos

Vamos a predecir los valores de w15l2, punto de marchitez permanente, en función del resto. Primero, dividiremos los datos del set en entrenamiento (80%,p=0.8) y testeo (20%). Luego seleccionaremos las variables para el testeo.

library(caret)
set.seed(1234)
dataset.part <- createDataPartition(dataset$w15l2, p = 0.8, list = FALSE)
training <- dataset[dataset.part,] #dataset para entrenar
testing <- dataset[-dataset.part,] #dataset para testear
x.testing <- testing[,1:9] # variables input para el testeo
y.testing <- testing[,10] # w15l2 en el set de testeo

Paso previo: Método de entrenamiento

El comando train nos permite buscar valores óptimos de los parámetros de ajuste, usando validación cruzada method="cv" u otras técnicas de remuestreo. Además, permite comparar modelos.

train.control<-trainControl(method="repeatedcv", number=10, repeats =5) #validación cruzada 10-fold

Desarrollo de los modelos

set.seed(1234)
lm.model<-train(w15l2 ~., data = training, #se selecciona los datos de training
                      method = "lmStepAIC", # método LM con Step-wise
                      trControl = train.control) #train.control definido anteriormente
set.seed(1234)
svm.model<-train(w15l2 ~., data = training,
                    method = "svmRadial", 
                    preProcess = c("center", "scale"), #paso preprocesamiento
                    trControl = train.control)
set.seed(1234)
rf.model<-train(w15l2 ~., data = training,
                    method = "rf", 
                    trControl = train.control)

Información de los modelos

Una vez que los modelos están listos, podemos sacar información detallada de su ajuste, parámetros, etc. Veremos un ejemplo para cada modelo.

Para el modelo lm.model podemos ver la ecuación y las variables seleccionadas

lm.model$finalModel
## 
## Call:
## lm(formula = .outcome ~ sand + silt + oc + p_nz + ph_h2o, data = dat)
## 
## Coefficients:
## (Intercept)         sand         silt           oc         p_nz  
##    -34.5352      -0.1584      -0.1162       1.7611       0.3032  
##      ph_h2o  
##      6.5767

Para el modelo svm.model podemos ver la importancia de las variables en la predicción

varImp(svm.model)
## loess r-squared variable importance
## 
##          Overall
## oc        100.00
## ph_h2o     65.53
## db_od      57.11
## w3cld      55.46
## p_nz       45.76
## silt       45.64
## base_sum   42.42
## sand       41.71
## clay        0.00

Para el modelo rf.model podemos ver como ajustaron los parámetros por defecto. Esto también se puede realizar de forma manual, lo cual es útil para optimizar los modelos.

plot(rf.model)

Ya están listos nuestros modelos predictivos, ahora podemos comparar los resultados entre ellos y evaluar su performance.

modelos<-resamples(list(LM = lm.model, SVM = svm.model, RF = rf.model))
summary(modelos) # resumen del cross-validation y resampleo
## 
## Call:
## summary.resamples(object = modelos)
## 
## Models: LM, SVM, RF 
## Number of resamples: 50 
## 
## MAE 
##         Min.  1st Qu.   Median     Mean  3rd Qu.     Max. NA's
## LM  1.399175 4.634382 6.543970 6.501128 8.173876 12.61716    0
## SVM 2.324841 3.987246 5.184039 5.952885 7.166975 14.09890    0
## RF  2.219670 4.030877 5.045302 5.633944 6.556685 13.92811    0
## 
## RMSE 
##         Min.  1st Qu.   Median     Mean  3rd Qu.     Max. NA's
## LM  1.802279 5.949179 7.897221 7.989510 9.204648 18.11070    0
## SVM 2.814433 5.091970 6.816799 7.910893 9.109436 22.74019    0
## RF  2.560375 5.361428 6.374673 7.470620 8.049140 20.72977    0
## 
## Rsquared 
##             Min.    1st Qu.    Median      Mean   3rd Qu.      Max. NA's
## LM  0.0003799249 0.08860291 0.5304505 0.4784419 0.8196291 0.9999250    0
## SVM 0.0014103830 0.15909263 0.5455629 0.4965747 0.7878340 0.9899102    0
## RF  0.0006877534 0.11645334 0.5322692 0.5134676 0.8516935 0.9999192    0
compare_models(lm.model,rf.model) #vemos si hay diferencias sig. entre dos modelos
## 
##  One Sample t-test
## 
## data:  x
## t = 1.733, df = 49, p-value = 0.08939
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.08282111  1.12060077
## sample estimates:
## mean of x 
## 0.5188898
bwplot(modelos, scales = list(relation = "free")) #graficamos estadistica para cada uno

Evaluación en el set de testeo (20%)

La función predict usa el modelo predictivo desarrollado y utiliza las variables x.testing para hacer una estimación

lm.pred<-predict(lm.model,x.testing) #hacemos la prediccion del 20%  
svm.pred<-predict(svm.model,x.testing)
rf.pred<-predict(rf.model,x.testing)

Ahora guardaremos los resultados de la predicción en el restante 20% y los almacenaremos ordenados como un data.frameen formato largo para graficar

lm.prediccion<-as.data.frame(cbind(modelo = rep("LM", nrow(x.testing)),
                                     predicho = lm.pred,
                                     observado = y.testing))
                             
svm.prediccion<-as.data.frame(cbind(modelo = rep("SVM", nrow(x.testing)),
                                     predicho = svm.pred,
                                     observado = y.testing))

rf.prediccion<-as.data.frame(cbind(modelo = rep("RF", nrow(x.testing)),
                                     predicho = rf.pred,
                                     observado = y.testing))

prediccion.modelos<-rbind(lm.prediccion,svm.prediccion,rf.prediccion)#unimos la informacion en un solo dataframe

prediccion.modelos$predicho <- as.numeric(as.character(prediccion.modelos$predicho))# nos aseguramos que sean valores numéricos
prediccion.modelos$observado <- as.numeric(as.character(prediccion.modelos$observado))

Gráficos para el set de testeo

Ahora usando el paquete ggplot2 graficaremos los resultados para el set de testeo. Hay un error ya que un punto no está en el rango del gráfico para el modelo LM.

library(ggplot2)
ggplot(prediccion.modelos, aes(x=observado, y=predicho, color=modelo, shape=modelo)) +
  geom_point() +  # grafico puntos scatter plot
  geom_abline()+ # agrego linea 1:1
  geom_smooth(method=lm)+  #agrego una regresion lineal con interval de confianza
  ggtitle("Comparación métodos") +
  xlim(0,50)+
  ylim(0,50)+
  ylab("% gravimetrico 1500 kPa predicho") +
  xlab("% gravimetrico 1500 kPa observado")

Vemos que todos los modelos subestiman de forme importante los valores de PMP sobre 40%. RF es el que muestra mejor performance en general.

Mapas en R

R tiene varias librerias para gráficas mapas. A continuación mostramos un ejemplo rápido para crear un mapa simple con los datos georeferenciados obtenidos al inicio del taller.Para mas info de mapas en R con ggplot:https://www.r-spatial.org/r/2018/10/25/ggplot2-sf.html

library(sf)
library(rnaturalearth)
mundo <- ne_countries(scale = "medium", returnclass = "sf") #cargo un mapa del mundo
ggplot(data = mundo) + #mapa tipo ggplot
  geom_sf()  + #datos geográficos
  geom_point(data = data.chile, aes(x = x, y = y), size = 1)+ #cargamos los puntos
  coord_sf(xlim = c(-75.64, -66.96), # zoom al mapa igual al bounding box inicial
           ylim = c(-55.61, -17.58), expand = FALSE)

Gracias por asistir! Cualquier consulta, pueden consultar al email: biofisica@ing.puc.cl

library(cowsay)
say("Gracias por venir \n Sígannos en @soilbiophysicslab", by="cat")
## Colors cannot be applied in this environment :( Try using a terminal or RStudio.
## 
##  -------------- 
## Gracias por venir 
##  Sígannos en @soilbiophysicslab 
##  --------------
##     \
##       \
##         \
##             |\___/|
##           ==) ^Y^ (==
##             \  ^  /
##              )=*=(
##             /     \
##             |     |
##            /| | | |\
##            \| | |_|/\
##       jgs  //_// ___/
##                \_)
##