4 ISLP Ch9 Q7
(ns assignment.islp-9-7
(:require
[calc-metric.patch]
[clojisr.v1.applications.plotting
:refer [plot->svg]]
[clojisr.v1.r :refer [bra r+ r- r->clj clj->r]]
[clojisr.v1.require :refer [require-r]]
[fastmath.stats :as stats]
[scicloj.kindly.v4.kind :as kind]
[scicloj.metamorph.core :as morph]
[scicloj.metamorph.ml.gridsearch :as grid]
[tablecloth.api :as tc]))(comment
(clojure.java.shell/sh "which" "R"))- In this problem, you will use support vector approaches in order to predict whether a given car gets high or low gas mileage based on the Auto data set. > (a) Create a binary variable that takes on a 1 for cars with gas mileage above the median, and a 0 for cars with gas mileage below the median.
4.1 Binary response
Load the required R libraries
(require-r '[base :refer [RNGkind set-seed summary plot $ expand-grid which-max subset
as-numeric factor levels as-character data-frame append]]
'[stats :refer [predict]]
'[ISLR :as islr]
'[caret :refer [createDataPartition trainControl modelLookup train
defaultSummary prSummary twoClassSummary mnLogLoss]]
'[kernlab]
'[e1071]
'[ggplot2 :refer [ggplot aes geom_point geom_line
facet_wrap theme_bw]])nilCall in datasets from R. This one comes from the islr library in a dataset called Auto.
(def auto
(-> (r->clj islr/Auto)
(tc/drop-columns :$row.names)))(stats/median (:mpg auto))22.754.1.1 Create binary response
(def auto-cat
(-> auto
(tc/map-columns :mpg-cat [:mpg]
#(if (>= % (stats/median (:mpg auto))) 1 0))
(tc/map-columns :mpg-cat str)))Later, working with the data, R doesn’t like :keywords like Clojure does. Create a R-compatible data.frame. Notice, I’m flipping through both Clojure and R data structures in R functions and vice cersa.
(def r-data
(tc/rename-columns auto-cat (fn [col]
(-> col
name
(clojure.string/replace #"-" ".")))))(keys auto-cat)(:mpg
:cylinders
:displacement
:horsepower
:weight
:acceleration
:year
:origin
:name
:mpg-cat)(keys r-data)("mpg"
"cylinders"
"displacement"
"horsepower"
"weight"
"acceleration"
"year"
"origin"
"name"
"mpg.cat")(summary auto-cat) mpg cylinders displacement horsepower weight
Min. : 9.00 Min. :3.000 Min. : 68.0 Min. : 46.0 Min. :1613
1st Qu.:17.00 1st Qu.:4.000 1st Qu.:105.0 1st Qu.: 75.0 1st Qu.:2225
Median :22.75 Median :4.000 Median :151.0 Median : 93.5 Median :2804
Mean :23.45 Mean :5.472 Mean :194.4 Mean :104.5 Mean :2978
3rd Qu.:29.00 3rd Qu.:8.000 3rd Qu.:275.8 3rd Qu.:126.0 3rd Qu.:3615
Max. :46.60 Max. :8.000 Max. :455.0 Max. :230.0 Max. :5140
acceleration year origin name
Min. : 8.00 Min. :70.00 Min. :1.000 amc matador : 5
1st Qu.:13.78 1st Qu.:73.00 1st Qu.:1.000 ford pinto : 5
Median :15.50 Median :76.00 Median :1.000 toyota corolla : 5
Mean :15.54 Mean :75.98 Mean :1.577 chevrolet impala: 4
3rd Qu.:17.02 3rd Qu.:79.00 3rd Qu.:2.000 amc hornet : 4
Max. :24.80 Max. :82.00 Max. :3.000 ford maverick : 4
(Other) :365
mpg-cat
Length:392
Class :character
Mode :character
- Fit a support vector classifier to the data with various values of C, in order to predict whether a car gets high or low gas mileage. Report the cross-validation errors associated with different values of this parameter. Comment on your results. Note you will need to fit the classifier without the gas mileage variable to produce sensible results
4.2 Fit SVMs
Partition data
(def index
(createDataPartition :y ($ r-data 'mpg.cat)
:p 0.7 :list false))Train and test data
(def training-data
(bra r-data index nil))(def test-data
(bra r-data (r- index) nil))Caret svmLinear
(RNGkind :sample.kind "Rounding")[1] "Mersenne-Twister" "Inversion" "Rounding"
(set-seed 0)NULL
Bootstrap cross-validation
(def train-control
(trainControl :method "boot" :number 20))(modelLookup "svmLinear") model parameter label forReg forClass probModel
1 svmLinear C Cost TRUE TRUE TRUE
Build model
(def svm-linear
(train '(tilde mpg.cat (- . mpg))
:data training-data :method "svmLinear"
:trControl train-control :metric "kappa"
:tuneGrid (expand-grid :C (range 0.01 0.125 0.025))))View final model
($ svm-linear 'finalModel)Support Vector Machine object of class "ksvm"
SV type: C-svc (classification)
parameter : cost C = 0.035
Linear (vanilla) kernel function.
Number of Support Vectors : 57
Objective Function Value : -1.7603
Training error : 0.072464
(comment
(def plot-svm-linear
(r.e1071/svm '(formula mpg.cat (- . mpg))
:data (tc/convert-types r-data "mpg.cat" :int32)
:kernel "linear" :cost 0.01)))4.2.1 Cross-validation errors
($ svm-linear 'results) C Accuracy Kappa AccuracySD KappaSD
1 0.010 0.9171101 0.8332584 0.01775988 0.03560476
2 0.035 0.9174666 0.8339890 0.02033831 0.04088029
3 0.060 0.9173710 0.8337879 0.02270234 0.04575695
4 0.085 0.9126572 0.8243661 0.02194878 0.04390802
5 0.110 0.9116561 0.8222342 0.02235704 0.04485107
I measured the goodness-of-fit versus errors. For each C, R built 20 bootstraped linear SVM models. The best Kappa per C is reported. Based on Kappa, \(C = 0.035\) is best.
- Now repeat (b), this time using SVMs with radial and polynomial basis kernels, with different values of gamma and degree and C. Comment on your results.
4.3 SVM radial and polynomial
4.3.1 svmRadial
(RNGkind :sample.kind "Rounding")[1] "Mersenne-Twister" "Inversion" "Rounding"
(set-seed 0)NULL
Hyperparameter check
(modelLookup "svmRadial") model parameter label forReg forClass probModel
1 svmRadial sigma Sigma TRUE TRUE TRUE
2 svmRadial C Cost TRUE TRUE TRUE
Build model
(def svm-radial
(train '(formula mpg.cat (- . mpg))
:data training-data :method "svmRadial" :trControl train-control
:metric "kappa"
:tuneGrid (tc/dataset
(grid/sobol-gridsearch
{:C (grid/linear 0.25 150 8)
:sigma (grid/linear 0.000001 0.00001 5)}))))View final model
($ svm-radial 'finalModel)Support Vector Machine object of class "ksvm"
SV type: C-svc (classification)
parameter : cost C = 150
Gaussian Radial Basis kernel function.
Hyperparameter : sigma = 1e-05
Number of Support Vectors : 64
Objective Function Value : -5873.092
Training error : 0.039855
(comment
(def plot-svm-radial
(r.e1071/svm '(formula mpg.cat (- . mpg))
:data (tc/convert-types r-data "mpg.cat" :int32)
:kernel "radial" :cost 107.214285714286 :sigma 3.25e-06)))4.3.2 Cross-validation errors
($ svm-radial 'results) C sigma Accuracy Kappa AccuracySD KappaSD
1 0.25000 1.00e-06 0.8833750 0.7656933 0.02474869 0.04959927
6 21.64286 1.00e-06 0.8901821 0.7794242 0.02189154 0.04388273
11 43.03571 1.00e-06 0.8920549 0.7831918 0.02473057 0.04941727
16 64.42857 1.00e-06 0.8935104 0.7860984 0.02644102 0.05270203
21 85.82143 1.00e-06 0.8944305 0.7879690 0.02549771 0.05071660
26 107.21429 1.00e-06 0.8959157 0.7909338 0.02456325 0.04906032
31 128.60714 1.00e-06 0.8969336 0.7929813 0.02504592 0.05000107
36 150.00000 1.00e-06 0.8984212 0.7959607 0.02488068 0.04971731
2 0.25000 3.25e-06 0.8847289 0.7683157 0.02244202 0.04541825
7 21.64286 3.25e-06 0.8974278 0.7940053 0.02429315 0.04837000
12 43.03571 3.25e-06 0.8941426 0.7874291 0.02359858 0.04706370
17 64.42857 3.25e-06 0.8965624 0.7922130 0.02282044 0.04542483
22 85.82143 3.25e-06 0.8974051 0.7938333 0.02347745 0.04706646
27 107.21429 3.25e-06 0.8964420 0.7919512 0.02409085 0.04837719
32 128.60714 3.25e-06 0.8949473 0.7889674 0.02409934 0.04849591
37 150.00000 3.25e-06 0.8964737 0.7920664 0.02448085 0.04912228
3 0.25000 5.50e-06 0.8852682 0.7694862 0.02442095 0.04904901
8 21.64286 5.50e-06 0.8941408 0.7874105 0.02386209 0.04764312
13 43.03571 5.50e-06 0.8949784 0.7889673 0.02704639 0.05416533
18 64.42857 5.50e-06 0.8929821 0.7850209 0.02624850 0.05269261
23 85.82143 5.50e-06 0.8944339 0.7879447 0.02578807 0.05175839
28 107.21429 5.50e-06 0.8934026 0.7858215 0.02374042 0.04766020
33 128.60714 5.50e-06 0.8968897 0.7927721 0.02607739 0.05246945
38 150.00000 5.50e-06 0.9003425 0.7997520 0.02857987 0.05747702
4 0.25000 7.75e-06 0.8856591 0.7703052 0.02281895 0.04614578
9 21.64286 7.75e-06 0.8930841 0.7852677 0.02436272 0.04873047
14 43.03571 7.75e-06 0.8929078 0.7849590 0.02838410 0.05704814
19 64.42857 7.75e-06 0.8954140 0.7900840 0.02544598 0.05083392
24 85.82143 7.75e-06 0.8969092 0.7930197 0.02787478 0.05573365
29 107.21429 7.75e-06 0.9029020 0.8049847 0.02936503 0.05871301
34 128.60714 7.75e-06 0.9064132 0.8120334 0.02916810 0.05820951
39 150.00000 7.75e-06 0.9060271 0.8112559 0.03078577 0.06134504
5 0.25000 1.00e-05 0.8875706 0.7741968 0.02488036 0.05014435
10 21.64286 1.00e-05 0.8924980 0.7840832 0.02836727 0.05711738
15 43.03571 1.00e-05 0.8939510 0.7871905 0.02735880 0.05474495
20 64.42857 1.00e-05 0.8980007 0.7952439 0.02905157 0.05803917
25 85.82143 1.00e-05 0.9010134 0.8012216 0.02900648 0.05780814
30 107.21429 1.00e-05 0.9058099 0.8106791 0.02781523 0.05558962
35 128.60714 1.00e-05 0.9078532 0.8147551 0.02702958 0.05401028
40 150.00000 1.00e-05 0.9093257 0.8176629 0.02677906 0.05357765
Much bigger grid to search through. A plot would make this easier.
4.3.3 svmRadial
(RNGkind :sample.kind "Rounding")[1] "Mersenne-Twister" "Inversion" "Rounding"
(set-seed 0)NULL
Hyperparameter check
(modelLookup "svmPoly") model parameter label forReg forClass probModel
1 svmPoly degree Polynomial Degree TRUE TRUE TRUE
2 svmPoly scale Scale TRUE TRUE TRUE
3 svmPoly C Cost TRUE TRUE TRUE
Build model
(comment
;too much time
(def svm-poly
(train '(tilde mpg.cat
(+ cylinders displacement horsepower weight
acceleration year origin name))
:data training-data :method "svmPoly"
:trControl train-control :metric "kappa"
:tuneGrid (tc/dataset
(take 6
(grid/sobol-gridsearch
{:C (grid/linear 0.25 2 5)
:scale (grid/linear 0.001 1 5)
:degree (grid/linear 1 2 2 :int16)})))))
;=> Support Vector Machine object of class "ksvm"
;
;SV type: C-svc (classification)
; parameter : cost C = 1.5625
;
;Polynomial kernel function.
; Hyperparameters : degree = 1 scale = 0.001 offset = 1
;
;Number of Support Vectors : 78
;
;Objective Function Value : -113.7188
;Training error : 0.094203
($ svm-poly 'finalModel))(def svm-poly
(train '(tilde mpg.cat
(+ cylinders displacement horsepower weight
acceleration year origin name))
:data training-data :method "svmPoly"
:trControl train-control :metric "kappa"
:tuneGrid (tc/dataset
{:C 1.5625
:scale 0.001
:degree 1})))View final model
($ svm-poly 'finalModel)Support Vector Machine object of class "ksvm"
SV type: C-svc (classification)
parameter : cost C = 1.5625
Polynomial kernel function.
Hyperparameters : degree = 1 scale = 0.001 offset = 1
Number of Support Vectors : 72
Objective Function Value : -97.6261
Training error : 0.068841
(def plot-svm-poly
(r.e1071/svm '(formula mpg.cat (- . mpg))
:data (tc/convert-types r-data "mpg.cat" :int32)
:kernel "polynomial" :cost 1.5625 :degree 1 :scale 0.001))(-> (plot plot-svm-poly
:data (-> (r->clj test-data)
(tc/drop-columns [:$row.names :mpg.cat])
clj->r)
:formula '(tilde displacement weight))
plot->svg)"<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"504pt\" height=\"504pt\" viewBox=\"0 0 504 504\" version=\"1.1\">\n<g id=\"surface91\">\n</g>\n</svg>\n"(predict plot-svm-poly test-data) 5 6 8 12 14 20
0.093884149 -0.127112610 -0.131835504 0.049035557 -0.032607662 0.924874267
21 26 27 31 32 33
0.834702676 -0.105369698 -0.032008226 0.821960241 0.936184818 0.466043634
36 37 41 42 44 49
0.379219163 0.386469818 0.013184390 -0.138358203 -0.175490237 0.851458132
50 54 57 64 68 69
0.906968346 1.022571797 0.930562366 0.015408963 -0.054924619 -0.051095613
70 74 75 77 80 89
-0.098763838 0.003993701 0.016675792 0.848765335 0.828472930 0.059221205
92 103 110 122 131 132
-0.034464354 -0.145905014 0.924478938 0.860436102 0.809709438 0.345093379
133 136 137 138 145 153
0.338051352 -0.061735951 -0.003350196 -0.014457039 0.922986143 0.378869456
154 161 163 169 170 171
0.410483645 0.332180080 0.456533639 0.811218804 0.901909615 0.946270965
182 186 188 190 191 203
0.927932948 0.050253581 0.068373214 0.448578887 0.411629347 1.016574970
207 210 211 212 215 218
0.064555419 0.477481333 -0.010333697 0.029146410 1.000284833 0.885044006
227 233 234 240 241 243
0.389889148 0.810420481 0.978950707 0.661027498 0.897052642 0.975836541
244 247 250 251 258 259
0.975490065 1.049862920 0.138400085 0.414194856 0.506510673 0.410251753
262 265 269 273 276 278
0.437350089 0.937877037 0.932893313 0.749171232 0.542348367 1.007347573
280 281 285 292 294 295
0.516670389 0.812230956 0.130105504 1.018648857 0.969676602 0.869772023
296 298 303 304 307 309
0.625985431 0.797638595 0.993188346 0.833773811 0.859628851 1.042122501
310 313 319 320 321 322
0.961846507 0.824304116 0.974988178 1.005300592 1.033585316 0.833026766
324 325 326 327 335 337
0.981389405 0.949713970 0.758088054 0.805468928 0.889552575 0.856759240
339 344 346 349 354 355
0.911943366 1.049921033 1.054456363 0.917688704 0.951089412 0.824438769
356 357 361 367 376 381
0.609667603 0.673222268 0.530197905 0.890472147 1.031906582 0.573524557
384 389
0.535200332 0.973548498
4.3.4 Cross-validation errors
($ svm-poly 'results) C scale degree Accuracy Kappa AccuracySD KappaSD
1 1.5625 0.001 1 0.9051167 0.809254 0.02522795 0.05062125
- Make some plots to back up your assertions in (b) and (c).
4.4 Plots
^kind/html
(-> (plot svm-linear :metric "Kappa")
plot->svg)^kind/html
(comment
(-> (r.kernlab/plot ($ svm-linear 'finalModel)
:data (-> (r->clj test-data)
(tc/drop-columns [:$row.names :mpg :mpg.cat])
clj->r)
:formula '(formula displacement weight)
:slice {'displacement 2 'weight 4})
plot->svg))According to our text, we could use > plot(svmfit, dat, x1 ∼ x4) to plot. I was unsuccessful.
^kind/html
(-> (plot plot-svm-poly
:data test-data
:formula '(formula displacement weight))
plot->svg)But predict works.
(predict plot-svm-poly test-data) 5 6 8 12 14 20
0.093884149 -0.127112610 -0.131835504 0.049035557 -0.032607662 0.924874267
21 26 27 31 32 33
0.834702676 -0.105369698 -0.032008226 0.821960241 0.936184818 0.466043634
36 37 41 42 44 49
0.379219163 0.386469818 0.013184390 -0.138358203 -0.175490237 0.851458132
50 54 57 64 68 69
0.906968346 1.022571797 0.930562366 0.015408963 -0.054924619 -0.051095613
70 74 75 77 80 89
-0.098763838 0.003993701 0.016675792 0.848765335 0.828472930 0.059221205
92 103 110 122 131 132
-0.034464354 -0.145905014 0.924478938 0.860436102 0.809709438 0.345093379
133 136 137 138 145 153
0.338051352 -0.061735951 -0.003350196 -0.014457039 0.922986143 0.378869456
154 161 163 169 170 171
0.410483645 0.332180080 0.456533639 0.811218804 0.901909615 0.946270965
182 186 188 190 191 203
0.927932948 0.050253581 0.068373214 0.448578887 0.411629347 1.016574970
207 210 211 212 215 218
0.064555419 0.477481333 -0.010333697 0.029146410 1.000284833 0.885044006
227 233 234 240 241 243
0.389889148 0.810420481 0.978950707 0.661027498 0.897052642 0.975836541
244 247 250 251 258 259
0.975490065 1.049862920 0.138400085 0.414194856 0.506510673 0.410251753
262 265 269 273 276 278
0.437350089 0.937877037 0.932893313 0.749171232 0.542348367 1.007347573
280 281 285 292 294 295
0.516670389 0.812230956 0.130105504 1.018648857 0.969676602 0.869772023
296 298 303 304 307 309
0.625985431 0.797638595 0.993188346 0.833773811 0.859628851 1.042122501
310 313 319 320 321 322
0.961846507 0.824304116 0.974988178 1.005300592 1.033585316 0.833026766
324 325 326 327 335 337
0.981389405 0.949713970 0.758088054 0.805468928 0.889552575 0.856759240
339 344 346 349 354 355
0.911943366 1.049921033 1.054456363 0.917688704 0.951089412 0.824438769
356 357 361 367 376 381
0.609667603 0.673222268 0.530197905 0.890472147 1.031906582 0.573524557
384 389
0.535200332 0.973548498
(comment
(-> (r->clj test-data)
(tc/drop-columns [:$row.names :mpg.cat])
clj->r))4.4.0.1 SVM Radial
^kind/html
(-> (plot svm-radial :metric "Kappa")
plot->svg)^kind/html
(-> (plot svm-radial :metric "Kappa" :plotType "level")
plot->svg)4.4.0.2 SVM Polynomial
This model took the longest to train. Having 5 sets of hyperparameters added two minutes to rendering. If I ran the full commented out :tuneGrid, the following three plots will view accuracy measures, like the above plots. But with one data we have:
^kind/html
(-> (ggplot :data ($ svm-poly 'results)
(aes :x 'C :y 'Kappa :color '(factor scale)))
(r+ (geom_point)
(geom_line)
(facet_wrap '(formula nil degree))
(theme_bw))
plot->svg)(comment
^kind/html
(-> (plot svm-poly :metric "Kappa")
plot->svg)
^kind/html
(-> (plot svm-poly :metric "Kappa" :plotType "level")
plot->svg))4.5 Evaluate model
(defn eval-list [model]
(let [pred (predict model test-data)
obs (factor ($ test-data 'mpg.cat))
df (data-frame :pred pred :obs obs)
ds (defaultSummary df :lev (levels (factor ($ test-data 'mpg.cat))))
tcs (twoClassSummary df :lev (levels (factor ($ test-data 'mpg.cat))))]
(append ds tcs)))(eval-list svm-linear) Accuracy Kappa ROC Sens Spec
0.8620690 0.7241379 NA 0.8103448 0.9137931
(eval-list svm-radial) Accuracy Kappa ROC Sens Spec
0.8706897 0.7413793 NA 0.8275862 0.9137931
(eval-list svm-poly) Accuracy Kappa ROC Sens Spec
0.8706897 0.7413793 NA 0.8448276 0.8965517
Polynomial model preformed best. Turns out, the best was degree 1, so linear. Then linear itself.
source: src/assignment/islp_9_7.clj