`library(biclustermd)`

Suppose you are a data scientist who is interested in modeling delays in flight arrivals at airports around the world. You think that arrivals could be correlated by month. Hence, since you’re looking at delays in arrivals you cluster average delay (in minutes) by month and destination airport.

To illustrate this problem I’m using the dataset `flights`

from Hadley Wickham’s package `nycflights13`

:

```
# install.packages("nycflights13")
library(nycflights13)
data("flights")
```

Per the documentation, `flights`

contains data on all flights in 2013 that departed NYC via JFK, LaGuardia, or Newark. The variables we’re interested in here are `month`

, `dest`

, and `arr_delay`

.

```
library(dplyr)
flights <- flights %>%
select(month, dest, arr_delay)
```

`bicluster()`

requires data to be in table format, which is what we will do here using `spread()`

from `tidyr`

after using `dplyr`

to summarize the data since we are analyzing average arrival delay.

```
library(tidyr)
flights <- flights %>%
group_by(month, dest) %>%
summarise(mean_arr_delay = mean(arr_delay, na.rm = TRUE)) %>%
spread(dest, mean_arr_delay) %>%
as.data.frame()
```

Now we name the rows using the text version of `month`

followed by removing `month`

and converting the data to a matrix.

```
rownames(flights) <- month.name[flights$month]
flights <- as.matrix(flights[, -1])
flights[1:5, 1:5]
```

```
## ABQ ACK ALB ANC ATL
## January NA NA 35.17460 NA 4.152047
## February NA NA 17.38889 NA 5.174092
## March NA NA 17.16667 NA 7.029286
## April 12.222222 NA 18.00000 NA 11.724280
## May -6.516129 3.904762 10.19643 NA 8.187036
```

Now that our data is in correct form we can bicluster it.

We first need to determine how many groups of months (*r*) and how many groups of destination airports (*c*) to define. Since rows correspond to months, in this analysis, setting *r* = 4 may not be a bad idea: create a group for each season/quarter of the year. For the columns (airports) the number of groups is a bit more ambiguous, just like with k-means clustering. Since each continent of the world may take a different amount of time and may have different policies for flights from the US, let’s set *c* = 6.

Now we’re ready to bicluster. I’ll put the code first and give brief descriptions to selected arguments. This will mostly be a repeat of the documentation for `biclustermd()`

.

```
library(biclustermd)
bc <- biclustermd(data = flights, col_clusters = 6, row_clusters = 4)
bc
```

```
##
## Data has 1260 values, 11.75% of which are missing
## 11 Iterations
## Initial SSE = 195467; Final SSE = 80940, a 58.6% reduction
## Rand similarity used; Indices: Columns (P) = 1, Rows (Q) = 1
```

`biclustermd`

outputs a list of class `biclustermd`

with useful information. Execute `?biclustermd`

to view the output in detail.

This list is used to make plots and training data, for example.

To plot SSE by iteration, use `autoplot.biclustermd_sse()`

, which outputs points by default:

```
library(ggplot2)
?autoplot.biclustermd_sse
autoplot(bc$SSE)
```

Next, cluster similarity measures. If a similarity measure for rows equals one, that means that the row shuffling from the last iteration and second to last iteration are identical. If zero, the two shufflings have nothing in common. The same goes for column similarities. `autoplot.biclustermd_sim()`

plots the RIs by iteration.

```
?autoplot.biclustermd_sim
autoplot(bc$Similarities)
```

`autoplot.biclustermd()`

makes visual analysis of biclustering results easy by making a heat map of the clustered data. The algorithm does use randomness, so my results may look different from yours.

```
?autoplot.biclustermd
autoplot(bc) +
scale_fill_viridis_c(na.value = "white") +
labs(x = "Destination Airport", y = "Month", fill = "Average Delay")
```

Often times it is helpful to run the data through an S-shaped function before plotting. This is easy with `autoplot.biclustermd()`

: set `transform_colors = TRUE`

and specify by what constant you’d like to scale your data by (with `c`

) before running it through a standard normal CDF:

```
autoplot(bc, transform_colors = TRUE, c = 1/10) +
scale_fill_viridis_c(na.value = "white", limits = c(0, 1)) +
labs(x = "Destination Airport", y = "Month", fill = "Average Delay")
```

To make the results really stand out we can reorder the row and column groups with the `reorder`

argument:

```
autoplot(bc, transform_colors = TRUE, c = 1/10, reorder = TRUE) +
scale_fill_viridis_c(na.value = "white", limits = c(0, 1)) +
labs(x = "Destination Airport", y = "Month", fill = "Average Delay")
```

Since the algorithm uses purposeful randomness, it is recommend that analysts run the biclustering multiple times and keep the one with the smallest SSE. The function `rep_biclustermd()`

does exactly that. Arguments to `rep_biclustermd()`

are the same as those for `biclustermd()`

with 3 additional arguments: 1. `nrep`

: the number of times to repeat the biclustering. 2. `parallel`

: a logical indicating if parallelization should be used or not. By default this is `FALSE`

. 3. `ncores`

: the number of cores to parallelize over. Ignored if `parallel = FALSE`

.

Below we run biclustering 100 times without using parallelization.

```
repeated_bc <- rep_biclustermd(flights, nrep = 100, col_clusters = 6, row_clusters = 4)
repeated_bc
```

```
## $best_bc
##
## Data has 1260 values, 11.75% of which are missing
## 11 Iterations
## Initial SSE = 199530; Final SSE = 78972, a 60.4% reduction
## Rand similarity used; Indices: Columns (P) = 1, Rows (Q) = 1
##
## $rep_sse
## [1] 87119.92 82594.24 81608.50 86624.23 83149.45 83912.02 82932.12
## [8] 81893.68 83342.92 85642.50 90324.74 84626.87 83925.14 81583.22
## [15] 84638.64 86700.00 81827.78 80057.36 90861.00 81608.50 88219.46
## [22] 90478.87 81608.50 86778.37 83720.40 83415.29 81453.21 88219.02
## [29] 82545.30 90926.11 87243.40 81955.76 84772.49 81878.69 86082.81
## [36] 81241.56 84158.32 86334.53 84557.68 82806.47 81662.06 84073.19
## [43] 83063.22 83522.82 79713.51 81617.01 81879.69 83390.13 80727.01
## [50] 81597.20 81608.50 85248.55 83242.92 81523.99 87534.30 81617.01
## [57] 82454.75 80906.94 86523.78 92567.57 86097.43 87132.13 83895.40
## [64] 85362.92 80027.62 84746.26 78971.82 79788.19 81622.42 90151.03
## [71] 83492.97 83106.76 81608.50 85007.90 85944.64 81650.03 80081.65
## [78] 85674.39 88562.24 87697.29 88377.38 82451.65 81995.61 80057.36
## [85] 82361.22 82699.95 81567.19 83525.05 81739.17 88726.92 80965.93
## [92] 85214.42 86129.56 81850.54 82674.72 84963.55 82940.31 81608.50
## [99] 84227.91 82177.09
##
## $runtime
## user system elapsed
## 15.459 0.222 15.884
```

`repeated_bc$best_bc`

returns the best biclustering, and has the same structure as the output of`biclustermd()`

.`repeated_bc$rep_sse`

is the SSE for each repeat, in order.`repeated_bc$runtime`

provides the CPU runtime for the user and system as well as elapsed time in seconds.

The 100 repeats take approximately 20 seconds to complete. Parallelization is quicker given sufficient complexity. However, as `nrep`

gets large, parallelization gets slow, since all `nrep`

objects have to be written to memory versus the minimum SSE object when processing serially. These remarks are current as of July 17, 2019.

`repeated_bc$runtime`

```
## user system elapsed
## 15.459 0.222 15.884
```

We can plot the best SSE at each repeat, that is, we can take the cumulative minimum of `rep_sse`

and plot:

`plot(cummin(repeated_bc$rep_sse), type = 'o', ylab = 'Cumulative Minimum', xlab = 'Repeat Number')`

We’ll end this tutorial by making a training dataset from the clustered data using `gather.biclustermd()`

, which gives the name of the row and column the data point comes from as well as the row and column group it belongs to.

```
training <- gather(repeated_bc$best_bc)
training %>% head()
```

```
## row_name col_name row_cluster col_cluster bicluster_no value
## 1 January ABQ 1 1 1 NA
## 2 April ABQ 1 1 1 12.222222
## 3 January ACK 1 1 1 NA
## 4 April ACK 1 1 1 NA
## 5 January BOS 1 1 1 -2.537891
## 6 April BOS 1 1 1 3.512500
```

Now we’ll plot cell (1, 4) and fit a linear model to it.

```
autoplot(repeated_bc$best_bc, row_clusts = 1, col_clusts = 4) +
scale_fill_viridis_c(na.value = "white") +
labs(x = "Destination Airport", y = "Month", fill = "Average Delay")
```

```
model <- training %>%
filter(row_cluster == 1, col_cluster == 4) %>%
lm(value ~ row_name + col_name, data = .)
summary(model)
```

```
##
## Call:
## lm(formula = value ~ row_name + col_name, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.725 -2.035 0.000 2.035 26.725
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.663093 5.094234 1.897 0.0649 .
## row_nameJanuary -3.449859 1.553727 -2.220 0.0320 *
## col_nameAUS 3.203479 7.120071 0.450 0.6551
## col_nameBDL 3.451832 7.120071 0.485 0.6304
## col_nameBGR 13.543804 8.754806 1.547 0.1295
## col_nameBNA 7.982674 7.120071 1.121 0.2687
## col_nameBQN 3.951582 7.120071 0.555 0.5819
## col_nameBTV 1.311083 7.120071 0.184 0.8548
## col_nameBUF 1.597314 7.120071 0.224 0.8236
## col_nameBUR -4.488163 7.120071 -0.630 0.5320
## col_nameBWI 7.540385 7.120071 1.059 0.2958
## col_nameBZN 3.286766 8.754806 0.375 0.7093
## col_nameCHO -7.463093 8.754806 -0.852 0.3989
## col_nameCHS 8.933172 7.120071 1.255 0.2167
## col_nameCLE -0.135876 7.120071 -0.019 0.9849
## col_nameCLT 0.726672 7.120071 0.102 0.9192
## col_nameCMH 2.956728 7.120071 0.415 0.6801
## col_nameDCA 2.652071 7.120071 0.372 0.7115
## col_nameDEN 3.679914 7.120071 0.517 0.6080
## col_nameDTW -1.706056 7.120071 -0.240 0.8118
## col_nameEYW 12.061837 7.120071 1.694 0.0978 .
## col_nameFLL 0.003447 7.120071 0.000 0.9996
## col_nameHOU 2.713753 7.120071 0.381 0.7051
## col_nameIAD 6.789118 7.120071 0.954 0.3459
## col_nameIND 8.477957 7.120071 1.191 0.2406
## col_nameJAX 6.480525 7.120071 0.910 0.3680
## col_nameMCO -1.909092 7.120071 -0.268 0.7899
## col_nameMDW 7.368741 7.120071 1.035 0.3068
## col_nameMEM 5.682258 7.120071 0.798 0.4294
## col_nameMSP 4.185336 7.120071 0.588 0.5599
## col_nameMSY 0.401768 7.120071 0.056 0.9553
## col_nameORD 5.295359 7.120071 0.744 0.4613
## col_nameORF 0.111919 7.120071 0.016 0.9875
## col_namePBI 1.245718 7.120071 0.175 0.8620
## col_namePHL 5.274161 7.120071 0.741 0.4631
## col_namePIT -2.170265 7.120071 -0.305 0.7621
## col_namePSE -10.462357 7.120071 -1.469 0.1494
## col_namePWM 5.790581 7.120071 0.813 0.4208
## col_nameRDU 3.061238 7.120071 0.430 0.6695
## col_nameROC 5.613774 7.120071 0.788 0.4350
## col_nameSAT -8.775317 7.120071 -1.232 0.2248
## col_nameSMF -0.193163 7.120071 -0.027 0.9785
## col_nameSTL 7.152415 7.120071 1.005 0.3210
## col_nameSYR -0.593272 7.120071 -0.083 0.9340
## col_nameTPA 1.469188 7.120071 0.206 0.8375
## col_nameXNA 5.575974 7.120071 0.783 0.4380
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.12 on 41 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.5058, Adjusted R-squared: -0.03653
## F-statistic: 0.9326 on 45 and 41 DF, p-value: 0.5917
```

`sqrt(mean(resid(model) ^ 2))`

`## [1] 4.887836`

If you have questions you can email me at johntreisner@gmail.com. In the case you’ve found an error please open up an issue.