<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/nyc-buses-simple-cubist-regression/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://datascienceplus.com/nyc-buses-simple-cubist-regression/"> R Programming – DataScience+</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<ul></ul>
<p><strong>
<p>Categories</p>
<p></p></strong>
<ol>
<li><a rel="nofollow" href="https://datascienceplus.com/category/advanced-modeling/" title="Permanent Link toAdvanced Modeling" target="_blank">Advanced Modeling</a></li>
</ol>
<p><strong>
<p>Tags</p>
<p></p></strong>
<ol>
<li><a rel="nofollow" href="https://datascienceplus.com/tag/data-manipulation/" title="Permanent Link toData Manipulation" target="_blank">Data Manipulation</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/tag/machine-learning/" title="Permanent Link toMachine Learning" target="_blank">Machine Learning</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/tag/rstats/" title="Permanent Link toR Programming" target="_blank">R Programming</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/tag/tips-tricks/" title="Permanent Link toTips &amp; Tricks" target="_blank">Tips &amp; Tricks</a></li>
</ol>
<p>We continue on from our <a rel="nofollow" href="https://datascienceplus.com/nyc-bus-delays/" target="_blank">previously loaded data of NYC bus delays</a>. Although we are fairly sure that our target variable (time_delayed) is not accurate, it is all we have. We can thank the bus drivers of New York that we have anything at all, I suppose, but the inaccuracy of this variable is probably one of the reasons not very many people have touched this data set.</p>
<h2>Data preparation</h2>
<p>Some data preparation first. We have initially included the variable 'time_diff_report' in the ii main data set. This is the time between the incident and when the incident is called in to operations. The difficult with this number is that it is obviously included in the time estimate given for the breakdown. The time difference is only ever less than the breakdown time for cases which are reported after they are resolved. Obviously, this number can be used as a strong predictor, while also providing no useful predictive information for our real problem: working out how much longer the bus will be delayed. We correct all of this by subtracting this value from the delayed time, so that 'time_delayed' becomes the amount of time from after operations is called that the bus is delayed. We also remove cases which are resolved before they are reported.</p>
<pre><code>#remove zero variance columns
zero_var &lt;- function(dat) {
  out &lt;- lapply(dat, function(x) length(unique(x)))
  keep &lt;- rownames(data.frame(which(out > 1)))
  dat %>%
    select(keep)
}

in_csv &lt;- &quot;../output/intermediate/ii_spread.csv&quot;
ii_spread &lt;- read_csv(in_csv) %>%
  mutate(
    time_delayed = as.numeric(time_delayed - time_diff_report),
    time_diff_report = NULL
  ) %>%
  filter(reported_before_resolved == 1) %>%
  select(-School_Year, -Route_Number, -Schools_Serviced, -Bus_Company_Name, -Occurred_On) %>%
  zero_var

halfFold &lt;- createFolds(ii_spread$time_delayed, k = 5, list=TRUE)
ii_small &lt;- ii_spread[halfFold[[2]],]
ii_test &lt;-ii_spread[halfFold[[3]],]
ii_train &lt;- ii_spread %>%
  anti_join(ii_test, by = &quot;Busbreakdown_ID&quot;)

rctrl_manual &lt;- trainControl(method = &quot;none&quot;, returnResamp = &quot;all&quot;)
rctrl_repcv &lt;- trainControl(method = &quot;repeatedcv&quot;, number = 2, repeats = 5, returnResamp = &quot;all&quot;)
</code></pre>
<h2>Cross fold validation testing</h2>
<p>Since I am running all of this locally on a laptop which 'benchmarkme' tells me is the slowest system to have been benchmarked in 2018, I build data set 'ii_small' to do cross-fold validation testing, with root mean square error as my target. During my run, I test different numbers of committees and rules. The numbers presented here are what I ended up with. Neighbours always degrade the solution in this case are fixed to zero. We can see the full rule set of the model by calling 'summary', but this results in a lot of spam. A more succinct presentation is available by looking at the 'finalModel' components of the derived model. 'splits' gives a list of the derived rules, where multiple entries per rule indicates an 'AND' condition (both conditions satisfied). 'coefficients' gives the resulting slope and intercept of the application of that rule. Finally, we can get a full listing of how many times a predictor appears in the ruleset and the coefficient set by looking at 'usage': this is somewhat equivalent to a variable importance call.</p>
<pre><code>#---cubist
cubist_grd &lt;- expand.grid(
  committees = 1,
  neighbors = 0
)

set.seed(849)
cubist_cv &lt;- train(time_delayed ~ . -Busbreakdown_ID,
  data = ii_small,
  method = &quot;cubist&quot;,
  metric=&quot;RMSE&quot;,
  na.action = na.pass,
  trControl = rctrl_repcv,
  tuneGrid = cubist_grd,
  control = Cubist::cubistControl(unbiased = T, rules = 20, sample = 0)
)
cubist_cv </code>
<em>## Cubist ## ## 40908 samples
## 25 predictor
## ## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times) ## Summary of sample sizes: 20454, 20454, 20453, 20455, 20455, 20453, ... ## Resampling results:
## ## RMSE Rsquared MAE ## 13.52728 0.1336734 9.883126
## ## Tuning parameter 'committees' was held constant at a value of 1
## ## Tuning parameter 'neighbors' was held constant at a value of 0
</em></pre>
<pre><code>#summary(cubist_cv)
cubist_cv$finalModel$splits %>% glimpse</code>
<em>## Observations: 89
## Variables: 8
## $ committee </em></pre>