Arbres de regressió i de classificació (CART)

Notes d’estadística aplicada a les ciències ambientals: Nota Núm. 7

library(tibble)
library(dplyr)
library(ggplot2)

Introducció

Els Mètodes aquí estudiats solen venir englobats habitualment en el que s’anomena Mineria de Dades o Data Mining perquè solen utilitzar-se amb dades de gran dimensió (número p de variables molt alt) i/o enorme mida mostral (\(n\) molt gran) i, en ocasions, amb \(p >> n\) el que crea greus problemes d’aplicació. Una qüestió dʻinterès que ens agradaria ressaltar és que, en contra del que comúnment es creu, aquest tipus de tècniques solen ser poc robustes, és a dir, solen ser sensibles a la presència de dades anòmales a la mostra.

Referències

• L. Cayuela, M. De la Cruz, Análisis de datos ecológicos en R, 2022.
• A. García Pérez, Cuadernos de estadística aplicada: biología y ciencias ambientales, 2023.
• A. García pérez, Cuadernos de estadística aplcada: área de la salud, 2021.

Base teòrica

Els Arbres de Classificació i Regressió (Classification and Regression Trees), habitualment coneguts pel seu acrònim anglosaxó CARTs, són una tècnica consistent a descobrir relacions (condicionals) entre un gran nombre de covariables explicatives \(X_i\) i una dependent \(Y\), qualitativa o contínua.

S’anomenen Arbres de Classificació quan s’apliquen a variables dependents qualitatives i Arbres de Regressió quan la variable dependent és de tipus continu.

Ambdues tècniques, degudes a Breiman et al. (1993), suposen l’aplicació d’un algorisme que va dividint el conjunt d’individus de la mostra en subgrups, de manera que es minimitzi l’heterogeneïtat (anomenada impuresa del node) dins dels nous grups formats.

Aquesta tècnica ordena per importància a les covariables \(X_i\), la qual cosa implica que podem quedar-nos només unes quantes covariables explicatives abans de fer una regressió de models GLM o GAM.

Construcció de l’arbre

Les dades de les quals disposarem seran observacions d’una variable dependent \(Y\) i \(p\) variables independents \(X_1, ..., X_p\) les quals pensem serveixen per predir la variable dependent.

Un Arbre es construeix determinant primer la variable \(X_j\) més predictiva de \(Y\) en el sentit que veurem més avall. Suposem per començar que aquesta variable \(X_j\) prengués només dos valors; els individus de la mostra, els quals inicialment estan tots en un conjunt, anomenat node arrel, o node pare i que representarem per \(Ω\), es dividiran en dos subconjunts o nodes filla, \(Ω_1\) i \(Ω_2\) segons els valors d’aquesta variable \(X_j\) . Si prengués més valors –per exemple, fora de tipus continu–, els dos grups es formarien depenent de si \(X_j < c\), \(X_j ≥ c\), essent algun valor possible de \(X_j\).

S’escull a continuació la segona variable més predictiva de \(Y\) en cadascun dels nodes filla i s’aplica de nou una regla similar a cadascun dels dos nodes filla; i així se segueix particionant la mostra fins a un determinat moment fixat per una regla de parada (per exemple que el node tingui menys de tres individus). Advertim que el mètode pot conduir a arbres asimètrics.

Quan hàgim construït l’Arbre haurem seleccionat unes quantes covariables, les més influents en la variable dependent, i a més per ordre d’importància. L’elecció de la variable més predictiva i la regla de classificació a partir d’ella, es basa en el que s’anomena impuresa del node (és a dir, la seva heterogeneïtat) \(I(Ω)\) per a la qual hi ha diverses opcions; no obstant, sol dir-se que totes condueixen bàsicament al mateix arbre. La variable predictiva (i la regla de classificació basada en ella) s’escull com aquella que maximitzi \[ I(\Omega)-I(\Omega_D)-I(\Omega_I) \]

sent \(Ω_D\) i \(Ω_I\) els dos nodes filla (Dreta i Esquerra) del node \(Ω\) obtinguts després d’aplicar la regla considerada.

Variable dependent dicotàmica

De les diverses mesures d’impuresa del node, la més habitual per al cas que la variable dependent \(Y\) sigui dicotàmica, és l’Índex de Gini definit per \(I(Ω) = 2pΩ (1−pΩ)\) on \(pΩ=P(Y = 1|Ω)\) és la probabilitat que \(Y\) sigui 1, condicionat a pertànyer al node \(Ω\).

La llibreria rpart

Una de les eines més utilitzades per construir arbres de classificació i regressió és rpart.

Presentació de l’exemple

En aquest capítol farem servir el data frame defo de la llibreria ADER que recull dades de la defoliació de pinars a la Sierra de Filabres. La taula de dades recull informació de pinars de 2 espècie de pi, Pinus nigra i Pinus sylvestris, on s’ha estimat visualment el percentatge de defoliació, així com altres variables: l’elevació, la pendent, l’orientació, la insolació i el potencial hídric:

library(ADER)
data('defo')
str(defo)

## 'data.frame':    76 obs. of  12 variables:
##  $ x                : num  543726 543486 542903 535211 539086 ...
##  $ y                : num  4124164 4124224 4122417 4118071 4121908 ...
##  $ Especie          : Factor w/ 2 levels "Pinus nigra",..: 1 1 1 1 1 1 1 2 1 1 ...
##  $ Defoliacion      : int  55 55 30 55 30 55 55 60 45 55 ...
##  $ Area_basimetrica : num  25.8 28.1 24.2 22 18 ...
##  $ Altura_media     : num  5.2 6.5 6.4 8 8.2 4 6.8 7 5.5 3.8 ...
##  $ Densidad_pinos   : int  286 286 1592 923 1974 223 103 1293 1401 1496 ...
##  $ Elevacion        : num  1807 1821 1877 1909 1688 ...
##  $ Pendiente        : num  17.7 16.5 18.2 12.1 14.8 ...
##  $ Orientacion      : num  220.2 189 248.8 174.7 98.4 ...
##  $ Insolacion       : num  8526 8560 8480 8680 8481 ...
##  $ Potencial_hidrico: num  6.7 5.01 3.08 6.66 7.8 ...

En aquest exemple interessa saber com influeixen les diferents variables biofísiques sobre la defoliació de les masses forestals. Incloem com a variable explicativa l’espècie de pi, la densitat dels pins, l’elevació, l’orientació de la parcel·la, la insolació i el potencial hídric.

Ajustament de del model

Amb la funció rpart() l’ajustament es duu a terme amb una fórmula semblant a la del models lineals, amb l’excepció que no es poden especificar interaccions entre les variables explicatives.

set.seed(2)
rpart.def <- rpart::rpart(Defoliacion ~ Especie+Densidad_pinos+Elevacion+Orientacion+
               Insolacion+Potencial_hidrico, data = defo)

La interpretació dels resultats de rpart

rpart.def

## n= 76 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 76 30100.990 45.98684  
##    2) Elevacion>=1955.836 8     0.000 20.00000 *
##    3) Elevacion< 1955.836 68 24062.870 49.04412  
##      6) Orientacion< 153.2072 40 13294.380 42.12500  
##       12) Elevacion< 1661.86 7  1421.429 25.71429 *
##       13) Elevacion>=1661.86 33  9587.879 45.60606  
##         26) Orientacion>=73.70248 20  6655.000 41.50000  
##           52) Insolacion>=8397.123 11  2772.727 35.45455 *
##           53) Insolacion< 8397.123 9  2988.889 48.88889 *
##         27) Orientacion< 73.70248 13  2076.923 51.92308 *
##      7) Orientacion>=153.2072 28  6117.857 58.92857  
##       14) Especie=Pinus nigra 15  1573.333 48.66667 *
##       15) Especie=Pinus sylvestris 13  1142.308 70.76923 *

Cada fila representa una branca de l’arbre (split).

• El criteri que separa una branca de la resta apareix en primer lloc: per exemple, a la línia 2, Elevacion>=1955.836
• seguit pel nombre de casos inclosos en aquesta branca, 8.
• Seguidament ve la quantitat d’error residual que queda per explicar (la suma dels quadrats de l’error de les observacions al nus), 0.000.
• finalment, s’indica el valor estimat de la variable de resposta per a aquell nus, és a dir 20.0000 (és el % de defoliació). Si es tracta d’un nus terminal (fulla), això s’indica amb asterisc.

Per exemple, la primera partició es fa tenint en compte que l’elevació sigui \(\geq 1955.836\), amb 8 de 76 observacions complint aquest criteri, i veiem que el % de defoliació per a aquest nus serà del 20%, mentre que, de la variabilitat total de 30100 en la variable de resposta, resta sense explicar 0+24062 (veure branca de l’esquerra i de la dreta). La branca de elevació \(<1955.836\) seguirà bifurcant-se múltiples cops.

Representació gràfica

Tot i que es pot utilitzar la funció plot amb rpart, utilitzarem la funció rpart.plot() del paquet rpart.plot, ja que té més qualitat.

rpart.plot::rpart.plot(rpart.def, type = 2)

La funció fa un ombrejat proporcional al valor de la predicció (més alt, més fosc). Cal notar es fan servir només 4 de les variables originals per construir l’arbre

Avaluació de l’ajustament

Els CART són mètodes no paramètrics on no es fan supòsits específics sobre l’estructura del model i en particular sobre l’estructura dels errors. Tot i això, cal tenir en compte que els arbres són sensibles als valors atípics i a l’heteroscedasticitat, ja que es basen en mitjanes locals (estimacions de la mitjana d’una variable calculades utilitzant únicament un subconjunt de dades properes o veïnes, en lloc de tota la mostra).

Per tant és sempre recomanable representar el gràfic de residuals / valors esperats i el gràfic Q-Q per identificar valors atípics. Si la variança no fos constant, es podria considerar transformar la resposta.

q1 <- ggplot(data = data.frame(x=predict(rpart.def), y=resid(rpart.def)),
       aes(x,y)) +
  geom_point() +
  labs(x = "Valors esperats", 
       y = "Valors residuals",
       title = "Residuals vs Fitted") +
  geom_hline(yintercept = 0, colour= "red") +
  theme_bw()
q2 <- qplot(sample = resid(rpart.def)) + 
  stat_qq_line(colour= "red") +
  labs(title = "Q-Q Residuals",
       x= "Quantils teòrics",
       y= "Valors residuals") +
  theme_bw()

library(patchwork)
q1 + q2

Si la variança no fos constant es podria considerar transformar la variable de resposta, encara que aquí no sembla ser el cas.

Validació creuada i poda

Encara que els CARTs no assumeixen cap model de partida, quan ajustem l’algoritme obtenim un model final. Tot i així, ens podem preguntar si aquest realment representa el millor model final. Per determinar quina és la grandària òptima de l’arbre que construïm, podem seguir un procés de validació creuada.

La funció rpart() no només ajusta l’arbre, sinó que a més a més du a terme la validació creuada i poda.

Resulta que, a més de fer créixer l’arbre, entre bastidors rpart:

• ja ha utilitzat la “poda de complexitat de costos” (cost complexity pruning) per obtenir la seqüència anidada d’arbres
• la ha aplicat la validació creuada de 10 vegades per calcular les estimacions de validació creuada i els errors estàndard per a cada valor de α.

Només hem de executar la funció printcp() per obtenir un resum de tota aquesta informació:

rpart::printcp(rpart.def)

## 
## Regression tree:
## rpart::rpart(formula = Defoliacion ~ Especie + Densidad_pinos + 
##     Elevacion + Orientacion + Insolacion + Potencial_hidrico, 
##     data = defo)
## 
## Variables actually used in tree construction:
## [1] Elevacion   Especie     Insolacion  Orientacion
## 
## Root node error: 30101/76 = 396.07
## 
## n= 76 
## 
##         CP nsplit rel error  xerror    xstd
## 1 0.200595      0   1.00000 1.02876 0.10561
## 2 0.154501      1   0.79940 0.97043 0.12775
## 3 0.113027      2   0.64490 0.95057 0.12000
## 4 0.075913      3   0.53188 0.85548 0.11688
## 5 0.029058      4   0.45596 0.76829 0.11038
## 6 0.010000      6   0.39785 0.77829 0.11781

Centrem-nos en la taula a dalt. Cada fila correspon a un arbre de la seqüència obtinguda per poda.

Analitzem cada columna per parts:

• La columna CP és el “paràmetre de complexitat”. Està relacionat amb el paràmetre \(α\) del cost de complexitat (\(\small CC=R(T)+\alpha|T|\)). Compte que la terminologia “paràmetre de complexitat” és una mica enganyosa perquè els paràmetres de complexitat més alta corresponen a models menys complexos (igual que lambda en la regressió penalitzada).
• nsplit és el nombre de particions de l’arbre. Tingueu en compte que 1+nsplit és el nombre de nodes terminals de l’arbre.
• rel error és l’error d’entrenament RSS de l’arbre, normalitzat per la variància total de la resposta; equivalentment, és a dir \(1 − R^2\), com en regressió lineal. L’error d’entrenament disminueix a mesura que augmenta la complexitat. Si volguéssim saber l’error absolut, hauríem de multiplicar l’error relatiu per la variança del nus arrel que hem vist que és 30100.
• xerror és l’estimació de l’error de validació creuada, l’error mitjà de prediccions (MSE), en la mateixa escala que l’error relatiu
• xstd és l’error estàndard de validació creuada.

Si representem gràficament podrem visualitzar com es realitza la poda: se selecciona com a màxim de particions la primera que queda per sota de la suma entre MSE (xerror) mínim i 1 corresponent desviació estàndard (xstd): en aquest cas, 0.7682870 +0.1103756= 0.8786626, que és on es col·locarà la línia horitzontal del gràfic: és a dir a la tercera partició.

par(mar=c(5,5,7,2))
rpart::plotcp(rpart.def, upper = "splits")
title("Validació creuada de l'arbre de regressió", line = 5, adj=0, cex.main=1.2)

Important entendre la diferència entre els valors de CP de la taula i els de CP del gràfic: la taula aporta els valors mínims, allà per on s’ha de fer la poda; el gràfic dona les mitjanes geomètriques dels intervals de valors. Per exemple, a la partició 3 on s’ha de fer la poda, el valor de poda serà \(0.075913\) (el mínim de l’interval), mentre el valor del gràfic serà la mitjana geomètrica de l’interval entre les particions 2 i 3: \(\small\sqrt{0.113027·0.075913}\approx0.093\)

Arbre podat

Triarem per tant l’arbre amb tres particions, podant l’arbre amb el valor de poda que hem vist.

poda.rpart.def <- rpart::prune.rpart(rpart.def, cp=0.075913)

rpart.plot::rpart.plot(poda.rpart.def,
                       main= "Arbre de regressió de dades de defoliació i podat a un valor de cp=0.0759", cex.main=1)

Arbre de classificació

Els arbre de classificació s’ajusten i es poden de la mateixa manera. Per exemple podríem posar com a variable resposta, en compte d’una amb valors continus, una amb valors categòrics, amb nivells de defoliació. Establint tres nivells de defoliació: Baixa (\(\small<25\)), Mitjana (\(\small\geq25<50\)), Alta (\(\small\geq50\leq10\)) de defoliació, generarem una variable categòrica:

defo$Defo.cat <- cut(defo$Defoliacion, 
                     breaks = c(0,25,50,100),
                     labels = c("Baixa","Mitjana", "Alta"))

R, d’aquesta manera genera una variable factorial de 3 nivells i, automàticament, rpart() ajusta l’arbre de classificació.

set.seed(2)
rpart.cat <- rpart::rpart(Defo.cat ~ Especie+Densidad_pinos+Elevacion+Orientacion+
               Insolacion+Potencial_hidrico, data = defo)

rpart.plot::rpart.plot(rpart.cat,
                       main= "Arbre de classificació de dades de defoliació per categories", cex.main=1)

Aquí els colors indiquen la pertinència a una de les 3 categories de la variables de resposta.

A la línia central de cada fulla s’indiquen les proporcions de com s’agrupen per categories les observacions reunides en aquest nus. Com hi ha tres categories, hi ha 3 valors que sempre han de sumar 1.

Parèntesi sobre la intensitat dels colors

Amb rpart.plot, la intensitat del color de la caixa d’un node representa la puresa o la confiança de la resposta categòrica predita. Els colors més foscs i saturats indiquen una confiança predictiva més alta (alta puresa), mentre que els colors més pàl·lids i esvaïts indiquen que un node està més barrejat o més a prop d’una divisió 50/50 (menor puresa).

És a dir:

1. Els tons de color representen la classe predita. El color específic (per exemple, verd, gris, taronja) s’assigna en funció de la classe que el model prediu per a aquest node específic (per exemple, “Baixa”, “Mitjana”, o “Alta”). Aquests colors assignats corresponen als nivells de factor de la variable de resposta.
2. La intensitat del color representa la “puresa” del node. La intensitat (foscor o brillantor) d’aquest to és directament proporcional a la probabilitat de la classe predita en aquest node.

• Colors foscos/brillants: La classe majoritària constitueix un percentatge molt gran de les observacions en aquest node (per exemple, \(\small\geq 85\%\)). El node és altament “pur” i la predicció del model és molt fiable.
• Colors tènues/pàl·lids: Les observacions del node són una barreja (per exemple, una divisió de \(\small55\%/45\%\) de dues classes). El node és “impur” i la predicció del model és menys fiable.

Exemple: pesticides a un riu d’Oregon

El conjunt de dades de contaminació de Willamette d’Anderson de 1997, Distribution of Dissolved Pesticides and Other Water Quality Constituents in Small Streams, and their Relation to Land Use, in the Willamette River Basin, Oregon, 1996 és un estudi històric de la qualitat de l’aigua publicat pel Servei Geològic dels Estats Units (USGS). Proporciona nivells de concentració de pesticides agrícoles i urbans, inclòs el diuró (codi de paràmetre USGS P49300 assignat), en petits rierols de la conca del riu Willamette, Oregon.

En l’estudi de 95 mostres per detectar 36 pesticides (29 herbicides i 7 insecticides) i diverses covariables potencialment predictores (veure la matriu de dades Willamette.csv) com a explicatives de la variable dependent concentració de pesticides P49300 (=diuron).

pesticides <- readRDS(url("https://analysis.cat/r-data/pesticides.rds"))
str(pesticides)

## 'data.frame':    95 obs. of  70 variables:
##  $ X        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ STATION  : num  14190970 14190970 14190970 14190970 14190970 ...
##  $ DATE     : int  19960418 19960515 19960723 19961014 19961114 19960419 19960514 19960725 19961019 19961118 ...
##  $ TIME     : int  1530 1330 1300 1740 1320 1430 1120 1220 1300 1010 ...
##  $ R        : chr  NA NA "<" NA ...
##  $ NH4      : num  NA NA 0.02 0.114 0.01 NA 0.03 0.036 0.027 0.05 ...
##  $ R.1      : chr  NA NA "<" "E" ...
##  $ NO2      : num  NA NA 0.005 0.01 0.01 NA NA 0.005 0.005 NA ...
##  $ R.2      : chr  NA NA NA NA ...
##  $ TKN      : num  NA NA 0.492 1.19 0.256 NA NA 0.283 0.291 0.4 ...
##  $ R.3      : chr  NA NA NA NA ...
##  $ N2.3     : num  1.1 NA 0.617 0.329 0.68 1.2 0.82 0.368 0.527 1.1 ...
##  $ R.4      : chr  NA NA NA NA ...
##  $ TOTP     : num  0.05 NA 0.092 0.273 0.051 0.05 0.07 0.137 0.086 0.05 ...
##  $ R.5      : chr  NA NA "<" NA ...
##  $ SRP      : num  NA NA 0.01 0.058 0.01 NA 0.056 0.047 0.025 0.025 ...
##  $ AGY      : chr  "USGS" NA "USA" "USA" ...
##  $ R.6      : chr  NA NA "<" NA ...
##  $ BOD      : num  NA NA 2 8.5 2 NA 0.7 2 2 NA ...
##  $ AGY.1    : chr  NA NA "POR" "POR" ...
##  $ R.7      : chr  NA NA NA NA ...
##  $ ECOL     : int  NA NA 4300 1000 NA NA 385 1300 NA 150 ...
##  $ R.8      : chr  NA NA NA NA ...
##  $ FECAL    : int  NA NA 5000 7700 NA NA 460 2600 NA 190 ...
##  $ AGY.2    : chr  NA NA "EUG" "EUG" ...
##  $ R.9      : logi  NA NA NA NA NA NA ...
##  $ wtemp    : num  12.2 15.2 23.5 13.9 11.9 ...
##  $ R.10     : logi  NA NA NA NA NA NA ...
##  $ bpres    : num  757 760 759 761 765 762 NA NA 759 751 ...
##  $ R.11     : chr  NA NA NA NA ...
##  $ flow     : num  38.3 16.9 11.2 11.8 24.2 ...
##  $ R.12     : chr  NA NA NA NA ...
##  $ cond     : num  70 71.1 75.5 67.7 87 ...
##  $ R.13     : chr  NA NA NA NA ...
##  $ DO       : num  10.3 10.8 11.16 6.64 10.2 ...
##  $ R.14     : logi  NA NA NA NA NA NA ...
##  $ DOPCT    : num  96.7 108 131.9 64.3 94 ...
##  $ R.15     : logi  NA NA NA NA NA NA ...
##  $ pH       : num  7.46 7.78 7.34 7.35 7.21 6.62 6.83 6.9 6.72 6.6 ...
##  $ R.16     : logi  NA NA NA NA NA NA ...
##  $ SSC      : int  NA 14 51 52 9 NA 27 11 27 12 ...
##  $ R.17     : logi  NA NA NA NA NA NA ...
##  $ SSF      : int  NA NA 96 NA 100 NA NA 88 80 NA ...
##  $ R.21     : chr  NA NA NA NA ...
##  $ P04035   : num  0.058 0.019 0.043 0.221 0.0388 0.49 0.065 0.013 0.0124 0.0102 ...
##  $ R.23     : chr  "E" "E" "E" "<" ...
##  $ P04040   : num  0.007 0.012 0.006 0.002 0.0075 0.01 0.008 0.005 0.0038 0.0046 ...
##  $ R.41     : chr  NA NA NA NA ...
##  $ P39415   : num  0.005 0.006 0.12 0.106 0.0092 0.002 0.002 0.002 0.0037 0.002 ...
##  $ R.44     : chr  NA NA NA NA ...
##  $ P39572   : num  0.009 0.031 0.007 0.0366 0.002 0.002 0.002 0.002 0.002 0.002 ...
##  $ R.45     : chr  NA NA NA NA ...
##  $ P39632   : num  0.027 0.073 0.071 0.0444 0.0215 0.077 0.078 0.015 0.0133 0.0124 ...
##  $ R.62     : chr  NA NA "E" "E" ...
##  $ P49300   : num  0.66 0.35 0.02 0.03 0.44 0.02 0.02 0.02 0.02 0.02 ...
##  $ R.88     : chr  NA NA NA "<" ...
##  $ P82670   : num  0.063 0.16 0.019 0.04 0.0905 0.01 0.01 0.01 0.01 0.01 ...
##  $ R.94     : chr  "<" "<" "<" "<" ...
##  $ P82676   : num  0.003 0.003 0.003 0.003 0.003 0.003 0.003 0.003 0.003 0.003 ...
##  $ MapNumber: chr  "U3" "U3" "U3" "U3" ...
##  $ Longitude: num  -123 -123 -123 -123 -123 ...
##  $ Latitude : num  44.9 44.9 44.9 44.9 44.9 ...
##  $ Class    : chr  "URBAN" "URBAN" "URBAN" "URBAN" ...
##  $ Size     : int  5484 5484 5484 5484 5484 6173 6173 6173 6173 6173 ...
##  $ LU.Ag    : int  20 20 20 20 20 40 40 40 40 40 ...
##  $ LU.For   : int  0 0 0 0 0 49 49 49 49 49 ...
##  $ LU.Resid : int  80 80 80 80 80 7 7 7 7 7 ...
##  $ LU.Other : int  0 0 0 0 0 4 4 4 4 4 ...
##  $ Month    : chr  "Apr" "May" "Jul" "Oct" ...
##  $ NumCrops : int  NA NA NA NA NA 21 21 21 21 21 ...

Ajustament amb rpart()

Com que realitzar una regressió lineal múltiple no és un mètode viable per les característiques de les dades (rara vegada les covariables actuen de forma lineal en ciències ambientals) i el nombre de potencials models no lineals fa aquesta tècnica també inviable. Els autors del treball consideraren doncs com a variables predictives les que figuren en la funció rpart, predictives sobre el logaritme de la variable dependent de tipus continu.

set.seed(18061962)
pest.rpart <- rpart::rpart(log(P49300)~NH4+NO2+TKN+N2.3+TOTP+SRP+BOD+ECOL+FECAL+Longitude+Latitude+
                  Size+LU.Ag+LU.For+LU.Resid+LU.Other+NumCrops+Month, data=pesticides,
                method="anova")

rpart.plot::rpart.plot(pest.rpart, type = 2, # És el tipus que surt per defecte
                       main= "Arbre de regressió de les dades de contaminació del riu Willamette",
                       cex.main=1)

Cal recordar que executant el type=2 de la funció, a cada nus apareix, a més del criteri de partició, una etiqueta amb dos valors: el primer és el valor que agafaria la variable de resposta en aquell nus, i el segon és el % d’observacions en aquell nus (si és l’últim, serà una fulla, com sabem). Per exemple, mirant el nus més a baix i a l’esquerra, veiem que en aquell nus el valor de de la resposta és -3.8 i que agrupa un 26% de les observacions. Cal recordar també que l’ombrejat és proporcional al valor estimat de la resposta.

Veiem que, de les 18 variables explicatives introduïdes, rpart en descarta 11, deixant-ne tan sol 5. Al node arrel que la variable més predictiva és LU.Ag i els valors d’aquesta variable divideixen els individus de la mostra en dos grups: aquells per als quals LU.Ag < 74,5 (percentatge de terra utilitzada per a treballs agrícoles regada pel riu en qüestió), formant el node fill esquerre, i aquells per als quals LU.Ag ≥ 74,5, formant el node fill dret.

Per obtenir més informació sobre l’arbre resultant imprimirem l’objecte creat:

print(pest.rpart)

## n=94 (1 observation deleted due to missingness)
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 94 450.947300 -1.3326380  
##    2) LU.Ag< 74.5 63 183.558600 -2.2842210  
##      4) N2.3< 1.735 49 101.668000 -2.8005290  
##        8) LU.For>=5.5 32  18.751230 -3.4943280  
##         16) LU.Resid>=6.5 24   1.044089 -3.8400550 *
##         17) LU.Resid< 6.5 8   6.232528 -2.4571470 *
##        9) LU.For< 5.5 17  38.518650 -1.4945550 *
##      5) N2.3>=1.735 14  23.110940 -0.4771406 *
##    3) LU.Ag>=74.5 31  94.407100  0.6012235  
##      6) NH4< 0.0835 10  16.350190 -0.8724933 *
##      7) NH4>=0.0835 21  45.996400  1.3029930  
##       14) LU.Ag< 83.5 10  18.672770  0.3282101 *
##       15) LU.Ag>=83.5 11   9.183402  2.1891600 *

D’aquesta informació podem extreure com a observacions principals:

• dels 95 individus, 1 és no considerat,
• en els dos grups de què parlàvem més amunt, hi ha 63 individus al node filla de l’esquerra, i 31 al node filla de la dreta (recordar: 1. criteri de separació, 2. nombre de casos, 3. error residual, 4. valor estimat de la resposta, les fulles es marquen amb asterisc).

Així anirem comentant l’arbre de regressió, però el més interessant per a la predicció és en la informació subministrada en els nodes terminals

Per exemple, si el percentatge de terreny utilitzat per a labors agrícoles regat pel riu en qüestió és superior o igual al 74′5 %, el logaritme de la concentració de diüron és −0′87, si la concentració de NH4 és menor que 0′0835. De manera que (prenent antilogaritmes), en aquest cas, la concentració de diüron seria \(\small \exp(−0′87) = 0′4189515~ \mu/L\) contra \(\small3.6695~ \mu/L\) si el contingut de NH4 és superior.

Avaluació de l’ajustament

pest1 <- ggplot(data = data.frame(x=predict(pest.rpart), y=resid(pest.rpart)),
       aes(x,y)) +
  geom_point() +
  labs(x = "Valors esperats", 
       y = "Valors residuals",
       title = "Residuals vs Fitted") +
  geom_hline(yintercept = 0, colour= "red") +
  theme_bw()
pest2 <- qplot(sample = resid(pest.rpart)) + 
  stat_qq_line(colour= "red") +
  labs(title = "Q-Q Residuals",
       x= "Quantils teòrics",
       y= "Valors residuals") +
  theme_bw()

library(patchwork)
pest1 + pest2

Els valors residuals mostren un nivell de variança acceptable. Cal tenir en compte que hem transformat logarítmicament la variable de resposta.

Validació creuada i poda

Per veure si aquest representa el millor model final i saber quina és la grandària òptima de l’arbre seguirem també aquí un procés de validació creuada.

rpart::printcp(pest.rpart)

## 
## Regression tree:
## rpart::rpart(formula = log(P49300) ~ NH4 + NO2 + TKN + N2.3 + 
##     TOTP + SRP + BOD + ECOL + FECAL + Longitude + Latitude + 
##     Size + LU.Ag + LU.For + LU.Resid + LU.Other + NumCrops + 
##     Month, data = pesticides, method = "anova")
## 
## Variables actually used in tree construction:
## [1] LU.Ag    LU.For   LU.Resid N2.3     NH4     
## 
## Root node error: 450.95/94 = 4.7973
## 
## n=94 (1 observation deleted due to missingness)
## 
##         CP nsplit rel error  xerror    xstd
## 1 0.383596      0   1.00000 1.04534 0.10602
## 2 0.130347      1   0.61640 0.83603 0.12760
## 3 0.098455      2   0.48606 0.73806 0.11769
## 4 0.071096      3   0.38760 0.72092 0.12650
## 5 0.040227      4   0.31651 0.65376 0.12429
## 6 0.025446      5   0.27628 0.60900 0.11057
## 7 0.010000      6   0.25083 0.55094 0.10073

par(mar=c(5,5,7,2))
rpart::plotcp(pest.rpart, upper = "splits")
title("Validació creuada de l'arbre de regressió", line = 5, adj=0, cex.main=1.2)

En aquest cas, el nivell de poda serà amb un MSE de 0.55094 + 0.10073 = 0.65167, just després de la 4ª partició (0.65376).

Arbre podat

poda.pest.rpart <- rpart::prune.rpart(pest.rpart, cp=0.025446)

rpart.plot::rpart.plot(poda.pest.rpart,
                       main= "Arbre de regressió de dades de pesticides i podat a un valor de cp=0.025446", cex.main=1)