<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/chaid-vs-ranger-vs-xgboost-a-comparison/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://datascienceplus.com/chaid-vs-ranger-vs-xgboost-a-comparison/"> R Programming – DataScience+</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>In an earlier <a rel="nofollow" href="https://ibecav.github.io/chaidtutor1/" target="_blank">post</a>, I focused on an in-depth visit with CHAID (Chi-square automatic interaction detection). Quoting myself, I said “As the name implies it is fundamentally based on the venerable Chi-square test – and while not the most powerful (in terms of detecting the smallest possible differences) or the fastest, it really is easy to manage and more importantly to tell the story after using it”. In this post I'll spend a little time comparing CHAID with a random forest algorithm in the <code>ranger</code> library and with a gradient boosting algorithm via the <code>xgboost</code> library. I'll use the exact same data set for all three so we can draw some easy comparisons about their speed and their accuracy.</p>
<p>I do believe CHAID is a great choice for some sets of data and some circumstances but I'm interested in some empirical information, so off we go. </p>
<h2>Setup and library loading</h2>
<p>If you've never used <code>CHAID</code> before you may also not have <code>partykit</code>. <code>CHAID</code> isn't on <code>CRAN</code> but I have provided the commented out install command below. <code>ranger</code> and <code>xgboost</code> are available from CRAN and are straightforward to install. You'll also get a variety of messages, none of which is relevant to this example so I've suppressed them.</p>
<pre><code># install.packages(&quot;partykit&quot;)
# install.packages(&quot;CHAID&quot;, repos=&quot;http://R-Forge.R-project.org&quot;)
# install.packages(&quot;ranger&quot;)
# install.packages(&quot;xgboost&quot;)
require(dplyr)
require(tidyr)
require(ggplot2)
require(CHAID)
require(purrr)
require(caret)
require(ranger)
require(xgboost)
require(kableExtra) # just to make the output nicer
theme_set(theme_bw()) # set theme for ggplot2
</code></pre>
<h2>Predicting customer churn for a fictional TELCO company</h2>
<p>We're going to use a dataset that comes to us from the <a rel="nofollow" href="https://www.ibm.com/communities/analytics/watson-analytics-blog/using-customer-behavior-data-to-improve-customer-retention/" target="_blank">IBM Watson Project</a>. It's a very practical example and an understandable dataset. A great use case for the algorithms we'll be using. Imagine yourself in a fictional company faced with the task of trying to predict which customers are going to leave your business for another provider a.k.a. churn. Obviously we'd like to be able to predict this phenomenon and potentially target these customers for retention or just better project our revenue. Being able to predict churn even a little bit better could save us lots of money, especially if we can identify the key indicators and influence them.</p>
<p>In the original <a rel="nofollow" href="https://ibecav.github.io/chaidtutor1/" target="_blank">posting</a> I spent a great deal of time explaining the mechanics of loading and prepping the data. This time we'll do that quickly and efficiently and if you need an explanation of what's going on please refer back. I've embedded some comments in the code where I think they'll be most helpful. First we'll grab the data from the IBM site using <code>read.csv</code>, in this case I'm happy to let it tag most of our variables as <code>factors</code> since that's what we'll want for our <code>CHAID</code> work. </p>
<pre><code>set.seed(2018)
churn &lt;- read.csv(&quot;https://community.watsonanalytics.com/wp-content/uploads/2015/03/WA_Fn-UseC_-Telco-Customer-Churn.csv&quot;)
str(churn)</code>
<em>## 'data.frame': 7043 obs. of 21 variables:
## $ customerID : Factor w/ 7043 levels &quot;0002-ORFBO&quot;,&quot;0003-MKNFE&quot;,..: 5376 3963 2565 5536 6512 6552 1003 4771 5605 4535 ...
## $ gender : Factor w/ 2 levels &quot;Female&quot;,&quot;Male&quot;: 1 2 2 2 1 1 2 1 1 2 ...
## $ SeniorCitizen : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Partner : Factor w/ 2 levels &quot;No&quot;,&quot;Yes&quot;: 2 1 1 1 1 1 1 1 2 1 ...
## $ Dependents : Factor w/ 2 levels &quot;No&quot;,&quot;Yes&quot;: 1 1 1 1 1 1 2 1 1 2 ...
## $ tenure : int 1 34 2 45 2 8 22 10 28 62 ...
## $ PhoneService : Factor w/ 2 levels &quot;No&quot;,&quot;Yes&quot;: 1 2 2 1 2 2 2 1 2 2 ...
## $ MultipleLines : Factor w/ 3 levels &quot;No&quot;,&quot;No phone service&quot;,..: 2 1 1 2 1 3 3 2 3 1 ...
## $ InternetService : Factor w/ 3 levels &quot;DSL&quot;,&quot;Fiber optic&quot;,..: 1 1 1 1 2 2 2 1 2 1 ...
## $ OnlineSecurity : Factor w/ 3 levels &quot;No&quot;,&quot;No internet service&quot;,..: 1 3 3 3 1 1 1 3 1 3 ...
## $ OnlineBackup : Factor w/ 3 levels &quot;No&quot;,&quot;No internet service&quot;,..: 3 1 3 1 1 1 3 1 1 3 ...
## $ DeviceProtection: Factor w/ 3 levels &quot;No&quot;,&quot;No internet service&quot;,..: 1 3 1 3 1 3 1 1 3 1 ...
## $ TechSupport : Factor w/ 3 levels &quot;No&quot;,&quot;No internet service&quot;,..: 1 1 1 3 1 1 1 1 3 1 ...
## $ StreamingTV : Factor w/ 3 levels &quot;No&quot;,&quot;No internet service&quot;,..: 1 1 1 1 1 3 3 1 3 1 ...
## $ StreamingMovies : Factor w/ 3 levels &quot;No&quot;,&quot;No internet service&quot;,..: 1 1 1 1 1 3 1 1 3 1 ...
## $ Contract : Factor w/ 3 levels &quot;Month-to-month&quot;,..: 1 2 1 2 1 1 1 1 1 2 ...
## $ PaperlessBilling: Factor w/ 2 levels &quot;No&quot;,&quot;Yes&quot;: 2 1 2 1 2 2 2 1 2 1 ...
## $ PaymentMethod : Factor w/ 4 levels &quot;Bank transfer (automatic)&quot;,..: 3 4 4 1 3 3 2 4 3 1 ...
## $ MonthlyCharges : num 29.9 57 53.9 42.3 70.7 ...
## $ TotalCharges : num 29.9 1889.5 108.2 1840.8 151.7 ...
## $ Churn : Factor w/ 2 levels &quot;No&quot;,&quot;Yes&quot;: 1 1 2 1 2 2 1 1 2 1 ...
</em></pre>
<p>We have data on 7,043 customers across 21 variables. <code>customerID</code> can't really be a predictor but we will use it in a little bit. <code>Churn</code> is what we want to predict so we have 19 potential predictor variables to work with. Four of them were not automatically converted to factors so we'll have to look into them for <code>CHAID</code>. For a review of what the output means and how CHAID works please refer <a rel="nofollow" href="https://ibecav.github.io/chaidtutor1/" target="_blank">back</a>.</p>
<p>Let's address the easiest thing first. <code>SeniorCitizen</code> is coded zero and one instead of yes/no so let's recode that in a nice conservative fashion and see what the breakdown is.</p>
<pre><code># Fix senior citizen status
churn$SeniorCitizen &lt;- recode_factor(
  churn$SeniorCitizen,
  `0` = &quot;No&quot;,
  `1` = &quot;Yes&quot;,
  .default = &quot;Should not happen&quot;
)
summary(churn$SeniorCitizen)</code>
<em>## No Yes ## 5901 1142
</em></pre>
<p>We have three variables left that are numeric, now that we have addressed senior citizen status. Let's use a combination of <code>dplyr</code> and <code>ggplot2</code> to see what the distribution looks like using a density plot.</p>
<pre><code>churn %>%
   select_if(is.numeric) %>%
   gather(metric, value) %>%
   ggplot(aes(value, fill = metric)) +
   geom_density(show.legend = FALSE) +
   facet_wrap( ~ metric, scales = &quot;free&quot;)</code>
<em>## Warning: Removed 11 rows containing non-finite values (stat_density).
</em></pre>
<p>The plot:<br>
<a rel="nofollow" href="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx3-1.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx3-1-490x218.png?w=450&ssl=1" alt="" class="wp-image-15722" data-recalc-dims="1"></a></p>
<p>Well those aren't the most normal looking distributions and we have this message <code>## Warning: Removed 11 rows containing non-finite values (stat_density).</code> which alerts us to the fact that there are some missing values in our data. Let's first figure out where the missing data is:</p>
<pre><code>churn %>%
  select_if(anyNA) %>% summary</code>
<em>## TotalCharges ## Min. : 18.8 ## 1st Qu.: 401.4 ## Median :1397.5 ## Mean :2283.3 ## 3rd Qu.:3794.7 ## Max. :8684.8 ## NA's :11
</em></pre>
<p>Now we know that total customer charges are missing 11 entries. Our three algorithms vary as to how gracefully they handle missing values but at this point, we have several options including:</p>
<ul>
<li>Eliminate the entire customer record if anything is missing</li>
<li>Impute or substitute in some reasonable value like the mean or the median for missing values</li>
<li>Do some fancier imputation to make sure we substitute in the most plausible value for <code>TotalCharges</code> </li>
</ul>
<p>Elimination is easy, efficient, and conservative and since it is a very small percentage of our total data set unlikely to cost us a lot of information for the models that don't handle missing values well. But for purposes of this blog post and to help demonstrate some of the capabilities within <code>caret</code> (since we're going to use it anyway) we'll try median and knn (k nearest neighbor) imputation.</p>
<p>First let's make a vector that contains the <code>customerID</code> numbers of the eleven cases in question.</p>
<pre><code>xxx &lt;- churn %>%
         filter_all(any_vars(is.na(.))) %>% 
         select(customerID)
xxx &lt;- as.vector(xxx$customerID)
xxx</code>
<em>## [1] &quot;4472-LVYGI&quot; &quot;3115-CZMZD&quot; &quot;5709-LVOEQ&quot; &quot;4367-NUYAO&quot; &quot;1371-DWPAZ&quot; &quot;7644-OMVMY&quot; &quot;3213-VVOLG&quot; &quot;2520-SGTTA&quot; &quot;2923-ARZLG&quot; &quot;4075-WKNIU&quot; &quot;2775-SEFEE&quot;
</em></pre>
<pre><code>churn %>% filter(customerID %in% xxx)</code>
<em>## customerID gender SeniorCitizen Partner Dependents tenure PhoneService MultipleLines InternetService OnlineSecurity OnlineBackup DeviceProtection TechSupport StreamingTV StreamingMovies Contract PaperlessBilling PaymentMethod MonthlyCharges TotalCharges Churn
## 1 4472-LVYGI Female No Yes Yes 0 No No phone service DSL Yes No Yes Yes Yes No Two year Yes Bank transfer (automatic) 52.55 NA No
## 2 3115-CZMZD Male No No Yes 0 Yes No No No internet service No internet service No internet service No internet service No internet service No internet service Two year No Mailed check 20.25 NA No
## 3 5709-LVOEQ Female No Yes Yes 0 Yes No DSL Yes Yes Yes No Yes Yes Two year No Mailed check 80.85 NA No
## 4 4367-NUYAO Male No Yes Yes 0 Yes Yes No No internet service No internet service No internet service No internet service No internet service No internet service Two year No Mailed check 25.75 NA No
## 5 1371-DWPAZ Female No Yes Yes 0 No No phone service DSL Yes Yes Yes Yes Yes No Two year No Credit card (automatic) 56.05 NA No
## 6 7644-OMVMY Male No Yes Yes 0 Yes No No No internet service No internet service No internet service No internet service No internet service No internet service Two year No Mailed check 19.85 NA No
## 7 3213-VVOLG Male No Yes Yes 0 Yes Yes No No internet service No internet service No internet service No internet service No internet service No internet service Two year No Mailed check 25.35 NA No
## 8 2520-SGTTA Female No Yes Yes 0 Yes No No No internet service No internet service No internet service No internet service No internet service No internet service Two year No Mailed check 20.00 NA No
## 9 2923-ARZLG Male No Yes Yes 0 Yes No No No internet service No internet service No internet service No internet service No internet service No internet service One year Yes Mailed check 19.70 NA No
## 10 4075-WKNIU Female No Yes Yes 0 Yes Yes DSL No Yes Yes Yes Yes No Two year No Mailed check 73.35 NA No
## 11 2775-SEFEE Male No No Yes 0 Yes Yes DSL Yes Yes No Yes No No Two year Yes Bank transfer (automatic) 61.90 NA No
</em></pre>
<p>As you look at those eleven records it doesn't appear they are “average”! In particular, I'm worried that the <code>MonthlyCharges</code> look small and they have 0 <code>tenure</code> for this group. No way of knowing for certain but it could be that these are just the newest customers with very little time using our service. Let's use our list to do some comparing of these eleven versus the total population, that will help us decide what to do about the missing cases. Replacing with the median value is simple and easy but it may well not be the most accurate choice.</p>
<pre><code>churn %>% 
   filter(customerID %in% xxx) %>% 
   summarise(median(MonthlyCharges))</code>
<em>## median(MonthlyCharges)
## 1 25.75
</em></pre>
<pre><code>median(churn$MonthlyCharges, na.rm = TRUE)</code>
<em>## [1] 70.35
</em></pre>
<pre><code>churn %>% 
   filter(customerID %in% xxx) %>% 
   summarise(median(tenure))</code>
<em>## median(tenure)
## 1 0
</em></pre>
<pre><code>median(churn$tenure, na.rm = TRUE)</code>
<em>## [1] 29
</em></pre>
<p>The median <code>MonthlyCharges</code> are much lower and instead of two years or so of median tenure this group has none. Let's use the <code>preProcess</code> function in <code>caret</code> to accomplish several goals. We'll ask it to impute the missing values for us using both <code>knnImpute</code> (k nearest neighbors) and a pure median <code>medianImpute</code>. From the <code>?preProcess</code> help pages:</p>
<p>k-nearest neighbor imputation is carried out by finding the k closest samples (Euclidian distance) in the training set. Imputation via bagging fits a bagged tree model for each predictor (as a function of all the others). This method is simple, accurate and accepts missing values, but it has much higher computational cost. Imputation via medians takes the median of each predictor in the training set, and uses them to fill missing values. This method is simple, fast, and accepts missing values, but treats each predictor independently, and may be inaccurate.</p>
<p>We'll also have it transform our numeric variables using <code>YeoJohnson</code> and identify any predictor variables that have near zero variance <code>nzv</code>.</p>
<pre><code># using k nearest neighbors
pp_knn &lt;- preProcess(churn, method = c(&quot;knnImpute&quot;, &quot;YeoJohnson&quot;, &quot;nzv&quot;))
# simple output
pp_knn</code>
<em>## Created from 7032 samples and 21 variables
## ## Pre-processing:
## - centered (3)
## - ignored (18)
## - 5 nearest neighbor imputation (3)
## - scaled (3)
## - Yeo-Johnson transformation (3)
## ## Lambda estimates for Yeo-Johnson transformation:
## 0.45, 0.93, 0.25
</em></pre>
<pre><code># more verbose
pp_knn$method</code>
<em>## $knnImpute
## [1] &quot;tenure&quot; &quot;MonthlyCharges&quot; &quot;TotalCharges&quot; ## ## $YeoJohnson
## [1] &quot;tenure&quot; &quot;MonthlyCharges&quot; &quot;TotalCharges&quot; ## ## $ignore
## [1] &quot;customerID&quot; &quot;gender&quot; &quot;SeniorCitizen&quot; &quot;Partner&quot; &quot;Dependents&quot; &quot;PhoneService&quot; &quot;MultipleLines&quot; &quot;InternetService&quot; &quot;OnlineSecurity&quot; &quot;OnlineBackup&quot; &quot;DeviceProtection&quot; &quot;TechSupport&quot; &quot;StreamingTV&quot; &quot;StreamingMovies&quot; &quot;Contract&quot; &quot;PaperlessBilling&quot; &quot;PaymentMethod&quot; &quot;Churn&quot; ## ## $center
## [1] &quot;tenure&quot; &quot;MonthlyCharges&quot; &quot;TotalCharges&quot; ## ## $scale
## [1] &quot;tenure&quot; &quot;MonthlyCharges&quot; &quot;TotalCharges&quot;
</em></pre>
<pre><code># using medians
pp_median &lt;- preProcess(churn, method = c(&quot;medianImpute&quot;, &quot;YeoJohnson&quot;, &quot;nzv&quot;))
pp_median</code>
<em>## Created from 7032 samples and 21 variables
## ## Pre-processing:
## - ignored (18)
## - median imputation (3)
## - Yeo-Johnson transformation (3)
## ## Lambda estimates for Yeo-Johnson transformation:
## 0.45, 0.93, 0.25
</em></pre>
<pre><code>pp_median$method</code>
<em>## $medianImpute
## [1] &quot;tenure&quot; &quot;MonthlyCharges&quot; &quot;TotalCharges&quot; ## ## $YeoJohnson
## [1] &quot;tenure&quot; &quot;MonthlyCharges&quot; &quot;TotalCharges&quot; ## ## $ignore
## [1] &quot;customerID&quot; &quot;gender&quot; &quot;SeniorCitizen&quot; &quot;Partner&quot; &quot;Dependents&quot; &quot;PhoneService&quot; &quot;MultipleLines&quot; &quot;InternetService&quot; &quot;OnlineSecurity&quot; &quot;OnlineBackup&quot; &quot;DeviceProtection&quot; &quot;TechSupport&quot; &quot;StreamingTV&quot; &quot;StreamingMovies&quot; &quot;Contract&quot; &quot;PaperlessBilling&quot; &quot;PaymentMethod&quot; &quot;Churn&quot;
</em></pre>
<p>The <code>preProcess</code> function creates a list object of <code>class</code> preProcess that contains information about what needs to be done and what the results of the transformations will be, but we need to apply the <code>predict</code> function to actually make the changes proposed. So at this point let's create two new dataframes <code>nchurn1</code> and <code>nchurn2</code> that contain the data after the pre-processing has occurred. Then we can see how the results compare.</p>
<pre><code>nchurn1 &lt;- predict(pp_knn,churn)
nchurn2 &lt;- predict(pp_median,churn)
nchurn2 %>% 
   filter(customerID %in% xxx) %>% 
   summarise(median(TotalCharges))</code>
<em>## median(TotalCharges)
## 1 20.79526
</em></pre>
<pre><code>median(nchurn2$TotalCharges, na.rm = TRUE)</code>
<em>## [1] 20.79526
</em></pre>
<pre><code>nchurn1 %>% 
   filter(customerID %in% xxx) %>% 
   summarise(median(TotalCharges))</code>
<em>## median(TotalCharges)
## 1 -1.849681
</em></pre>
<pre><code>median(nchurn1$TotalCharges, na.rm = TRUE)</code>
<em>## [1] 0.01820494
</em></pre>
<p>May also be useful to visualize the data as we did earlier to see how the transformations have changed the density plots.</p>
<pre><code>nchurn1 %>%
  select_if(is.numeric) %>%
  gather(metric, value) %>%
  ggplot(aes(value, fill = metric)) +
      geom_density(show.legend = FALSE) +
      facet_wrap( ~ metric, scales = &quot;free&quot;)
</code></pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx9-1.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx9-1-490x218.png?w=450&ssl=1" alt="" class="wp-image-15723" data-recalc-dims="1"></a></p>
<pre><code>nchurn2 %>%
  select_if(is.numeric) %>%
  gather(metric, value) %>%
  ggplot(aes(value, fill = metric)) +
      geom_density(show.legend = FALSE) +
      facet_wrap( ~ metric, scales = &quot;free&quot;)
</code></pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx9-2.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx9-2-490x218.png?w=450&ssl=1" alt="" class="wp-image-15724" data-recalc-dims="1"></a></p>
<p>If you compare the two plots you can see that they vary imperceptibly except for the y axis scale. There is no warning about missing values and if you scroll back and compare with the original plots of the raw variables the shape of <code>tenure</code> and <code>TotalCharges</code> have changed significantly because of the transformation.</p>
<p>I'm pretty convinced that knn provides a much better approximation of those eleven missing values than a mere median substitution so let's make those changes and move on to comparing models. While we're at it, let's go ahead and remove the unique customer ID number as well. We really only needed it to compare a few specific cases.</p>
<pre><code>churn &lt;- predict(pp_knn,churn)
churn$customerID &lt;- NULL
str(churn)</code>
<em>## 'data.frame': 7043 obs. of 20 variables:
## $ gender : Factor w/ 2 levels &quot;Female&quot;,&quot;Male&quot;: 1 2 2 2 1 1 2 1 1 2 ...
## $ SeniorCitizen : Factor w/ 2 levels &quot;No&quot;,&quot;Yes&quot;: 1 1 1 1 1 1 1 1 1 1 ...
## $ Partner : Factor w/ 2 levels &quot;No&quot;,&quot;Yes&quot;: 2 1 1 1 1 1 1 1 2 1 ...
## $ Dependents : Factor w/ 2 levels &quot;No&quot;,&quot;Yes&quot;: 1 1 1 1 1 1 2 1 1 2 ...
## $ tenure : num -1.644 0.297 -1.495 0.646 -1.495 ...
## $ PhoneService : Factor w/ 2 levels &quot;No&quot;,&quot;Yes&quot;: 1 2 2 1 2 2 2 1 2 2 ...
## $ MultipleLines : Factor w/ 3 levels &quot;No&quot;,&quot;No phone service&quot;,..: 2 1 1 2 1 3 3 2 3 1 ...
## $ InternetService : Factor w/ 3 levels &quot;DSL&quot;,&quot;Fiber optic&quot;,..: 1 1 1 1 2 2 2 1 2 1 ...
## $ OnlineSecurity : Factor w/ 3 levels &quot;No&quot;,&quot;No internet service&quot;,..: 1 3 3 3 1 1 1 3 1 3 ...
## $ OnlineBackup : Factor w/ 3 levels &quot;No&quot;,&quot;No internet service&quot;,..: 3 1 3 1 1 1 3 1 1 3 ...
## $ DeviceProtection: Factor w/ 3 levels &quot;No&quot;,&quot;No internet service&quot;,..: 1 3 1 3 1 3 1 1 3 1 ...
## $ TechSupport : Factor w/ 3 levels &quot;No&quot;,&quot;No internet service&quot;,..: 1 1 1 3 1 1 1 1 3 1 ...
## $ StreamingTV : Factor w/ 3 levels &quot;No&quot;,&quot;No internet service&quot;,..: 1 1 1 1 1 3 3 1 3 1 ...
## $ StreamingMovies : Factor w/ 3 levels &quot;No&quot;,&quot;No internet service&quot;,..: 1 1 1 1 1 3 1 1 3 1 ...
## $ Contract : Factor w/ 3 levels &quot;Month-to-month&quot;,..: 1 2 1 2 1 1 1 1 1 2 ...
## $ PaperlessBilling: Factor w/ 2 levels &quot;No&quot;,&quot;Yes&quot;: 2 1 2 1 2 2 2 1 2 1 ...
## $ PaymentMethod : Factor w/ 4 levels &quot;Bank transfer (automatic)&quot;,..: 3 4 4 1 3 3 2 4 3 1 ...
## $ MonthlyCharges : num -1.158 -0.239 -0.343 -0.731 0.214 ...
## $ TotalCharges : num -1.81 0.254 -1.386 0.233 -1.249 ...
## $ Churn : Factor w/ 2 levels &quot;No&quot;,&quot;Yes&quot;: 1 1 2 1 2 2 1 1 2 1 ...
</em></pre>
<p>One more step before we start using <code>CHAID</code>, <code>ranger</code>, and <code>xgboost</code> and while we have the data in one frame. Let's take the 3 numeric variables and create 3 analogous variables as factors. This is necessary because <code>CHAID</code> requires categorical a.k.a. nominal data. If you'd like to review the options for how to “cut” the data please refer back to my earlier <a rel="nofollow" href="https://ibecav.github.io/chaidtutor1/" target="_blank">post</a>.</p>
<pre><code>churn &lt;- churn %>%
   mutate_if(is.numeric, 
             funs(factor = cut_number(., n=5, 
                                      labels = c(&quot;Lowest&quot;,&quot;Below Middle&quot;,&quot;Middle&quot;,&quot;Above Middle&quot;,&quot;Highest&quot;))))
summary(churn)</code>
<em>## gender SeniorCitizen Partner Dependents tenure PhoneService MultipleLines InternetService OnlineSecurity OnlineBackup DeviceProtection TechSupport StreamingTV StreamingMovies Contract PaperlessBilling PaymentMethod MonthlyCharges TotalCharges Churn tenure_factor MonthlyCharges_factor TotalCharges_factor
## Female:3488 No :5901 No :3641 No :4933 Min. :-1.8439 No : 682 No :3390 DSL :2421 No :3498 No :3088 No :3095 No :3473 No :2810 No :2785 Month-to-month:3875 No :2872 Bank transfer (automatic):1544 Min. :-1.5685 Min. :-1.929306 No :5174 Lowest :1481 Lowest :1420 Lowest :1409 ## Male :3555 Yes:1142 Yes:3402 Yes:2110 1st Qu.:-0.8555 Yes:6361 No phone service: 682 Fiber optic:3096 No internet service:1526 No internet service:1526 No internet service:1526 No internet service:1526 No internet service:1526 No internet service:1526 One year :1473 Yes:4171 Credit card (automatic) :1522 1st Qu.:-0.9632 1st Qu.:-0.783551 Yes:1869 Below Middle:1397 Below Middle:1397 Below Middle:1408 ## Median : 0.1183 Yes :2971 No :1526 Yes :2019 Yes :2429 Yes :2422 Yes :2044 Yes :2707 Yes :2732 Two year :1695 Electronic check :2365 Median : 0.2021 Median : 0.018205 Middle :1408 Middle :1411 Middle :1409 ## Mean : 0.0000 Mailed check :1612 Mean : 0.0000 Mean :-0.002732 Above Middle:1350 Above Middle:1407 Above Middle:1408 ## 3rd Qu.: 0.9252 3rd Qu.: 0.8341 3rd Qu.: 0.868066 Highest :1407 Highest :1408 Highest :1409 ## Max. : 1.3421 Max. : 1.7530 Max. : 1.758003
</em></pre>
<p>Okay now we have three additional variables that end in _factor, they're like their numeric equivalents only cut into more or less 5 equal bins.</p>
<h2>Training and testing our models</h2>
<p>We're going to use <code>caret</code> to train and test all three of the algorithms on our data. We could operate directly by invoking the individual model functions directly but <code>caret</code> will allow us to use some common steps. We'll employ <code>cross-validation</code> a.k.a. <code>cv</code> to mitigate the problem of over-fitting. This <a rel="nofollow" href="https://elitedatascience.com/overfitting-in-machine-learning#how-to-prevent" target="_blank">article</a> explains it well so I won't repeat that explanation here, I'll simply show you how to run the steps in <code>R</code>.</p>
<p>This is also a good time to point out that <code>caret</code> has extraordinarily comprehensive <a rel="nofollow" href="https://topepo.github.io/caret/" target="_blank">documentation</a> which I used extensively and I'm limiting myself to the basics.</p>
<p>As a first step, let's just take 30% of our data and put is aside as the testing data set. Why 30%? Doesn't have to be, could be as low as 20% or as high as 40% it really depends on how conservative you want to be, and how much data you have at hand. Since this is just a tutorial we'll simply use 30% as a representative number. I'm going to use <code>caret</code> syntax which is the line with <code>createDataPartition(churn$Churn, p=0.7, list=FALSE)</code> in it. That takes our data set <code>churn</code> makes a 70% split ensuring that we keep our outcome variable <code>Churn</code> as close to 70/30 as we can. <em>This is important because our data is already pretty lop-sided</em> for outcomes. The two subsequent lines serve to take the vector <code>intrain</code> and produce two separate dataframes, <code>testing</code> and <code>training</code>. They have 2112 and 4931 customers respectively. </p>
<pre><code>intrain &lt;- createDataPartition(churn$Churn, p=0.7, list=FALSE)
training &lt;- churn[intrain,]
testing &lt;- churn[-intrain,]
dim(training)</code>
<em>## [1] 4931 23
</em></pre>
<pre><code>dim(testing)</code>
<em>## [1] 2112 23
</em></pre>
<h3>CHAID</h3>
<p>Now that we have a training and testing dataset let's remove the numeric version of the variables <code>CHAID</code> can't use.</p>
<pre><code># first pass at CHAID
# remove numbers
training &lt;- training %>%
  select_if(is.factor)
dim(training)</code>
<em>## [1] 4931 20
</em></pre>
<pre><code>testing &lt;- testing %>%
  select_if(is.factor)
dim(testing)</code>
<em>## [1] 2112 20
</em></pre>
<p>The next step is a little counter-intuitive but quite practical. Turns out that many models do not perform well when you feed them a <code>formula</code> for the model even if they claim to support a formula interface (as CHAID does). Here's a Stack Overflow <a rel="nofollow" href="https://stackoverflow.com/questions/33088893/caret-random-forests-not-working-something-is-wrong-all-the-accuracy-metric" target="_blank">link</a> that discusses in detail but my suggestion to you is to always separate them and avoid the problem altogether. We're just taking our <code>predictors</code> or <code>features</code> and putting them in <code>x</code> while we put our outcome in <code>y</code>.</p>
<pre><code># create response and feature data
features &lt;- setdiff(names(training), &quot;Churn&quot;)
x &lt;- training[, features]
y &lt;- training$Churn
</code></pre>
<p><code>trainControl</code> is the next function within <code>caret</code> we need to use. Chapter 5 in the <code>caret</code> doco covers it in great detail. I'm simply going to pluck out a few sane and safe options. <code>method = &quot;cv&quot;</code> gets us cross-validation. <code>number = 5</code> is pretty obvious. I happen to like seeing the progress in case I want to go for coffee so <code>verboseIter = TRUE</code> (here I will turn it off since the static output is rather boring), and I play it safe and explicitly save my predictions <code>savePredictions = &quot;final&quot;</code>. We put everything in <code>train_control</code> which we'll use in a minute. <strong>We'll use this same <code>train_control</code> for all our models</strong></p>
<pre><code># set up 5-fold cross validation procedure
train_control &lt;- trainControl(method = &quot;cv&quot;,
                              number = 5,
#                              verboseIter = TRUE,
                              savePredictions = &quot;final&quot;)
</code></pre>
<p>By default <code>caret</code> allows us to adjust three parameters in our <code>chaid</code> model; <code>alpha2</code>, <code>alpha3</code>, and <code>alpha4</code>. As a matter of fact it will allow us to build a grid of those parameters and test all the permutations we like, using the same cross-validation process. I'm a bit worried that we're not being conservative enough. I'd like to train our model using p values for alpha that are not .05, .03, and .01 but instead the de facto levels in my discipline; .05, .01, and .001. The function in <code>caret</code> is <code>tuneGrid</code>. We'll use the base <code>R</code> function <code>expand.grid</code> to build a dataframe with all the combinations and then feed it to <code>caret</code> in our training via <code>tuneGrid = search_grid</code> in our call to <code>train</code>.</p>
<pre><code># set up tuning grid default
search_grid &lt;- expand.grid(
  alpha2 = c(.05, .01, .001),
  alpha4 = c(.05, .01, .001),
  alpha3 = -1
)
</code></pre>
<p>Now we can use the <code>train</code> function in <code>caret</code> to train our model! It wants to know what our <code>x</code> and <code>y</code>'s are, as well as our training control parameters which we've parked in <code>train_control</code>.</p>
<pre><code>chaid.model &lt;- train(
  x = x,
  y = y,
  method = &quot;chaid&quot;,
  trControl = train_control,
  tuneGrid = search_grid
)
chaid.model</code>
<em>## CHi-squared Automated Interaction Detection ## ## 4931 samples
## 19 predictor
## 2 classes: 'No', 'Yes' ## ## No pre-processing
## Resampling: Cross-Validated (5 fold) ## Summary of sample sizes: 3944, 3946, 3945, 3944, 3945 ## Resampling results across tuning parameters:
## ## alpha2 alpha4 Accuracy Kappa ## 0.001 0.001 0.7805723 0.3956853
## 0.001 0.010 0.7832088 0.3786619
## 0.001 0.050 0.7767165 0.3956067
## 0.010 0.001 0.7844263 0.3950757
## 0.010 0.010 0.7832088 0.3740510
## 0.010 0.050 0.7767167 0.3883397
## 0.050 0.001 0.7844263 0.3950757
## 0.050 0.010 0.7817889 0.3764934
## 0.050 0.050 0.7783394 0.3921129
## ## Tuning parameter 'alpha3' was held constant at a value of -1
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were alpha2 = 0.05, alpha3 = -1 and alpha4 = 0.001.
</em></pre>
<p>And after roughly two minutes it's done. Let's inspect what we have so far. The output gives us a nice concise summary. 4931 cases with 19 predictors. It gives us an idea of how many of the 4931 cases were used in the individual folds <code>Summary of sample sizes: 3944, 3946, 3945, 3944, 3945</code>. If you need a review of what <code>alpha2</code>, <code>alpha4</code>, and <code>alpha3</code> are please review the <code>?chaid</code> doco. </p>
<p>You'll notice that I stored the results in an object called <code>chaid.model</code>. That object has lots of useful information you can access (it's a <code>list</code> object of class “train”). As a matter of fact we will be creating one object per run and then using the stored information to build a nice comparison later. For now here are some useful examples of what's contained in the object…</p>
<ul>
<li> Produce the <code>confusionMatrix</code> across all folds <code>confusionMatrix(chaid.model)</code></li>
<li> Plot the effect of the tuning parameters on accuracy <code>plot(chaid.model)</code>. Note that the scaling deceives the eye and the results are close across the plot</li>
<li> Check on variable importance <code>varImp(chaid.model)</code></li>
<li> How long did it take? Look in <code>chaid.model$times</code></li>
</ul>
<p>If you need a refresher on what these represent please see the earlier <a rel="nofollow" href="https://ibecav.github.io/chaidtutor1/" target="_blank">post on CHAID</a>.</p>
<pre><code>confusionMatrix(chaid.model)</code>
<em>## Cross-Validated (5 fold) Confusion Matrix ## ## (entries are percentual average cell counts across resamples)
## ## Reference
## Prediction No Yes
## No 66.2 14.3
## Yes 7.2 12.2
## ## Accuracy (average) : 0.7844
</em></pre>
<pre><code>plot(chaid.model)
</code></pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx18-1.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx18-1-490x490.png?w=450&ssl=1" alt="" class="wp-image-15725" data-recalc-dims="1"></a></p>
<pre><code>varImp(chaid.model)</code>
<em>## ROC curve variable importance
## ## Importance
## Contract 100.0000
## tenure_factor 95.7457
## OnlineSecurity 81.2110
## TechSupport 74.1473
## TotalCharges_factor 57.0885
## MonthlyCharges_factor 53.3755
## OnlineBackup 48.3382
## DeviceProtection 45.3815
## PaperlessBilling 44.2540
## Partner 36.7387
## Dependents 35.4481
## PaymentMethod 29.7699
## SeniorCitizen 22.5287
## StreamingMovies 10.3811
## MultipleLines 9.1536
## InternetService 7.8258
## StreamingTV 6.7502
## gender 0.6105
## PhoneService 0.0000
</em></pre>
<pre><code>chaid.model$times</code>
<em>## $everything
## user system elapsed ## 146.121 2.658 149.002 ## ## $final
## user system elapsed ## 1.401 0.045 1.445 ## ## $prediction
## [1] NA NA NA
</em></pre>
<p>One of the nice aspects about CHAID as a method is that is relatively easy to “see”“ your model in either text or plot format. While there are packages that will help you &quot;see” a random forest; by definition (pardon the pun) it's hard to see the forest because of all the trees. Simply “printing” the final model with <code>chaid.model$finalModel</code> gives you the text representation while you can plot the final model with <code>plot(chaid.model$finalModel)</code>. As I explained in the earlier post it's nice being able to see where your model fits well and where it misses at a high level.</p>
<pre><code>chaid.model$finalModel</code>
<em>## ## Model formula:
## .outcome ~ gender + SeniorCitizen + Partner + Dependents + PhoneService + ## MultipleLines + InternetService + OnlineSecurity + OnlineBackup + ## DeviceProtection + TechSupport + StreamingTV + StreamingMovies + ## Contract + PaperlessBilling + PaymentMethod + tenure_factor + ## MonthlyCharges_factor + TotalCharges_factor
## ## Fitted party:
## [1] root
## | [2] Contract in Month-to-month
## | | [3] InternetService in DSL
## | | | [4] TotalCharges_factor in Lowest: No (n = 326, err = 49.7%)
## | | | [5] TotalCharges_factor in Below Middle: No (n = 255, err = 27.5%)
## | | | [6] TotalCharges_factor in Middle, Above Middle, Highest: No (n = 282, err = 15.2%)
## | | [7] InternetService in Fiber optic
## | | | [8] tenure_factor in Lowest
## | | | | [9] OnlineSecurity in No, No internet service: Yes (n = 414, err = 25.4%)
## | | | | [10] OnlineSecurity in Yes: No (n = 33, err = 42.4%)
## | | | [11] tenure_factor in Below Middle: Yes (n = 413, err = 41.4%)
## | | | [12] tenure_factor in Middle: No (n = 352, err = 44.3%)
## | | | [13] tenure_factor in Above Middle, Highest: No (n = 276, err = 31.9%)
## | | [14] InternetService in No: No (n = 378, err = 17.2%)
## | [15] Contract in One year
## | | [16] StreamingMovies in No: No (n = 303, err = 8.3%)
## | | [17] StreamingMovies in No internet service: No (n = 268, err = 2.2%)
## | | [18] StreamingMovies in Yes: No (n = 457, err = 20.6%)
## | [19] Contract in Two year
## | | [20] InternetService in DSL, No: No (n = 866, err = 1.2%)
## | | [21] InternetService in Fiber optic: No (n = 308, err = 8.1%)
## ## Number of inner nodes: 7
## Number of terminal nodes: 14
</em></pre>
<pre><code>plot(chaid.model$finalModel)
</code></pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx19-1.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx19-1-490x245.png?w=450&ssl=1" alt="" class="wp-image-15726" data-recalc-dims="1"></a></p>
<p>Finally, probably the most important step of all, we'll take our trained model and apply it to the <code>testing</code> data that we held back to see how well it fits this data it's never seen before. This is a key step because it reassures us that we have not overfit (if you want a fuller understanding please consider reading this post on EliteDataScience) <a rel="nofollow" href="https://elitedatascience.com/overfitting-in-machine-learning" target="_blank">our model</a>. We'll take our model we made with the training dataset <code>chaid.model</code> and have it predict against the testing dataset and see how we did with a <code>confusionMatrix</code></p>
<pre><code>confusionMatrix(predict(chaid.model, newdata = testing), testing$Churn)</code>
<em>## Confusion Matrix and Statistics
## ## Reference
## Prediction No Yes
## No 1443 331
## Yes 109 229
## ## Accuracy : 0.7917 ## 95% CI : (0.7737, 0.8088)
## No Information Rate : 0.7348 ## P-Value [Acc > NIR] : 7.793e-10 ## ## Kappa : 0.3878 ## Mcnemar's Test P-Value : &lt; 2.2e-16 ## ## Sensitivity : 0.9298 ## Specificity : 0.4089 ## Pos Pred Value : 0.8134 ## Neg Pred Value : 0.6775 ## Prevalence : 0.7348 ## Detection Rate : 0.6832 ## Detection Prevalence : 0.8400 ## Balanced Accuracy : 0.6693 ## ## 'Positive' Class : No ## </em></pre>
<p>Very nice! Our accuracy on <code>testing</code> actually exceeds the accuracy we achieved in training. </p>
<h3>Random Forest via <code>ranger</code></h3>
<p>One of the nicest things about using <code>caret</code> is that it is pretty straight-forward to move from one model to another. The amount of work we have to do while moving from <code>CHAID</code> to <code>ranger</code> and eventually <code>xgboost</code> is actually quite modest.</p>
<p><code>ranger</code> will accept a mix of factors and numeric variables so our first step will be to go back and recreate training and testing using the numeric versions of <code>tenure</code>, <code>MonthlyCharges</code>, and <code>TotalCharges</code> instead of the _factor versions. <code>intrain</code> still holds our list of rows that should be in <code>training</code> so we'll follow the exact same process just keep the numeric versions and arrive at <code>x</code> and <code>y</code> to feed to <code>caret</code> and <code>ranger</code>.</p>
<pre><code>##### using ranger
# intrain &lt;- createDataPartition(churn$Churn,p=0.7,list=FALSE)
training &lt;- churn[intrain,]
testing &lt;- churn[-intrain,]
dim(training)</code>
<em>## [1] 4931 23
</em></pre>
<pre><code>dim(testing)</code>
<em>## [1] 2112 23
</em></pre>
<pre><code>training &lt;- training %>%
  select(-ends_with(&quot;_factor&quot;))
dim(training)</code>
<em>## [1] 4931 20
</em></pre>
<pre><code># testing &lt;- testing %>%
#  select(-ends_with(&quot;_factor&quot;))
dim(testing)</code>
<em>## [1] 2112 23
</em></pre>
<pre><code># create response and feature data
features &lt;- setdiff(names(training), &quot;Churn&quot;)
x &lt;- training[, features]
y &lt;- training$Churn
</code></pre>
<p>As I mentioned earlier <code>train_control</code> doesn't have to change at all. So I'll just print it to remind you of what's in there.</p>
<p><code>search_grid</code> is almost always specific to the model and this is no exception. When we consult <a rel="nofollow" href="https://topepo.github.io/caret/train-models-by-tag.html#random-forest" target="_blank">the documentation</a> for ranger within caret we see that we can adjust <code>mtry</code>, <code>splitrule</code>, and <code>min.node.size</code>. We'll put in some reasonable values for those and then put the resulting grid into rf_grid. I tried to give ranger's search grid about the same amount of flexibility as I did for CHAID.</p>
<pre><code>##### reusing train_control
head(train_control)</code>
<em>## $method
## [1] &quot;cv&quot;
## ## $number
## [1] 5
## ## $repeats
## [1] NA
## ## $search
## [1] &quot;grid&quot;
## ## $p
## [1] 0.75
## ## $initialWindow
## NULL
</em></pre>
<pre><code># define a grid of parameter options to try with ranger
rf_grid &lt;- expand.grid(mtry = c(2:4),
                       splitrule = c(&quot;gini&quot;),
                       min.node.size = c(3, 5, 7))
rf_grid</code>
<em>## mtry splitrule min.node.size
## 1 2 gini 3
## 2 3 gini 3
## 3 4 gini 3
## 4 2 gini 5
## 5 3 gini 5
## 6 4 gini 5
## 7 2 gini 7
## 8 3 gini 7
## 9 4 gini 7
</em></pre>
<p>Okay, we're ready to train our model using ranger now. The only additional line we need (besides changing from <code>chaid</code> to <code>ranger</code> is to tell it what to use to capture variable importance e.g. “impurity”.</p>
<pre><code># re-fit the model with the parameter grid
rf.model &lt;- train(
                  x = x,
                  y = y,
                  method = &quot;ranger&quot;,
                  trControl = train_control,
                  tuneGrid = rf_grid,
                  importance = &quot;impurity&quot;)
rf.model</code>
<em>## Random Forest ## ## 4931 samples
## 19 predictor
## 2 classes: 'No', 'Yes' ## ## No pre-processing
## Resampling: Cross-Validated (5 fold) ## Summary of sample sizes: 3946, 3944, 3945, 3944, 3945 ## Resampling results across tuning parameters:
## ## mtry min.node.size Accuracy Kappa ## 2 3 0.7990299 0.4368130
## 2 5 0.8000437 0.4380797
## 2 7 0.7976106 0.4298171
## 3 3 0.7963928 0.4342273
## 3 5 0.7980151 0.4365761
## 3 7 0.7986242 0.4360132
## 4 3 0.7970009 0.4343468
## 4 5 0.7982171 0.4358659
## 4 7 0.7943648 0.4228518
## ## Tuning parameter 'splitrule' was held constant at a value of gini
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were mtry = 2, splitrule = gini and min.node.size = 5.
</em></pre>
<p>Now we can run the exact same set of commands as we did with <code>chaid.model</code> on <code>rf.model</code>.</p>
<pre><code>confusionMatrix(rf.model)</code>
<em>## Cross-Validated (5 fold) Confusion Matrix ## ## (entries are percentual average cell counts across resamples)
## ## Reference
## Prediction No Yes
## No 67.1 13.6
## Yes 6.4 12.9
## ## Accuracy (average) : 0.8
</em></pre>
<pre><code>plot(rf.model)
</code></pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx24-1.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx24-1-490x490.png?w=450&ssl=1" alt="" class="wp-image-15727" data-recalc-dims="1"></a></p>
<pre><code>varImp(rf.model)</code>
<em>## ranger variable importance
## ## Overall
## tenure 100.000
## TotalCharges 97.336
## MonthlyCharges 88.321
## Contract 66.882
## OnlineSecurity 43.533
## TechSupport 37.978
## PaymentMethod 34.685
## InternetService 29.435
## OnlineBackup 18.389
## DeviceProtection 14.376
## PaperlessBilling 13.363
## MultipleLines 10.114
## Partner 9.712
## StreamingTV 9.607
## StreamingMovies 8.183
## gender 7.998
## Dependents 7.679
## SeniorCitizen 7.314
## PhoneService 0.000
</em></pre>
<pre><code>rf.model$times</code>
<em>## $everything
## user system elapsed ## 108.364 1.834 18.655 ## ## $final
## user system elapsed ## 2.160 0.030 0.359 ## ## $prediction
## [1] NA NA NA
</em></pre>
<p>Now, the all-important prediction against the testing data set.</p>
<pre><code>confusionMatrix(predict(rf.model, newdata = testing), testing$Churn)</code>
<em>## Confusion Matrix and Statistics
## ## Reference
## Prediction No Yes
## No 1427 286
## Yes 125 274
## ## Accuracy : 0.8054 ## 95% CI : (0.7879, 0.8221)
## No Information Rate : 0.7348 ## P-Value [Acc > NIR] : 2.048e-14 ## ## Kappa : 0.4501 ## Mcnemar's Test P-Value : 2.969e-15 ## ## Sensitivity : 0.9195 ## Specificity : 0.4893 ## Pos Pred Value : 0.8330 ## Neg Pred Value : 0.6867 ## Prevalence : 0.7348 ## Detection Rate : 0.6757 ## Detection Prevalence : 0.8111 ## Balanced Accuracy : 0.7044 ## ## 'Positive' Class : No ## </em></pre>
<p>Very nice! Once again our accuracy on testing actually exceeds the accuracy we achieved in training. Looks like we were more accurate than CHAID but we'll come back to that after we finish <code>xgboost</code>.</p>
<h3>Extreme Gradient Boosting via <code>xgboost</code></h3>
<p>Moving from <code>ranger</code> to <code>xgboost</code> is even easier than it was from <code>CHAID</code>.</p>
<p><code>xgboost</code> like <code>ranger</code> will accept a mix of factors and numeric variables so there is no need to change our training and testing datasets at all. There's also no need to change our <code>train_control</code>. As far as tuning goes <code>caret</code> supports 7 of the <a rel="nofollow" href="http://xgboost.readthedocs.io/en/latest/parameter.html" target="_blank">many parameters</a> that you could feed to <code>?xgboost</code>. If you consult the <code>caret</code> documentation <a rel="nofollow" href="https://topepo.github.io/caret/train-models-by-tag.html#" target="_blank">here</a> under <code>xgbTree</code> you'll see them listed. If you don't provide any tuning guidance then it will provide a default set of pretty rational initial values. I initially ran it that way but below for purposes of this post have chosen only a few that seem to make the largest difference to accuracy and set the rest to a constant.</p>
<p>One final important note about the code below. Notice in the train command I am feeding a formula <code>Churn ~ .</code> to <code>train</code>. If you try to give it the same x = x &amp; y = y syntax I used with <code>ranger</code> it will fail. That's because as stated in the doco “xgb.train accepts only an xgb.DMatrix as the input. xgboost, in addition, also accepts matrix, dgCMatrix, or name of a local data file.” You could use commands like <code>xx &lt;- model.matrix(~. -1, data=x)[,-1]</code> &amp;<br>
<code>yy &lt;- as.numeric(y) -1</code> to convert them but since our dataset is small I'm just going to use the formula interface. </p>
<pre><code># reusing train_control
head(train_control)</code>
<em>## $method
## [1] &quot;cv&quot;
## ## $number
## [1] 5
## ## $repeats
## [1] NA
## ## $search
## [1] &quot;grid&quot;
## ## $p
## [1] 0.75
## ## $initialWindow
## NULL
</em></pre>
<pre><code># define a grid of parameter options to try with xgboost
xgb_grid &lt;- expand.grid(nrounds = c(100, 150, 200),
                       max_depth = 1,
                       min_child_weight = 1,
                       subsample = 1,
                       gamma = 0,
                       colsample_bytree = 0.8,
                       eta = c(.2, .3, .4))
xgb_grid</code>
<em>## nrounds max_depth min_child_weight subsample gamma colsample_bytree eta
## 1 100 1 1 1 0 0.8 0.2
## 2 150 1 1 1 0 0.8 0.2
## 3 200 1 1 1 0 0.8 0.2
## 4 100 1 1 1 0 0.8 0.3
## 5 150 1 1 1 0 0.8 0.3
## 6 200 1 1 1 0 0.8 0.3
## 7 100 1 1 1 0 0.8 0.4
## 8 150 1 1 1 0 0.8 0.4
## 9 200 1 1 1 0 0.8 0.4
</em></pre>
<pre><code># Fit the model with the parameter grid
xgboost.model &lt;- train(Churn ~ ., 
                       training , 
                       method = &quot;xgbTree&quot;, 
                       tuneGrid = xgb_grid,
                       trControl = train_control)
xgboost.model</code>
<em>## eXtreme Gradient Boosting ## ## 4931 samples
## 19 predictor
## 2 classes: 'No', 'Yes' ## ## No pre-processing
## Resampling: Cross-Validated (5 fold) ## Summary of sample sizes: 3945, 3944, 3945, 3944, 3946 ## Resampling results across tuning parameters:
## ## eta nrounds Accuracy Kappa ## 0.2 100 0.7974010 0.4272991
## 0.2 150 0.7998353 0.4414212
## 0.2 200 0.8002400 0.4465105
## 0.3 100 0.8002410 0.4449462
## 0.3 150 0.8020667 0.4532478
## 0.3 200 0.8028779 0.4565789
## 0.4 100 0.8010528 0.4478758
## 0.4 150 0.8020655 0.4530151
## 0.4 200 0.8024716 0.4546512
## ## Tuning parameter 'max_depth' was held constant at a value of 1
## Tuning parameter 'gamma' was held constant at a value of 0
## Tuning parameter 'colsample_bytree' was held constant at a value of 0.8
## Tuning parameter 'min_child_weight' was held constant at a value of 1
## Tuning parameter 'subsample' was held constant at a value of 1
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were nrounds = 200, max_depth = 1, eta = 0.3, gamma = 0, colsample_bytree = 0.8, min_child_weight = 1 and subsample = 1.
</em></pre>
<p>After a (relatively) brief moment the results are back. Average accuracy on the <code>training</code> is .8029 which is better than <code>CHAID</code> or <code>ranger</code>. We can run the same additional commands simply by listing <code>xgboost.model</code>.</p>
<pre><code>confusionMatrix(xgboost.model)</code>
<em>## Cross-Validated (5 fold) Confusion Matrix ## ## (entries are percentual average cell counts across resamples)
## ## Reference
## Prediction No Yes
## No 66.5 12.7
## Yes 7.0 13.8
## ## Accuracy (average) : 0.8029
</em></pre>
<pre><code>plot(xgboost.model)
</code></pre>
<p>The plot:<br>
<a rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx27-1.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2018/07/crx27-1-490x490.png?w=450&ssl=1" alt="" class="wp-image-15728" data-recalc-dims="1"></a></p>
<pre><code>varImp(xgboost.model)</code>
<em>## xgbTree variable importance
## ## only 20 most important variables shown (out of 30)
## ## Overall
## tenure 100.00000
## InternetServiceFiber optic 64.17202
## ContractTwo year 48.23551
## InternetServiceNo 26.81070
## ContractOne year 10.38092
## OnlineSecurityYes 7.23882
## StreamingTVYes 5.04154
## MultipleLinesYes 2.87755
## StreamingMoviesYes 2.74796
## TechSupportYes 2.14009
## PhoneServiceYes 1.52273
## DependentsYes 1.33244
## SeniorCitizenYes 0.50710
## genderMale 0.23151
## DeviceProtectionYes 0.03949
## MonthlyCharges 0.00000
## MultipleLinesNo phone service 0.00000
## StreamingMoviesNo internet service 0.00000
## OnlineBackupNo internet service 0.00000
## PartnerYes 0.00000
</em></pre>
<pre><code>xgboost.model$times</code>
<em>## $everything
## user system elapsed ## 4.398 0.096 4.497 ## ## $final
## user system elapsed ## 0.266 0.002 0.269 ## ## $prediction
## [1] NA NA NA
</em></pre>
<p>Now, the all-important prediction against the testing data set.</p>
<pre><code>confusionMatrix(predict(xgboost.model, newdata = testing), testing$Churn)</code>
<em>## Confusion Matrix and Statistics
## ## Reference
## Prediction No Yes
## No 1405 262
## Yes 147 298
## ## Accuracy : 0.8063 ## 95% CI : (0.7888, 0.823)
## No Information Rate : 0.7348 ## P-Value [Acc > NIR] : 9.038e-15 ## ## Kappa : 0.4682 ## Mcnemar's Test P-Value : 1.731e-08 ## ## Sensitivity : 0.9053 ## Specificity : 0.5321 ## Pos Pred Value : 0.8428 ## Neg Pred Value : 0.6697 ## Prevalence : 0.7348 ## Detection Rate : 0.6652 ## Detection Prevalence : 0.7893 ## Balanced Accuracy : 0.7187 ## ## 'Positive' Class : No ## </em></pre>
<p>Very nice! Once again our accuracy on testing .8063 actually exceeds the accuracy we achieved in training. Looks like we were more accurate than either CHAID or ranger and we'll focus on the comparison in the next section. </p>
<h2>Comparing Models</h2>
<p>At this juncture we're faced with a problem I've had before. We're drowning in data from the individual <code>confusionMatrix</code> results. We'll resort to the same <code>purrr</code> solution to give us a far more legible table of results focusing on the metrics I'm most interested in. To do that we need to:</p>
<ul>
<li> Make a <code>named list</code> called <code>modellist</code> that contains our 3 models with a descriptive name for each</li>
<li> Use <code>map</code> from <code>purrr</code> to apply the <code>predict</code> command to each model in turn to our <code>testing</code> dataset</li>
<li> Pipe those results to a second <code>map</code> command to generate a confusion matrix comparing our predictions to <code>testing$Churn</code> which are the actual outcomes.</li>
<li> Pipe those results to a complex <code>map_dfr</code> that creates a dataframe of all the results with each model as a row.</li>
<li> Separately grab the elapsed times for training with commands like <code>chaid.model$times$everything[[3]]</code></li>
<li> Separately grab the best accuracy for training with commands like <code>max(chaid.model$results$Accuracy)</code></li>
<li> Then use <code>kable</code> to make a pretty table that is much easier to understand.</li>
</ul>
<pre><code>modellist &lt;- list(&quot;CHAID&quot; = chaid.model,
                  &quot;ranger&quot; = rf.model,
                  &quot;xgboost&quot; = xgboost.model)

CompareResults &lt;- map(modellist, ~ predict(.x, newdata = testing)) %>%
                  map(~ confusionMatrix(testing$Churn, .x)) %>%
                  map_dfr(~ cbind(as.data.frame(t(.x$overall)), 
                    as.data.frame(t(.x$byClass))), 
                    .id = &quot;Model&quot;)

CompareResults[1,&quot;ETime&quot;] &lt;- chaid.model$times$everything[[3]]
CompareResults[2,&quot;ETime&quot;] &lt;- rf.model$times$everything[[3]]
CompareResults[3,&quot;ETime&quot;] &lt;- xgboost.model$times$everything[[3]]
CompareResults[1,&quot;BestTrain&quot;] &lt;- max(chaid.model$results$Accuracy)
CompareResults[2,&quot;BestTrain&quot;] &lt;- max(rf.model$results$Accuracy)
CompareResults[3,&quot;BestTrain&quot;] &lt;- max(xgboost.model$results$Accuracy)

kable(CompareResults, &quot;html&quot;) %>% 
   kable_styling(bootstrap_options = c(&quot;striped&quot;, &quot;hover&quot;, &quot;condensed&quot;, &quot;responsive&quot;))
</code></pre>
<h2>What do we know?</h2>
<p>Well our table looks very nice but there's probably still too much information. What data should we focus on and what conclusions can we draw from our little exercise in comparative modeling? I will draw your attention back to this <a rel="nofollow" href="https://www.dataschool.io/simple-guide-to-confusion-matrix-terminology/" target="_blank">webpage</a> to review the terminology for classification models and how to interpret a confusion matrix.</p>
<p>So <code>Accuracy</code>, <code>Kappa</code>, and <code>F1</code> are all measures of overall accuracy. There are merits to each. <code>Pos Pred Value</code>, and <code>Neg Pred Value</code> are related but different nuanced ideas we'll discuss in a minute. We'll also want to talk about time to complete training our model with <code>ETime</code> and training accuracy with <code>BestTrain</code>.</p>
<p>Let's use <code>dplyr</code> to select just these columns we want and see what we can glean from this reduced table.</p>
<pre><code>CompareResults %>%
  select(Model, ETime, BestTrain, Accuracy, Kappa, F1, 'Pos Pred Value', 'Neg Pred Value') %>%
  kable(&quot;html&quot;) %>% 
    kable_styling(bootstrap_options = c(&quot;striped&quot;, &quot;hover&quot;, &quot;condensed&quot;, &quot;responsive&quot;))
</code></pre>
<table>
<thead>
<tr>
<th>
Model
</th>
<th>
ETime
</th>
<th>
BestTrain
</th>
<th>
Accuracy
</th>
<th>
Kappa
</th>
<th>
F1
</th>
<th>
Pos Pred Value
</th>
<th>
Neg Pred Value
</th>
</tr>
</thead>
<tbody>
<tr>
<td>
CHAID
</td>
<td>
162.493
</td>
<td>
0.7844263
</td>
<td>
0.7916667
</td>
<td>
0.3878325
</td>
<td>
0.8677090
</td>
<td>
0.9297680
</td>
<td>
0.4089286
</td>
</tr>
<tr>
<td>
ranger
</td>
<td>
19.375
</td>
<td>
0.8000437
</td>
<td>
0.8053977
</td>
<td>
0.4501003
</td>
<td>
0.8741194
</td>
<td>
0.9194588
</td>
<td>
0.4892857
</td>
</tr>
<tr>
<td>
xgboost
</td>
<td>
4.650
</td>
<td>
0.8028779
</td>
<td>
0.8063447
</td>
<td>
0.4681509
</td>
<td>
0.8729419
</td>
<td>
0.9052835
</td>
<td>
0.5321429
</td>
</tr>
</tbody>
</table>
<p>Clearly <code>xgboost</code> is the fastest to train a model, more than 30 times faster than <code>CHAID</code>, and 3 times faster than ranger for this data. Not really surprising since <a rel="nofollow" href="https://github.com/dmlc/xgboost" target="_blank">xgboost</a> is a very modern set of code designed from the ground up to be fast and efficient.</p>
<p>One interesting fact you can glean from all 3 models is that they all did better on <code>testing</code> than they did on <code>training</code>. This is slightly unusual since one would expect some differences to be missed but is likely simply due to a lucky split in our data with more of the difficult to predict cases falling in <code>training</code> than <code>testing</code>. The good news is it leaves us feeling comfortable that we did not overfit our model to the training data, which is why we were conservative in our fitting and cross-validated the training data.</p>
<p>No matter which “accuracy measure” we look at <code>Accuracy</code>, <code>F1</code> or <code>Kappa</code> the answer is pretty consistent, <code>xgboost</code> “wins” or is the most accurate. The exception is <code>F1</code> where <code>ranger</code> edges is out by 0.11775% which means it was correct on about 3 more cases out of 2112 cases in the <code>testing</code> set. </p>
<p>Notice that the differences in accuracy are not large as percentages <code>xgboost</code> is 1.4678% more accurate than <code>CHAID</code> or it correctly predicted 31 more customers. While more accurate is always “better” the practical significance is also a matter of what the stakes are. If a wrong prediction costs you $1,000.00 dollars that additional accuracy is more concerning than a lesser dollar amount.</p>
<p>I also deliberately included <a rel="nofollow" href="https://en.wikipedia.org/wiki/Positive_and_negative_predictive_values" target="_blank">Positive and Negative Predictive Values</a> the columns labelled <code>Pos Pred Value</code> and <code>Neg Pred Value</code> for a very specific reason. Notice that <code>CHAID</code> has the highest <code>Pos Pred Value</code> that means is is the most accurate at predicting customers who did not “churn”. Of the 1,552 customers who did not leave us is correctly predicted 1,443 of them. <code>xgboost</code> on the other hand was much much better at <code>Neg Pred Value</code> correctly predicting 298 out of 560 customers who left us. While <code>Accuracy</code>, <code>Kappa</code> and <code>F1</code> take different approaches to finding “balanced” accuracy sometimes one case negative or positive has more important implications for your business and you should choose those measures.</p>
<p>At least at this point after a possible tl;dr journey we have some empirical data to inform my original statement about CHAID: “As the name implies it is fundamentally based on the venerable Chi-square test – and while not the most powerful (in terms of detecting the smallest possible differences) or the fastest, it really is easy to manage and more importantly to tell the story after using it”.</p>
<h2>What don't we know?</h2>
<ul>
<li>That this example would apply to other types of datasets. Absolutely not! This sort of data is almost ideal for CHAID since it involves a lot of nominal/categorical and/or ordinal data. CHAID will get much slower faster as we add more columns. More generally this was one example relatively small dataset more about learning something about <code>caret</code> and process than a true comparison of accuracy across a wide range of cases.</li>
<li>This is the “best” these models can do with this data Absolutely not! I made no attempt to seriously tune any of them. Tried some mild comparability. Also made no effort to feature engineer or adjust. I'm pretty certain if you tried you can squeeze a little more out of all three. Even wth <code>CHAID</code> there's more we could do very easily. I arbitrarily divided <code>tenure</code> into 5 equal sized bins. Why not 10? Why not equidistant instead of equal sized?</li>
</ul>
<p>I hope you've found this useful. I am always open to comments, corrections and suggestions.</p>
<ul></ul>
<p><strong>
<p>Related Post</p>
<p></p></strong>
<ol>
<li><a rel="nofollow" href="https://datascienceplus.com/common-mistakes-to-avoid-when-learning-to-code-in-python/" title="Permanent Link toCommon Mistakes to Avoid When Learning to Code in Python" target="_blank">Common Mistakes to Avoid When Learning to Code in Python</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/top-python-libraries-for-machine-learning/" title="Permanent Link toTop Python Libraries for Machine Learning" target="_blank">Top Python Libraries for Machine Learning</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/get-popular-programming-languages-from-tiobe-index-using-r/" title="Permanent Link toGet Popular Programming Languages from TIOBE Index using R" target="_blank">Get Popular Programming Languages from TIOBE Index using R</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/how-to-make-seaborn-pairplot-and-heatmap-in-r-write-python-in-r/" title="Permanent Link toHow to make Seaborn Pairplot and Heatmap in R (Write Python in R)" target="_blank">How to make Seaborn Pairplot and Heatmap in R (Write Python in R)</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/topic-modeling-in-python-with-nltk-and-gensim/" title="Permanent Link toTopic Modeling in Python with NLTK and Gensim" target="_blank">Topic Modeling in Python with NLTK and Gensim</a></li>
</ol>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://datascienceplus.com/chaid-vs-ranger-vs-xgboost-a-comparison/"> R Programming – DataScience+</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/F5OKj3X8EI0" height="1"></a></div>