4  Caret with R Interop

(ns assignment.r-caret
  (:require
    [assignment.eda :refer [concrete-data]]
    [calc-metric.patch]
    [clojisr.v1.applications.plotting
     :refer [plot->svg]]
    [clojisr.v1.r :refer [bra r+ r-]]
    [clojisr.v1.require :refer [require-r]]
    [scicloj.kindly.v4.kind :as kind]
    [scicloj.ml.dataset :as ds]))
(comment
  (clojure.java.shell/sh "which" "R"))

In this chapter, I will graph a tuned decision tree. The issue with Smile’s implementation is that there is no .dot method of class gradient.tree like there is with the classification tree in the scicloj models tutorial.

In native Clojure, I tried the scicloj Kroki method, but failed. So I turned to my trusty R… in Clojure. While R has a nice plotting function for tree models, viz. prp and rpart.plot, neither of them returned a plottable output. Instead they returned the data structure that would, I suppose, pass to a graphing library.

Because I find prp and rpart.plot outputs beautiful, I tried to plot them in Clojure with R’s Rgraphviz, igraph, and tree.data libraries, all to no success.

4.1 Load the required R libraries

(require-r '[base :refer [RNGkind set-seed summary plot $ which-min]]
           '[stats :refer [predict]]
           '[caret :refer [createDataPartition trainControl modelLookup
                           train defaultSummary postResample]]
           '[rpart :refer [rpart]]
           '[rpart.plot :refer [rpart-plot prp]]
           '[ggplot2 :refer [ggplot aes geom_segment geom_text scale_size]]
           '[ggdendro :refer [dendro_data theme_dendro segment
                              label leaf_label]])
nil
(summary concrete-data)
     cement      blast-furnace-slag    fly-ash           water      
 Min.   :102.0   Min.   :  0.0      Min.   :  0.00   Min.   :121.8  
 1st Qu.:192.4   1st Qu.:  0.0      1st Qu.:  0.00   1st Qu.:164.9  
 Median :272.9   Median : 22.0      Median :  0.00   Median :185.0  
 Mean   :281.2   Mean   : 73.9      Mean   : 54.19   Mean   :181.6  
 3rd Qu.:350.0   3rd Qu.:142.9      3rd Qu.:118.30   3rd Qu.:192.0  
 Max.   :540.0   Max.   :359.4      Max.   :200.10   Max.   :247.0  
 superplasticizer coarse-aggregate fine-aggregate       age        
 Min.   : 0.000   Min.   : 801.0   Min.   :594.0   Min.   :  1.00  
 1st Qu.: 0.000   1st Qu.: 932.0   1st Qu.:731.0   1st Qu.:  7.00  
 Median : 6.400   Median : 968.0   Median :779.5   Median : 28.00  
 Mean   : 6.205   Mean   : 972.9   Mean   :773.6   Mean   : 45.66  
 3rd Qu.:10.200   3rd Qu.:1029.4   3rd Qu.:824.0   3rd Qu.: 56.00  
 Max.   :32.200   Max.   :1145.0   Max.   :992.6   Max.   :365.00  
 concrete-compressive-strength
 Min.   : 2.33                
 1st Qu.:23.71                
 Median :34.45                
 Mean   :35.82                
 3rd Qu.:46.13                
 Max.   :82.60                

4.2 Setup dataset

The R interop does not like Clojure’s convention of columns as keyword (with their prefixed colons) nor the kabab-casing. Below is my assisting Clojure-R op to work with the data.

(def r-data
  (ds/rename-columns concrete-data (fn [col]
                                     (-> col
                                         name               ; removes ":" from type keyword
                                         (clojure.string/replace #"-" ".")))))

4.3 Partition data

(def index
  (createDataPartition :y (:concrete-compressive-strength concrete-data)
                       :p 0.7 :list false))
(def training-data
  (bra r-data index nil))
(def test-data
  (bra r-data (r- index) nil))

4.4 Caret decision tree

(RNGkind :sample.kind "Rounding")
[1] "Mersenne-Twister" "Inversion"        "Rounding"        
(set-seed 0)
NULL

4.4.1 Bootstrap cross-validation

(def train-control
  (trainControl :method "boot" :number 25))
(modelLookup "rpart")
  model parameter                label forReg forClass probModel
1 rpart        cp Complexity Parameter   TRUE     TRUE      TRUE

4.4.2 Build model

(def decision-tree
  (train '(tilde concrete.compressive.strength
                 (+ cement blast.furnace.slag fly.ash water
                    superplasticizer coarse.aggregate fine.aggregate
                    age concrete.compressive.strength)) :data training-data
         :method "rpart" :trControl train-control :metric "MAE" :tuneLength 20))
decision-tree
CART 

722 samples
  8 predictor

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 722, 722, 722, 722, 722, 722, ... 
Resampling results across tuning parameters:

  cp           RMSE       Rsquared   MAE      
  0.005851617   8.682246  0.7356957   6.673698
  0.006316008   8.747360  0.7316275   6.719544
  0.006414432   8.744562  0.7316855   6.720624
  0.006681596   8.777276  0.7293616   6.747759
  0.006786899   8.787395  0.7288109   6.751666
  0.007753180   8.878819  0.7230280   6.827006
  0.008293984   8.921629  0.7203711   6.864426
  0.011267153   9.229526  0.7001990   7.143066
  0.011320380   9.233150  0.7001072   7.150360
  0.014689720   9.428913  0.6868843   7.317053
  0.016414669   9.510373  0.6812937   7.398606
  0.022825559  10.045104  0.6439139   7.859829
  0.023449328  10.080486  0.6413338   7.885305
  0.025516891  10.161492  0.6356217   7.972034
  0.035133607  10.564518  0.6067419   8.298752
  0.043982862  11.087541  0.5662961   8.701584
  0.061654418  11.778771  0.5118749   9.262753
  0.075570047  12.230407  0.4726947   9.670118
  0.196354736  14.923635  0.2161278  12.081176
  0.239627405  15.449365  0.2005148  12.519071

MAE was used to select the optimal model using the smallest value.
The final value used for the model was cp = 0.005851617.
^kind/hiccup
(-> (plot decision-tree)
    plot->svg)
($ decision-tree 'results)
            cp      RMSE  Rsquared       MAE    RMSESD RsquaredSD     MAESD
1  0.005851617  8.682246 0.7356957  6.673698 0.3108629 0.02138865 0.2447274
2  0.006316008  8.747360 0.7316275  6.719544 0.3280370 0.02179039 0.2834077
3  0.006414432  8.744562 0.7316855  6.720624 0.3315382 0.02196734 0.2813194
4  0.006681596  8.777276 0.7293616  6.747759 0.3495973 0.02301042 0.2965770
5  0.006786899  8.787395 0.7288109  6.751666 0.3588749 0.02342017 0.2968638
6  0.007753180  8.878819 0.7230280  6.827006 0.3976424 0.02328981 0.3191092
7  0.008293984  8.921629 0.7203711  6.864426 0.3919333 0.02272673 0.3247436
8  0.011267153  9.229526 0.7001990  7.143066 0.4663804 0.02581972 0.3823203
9  0.011320380  9.233150 0.7001072  7.150360 0.4650576 0.02606384 0.3746424
10 0.014689720  9.428913 0.6868843  7.317053 0.5124020 0.03064604 0.4308113
11 0.016414669  9.510373 0.6812937  7.398606 0.5982501 0.03694672 0.5153925
12 0.022825559 10.045104 0.6439139  7.859829 0.8174476 0.05387524 0.6745357
13 0.023449328 10.080486 0.6413338  7.885305 0.8120077 0.05418223 0.6694464
14 0.025516891 10.161492 0.6356217  7.972034 0.7841843 0.05219744 0.6575686
15 0.035133607 10.564518 0.6067419  8.298752 0.7355509 0.04995815 0.5921249
16 0.043982862 11.087541 0.5662961  8.701584 0.7118188 0.05087474 0.5950899
17 0.061654418 11.778771 0.5118749  9.262753 0.7743284 0.05178480 0.7723233
18 0.075570047 12.230407 0.4726947  9.670118 0.9111771 0.06863740 0.8836898
19 0.196354736 14.923635 0.2161278 12.081176 0.7459286 0.07226630 0.6810511
20 0.239627405 15.449365 0.2005148 12.519071 0.9548339 0.03016883 0.7238487

4.4.3 Best hyperparameter

(-> ($ decision-tree 'results)
    (bra (which-min
           (bra ($ decision-tree 'results) nil 4))
         nil))
           cp     RMSE  Rsquared      MAE    RMSESD RsquaredSD     MAESD
1 0.005851617 8.682246 0.7356957 6.673698 0.3108629 0.02138865 0.2447274

4.4.4 View tree

($ decision-tree 'finalModel)
n= 722 

node), split, n, deviance, yval
      * denotes terminal node

  1) root 722 205430.5000 35.85429  
    2) age< 21 228  36002.1700 23.70004  
      4) cement< 354.5 159  11463.3300 18.26421  
        8) age< 10.5 116   4186.1150 14.76836 *
        9) age>=10.5 43   2035.2640 27.69488 *
      5) cement>=354.5 69   9014.4550 36.22609  
       10) blast.furnace.slag< 101.7 46   4662.6590 32.13065 *
       11) blast.furnace.slag>=101.7 23   2037.1790 44.41696  
         22) age< 5 11    104.0818 35.42727 *
         23) age>=5 12    229.2604 52.65750 *
    3) age>=21 494 120201.5000 41.46395  
      6) cement< 355.95 379  60468.7600 36.48636  
       12) cement< 164.8 83   7934.0720 25.56940  
         24) blast.furnace.slag< 115.5 27    321.9453 14.59778 *
         25) blast.furnace.slag>=115.5 56   2794.9200 30.85929 *
       13) cement>=164.8 296  39868.9900 39.54753  
         26) blast.furnace.slag< 14.3 164  15963.1300 35.11744  
           52) age< 42 81   4589.5770 29.70469  
            104) coarse.aggregate>=1086.7 13    204.2960 20.48000 *
            105) coarse.aggregate< 1086.7 68   3067.5610 31.46824 *
           53) age>=42 83   6684.4910 40.39976  
            106) superplasticizer< 8.25 57   1718.0680 36.32737 *
            107) superplasticizer>=8.25 26   1948.7070 49.32769 *
         27) blast.furnace.slag>=14.3 132  16688.3400 45.05159  
           54) superplasticizer< 7.3 78   6466.7660 40.84615  
            108) age< 42 43   3173.6050 36.76930  
              216) cement< 263.5 25    898.6939 31.93760 *
              217) cement>=263.5 18    880.6752 43.48000 *
            109) age>=42 35   1700.4210 45.85486 *
           55) superplasticizer>=7.3 54   6849.4970 51.12611  
            110) blast.furnace.slag< 157.7 42   2774.6890 47.61833 *
            111) blast.furnace.slag>=157.7 12   1749.2570 63.40333 *
      7) cement>=355.95 115  19395.5400 57.86835  
       14) water>=183.05 45   4939.5690 46.81311  
         28) superplasticizer< 3.45 36   2420.7300 44.05167 *
         29) superplasticizer>=3.45 9   1146.2350 57.85889 *
       15) water< 183.05 70   5420.5510 64.97529  
         30) blast.furnace.slag< 170.1 55   3594.8440 62.72691 *
         31) blast.furnace.slag>=170.1 15    528.2065 73.21933 *

4.5 Plot tree

(def plot-tree
  (rpart-plot ($ decision-tree 'finalModel)
              :type 1 :extra 1 :under true :split.font 2 :varlen -10))
(def ddata
  (dendro_data ($ decision-tree 'finalModel)))
^kind/hiccup
(-> (ggplot)
    (r+ (geom_segment :data ($ ddata 'segments)
                      (aes :x 'x :y 'y :xend 'xend :yend 'yend)))
    (r+ (geom_text :data ($ ddata 'labels)
                   (aes :x 'x :y 'y :label 'label)
                   :size 3 :vjust -0.5))
    (r+ (geom_text :data ($ ddata 'leaf_labels)
                   (aes :x 'x :y 'y :label 'label)
                   :size 3 :vjust 1))
    (r+ (theme_dendro))
    plot->svg)

Truly an ugly tree. As I said above, prp and rpart.plot do a much better job, but if you see the plot-tree variable defined above, you’d see that it is only a data structure, not a plot as it would be in R.

^kind/hiccup
(-> (ggplot (segment ddata))
    (r+ (geom_segment (aes :x 'x :y 'y :size 'n :xend 'xend :yend 'yend)
                      :color "blue" :alpha 0.5))
    (r+ (geom_text :data (label ddata)
                   (aes :x 'x :y 'y :label 'label)
                   :size 3 :vjust -0.5))
    (r+ (geom_text :data (leaf_label ddata)
                   (aes :x 'x :y 'y :label 'label)
                   :size 3 :vjust 1))
    (r+ (theme_dendro))
    plot->svg)

This plot is almost worse in terms of good looks. But there was a vignette of it and I tried sooo many other ways of plotting a tree, I had to include it.

4.6 Evaluate model

(def pred
  (predict ($ decision-tree 'finalModel) test-data))
(postResample pred ($ test-data 'concrete.compressive.strength))
     RMSE  Rsquared       MAE 
9.4059631 0.6778407 7.1372076 

Somewhat surprisingly, Smile did much better in terms of R\(^2\). However, the amount of hyperparameters a gradient tree can take is much more than the rpart tree. Not to mention, how wide I let the grid be. Realistically, I would higher the node-size as to not allow more depth as to lower variance. The MAE and RMSE metrics are almost the same.