-
Notifications
You must be signed in to change notification settings - Fork 1
/
classification.html
1206 lines (1081 loc) · 52.5 KB
/
classification.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8" />
<meta name="generator" content="pandoc" />
<meta http-equiv="X-UA-Compatible" content="IE=EDGE" />
<title>Classification with Logistic Regression</title>
<script src="site_libs/header-attrs-2.25/header-attrs.js"></script>
<script src="site_libs/jquery-3.6.0/jquery-3.6.0.min.js"></script>
<meta name="viewport" content="width=device-width, initial-scale=1" />
<link href="site_libs/bootstrap-3.3.5/css/flatly.min.css" rel="stylesheet" />
<script src="site_libs/bootstrap-3.3.5/js/bootstrap.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/html5shiv.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/respond.min.js"></script>
<style>h1 {font-size: 34px;}
h1.title {font-size: 38px;}
h2 {font-size: 30px;}
h3 {font-size: 24px;}
h4 {font-size: 18px;}
h5 {font-size: 16px;}
h6 {font-size: 12px;}
code {color: inherit; background-color: rgba(0, 0, 0, 0.04);}
pre:not([class]) { background-color: white }</style>
<script src="site_libs/navigation-1.1/tabsets.js"></script>
<link href="site_libs/highlightjs-9.12.0/default.css" rel="stylesheet" />
<script src="site_libs/highlightjs-9.12.0/highlight.js"></script>
<link href="site_libs/vembedr-0.1.5/css/vembedr.css" rel="stylesheet" />
<link href="site_libs/font-awesome-6.4.2/css/all.min.css" rel="stylesheet" />
<link href="site_libs/font-awesome-6.4.2/css/v4-shims.min.css" rel="stylesheet" />
<style type="text/css">
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
span.underline{text-decoration: underline;}
div.column{display: inline-block; vertical-align: top; width: 50%;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
</style>
<style type="text/css">code{white-space: pre;}</style>
<script type="text/javascript">
if (window.hljs) {
hljs.configure({languages: []});
hljs.initHighlightingOnLoad();
if (document.readyState && document.readyState === "complete") {
window.setTimeout(function() { hljs.initHighlighting(); }, 0);
}
}
</script>
<style type = "text/css">
.main-container {
max-width: 940px;
margin-left: auto;
margin-right: auto;
}
img {
max-width:100%;
}
.tabbed-pane {
padding-top: 12px;
}
.html-widget {
margin-bottom: 20px;
}
button.code-folding-btn:focus {
outline: none;
}
summary {
display: list-item;
}
details > summary > p:only-child {
display: inline;
}
pre code {
padding: 0;
}
</style>
<style type="text/css">
.dropdown-submenu {
position: relative;
}
.dropdown-submenu>.dropdown-menu {
top: 0;
left: 100%;
margin-top: -6px;
margin-left: -1px;
border-radius: 0 6px 6px 6px;
}
.dropdown-submenu:hover>.dropdown-menu {
display: block;
}
.dropdown-submenu>a:after {
display: block;
content: " ";
float: right;
width: 0;
height: 0;
border-color: transparent;
border-style: solid;
border-width: 5px 0 5px 5px;
border-left-color: #cccccc;
margin-top: 5px;
margin-right: -10px;
}
.dropdown-submenu:hover>a:after {
border-left-color: #adb5bd;
}
.dropdown-submenu.pull-left {
float: none;
}
.dropdown-submenu.pull-left>.dropdown-menu {
left: -100%;
margin-left: 10px;
border-radius: 6px 0 6px 6px;
}
</style>
<script type="text/javascript">
// manage active state of menu based on current page
$(document).ready(function () {
// active menu anchor
href = window.location.pathname
href = href.substr(href.lastIndexOf('/') + 1)
if (href === "")
href = "index.html";
var menuAnchor = $('a[href="' + href + '"]');
// mark the anchor link active (and if it's in a dropdown, also mark that active)
var dropdown = menuAnchor.closest('li.dropdown');
if (window.bootstrap) { // Bootstrap 4+
menuAnchor.addClass('active');
dropdown.find('> .dropdown-toggle').addClass('active');
} else { // Bootstrap 3
menuAnchor.parent().addClass('active');
dropdown.addClass('active');
}
// Navbar adjustments
var navHeight = $(".navbar").first().height() + 15;
var style = document.createElement('style');
var pt = "padding-top: " + navHeight + "px; ";
var mt = "margin-top: -" + navHeight + "px; ";
var css = "";
// offset scroll position for anchor links (for fixed navbar)
for (var i = 1; i <= 6; i++) {
css += ".section h" + i + "{ " + pt + mt + "}\n";
}
style.innerHTML = "body {" + pt + "padding-bottom: 40px; }\n" + css;
document.head.appendChild(style);
});
</script>
<!-- tabsets -->
<style type="text/css">
.tabset-dropdown > .nav-tabs {
display: inline-table;
max-height: 500px;
min-height: 44px;
overflow-y: auto;
border: 1px solid #ddd;
border-radius: 4px;
}
.tabset-dropdown > .nav-tabs > li.active:before, .tabset-dropdown > .nav-tabs.nav-tabs-open:before {
content: "\e259";
font-family: 'Glyphicons Halflings';
display: inline-block;
padding: 10px;
border-right: 1px solid #ddd;
}
.tabset-dropdown > .nav-tabs.nav-tabs-open > li.active:before {
content: "\e258";
font-family: 'Glyphicons Halflings';
border: none;
}
.tabset-dropdown > .nav-tabs > li.active {
display: block;
}
.tabset-dropdown > .nav-tabs > li > a,
.tabset-dropdown > .nav-tabs > li > a:focus,
.tabset-dropdown > .nav-tabs > li > a:hover {
border: none;
display: inline-block;
border-radius: 4px;
background-color: transparent;
}
.tabset-dropdown > .nav-tabs.nav-tabs-open > li {
display: block;
float: none;
}
.tabset-dropdown > .nav-tabs > li {
display: none;
}
</style>
<!-- code folding -->
</head>
<body>
<div class="container-fluid main-container">
<div class="navbar navbar-default navbar-fixed-top" role="navigation">
<div class="container">
<div class="navbar-header">
<button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-bs-toggle="collapse" data-target="#navbar" data-bs-target="#navbar">
<span class="icon-bar"></span>
<span class="icon-bar"></span>
<span class="icon-bar"></span>
</button>
<a class="navbar-brand" href="index.html">Machine Learning for Public Policy</a>
</div>
<div id="navbar" class="navbar-collapse collapse">
<ul class="nav navbar-nav">
</ul>
<ul class="nav navbar-nav navbar-right">
<li>
<a href="index.html">
<span class="fa fa-home"></span>
Home
</a>
</li>
<li>
<a href="intro.html">
<span class="fa fa-duotone fa-robot"></span>
Introduction
</a>
</li>
<li>
<a href="predictionpolicy.html">
<span class="fa fa-line-chart"></span>
Prediction Policy Problems
</a>
</li>
<li>
<a href="classification.html">
<span class="fa fa-solid fa-gears"></span>
Classification:Logistic
</a>
</li>
<li>
<a href="treebasedmodels.html">
<span class="fa fa-tree"></span>
TreeModels:RandomForests
</a>
</li>
<li>
<a href="fairml.html">
<span class="fa fa-graduation-cap"></span>
Fair ML/Data Ethics
</a>
</li>
<li>
<a href="NeuralNets.html">
<span class="fa fa-superpowers"></span>
Neural Networks
</a>
</li>
<li>
<a href="discussionboard.html">
<span class="fa fa-solid fa-comments"></span>
Discussion Board
</a>
</li>
</ul>
</div><!--/.nav-collapse -->
</div><!--/.container -->
</div><!--/.navbar -->
<div id="header">
<h1 class="title toc-ignore">Classification with Logistic
Regression</h1>
</div>
<style>
body {
text-align: justify}
</style>
<div
id="prediction-policy-problems-classification-with-logistic-regression"
class="section level2 tabset tabset-fade tabset-pills">
<h2 class="tabset tabset-fade tabset-pills"><strong>Prediction Policy
Problems: Classification with Logistic Regression</strong></h2>
<p>Have you heard the English proverb, “Birds of a feather flock
together”? It references and old saying that indicates that people with
similar characteristics tend to group and stay together. In Machine
Learning, Classification problems deal with the evaluation of models of
categorical response, such as:</p>
<ul>
<li><p>Predictive classification: E.g. is this spam or not? Predictive
classification concerns itself with unlabeled data, and groups them by
the proportion of characteristics they commonly share. After which, it
classifies them into some predetermined category. A common, ‘lazy’
method is kNearest Neighbors.</p></li>
<li><p><strong>Binary classification:</strong> You may already be
familiar with probit or logistic regression models. You obtain two types
of predictions from such models: proportions, and the generation of a
predicted discrete choice. For Policy purposes, we are interested in the
discrete choice. E.g. filtering low-income individuals to select those
who will receive social assistance and those who will not, based on some
income/expenditure threshold. But, we still need the probability
estimates of each of the two categories. They are relevant when working
out the model’s confidence about the predicted discrete choice.</p></li>
<li><p>Multi-label classification: Not unlike binary classification, it
is a labeled data model that relies on techniques such as multinomial
logistic regression. It deals with data with more than two categories,
and generates discrete choices, which policymakers then rely on to make
decisions.</p></li>
</ul>
<p>In the video-lecture below you’ll get an intuitive explanation of
what a logistic regression model is, and how we can use it in the
context of a prediction policy framework.</p>
<center>
<div class="vembedr">
<div>
<iframe src="https://www.youtube.com/embed/A9qVrFhlRMY" width="533" height="300" frameborder="0" allowfullscreen="" data-external="1"></iframe>
</div>
</div>
</center>
<p>After watching the video, below you’ll find a continuation of our
previous exercise. Previously, we were working on predicting per capita
monthly expenditures of a sample of individuals from Malawi. Our
assumption is that by predicting how much a person spends per month, we
can infer whether they are in poverty (or not) by contrasting that value
to other relevant information, such as the cost of food and rent in the
country. Another way to go about this is to use the estimated poverty
line, and generate a variable that takes on the value <span
class="math inline">\(1\)</span> if the person’s expenditure is below
the poverty line (they are poor) and <span
class="math inline">\(0\)</span> otherwise (not poor). Thus, our policy
problem becomes one of classification.</p>
<div id="r-practical" class="section level3">
<h3><strong>R practical</strong></h3>
<p>We will continue to work with the Malawi dataset, which can be
downloaded in the (Prediction Policy Problems)[<a
href="https://www.ml4publicpolicy.com/predictionpolicy.html"
class="uri">https://www.ml4publicpolicy.com/predictionpolicy.html</a>]
tab of this website.</p>
<h3>
<ol style="list-style-type: decimal">
<li>Preliminaries: working directory, libraries, data upload
</h3>
<br></li>
</ol>
<pre class="r"><code>rm(list = ls()) # this line cleans your Global Environment.
setwd("/Users/michellegonzalez/Documents/GitHub/Machine-Learning-for-Public-Policy") # set your working directory
# Do not forget to install a package with the install.packages() function if it's the first time you use it!
# install.packages(caTools, plotROC) # these guys are new for us
library(dplyr) # core package for dataframe manipulation. Usually installed and loaded with the tidyverse, but sometimes needs to be loaded in conjunction to avoid warnings.
library(tidyverse) # a large collection of packages for data manipulation and visualisation.
library(caret) # a library with key functions that streamline the process for predictive modelling
library(skimr) # a package with a set of functions to describe dataframes and more
library(plyr) # a package for data wrangling
library(caTools) # a library with several basic utility functions (e.g. ROC curves, LogitBoos classifier, etc.)
library(plotROC) # a companion to ggplot2 (loaded with the tidyverse) for plotting ROC curves
data_malawi <- read_csv("malawi.csv") # the file is directly read from the working directory/folder previously set</code></pre>
<h3>
<ol start="2" style="list-style-type: decimal">
<li>Data pre-processing
</h3>
<br></li>
</ol>
<p>This section will not be a thorough step-by-step of the
pre-processing and visualisation of our data because we have already
done that. However, we have to do something very important: recover a
static variable from the original dataset that contains a single number:
the poverty line in Malawi.</p>
<p><strong>Feature selection: subsetting the dataset </strong></p>
<p>The variable that we’re interested in recovering is
<strong>lnzline</strong>. The code below reproduces the dataframe
subsetting from our previous exercise. Except, this time we will NOT
delete de static vector lnzline.</p>
<pre class="r"><code># object:vector that contains the names of the variables that we want to get rid of (notice this time lnzline is still there)
cols <- c("ea", "EA", "psu","hhwght", "strataid", "case_id","eatype")
# subset of the data_malawi object:datframe
data_malawi <- data_malawi[,-which(colnames(data_malawi) %in% cols)] # the minus sign indicates deletion of cols
colnames(data_malawi) # print the names of the remaining vectors in our dataframe</code></pre>
<pre><code>## [1] "lnexp_pc_month" "hhsize" "hhsize2" "agehead"
## [5] "agehead2" "north" "central" "rural"
## [9] "nevermarried" "sharenoedu" "shareread" "nrooms"
## [13] "floor_cement" "electricity" "flushtoilet" "soap"
## [17] "bed" "bike" "musicplayer" "coffeetable"
## [21] "iron" "dimbagarden" "goats" "dependratio"
## [25] "hfem" "grassroof" "mortarpestle" "table"
## [29] "clock" "region" "lnzline"</code></pre>
<p><br></p>
<p>At this point, we still need to do two more pre-processing steps:
correctly define the vector/variable class in the dataframe, and create
the binary outcome/target variable. We will repeat the
class-transformation code chunk below so that you have all that is
needed in one section. However, we won’t spend time explaining it in
detail as that was done in the previous session.</p>
<pre class="r"><code># transform all binary/categorical data into factor class
min_count <- 3 # vector: 3 categories is our max number of categories found
# store boolean (true/false) if the number of unique values is lower or equal to the min_count vector
n_distinct2 <- apply(data_malawi, 2, function(x) length(unique(x))) <= min_count
# select the identified categorical variables and transform them into factors
data_malawi[n_distinct2] <- lapply(data_malawi[n_distinct2], factor)
# recall poverty line contains 1 unique value (it is static), let's transform the variable into numeric again
data_malawi$lnzline <- as.numeric(as.character(data_malawi$lnzline))
# you can use ``skim(data_malawi)'' to check that the dataframe is in working order</code></pre>
<p><br></p>
<p><strong>Feature creation: create a binary variable</strong></p>
<p><br></p>
<pre class="r"><code># print summary statistics of target variable
summary(data_malawi$lnexp_pc_month)</code></pre>
<pre><code>## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.777 6.893 7.305 7.359 7.758 11.064</code></pre>
<pre class="r"><code># if the log of per capita expenditure is below the estimated poverty line, classify individual as poor, else classify individual as not poor. Store as factor (default with text is class character)
data_malawi$poor <- as.factor(ifelse(data_malawi$lnexp_pc_month<= data_malawi$lnzline,"Y","N")) # Y(es) N(o)
# make sure that the factor target variable has poor = Y as reference category (this step is important when running the logistic regression)
data_malawi$poor <- relevel(data_malawi$poor, ref="Y") # make Y reference category
# print a proportions table to get a first impression of the state of poverty in Malawi
prop.table(table(data_malawi$poor))</code></pre>
<pre><code>##
## Y N
## 0.65 0.35</code></pre>
<p>According to our sample, about 65% of Malawians are considered poor.
This number is not unreasonable. According to The World Bank’s <a
href="https://databankfiles.worldbank.org/public/ddpext_download/poverty/987B9C90-CB9F-4D93-AE8C-750588BF00QA/current/Global_POVEQ_MWI.pdf">Country
Report</a> for Malawi, ca. <span class="math inline">\(70\%\)</span> of
the population lives with under <span
class="math inline">\(\$2.15\)</span> a day, and the poverty rate is
estimated to be at <span class="math inline">\(50\%\)</span>. About half
of their population is labelled as poor. These estimates were done with
<span class="math inline">\(2019\)</span> data (so, a bit more recent
than our sample).</p>
<p><br></p>
<pre class="r"><code># Final data pre-processing: delete static variable (poverty line)
# and along with it: remove the continuous target (as it perfectly predicts the binary target)
which(colnames(data_malawi)=="lnzline") # returns column number 31</code></pre>
<pre><code>## [1] 31</code></pre>
<pre class="r"><code>which(colnames(data_malawi)=="lnexp_pc_month") # returns column number 1</code></pre>
<pre><code>## [1] 1</code></pre>
<pre class="r"><code>data_malawi <- data_malawi[,-c(1,31)] # delete columns no. 1 and 31 from the dataset</code></pre>
<br>
<h3>
<ol start="3" style="list-style-type: decimal">
<li>Model Validation
</h3>
<br></li>
</ol>
<p>Let’s use a simple 80:20 split of our data. We will use the caret
package again.</p>
<pre class="r"><code>set.seed(1234) # ensures reproducibility of our data split
# data partitioning: train and test datasets
train_idx <- createDataPartition(data_malawi$poor, p = .8, list = FALSE, times = 1)
Train_df <- data_malawi[ train_idx,]
Test_df <- data_malawi[-train_idx,]</code></pre>
<p><br> Now, let’s fit a logistic model: <br></p>
<pre class="r"><code># Step 1: create trainControl object
TrControl <- trainControl(
method = "cv",
number = 5,
summaryFunction = twoClassSummary,
classProbs = TRUE, # IMPORTANT!
verboseIter = FALSE
)</code></pre>
<p>We’re going to pass the TrControl object onto the caret model
estimation to ask for the following: - cross-validate with 5 folds -
show model summary: performance metrics for when we have two distinct
classes (binary outcome), including the area under the ROC curve, the
sensitivity and specificity. - the ROC curve is based on the predicted
class probabilities, so the classProbs = TRUE parameter must accompany a
twoClassSummary setup. - veboseIter = TRUE shows you the output for each
iteration (but we don’t want to display all the details atm).</p>
<pre class="r"><code># Step 2: train the model.
set.seed(12345)
m <- train(
poor ~ .,
Train_df,
method = "glm",
family="binomial",
trControl = TrControl,
preProcess=c("center", "scale")
)</code></pre>
<pre><code>## Warning in train.default(x, y, weights = w, ...): The metric "Accuracy" was not
## in the result set. ROC will be used instead.</code></pre>
<p><br> Notice the warning. If we want to report the “Accuracy” metric,
we should remove the twoClassSummary parameter specification in the
TrControl object. <br></p>
<pre class="r"><code># print the model's performance metrics
print(m) </code></pre>
<pre><code>## Generalized Linear Model
##
## 9025 samples
## 29 predictor
## 2 classes: 'Y', 'N'
##
## Pre-processing: centered (30), scaled (30)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 7219, 7220, 7220, 7221, 7220
## Resampling results:
##
## ROC Sens Spec
## 0.8825702 0.8943044 0.6774223</code></pre>
<p><strong>Performance metrics (on the train set!)</strong> <br></p>
<ul>
<li><p><strong>ROC:</strong> it is a probability curve plotted with the
True Positive Rate (y-axis) against the False Positive Rate (x-axis);
you can think of it as plotting the tradeoff between maximising the true
positive rate and minimising the false positive rate. The preferred area
under the curve is <span class="math inline">\(1\)</span>. Our estimated
<span class="math inline">\(0.88\)</span> score indicates that a
logistic classification is a good fit for our data (close to <span
class="math inline">\(1\)</span>).</p></li>
<li><p><strong>Sensitivity:</strong> it is a measure of the proportion
of the positive (<span class="math inline">\(1\)</span> = poor) values
that are correctly identified. Therefore, we have correctly identified
<span class="math inline">\(89\%\)</span> of the actual positives. The
formula is: <span class="math inline">\(\frac{tp}{tp + fn}\)</span>;
where tp = true positive and fn = false negative. In the video-lecture,
Stephan used the term <strong>Recall</strong>, where we now use
sensitivity. This means that our model does pretty well at
predicting/identifying people living below the poverty line in
Malawi!</p></li>
<li><p><strong>Specificity:</strong> measures the proportion of actual
negatives that are correctly identified by the model; i.e. the ability
of our model to predict if an observation doesn’t belong to a certain
category. The formula is: <span class="math inline">\(\frac{tn}{tn +
fp}\)</span>; where tn = true negative and fp = false positive. At <span
class="math inline">\(67\%\)</span>, we can trust a predicted negative
(<span class="math inline">\(0\)</span>) value to be real more than half
the time. Our model is not as good at predicting who doesn’t live below
the poverty line in Malawi.</p></li>
</ul>
<p><br> The performance metrics we have interpreted above are based on
the training dataset only. We are interested in our model’s ability to
make out-of-sample predictions. Therefore, we will use the definitions
above, but to take the scores on the test-dataset predictions to make
our final evaluation. <br></p>
<p><strong>Out-of-sample performance</strong> <br></p>
<p>Notice that we have used cross-validation in our training dataset. In
theory, our performance metrics have been validated in 5 different
folds. To a certain extent, that means that our performance metrics
above did reflect the model’s ability to extrapolate. Nevertheless, we
will still see how our trained model performs in our test dataset. You
can think of this step as predicting on a sixth fold. We know that the
performance of a logistic classification model on the train set is
relatively good, is it the same for the test dataset?</p>
<pre class="r"><code># First, use the logistic classification model (trained on the Train_df) to make predictions on the test dataset:
set.seed(12345)
pr1 <- predict(m, Test_df, type = "raw")
head(pr1) # Yes and No output</code></pre>
<pre><code>## [1] N Y Y Y Y Y
## Levels: Y N</code></pre>
<p>We have specified the type of prediction we want: raw. This will
return the predicted classification (<span
class="math inline">\(0\)</span> or <span
class="math inline">\(1\)</span>) as opposed to the individual’s
probability of falling into the selected category <span
class="math inline">\(1\)</span> (or the estimated probability of being
poor). There is a rule of thumb that says you will be categorised as
poor (or any chosen category) if your estimated probability is >= to
<span class="math inline">\(0.5\)</span>. With this information, we can
create a Confusion Matrix which will be accompanied by performance
metrics.</p>
<pre class="r"><code># Next, we call the caret package's confusionMatrix function, and select the two elements to be contrasted:
# the predicted classification vector, and the actual observed vector from the test dataframe.
confusionMatrix(pr1, Test_df[["poor"]], positive = "Y") # positive = "Y" indicates that our category of interest is Y (1)</code></pre>
<pre><code>## Confusion Matrix and Statistics
##
## Reference
## Prediction Y N
## Y 1294 270
## N 172 519
##
## Accuracy : 0.804
## 95% CI : (0.787, 0.8202)
## No Information Rate : 0.6501
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.5564
##
## Mcnemar's Test P-Value : 3.953e-06
##
## Sensitivity : 0.8827
## Specificity : 0.6578
## Pos Pred Value : 0.8274
## Neg Pred Value : 0.7511
## Prevalence : 0.6501
## Detection Rate : 0.5738
## Detection Prevalence : 0.6936
## Balanced Accuracy : 0.7702
##
## 'Positive' Class : Y
## </code></pre>
<p><br></p>
<p>The first element from the above function returns the confusion
matrix, a 2×2 table that shows the predicted values from the model
vs. the actual values from the test dataset. You may be acquainted with
this sort of table, but know it as a cross-tabulation. From the
confusion matrix, we obtain the information that we need to estimate
some performance metrics. If you need a reminder of what each cell in
the 2x2 matrix represents, recall that the structure of our target
variable is [<span class="math inline">\(Y(1),N(0)\)</span>]. Therefore,
the first cell would be the intersection of Predicted <span
class="math inline">\(Y\)</span> vs Observed <span
class="math inline">\(Y\)</span> (or True Positive) = <span
class="math inline">\(1294\)</span>, the fourth cell would be the
intersection of Predicted <span class="math inline">\(N\)</span> vs
Observed <span class="math inline">\(N\)</span> (or True Negative) =
<span class="math inline">\(519\)</span>. These guys are the predictions
that have hit the mark! On the other hand, the second cell would be the
intersection of Predicted <span class="math inline">\(Y\)</span> vs
Observed <span class="math inline">\(N\)</span> (or False Positive) =
<span class="math inline">\(270\)</span>, and the third cell Predicted
<span class="math inline">\(N\)</span> vs Observed <span
class="math inline">\(Y\)</span> (or False Negative) = <span
class="math inline">\(172\)</span>. These were incorrect predictions. We
use these counts (true positives, true negatives, false positives, false
negatives) to estimate performance metrics <br></p>
<p>Besides the performance metrics discussed previously, this function
also shows the Accuracy of our model (or <span
class="math inline">\(1\)</span> - the error rate) which, at <span
class="math inline">\(0.8\)</span>, indicates that our classification
algorithm is highly accurate.</p>
<p><br></p>
<pre><code>**Imbalanced data**
When you have a large number of zeros (or No, in this case), the Accuracy metric may not be the most reliable one. If we look at the formula: number of correct predictions / total number of predictions, we see why this might be an issue. It is a lot easier to correctly predict that of which there is plenty of (Yes), than the category for which we have less instances. </code></pre>
<p><br></p>
<p>Imbalance is not a problem for our target variable, as we have
roughly as many zeros as ones. In fact, we have more Yes (1) responses.
Nonetheless, this sets the stage for us to introduce the Kappa statistic
(<span class="math inline">\(0.55\)</span>), which is a measure of model
accuracy that is adjusted by accounting for the possibility of a correct
prediction by chance alone. It ranges from 0 to 1, and can be
interpreted using the following thresholds:</p>
<ul>
<li><p>Poor = Less than 0.20</p></li>
<li><p>Fair = 0.20 to 0.40</p></li>
<li><p>Moderate = 0.40 to 0.60</p></li>
<li><p>Good = 0.60 to 0.80</p></li>
<li><p>Very good = 0.80 to 1.00</p></li>
</ul>
<p>At <span class="math inline">\(0.55\)</span>, our classification
model performs moderately well. Finally, Sensitivity and Specificity
scores on the test dataset are very close to the ones obtained from the
train dataset. This is a good sign for the out-of-sample stability of
our model. <br></p>
<p><strong>Model Visualisation</strong> <br></p>
<p>We can visualise the performance of our classification model in
various ways. For now, we’ll focus on a simple ROC AUC.</p>
<pre class="r"><code># ROC AUC: Area Under the Curve
# colAUC function from the caTools library
# transform predicted values and observed values into class numeric (needed for the colAUC function)
pr_numeric <- as.numeric(as.factor(pr1))
# sanity check:
head(cbind(pr_numeric, pr1)) # the numeric values of both vectors seem to be the same. </code></pre>
<pre><code>## pr_numeric pr1
## [1,] 2 2
## [2,] 1 1
## [3,] 1 1
## [4,] 1 1
## [5,] 1 1
## [6,] 1 1</code></pre>
<pre class="r"><code>poor_numeric <- as.numeric(as.factor(Test_df$poor))
# sanity check
head(cbind(poor_numeric,Test_df$poor)) # all good </code></pre>
<pre><code>## poor_numeric
## [1,] 1 1
## [2,] 1 1
## [3,] 1 1
## [4,] 1 1
## [5,] 1 1
## [6,] 1 1</code></pre>
<pre class="r"><code># plot the ROC area under the curve
colAUC(pr_numeric, poor_numeric, plotROC = TRUE)</code></pre>
<p><img src="classification_files/figure-html/unnamed-chunk-13-1.png" width="672" /></p>
<pre><code>## [,1]
## 1 vs. 2 0.7702343</code></pre>
<pre class="r"><code># We can also plot the ROC AUC with ggplot
# First, we create a dataframe containing the observed and the predicted values (in numeric form)
roc_df <- data.frame(Observed = poor_numeric, Predicted = pr_numeric)
# Second, we add the geom_roc() layer to a ggplot object
roc_gg <- ggplot(roc_df, aes (d = Observed, m = Predicted)) +
geom_roc(labels = FALSE, color='orange') +
style_roc(theme = theme_bw, guide = TRUE) # guide=TRUE adds a diagonal guideline, style_roc() adds minor grid lines, and optionally direct labels to ggplot objects containing a geom_roc layer
direct_label(roc_gg, nudge_y = 0.2) # direct_label tells you what the plotted line represents, nudge_y option places the label (you can play around with that number to see where different values place the label)</code></pre>
<pre><code>## Warning in verify_d(data$d): D not labeled 0/1, assuming 1 = 0 and 2 = 1!
## Warning in verify_d(data$d): D not labeled 0/1, assuming 1 = 0 and 2 = 1!</code></pre>
<p><img src="classification_files/figure-html/unnamed-chunk-13-2.png" width="672" /></p>
<p>Notice the warning on the ggplot2 ROC AUC plot. The assumption that
they are making is correct, so we do not need to do anything else at
this moment. You can check this by contrasting the values of the
labelled vs. the numeric vectors (use the head() function).</p>
<p><strong>Conclusion</strong> <br></p>
<p>How did our classification model do? Is a logistic regression the
right algorithm? I trust you can form your own judgement based on the
performance metrics above. Personally, I think we have improved from a
linear regression, but perhaps we can do better with Ensemble Learning
techniques!</p>
</div>
<div id="python-practical" class="section level3">
<h3><strong>Python practical</strong></h3>
<p>We will continue to work with the Malawi dataset, which can be
downloaded in the (Prediction Policy Problems)[<a
href="https://www.ml4publicpolicy.com/predictionpolicy.html"
class="uri">https://www.ml4publicpolicy.com/predictionpolicy.html</a>]
tab of this website.</p>
<h3>
<ol style="list-style-type: decimal">
<li>Preliminaries: libraries and data upload
</h3>
<br></li>
</ol>
<pre class="python"><code>#==== Python version: 3.10.12 ====#
# Opening libraries
import sklearn as sk # our trusted Machine Learning library
from sklearn.model_selection import train_test_split # split the dataset into train and test
from sklearn.model_selection import cross_val_score # to obtain the cross-validation score
from sklearn.model_selection import cross_validate # to perform cross-validation
from sklearn.linear_model import LogisticRegression # computation of logistic model for classification
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import confusion_matrix # returns a confusion matrix
from sklearn.metrics import roc_curve, roc_auc_score, auc # roc area under the curve
from sklearn.metrics import accuracy_score # performance metric for classification models
from sklearn.metrics import classification_report # general report for classification model
# Non-ML libraries
import random # for random state
import csv # a library to read and write csv files
import numpy as np # a library for handling
import pandas as pd # a library to help us easily navigate and manipulate dataframes
import seaborn as sns # a data visualisation library
import matplotlib.pyplot as plt # a data visualisation library
import statsmodels.api as sm # computation of linear and logistic regressions
# Uploading data
malawi = pd.read_csv('/Users/michellegonzalez/Documents/GitHub/Machine-Learning-for-Public-Policy/malawi.csv')</code></pre>
<br>
<h3>
<ol start="2" style="list-style-type: decimal">
<li>Data pre-processing
</h3>
<br></li>
</ol>
<p>This section will not be a thorough step-by-step of the
pre-processing and visualisation of our data because we have already
done that in the previous session. However, we have to do something very
important: recover a static variable from the original dataset that
contains a single number: the poverty line in Malawi. <br></p>
<p><strong>Feature selection: subsetting the dataset </strong> <br></p>
<p>The variable that we’re interested in recovering is
<strong>lnzline</strong>. The code below reproduces the dataframe
subsetting from our previous exercise. Except, this time we will make
sure to store the lnzline vector for future use.</p>
<p><br></p>
<pre class="python"><code># Retrieve information from the poverty line static vector:
malawi['lnzline'].describe()</code></pre>
<pre><code>## count 1.128000e+04
## mean 7.555000e+00
## std 1.776436e-15
## min 7.555000e+00
## 25% 7.555000e+00
## 50% 7.555000e+00
## 75% 7.555000e+00
## max 7.555000e+00
## Name: lnzline, dtype: float64</code></pre>
<pre class="python"><code># Let's store this in an object outside of the dataframe and get rid of it in the dataframe
# Static variables have zero-variance, and we already know zero-variance predictors are troublesome..
lnzline = malawi['lnzline']
# Instead of deleting case_id, we will set it as an index (we did not do this last time!).
# This is essentially using the case_id variable as row names (and won't be included in your ML model)
malawi = malawi.set_index('case_id')
# sanity check: notice that case_id is there, but doesn't have a corresponding column (it is now only considered a row name)
malawi.head()</code></pre>
<pre><code>## lnexp_pc_month hhsize hhsize2 ... psu strataid lnzline
## case_id ...
## 10101002025 6.900896 7 49 ... 10101002 1 7.555
## 10101002051 7.064378 3 9 ... 10101002 1 7.555
## 10101002072 6.823851 6 36 ... 10101002 1 7.555
## 10101002079 6.894722 6 36 ... 10101002 1 7.555
## 10101002095 6.465989 6 36 ... 10101002 1 7.555
##
## [5 rows x 37 columns]</code></pre>
<pre class="python"><code># deleting variables from pandas dataframe
cols2delete = ['ea', 'EA', 'hhwght', 'psu', 'strataid', 'lnzline', 'eatype', 'region']
malawi = malawi.drop(cols2delete,axis=1) # axis=0 means delete rows and axis=1 means delete columns
# check if we have deleted the columns: we originally had 37 variables, we now should have 29
print(malawi.shape)</code></pre>
<pre><code>## (11280, 29)</code></pre>
<p>At this point, we still need to do two more pre-processing steps:
correctly define the vector/variable class in the dataframe, and create
the binary outcome/target variable. We will repeat the
class-transformation code chunk below so that you have all that is
needed in one section. However, we won’t spend time explaining it in
detail as that was done in the previous session.</p>
<p><br></p>
<p><strong>Feature creation: create a binary variable</strong></p>
<pre class="python"><code>#==== Correctly identify each vector type: ====#
# for-loop that iterates over variables in dataframe, if they have 2 unique values, transform vector into categorical
for column in malawi:
if malawi[column].nunique() == 2:
malawi[column] = pd.Categorical(malawi[column])
#==== Create a binary target variable: ====#
# print summary statistics of target variable
malawi['lnexp_pc_month'].describe()</code></pre>
<pre><code>## count 11280.000000
## mean 7.358888
## std 0.675346
## min 4.776855
## 25% 6.892941
## 50% 7.305191
## 75% 7.757587
## max 11.063562
## Name: lnexp_pc_month, dtype: float64</code></pre>
<pre class="python"><code># if the log of per capita expenditure is below the estimated poverty line, classify individual as poor, else classify individual as not poor. Store as factor (default with text is class character)
print(lnzline[0]) # copy the static number 7.555 to use as threshold / we're printing the first row of a static vector --- i.e. all rows contain the same number</code></pre>
<pre><code>## 7.5549998</code></pre>
<pre class="python"><code># use numPy to create a binary vector (notice we have rounded up the threshold)
malawi['Poor'] = np.where(malawi['lnexp_pc_month'] <= 7.555, 1, 0)
# sanity check
malawi['Poor'].describe() # returns binary [0,1] float vector, let's turn this into a categorical vector</code></pre>
<pre><code>## count 11280.000000
## mean 0.650000
## std 0.476991
## min 0.000000
## 25% 0.000000
## 50% 1.000000
## 75% 1.000000
## max 1.000000
## Name: Poor, dtype: float64</code></pre>
<pre class="python"><code>malawi['Poor'] = pd.Categorical(malawi['Poor']) # use malawi['Poor'].info() if you want to see the transformation
# alternatively...
# malawi['Poor'] = (malawi['lnexp_pc_month'] <= 7.555).astype(bool) # anlso numpy, but directly specifying boolean type (true (poor) /false (not poor))
# print a proportions table to get a first impression of the state of poverty in Malawi
malawi['Poor'].value_counts(normalize=True)</code></pre>
<pre><code>## 1 0.65
## 0 0.35
## Name: Poor, dtype: float64</code></pre>
<p><br></p>
<p>According to our sample, about 65% of Malawians are considered poor.
This number is not unreasonable. According to The World Bank’s <a
href="https://databankfiles.worldbank.org/public/ddpext_download/poverty/987B9C90-CB9F-4D93-AE8C-750588BF00QA/current/Global_POVEQ_MWI.pdf">Country
Report</a> for Malawi, ca. <span class="math inline">\(70\%\)</span> of
the population lives with under <span
class="math inline">\(\$2.15\)</span> a day, and the poverty rate is
estimated to be at <span class="math inline">\(50\%\)</span>. About half
of their population is labelled as poor. These estimates were done with
<span class="math inline">\(2019\)</span> data (so, a bit more recent
than our sample).</p>
<p><br></p>
<pre class="python"><code>
# Final data pre-processing: remove the continuous target (as it perfectly predicts the binary target in a non-informative way)
cont_target = ['lnexp_pc_month']
malawi = malawi.drop(cont_target , axis=1)
</code></pre>
<br>
<h3>
<ol start="3" style="list-style-type: decimal">
<li>Model Validation
</h3>
<br></li>
</ol>
<p>Let’s use a simple 80:20 split of our data.</p>
<pre class="python"><code># First, recall the df structure
malawi.info() # returns the column number, e.g. hhsize = column number 0, hhsize2 = 1... etc.</code></pre>
<pre><code>## <class 'pandas.core.frame.DataFrame'>
## Int64Index: 11280 entries, 10101002025 to 31202086374
## Data columns (total 29 columns):
## # Column Non-Null Count Dtype
## --- ------ -------------- -----
## 0 hhsize 11280 non-null int64
## 1 hhsize2 11280 non-null int64
## 2 agehead 11280 non-null int64
## 3 agehead2 11280 non-null int64
## 4 north 11280 non-null category
## 5 central 11280 non-null category
## 6 rural 11280 non-null category
## 7 nevermarried 11280 non-null category
## 8 sharenoedu 11280 non-null float64
## 9 shareread 11280 non-null float64
## 10 nrooms 11280 non-null int64
## 11 floor_cement 11280 non-null category
## 12 electricity 11280 non-null category
## 13 flushtoilet 11280 non-null category
## 14 soap 11280 non-null category
## 15 bed 11280 non-null category
## 16 bike 11280 non-null category
## 17 musicplayer 11280 non-null category
## 18 coffeetable 11280 non-null category
## 19 iron 11280 non-null category
## 20 dimbagarden 11280 non-null category
## 21 goats 11280 non-null category
## 22 dependratio 11280 non-null float64
## 23 hfem 11280 non-null category
## 24 grassroof 11280 non-null category
## 25 mortarpestle 11280 non-null category
## 26 table 11280 non-null category
## 27 clock 11280 non-null category
## 28 Poor 11280 non-null category
## dtypes: category(21), float64(3), int64(5)
## memory usage: 1.0 MB</code></pre>
<pre class="python"><code>
# Then, split!
X = malawi.iloc[:, 0:27] # x is a matrix containing all variables except the last one, which conveniently is our binary target variable
y = malawi.iloc[:, 28] # y is a vector containing our target variable
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=12345) # random_state is for reproducibility purposes</code></pre>
<p><br></p>
<p><strong>Fit a logistic model</strong> <br></p>
<pre class="python"><code>
#==== Create a Logistic Regression Object and Fit the Model ====#
# We are going to add a random_state (reproducibility) and increase the number of iterations from the default (100) to 1000
# We have also stated that we want to fit the intercept and have no penalty (remember the penalisation parameter from the Lasso regression? sklearn has the option of including it in the Logistic Regression as well).
m = LogisticRegression(random_state=12345, max_iter=1000, fit_intercept = True, penalty
= None).fit(X_train, y_train) </code></pre>
<p><br></p>
<p>We have successfully fit a logistic classification model. A
limitation of the sklearn python library is that we cannot easily access
the output of the model. Instead, we immediately estimate performance