Utilizaremos el algoritmo Random Forest para determinar variables relevantes obtenidas de un análisis de regresión entrenado asimismo con RF.

1 Cargamos nuestras librerías:

library("readxl")
library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
library(reshape)

2 Leemos la data y desplegamos las 3 primeras filas:

my_data <- read_excel("data_e.xlsx", sheet = 2, col_names = TRUE, col_types = NULL, na = "", skip = 0)
head(my_data, 3)
## # A tibble: 3 x 32
##   SOC_TCH SOC_Modelado Aspect    BSI   EVI Ecosystm Elevation    IBI Litologia
##     <dbl>        <dbl>  <dbl>  <dbl> <dbl>    <dbl>     <dbl>  <dbl>     <dbl>
## 1    395.         310.   328. NA     NA          NA      2.25 NA            NA
## 2    120.         295.   354. NA     NA          NA   1171.   NA            NA
## 3    548.         438.   302. -0.272  2.41        0      7.83 -0.251         0
## # ... with 23 more variables: MNDVI <dbl>, MeanCurvature <dbl>, NDFI <dbl>,
## #   NDSI <dbl>, NDVI <dbl>, NDWBI <dbl>, NDWI <dbl>, NIRv <dbl>,
## #   Precipitacionmean <dbl>, SAVI <dbl>, Slope <dbl>, Suelo_Map <dbl>,
## #   Temp_Cmean <dbl>, VH <dbl>, VV <dbl>, claymean <dbl>, nitrogenmean <dbl>,
## #   ocdmean <dbl>, ocsmean <dbl>, phh2omean <dbl>, sandmean <dbl>,
## #   siltmean <dbl>, socmean <dbl>

3 Excluímos del análisis la columna SOC_Modelado y asignamos NA a los valores 0:

my_data <- my_data[,-(2)] 
my_data[my_data == 0] <- NA
my_data <- my_data[,-(11:14)] 
head(my_data, 10)
## # A tibble: 10 x 27
##    SOC_TCH Aspect    BSI   EVI Ecosystm Elevation    IBI Litologia  MNDVI
##      <dbl>  <dbl>  <dbl> <dbl>    <dbl>     <dbl>  <dbl>     <dbl>  <dbl>
##  1   395.   328.  NA     NA          NA      2.25 NA            NA NA    
##  2   120.   354.  NA     NA          NA   1171.   NA            NA NA    
##  3   548.   302.  -0.272  2.41       NA      7.83 -0.251        NA  0.652
##  4    57.2  353.  -0.367  2.62       NA      1.04 -0.326        NA  0.731
##  5   226.    69.9 -0.534  2.56       NA     11.7  -0.498        NA  0.864
##  6   448.   160.  -0.464  2.64       NA     17.1  -0.441        30  0.809
##  7   437.    NA   -0.492  2.32      398      3.46 -0.473        50  0.850
##  8   196.   205.  -0.444  2.43       NA      5.34 -0.471        NA  0.790
##  9   245.   214.  -0.410  2.22       NA      4.82 -0.414        52  0.781
## 10   500.   230.  -0.488  2.38       NA      4.52 -0.491         5  0.815
## # ... with 18 more variables: MeanCurvature <dbl>, NDWI <dbl>, NIRv <dbl>,
## #   Precipitacionmean <dbl>, SAVI <dbl>, Slope <dbl>, Suelo_Map <dbl>,
## #   Temp_Cmean <dbl>, VH <dbl>, VV <dbl>, claymean <dbl>, nitrogenmean <dbl>,
## #   ocdmean <dbl>, ocsmean <dbl>, phh2omean <dbl>, sandmean <dbl>,
## #   siltmean <dbl>, socmean <dbl>

4 El modelo

Creamos el objeto modelo con la función randomForest, siendo nuestra variable dependiente SOC_TCH y el resto las independientes.

Parámetros:

ntree nos permite especificar el número de árboles a realizar.

importance incluye en el objeto las “medidas de importancia”.

maxnodes indica el número máximo de nodos en nuestros árboles.

mtry indica el número máximo de variables en los modelos creados.

Para medir la importancia de las variables empleamos la función importance sobre el modelo creado. Esta función despliega dos columnas:

  1. el error cuadrático medio (%IncMSE) y

  2. la pureza del nodo (IncNodePurity)

modelo <- randomForest(SOC_TCH~., data = my_data,
              na.action = na.exclude, ntree = 500, importance = TRUE, maxnodes = 10, mtry = 20)

importancia = data.frame(importance(modelo))

importancia <- importancia[order(-importancia$X.IncMSE),] 
importancia
##                     X.IncMSE IncNodePurity
## Temp_Cmean         8.7674497     193676.56
## Suelo_Map          6.6064203     107035.32
## ocsmean            6.4017080     288936.17
## phh2omean          4.9833435      76824.34
## Slope              4.1426547     233423.43
## NDWI               3.5039528      76179.89
## Elevation          3.4775083      62754.79
## VV                 2.8973895      52420.60
## Aspect             2.6850337     169751.16
## IBI                2.6554081      50065.34
## SAVI               2.5556845      42414.61
## MNDVI              2.2320940      49930.40
## MeanCurvature      2.0648549     189900.68
## Precipitacionmean  1.9344040      54718.32
## nitrogenmean       1.8889987      50916.54
## siltmean           1.5865883      66796.83
## BSI                1.5049656      78008.55
## Ecosystm           1.1189054      36758.15
## socmean            0.9138285      47624.19
## ocdmean            0.9030106      38223.31
## sandmean           0.6922283      34036.76
## claymean           0.6588370      94749.56
## VH                 0.1062205      36122.25
## Litologia         -0.5319437      23187.75
## NIRv              -1.1126681      38007.67
## EVI               -1.8058776      34055.72
varImpPlot(modelo)

5 Gráfica del error cuadrático medio.

dotchart((importancia$X.IncMSE), xlim=c(-2,11), xlab="Error cuadratico medio", 
         labels = row.names(importancia),
         cex = 0.7)

6 Análisis del modelo.

El error cuadrático medio es la medida más sólida e informativa. Cuanto mayor es, más también lo es la variable a la que hace referencia.

IncNodePurity se relaciona con la función de pérdida elegida según las mejores divisiones. Las variables más útiles logran mayores aumentos en la pureza de los nodos, es decir, encuentran una división que tiene una alta “varianza” entre nodos y una pequeña “varianza” entre nodos. IncNodePurity está sesgado y sólo debe usarse si el tiempo de cálculo adicional para calcular el % IncMSE es inaceptable.

7 Conclusión.

Las variables más relevantes en el análisis son:

Temp_Cmean , Suelo_Map y ocsmean.