This vignette illustrates the bixplot function, which
visualizes a univariate dataset by combining an estimated density body,
a box-and-whisker summary, and a rug of individual data values. For
variables with more than one mode, the function automatically detects
clusters and displays each mode separately. The examples below start
with simple univariate comparisons and continue with more elaborate
multivariable layouts and rug coloring options.
## Warning: package 'vioplot' was built under R version 4.5.3
## Loading required package: sm
## Warning: package 'sm' was built under R version 4.5.3
## Package 'sm', version 2.2-6.0: type help(sm) for summary information
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.5.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:data.table':
##
## yearmon, yearqtr
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: perry
## Loading required package: parallel
## Loading required package: robustbase
##
## Attaching package: 'robustbase'
## The following object is masked from 'package:sm':
##
## aircraft
## The following object is masked from 'package:robCompositions':
##
## alcohol
We start with three simulated variables: one unimodal, one bimodal and one trimodal. This illustrates the core difference between a violin plot and a bixplot.
set.seed(1)
dat1 <- rnorm(120, 0, 2.5)
dat2 <- c(rnorm(80, -3, 1), rnorm(40, 3, 1))
dat3 <- c(rnorm(25, -4, 0.8), rnorm(50, 0, 0.8), rnorm(40, 4, 0.8))
xlist <- list(Unimodal = dat1, Bimodal = dat2, Multimodal = dat3)The default call already gives an informative result:
##
## bixplot(xlist, main = "bixplot")
##
## Variable 1 has length = 120:
## pvalue(dip.test) = 0.7848
## Selected k = 1.
##
## Variable 2 has length = 120:
## pvalue(dip.test) = 0.0022
## mykmax = 5
## Silhouette widths: [1] 0.792 0.630 0.605 0.537
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 80 40
##
## Variable 3 has length = 115:
## pvalue(dip.test) = 0.0013
## mykmax = 5
## Silhouette widths: [1] 0.635 0.742 0.639 0.593
## Selected k = 3.
## Cluster sizes:
## 1 2 3
## [1,] 25 50 40
The returned list bixout contains the cluster structure
and five-number summaries for each variable:
## $call
## bixplot(xlist, main = "bixplot")
##
## $p
## [1] 3
##
## $Unimodal
## $Unimodal$values
## [1] -5.536749718 -4.973379240 -4.512396572 -3.808917001 -3.676880960
## [6] -3.442648892 -3.191480521 -3.134083501 -3.061531537 -2.823407740
## [11] -2.610336566 -2.335244079 -2.277304121 -2.089071531 -2.051170960
## [16] -1.858183022 -1.774866077 -1.768737892 -1.721889236 -1.636461610
## [21] -1.626740883 -1.589341135 -1.566134527 -1.553101451 -1.550916693
## [26] -1.530065983 -1.473802365 -1.433163536 -1.421671832 -1.356300077
## [31] -1.195375138 -1.183501591 -1.154111826 -1.108229683 -1.037486408
## [36] -0.985724884 -0.982019824 -0.918053691 -0.799982171 -0.763470968
## [41] -0.760459809 -0.697783257 -0.633404200 -0.518451859 -0.443326206
## [46] -0.411308991 -0.389488767 -0.337946538 -0.337636510 -0.280865530
## [51] -0.256969318 -0.148283492 -0.140321849 -0.134512601 -0.112334023
## [56] -0.098100007 -0.040475658 0.002763379 0.070005397 0.105289683
## [61] 0.185853310 0.186412458 0.383133346 0.395071931 0.459108311
## [66] 0.471980749 0.667746977 0.728615589 0.823769430 0.832375928
## [71] 0.852799229 0.911454905 0.925047025 0.960463395 0.969179029
## [76] 0.974608091 0.995264701 1.044853900 1.188773822 1.218572631
## [81] 1.235470828 1.391657997 1.396216064 1.424299069 1.439453379
## [86] 1.484753303 1.484865469 1.526815884 1.549564370 1.724348406
## [91] 1.742408439 1.750534124 1.791768690 1.845811763 1.907939371
## [96] 1.921332311 1.955340752 2.053052988 2.202769316 2.275435574
## [101] 2.297443429 2.359590527 2.657749593 2.750063430 2.812327295
## [106] 2.901006539 2.945217491 3.019669515 3.396698879 3.580705596
## [111] 3.582559254 3.663887154 3.779452921 3.967083636 3.988202005
## [116] 4.205440201 4.418218173 4.950999746 5.431529176 6.004044401
##
## $Unimodal$fivenumbersummary
## [1] -5.5367497 -1.2758376 0.1455715 1.7333784 6.0040444
##
##
## $Bimodal
## $Bimodal$clustering
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [112] 2 2 2 2 2 2 2 2 2
##
## $Bimodal$cluster_1
## $Bimodal$cluster_1$members
## [1] -4.9143594 -4.6649724 -4.6406055 -4.5364498 -4.5183941 -4.4874603
## [7] -4.4672500 -4.4244947 -4.3844268 -4.2863005 -4.2313234 -4.2080828
## [13] -4.1159201 -4.0751923 -4.0479844 -4.0158475 -3.9293621 -3.9261095
## [19] -3.7660820 -3.7508190 -3.7317482 -3.6816605 -3.6520948 -3.6212667
## [25] -3.6192430 -3.5888945 -3.5282799 -3.5059575 -3.4635304 -3.4302118
## [31] -3.3810761 -3.3340008 -3.3242703 -3.3180684 -3.3009761 -3.2550270
## [37] -3.2386471 -3.2145794 -3.1795565 -3.1771040 -3.1587546 -3.1443996
## [43] -3.1001907 -3.0771529 -3.0735644 -3.0568968 -3.0376342 -3.0347260
## [49] -3.0185598 -2.9826044 -2.9398396 -2.8941976 -2.7924617 -2.7800752
## [55] -2.6934421 -2.5979882 -2.5880253 -2.5748996 -2.5498129 -2.5430012
## [61] -2.4789773 -2.4685038 -2.2873337 -2.2123604 -2.1696268 -2.1135773
## [67] -2.0161044 -1.9999712 -1.9726076 -1.9415170 -1.8234167 -1.7920916
## [73] -1.6569612 -1.5588423 -1.5354127 -1.1307094 -0.9247550 -0.9128335
## [79] -0.7938975 -0.6920216
##
## $Bimodal$cluster_1$fivenumbersummary
## [1] -4.9143594 -3.7412836 -3.1679293 -2.5109892 -0.6920216
##
##
## $Bimodal$cluster_2
## $Bimodal$cluster_2$members
## [1] 0.1110793 0.7147645 1.2667816 1.3944866 1.6297921 1.7467102 1.8135414
## [8] 1.8434276 2.0141733 2.3595183 2.3696997 2.5997533 2.6590314 2.6688680
## [15] 2.6690922 2.6912594 2.8356242 2.9018213 2.9402767 2.9552909 2.9866005
## [22] 2.9946560 3.0021319 3.1971934 3.2631756 3.4094018 3.4206946 3.5101084
## [29] 3.5413273 3.5608207 3.5705076 3.6422413 3.6670662 3.9878383 4.0967770
## [36] 4.5197450 4.5865884 4.6888733 4.8031419 5.4976616
##
## $Bimodal$cluster_2$fivenumbersummary
## [1] 0.1110793 2.3646090 2.9709457 3.5656642 5.4976616
##
##
##
## $Multimodal
## $Multimodal$clustering
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [75] 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [112] 3 3 3 3
##
## $Multimodal$cluster_1
## $Multimodal$cluster_1$members
## [1] -5.922477 -5.600132 -4.702966 -4.435833 -4.340214 -4.311390 -4.227465
## [8] -4.204537 -4.198131 -4.132897 -4.055724 -3.891022 -3.821216 -3.700220
## [15] -3.674266 -3.541808 -3.443559 -3.434151 -3.314072 -3.239190 -3.183629
## [22] -3.172714 -3.083017 -3.069628 -2.624298
##
## $Multimodal$cluster_1$fivenumbersummary
## [1] -5.922477 -4.227465 -3.821216 -3.314072 -2.624298
##
##
## $Multimodal$cluster_2
## $Multimodal$cluster_2$members
## [1] -1.83129918 -1.14100587 -1.05299613 -0.95129064 -0.83783852 -0.81963878
## [7] -0.80085774 -0.75186346 -0.68148567 -0.56098534 -0.53454289 -0.52462548
## [13] -0.46449144 -0.36330953 -0.33774721 -0.32602286 -0.31209493 -0.30690569
## [19] -0.26482638 -0.24465234 -0.20714607 -0.04045256 -0.02873794 0.06637259
## [25] 0.07926279 0.10755813 0.12480934 0.19533194 0.21604392 0.25840520
## [31] 0.30109623 0.31550333 0.31850412 0.34696172 0.45377673 0.59280093
## [37] 0.61247920 0.71493896 0.73584294 0.75614796 0.76410934 0.80412737
## [43] 0.83488997 0.90416581 1.05940690 1.20977016 1.32331624 1.42274343
## [49] 1.57706991 2.11933350
##
## $Multimodal$cluster_2$fivenumbersummary
## [1] -1.83129918 -0.46449144 0.09341046 0.71493896 2.11933350
##
##
## $Multimodal$cluster_3
## $Multimodal$cluster_3$members
## [1] 2.188089 2.504169 2.744734 2.872920 2.964688 3.149234 3.308771 3.576244
## [9] 3.612820 3.706039 3.717280 3.718142 3.732575 3.738319 3.768401 3.845761
## [17] 3.846977 3.903192 4.003519 4.136392 4.196969 4.303970 4.364908 4.385624
## [25] 4.395450 4.541815 4.543385 4.586200 4.591671 4.642627 4.651762 4.732815
## [33] 4.757269 4.855329 5.046321 5.091548 5.179105 5.197633 5.262313 5.509980
##
## $Multimodal$cluster_3$fivenumbersummary
## [1] 2.188089 3.711660 4.166680 4.647194 5.509980
For comparison, below we display the violin plot and the bixplot side by side. The bixplot separates the modes of the bimodal and trimodal variables and draws a dedicated box for each cluster.
ylim <- c(-8.5, 9)
par(las = 1, mfrow = c(1, 2))
par(mar = c(2.1, 2.2, 1.7, 2))
viocol <- adjustcolor("chocolate3", alpha.f = 0.5)
vioplot::vioplot(xlist, ylim = ylim, main = "", col = viocol)
title(main = "violin plot", line = 0.5, cex.main = 1)
par(mar = c(2.1, 2.2, 1.7, 0.2))
bixplot(xlist, ylim = ylim, main = "")##
## bixplot(xlist, ylim = ylim, main = "")
##
## Variable 1 has length = 120:
## pvalue(dip.test) = 0.7848
## Selected k = 1.
##
## Variable 2 has length = 120:
## pvalue(dip.test) = 0.0022
## mykmax = 5
## Silhouette widths: [1] 0.792 0.630 0.605 0.537
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 80 40
##
## Variable 3 has length = 115:
## pvalue(dip.test) = 0.0013
## mykmax = 5
## Silhouette widths: [1] 0.635 0.742 0.639 0.593
## Selected k = 3.
## Cluster sizes:
## 1 2 3
## [1,] 25 50 40
Next we apply bixplot to two real datasets: a latency
dataset with three variables and the penguin bill length variable
grouped by island.
## [1] 40 3
par(las = 1, mfrow = c(2, 2))
mar1 <- c(2.1, 2.4, 2.3, 2)
mar2 <- c(2.1, 2.2, 2.3, 0.2)
viocol <- adjustcolor("chocolate3", alpha.f = 0.5)
par(mar = mar1)
vioplot::vioplot(data_latenc, main = "", col = viocol)
title(main = "violin plot", cex.main = 1.2, line = 0.6)
par(mar = mar2)
mymodeCol <- c("cadetblue3", "hotpink2")
bixplot(data_latenc, main = "", cutmin = 0, cutmax = 300,
ylim = c(0, 300), modeCol = mymodeCol)##
## bixplot(data_latenc, modeCol = mymodeCol, cutmin = 0, cutmax = 300,
## ylim = c(0, 300), main = "")
##
## Variable 1 has length = 40:
## pvalue(dip.test) = 0.0018
## mykmax = 2
## Silhouette widths: [1] 0.691
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 27 13
##
## Variable 2 has length = 40:
## pvalue(dip.test) = 1e-04
## mykmax = 2
## Silhouette widths: [1] 0.733
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 20 20
##
## Variable 3 has length = 39:
## pvalue(dip.test) = 1e-04
## mykmax = 2
## Silhouette widths: [1] 0.752
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 24 15
title(main = "bixplot", cex.main = 1.2, line = 0.6)
# Score the islands in a meaningful order
islandscore <- rep(NA, length(penguins$island))
islandscore[penguins$island == "Torgersen"] <- 1
islandscore[penguins$island == "Biscoe"] <- 2
islandscore[penguins$island == "Dream"] <- 3
ylim <- c(29, 62)
par(mar = mar1)
vioplot::vioplot(bill_len ~ reorder(island, islandscore, mean),
data = penguins, main = "",
col = viocol, xlab = "", ylab = "",
ylim = ylim, cex.axis = 1)
title(main = "violin plot", cex.main = 1.2, line = 0.6)
par(mar = mar2)
bixplot(bill_len ~ reorder(island, islandscore, mean),
data = penguins, main = "", ylim = ylim,
bodyCol = "gray40", bodyOpaque = 0.3)##
## bixplot(bill_len ~ reorder(island, islandscore, mean), data = penguins,
## bodyCol = "gray40", bodyOpaque = 0.3, ylim = ylim, main = "")
##
## Variable 1 has length = 51:
## pvalue(dip.test) = 0.9909
## Selected k = 1.
##
## Variable 2 has length = 167:
## pvalue(dip.test) = 0.7905
## Selected k = 1.
##
## Variable 3 has length = 124:
## pvalue(dip.test) = 0.0055
## mykmax = 5
## Silhouette widths: [1] 0.683 0.619 0.622 0.572
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 62 62
For the latency variables, the bixplot correctly identifies a bimodal structure that is obscured by the smooth envelope of the violin plot. For the penguin bill lengths, the Biscoe island group is clearly bimodal, reflecting the presence of two penguin species.
The standardized iris dataset has four variables. The bixplot reveals that petal length and petal width have more than one mode (cluster).
par(mfrow = c(1, 2))
par(mar = c(4, 2, 2, 0.1))
diris <- data.frame(scale(iris[, 1:4]))
colnames(diris) <- c("Sepal.L", "Sepal.W", "Petal.L", "Petal.W")
mymodeCol <- c("cadetblue3", "hotpink2",
"cadetblue3", "hotpink2",
"lawngreen")
bixplot(diris, main = "", cut = 3, col = "gray75",
bodyOpaque = 0.6, rugW = 0.16,
rugoutCol = "red", curveLwd = 0.5,
modeCol = mymodeCol,
ylim = c(-3.6, 4.2), yaxs = "i",
xlab = "standardized variables")##
## bixplot(diris, col = "gray75", bodyOpaque = 0.6, modeCol = mymodeCol,
## curveLwd = 0.5, rugoutCol = "red", rugW = 0.16, cut = 3,
## ylim = c(-3.6, 4.2), main = "", xlab = "standardized variables",
## yaxs = "i")
##
## Variable 1 has length = 150:
## pvalue(dip.test) = 0.0789
## Selected k = 1.
##
## Variable 2 has length = 150:
## pvalue(dip.test) = 0.0177
## Selected k = 1.
##
## Variable 3 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.785 0.677 0.670 0.649
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 51 99
##
## Variable 4 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.725 0.727 0.700 0.693
## Selected k = 3.
## Cluster sizes:
## 1 2 3
## [1,] 50 52 48
We can also overlay marginal bixplots on top of a scatter plot, to show the marginal distributions of the axes variables. Here we add horizontal and vertical half-bixplots to a scatter plot of petal length versus petal width.
par(mar = c(4, 4, 2, 0.2))
x <- diris$Petal.L
y <- diris$Petal.W
xlim <- c(-2.3, 2.4)
ylim <- c(-2.2, 2.4)
xyratio <- (xlim[2] - xlim[1]) / (ylim[2] - ylim[1])
plot(x, y, xlim = xlim, ylim = ylim, pch = 16,
xlab = "", ylab = "", xaxs = "i", yaxs = "i")
title(xlab = "Petal.L", line = 2)
title(ylab = "Petal.W", line = 2)
title(main = "petal length versus petal width", cex.main = 1, line = 0.7)
bixplot(x, add = TRUE, horizontal = TRUE,
at = ylim[1] + 0.015, cutmin = xlim[1],
boxwex = 0.9, curveLwd = 0.5,
border = "black", side = "second",
bodyOpaque = 0.6)##
## bixplot(x, add = TRUE, at = ylim[1] + 0.015, horizontal = TRUE,
## bodyOpaque = 0.6, curveLwd = 0.5, border = "black", boxwex = 0.9,
## side = "second", cutmin = xlim[1])
##
## Variable 1 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.785 0.677 0.670 0.649
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 51 99
bixplot(y, add = TRUE, horizontal = FALSE,
at = xlim[1] + 0.015, cutmin = ylim[1],
boxwex = xyratio * 0.9,
modeCol = c("cadetblue3", "hotpink2", "lawngreen"),
curveLwd = 0.5, side = "second",
bodyOpaque = 0.6)##
## bixplot(y, add = TRUE, at = xlim[1] + 0.015, horizontal = FALSE,
## bodyOpaque = 0.6, modeCol = c("cadetblue3", "hotpink2", "lawngreen"),
## curveLwd = 0.5, boxwex = xyratio * 0.9, side = "second",
## cutmin = ylim[1])
##
## Variable 1 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.725 0.727 0.700 0.693
## Selected k = 3.
## Cluster sizes:
## 1 2 3
## [1,] 50 52 48
The argument xyratio scales the vertical bixplot so that
its density body has the same visual height as the horizontal one.
The bodysize argument controls how the density bodies of
individual modes within a variable are sized relative to each other. The
three options are illustrated below on the (bimodal) petal length
variable.
par(mfrow = c(1, 3))
par(mar = c(2.5, 2, 2, 1))
for (bs in c("width_is_constant", "area_is_constant", "area_from_count")) {
bixplot(diris[, 3], main = "", cut = 3, col = "gray75",
bodyOpaque = 0.5, bodysize = bs, curveLwd = 0.5,
ylim = c(-2.2, 2.4), yaxs = "i", names = "Petal.L")
title(main = switch(bs,
width_is_constant = "equal width",
area_is_constant = "equal area",
area_from_count = "area from count"),
cex.main = 1, line = 0.7)
}##
## bixplot(diris[, 3], names = "Petal.L", col = "gray75", bodyOpaque = 0.5,
## bodysize = bs, curveLwd = 0.5, cut = 3, ylim = c(-2.2, 2.4),
## main = "", yaxs = "i")
##
## Variable 1 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.785 0.677 0.670 0.649
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 51 99
##
## bixplot(diris[, 3], names = "Petal.L", col = "gray75", bodyOpaque = 0.5,
## bodysize = bs, curveLwd = 0.5, cut = 3, ylim = c(-2.2, 2.4),
## main = "", yaxs = "i")
##
## Variable 1 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.785 0.677 0.670 0.649
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 51 99
##
## bixplot(diris[, 3], names = "Petal.L", col = "gray75", bodyOpaque = 0.5,
## bodysize = bs, curveLwd = 0.5, cut = 3, ylim = c(-2.2, 2.4),
## main = "", yaxs = "i")
##
## Variable 1 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.785 0.677 0.670 0.649
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 51 99
With "area_from_count" (the default), the larger cluster
receives a proportionally higher area, making cluster sizes immediately
apparent.
Here we use the formula interface to plot bill length split by the
interaction of sex and island. We compare a standard multi-group bixplot
with the side = "both" layout, which pairs male and female
distributions on opposing sides of a shared axis.
par(las = 1, mfrow = c(1, 2))
islscore <- rep(NA, length(penguins$island))
islscore[penguins$island == "Torgersen"] <- 1
islscore[penguins$island == "Biscoe"] <- 2
islscore[penguins$island == "Dream"] <- 3
ylim <- c(27, 64)
mynames <- c("Torger.F", "Torger.M", "Biscoe.F",
"Biscoe.M", "Dream.F", "Dream.M")
mycol <- c("slateblue2", "orange")
mymodeCol <- c("darkorchid1", "slateblue3",
"goldenrod1", "darkorange2")
par(mar = c(2.9, 2, 0.8, 1))
bixplot(bill_len ~ sex + reorder(island, islscore, mean),
data = penguins, names = mynames, modeCol = mymodeCol,
bodyOpaque = 0.6, ylim = ylim, bodyW = 0.9,
col = mycol, main = "", rugCol = "black", las = 1)##
## bixplot(bill_len ~ sex + reorder(island, islscore, mean), data = penguins,
## names = mynames, col = mycol, bodyOpaque = 0.6, bodyW = 0.9,
## modeCol = mymodeCol, rugCol = "black", ylim = ylim, main = "",
## las = 1)
##
## Variable 1 has length = 24:
## Selected k = 1.
##
## Variable 2 has length = 23:
## Selected k = 1.
##
## Variable 3 has length = 80:
## pvalue(dip.test) = 0.7648
## Selected k = 1.
##
## Variable 4 has length = 83:
## pvalue(dip.test) = 0.7078
## Selected k = 1.
##
## Variable 5 has length = 61:
## pvalue(dip.test) = 3e-04
## mykmax = 4
## Silhouette widths: [1] 0.705 0.589 0.661
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 27 34
##
## Variable 6 has length = 62:
## pvalue(dip.test) = 0
## mykmax = 4
## Silhouette widths: [1] 0.829 0.619 0.513
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 28 34
legend("topleft", legend = c("female", "male"),
fill = c("slateblue3", "orange"), cex = 1.5)
par(mar = c(2.9, 3, 0.8, 0.1))
bixplot(bill_len ~ sex + reorder(island, islscore, mean),
data = penguins, main = "", rugCol = "black",
bodyOpaque = 0.6, ylim = ylim, bodyW = 0.9,
col = mycol, stickCol = "gray10", stickLwd = 1,
side = "both", modeCol = mymodeCol, las = 1,
names = c("Torgersen", "Biscoe", "Dream"))##
## bixplot(bill_len ~ sex + reorder(island, islscore, mean), data = penguins,
## names = c("Torgersen", "Biscoe", "Dream"), col = mycol, bodyOpaque = 0.6,
## bodyW = 0.9, modeCol = mymodeCol, rugCol = "black", stickCol = "gray10",
## stickLwd = 1, side = "both", ylim = ylim, main = "", las = 1)
##
## Variable 1 has length = 24:
## Selected k = 1.
##
## Variable 2 has length = 23:
## Selected k = 1.
##
## Variable 3 has length = 80:
## pvalue(dip.test) = 0.7648
## Selected k = 1.
##
## Variable 4 has length = 83:
## pvalue(dip.test) = 0.7078
## Selected k = 1.
##
## Variable 5 has length = 61:
## pvalue(dip.test) = 3e-04
## mykmax = 4
## Silhouette widths: [1] 0.705 0.589 0.661
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 27 34
##
## Variable 6 has length = 62:
## pvalue(dip.test) = 0
## mykmax = 4
## Silhouette widths: [1] 0.829 0.619 0.513
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 28 34
With side = "both" the number of axis labels is halved
and the male/female pair for each island shares one axis, making
comparisons within an island more direct. Here, male penguins tend to
have larger bill lengths than females.
The rugNumeric argument colors each rug tick mark
according to a continuous covariate via a smooth color palette. Here the
rug is colored by body mass, and a color bar legend is added
automatically.
par(las = 1, mfrow = c(1, 1))
par(mar = c(3.2, 3.2, 2, 2.6))
scpenguins <- scale(penguins[, c(4, 3, 5)])
varnames <- c("bill_depth", "bill_length", "flipper_length")
bixplot(scpenguins, rugNumeric = penguins$body_mass,
main = "", boxW = 0.30, rugW = 0.24,
names = varnames, boxLwd = 1.2, boxOpaque = 1,
curveLwd = 2, colorbarW = 0.16, cex.colorbar = 0.8,
bodysize = "width_is_constant",
bodyCol = "grey90", modeCol = "grey90",
xlab = "feature", ylab = "standardized value")##
## bixplot(scpenguins, names = varnames, bodyCol = "grey90", bodysize = "width_is_constant",
## modeCol = "grey90", curveLwd = 2, boxOpaque = 1, boxW = 0.3,
## boxLwd = 1.2, rugNumeric = penguins$body_mass, colorbarW = 0.16,
## rugW = 0.24, main = "", cex.colorbar = 0.8, xlab = "feature",
## ylab = "standardized value")
##
## Variable 1 has length = 342:
## pvalue(dip.test) = 0.1551
## Selected k = 1.
##
## Variable 2 has length = 342:
## pvalue(dip.test) = 0.0837
## Selected k = 1.
##
## Variable 3 has length = 342:
## pvalue(dip.test) = 0.0021
## mykmax = 5
## Silhouette widths: [1] 0.692 0.573 0.558 0.537
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 209 133
We see that heavier penguins (green) tend to have a larger flipper length, revealing a positive association between these two variables.
The rugFactor argument colors each rug tick mark by a
factor level. Below we display the same three penguin measurements with
rug color indicating species, this time using a horizontal layout.
par(mfrow = c(1, 1))
par(mar = c(3.2, 3.8, 2, 0.2))
rugFactorColors <- c("red", "blue", "forestgreen")
bixplot(scpenguins, rugFactor = penguins$species,
main = "", boxW = 0.30, rugW = 0.24,
ylim = c(-3, 4), names = varnames,
boxLwd = 1.2, boxOpaque = 1, curveLwd = 2,
horizontal = TRUE, bodysize = "width_is_constant",
bodyCol = "grey90", modeCol = "grey90",
xlab = "standardized value", ylab = "feature",
rugFactorColors = rugFactorColors, las = 0)##
## bixplot(scpenguins, names = varnames, horizontal = TRUE, bodyCol = "grey90",
## bodysize = "width_is_constant", modeCol = "grey90", curveLwd = 2,
## boxOpaque = 1, boxW = 0.3, boxLwd = 1.2, rugFactor = penguins$species,
## rugFactorColors = rugFactorColors, rugW = 0.24, ylim = c(-3,
## 4), main = "", xlab = "standardized value", ylab = "feature",
## las = 0)
## The variable rugFactor has the levels:
## [1] "Adelie" "Chinstrap" "Gentoo"
##
## Variable 1 has length = 342:
## pvalue(dip.test) = 0.1551
## Selected k = 1.
##
## Variable 2 has length = 342:
## pvalue(dip.test) = 0.0837
## Selected k = 1.
##
## Variable 3 has length = 342:
## pvalue(dip.test) = 0.0021
## mykmax = 5
## Silhouette widths: [1] 0.692 0.573 0.558 0.537
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 209 133
title(main = "penguins with rug color by species", line = 0.6)
legend(x = 2.36, y = 2.8,
legend = c("Adelie", "Chinstrap", "Gentoo"),
fill = rugFactorColors, cex = 0.88)The species labels in the rug indicate that Gentoo penguins tend to
have smaller bill_depth and longer
flipper_length.
We now illustrate bixplot on the Top Gear dataset from
the robustHD package, which contains performance and price
data for about three hundred cars.
data(TopGear, package = "robustHD")
scars <- TopGear[, c(14, 12, 5, 7, 8)]
scars[, 3] <- log(scars[, 3])
colnames(scars)[3] <- "log(Price)"
scars[, 1:4] <- scale(scars[, 1:4])An initial bixplot reveals a severe outlier in the Weight variable:
##
## bixplot(scars[, 1:4])
##
## Variable 1 has length = 264:
## pvalue(dip.test) = 0.9126
## Selected k = 1.
##
## Variable 2 has length = 293:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.606 0.489 0.562 0.572
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 185 108
##
## Variable 3 has length = 297:
## pvalue(dip.test) = 0.9912
## Selected k = 1.
##
## Variable 4 has length = 288:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.733 0.571 0.632 0.659
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 228 60
## [1] 199
## Maker Model Weight TopSpeed Price Displacement
## 199 Peugeot 107 210 100 9345 998
The Peugeot 107 has a recorded weight of 210, which is clearly a data
error. We set it to NA and replot:
##
## bixplot(scars[, 1:4])
##
## Variable 1 has length = 263:
## pvalue(dip.test) = 0.9111
## Selected k = 1.
##
## Variable 2 has length = 293:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.606 0.489 0.562 0.572
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 185 108
##
## Variable 3 has length = 297:
## pvalue(dip.test) = 0.9912
## Selected k = 1.
##
## Variable 4 has length = 288:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.733 0.571 0.632 0.659
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 228 60
There is also an outlier in TopSpeed: the value 50 of the Renault Twizy, that is a tiny city car. This value is correct, so it is retained.
## [1] 220
## Maker Model Weight TopSpeed Price Displacement
## 220 Renault Twizy 450 50 6950 NA
We now produce a final two-panel figure: a bixplot of the four standardized variables alongside a scatter plot of TopSpeed versus Displacement with marginal bixplots added on each axis.
par(las = 1, mfrow = c(1, 2))
par(mar = c(4, 2, 2, 0.1))
bixplot(scars[, 1:4], ylim = c(-3.8, 5.2), main = "",
col = "darkgoldenrod2", yaxs = "i",
modeCol = c("cadetblue3", "hotpink2"),
bodyOpaque = 0.6, las = 1)##
## bixplot(scars[, 1:4], col = "darkgoldenrod2", bodyOpaque = 0.6,
## modeCol = c("cadetblue3", "hotpink2"), ylim = c(-3.8, 5.2),
## main = "", yaxs = "i", las = 1)
##
## Variable 1 has length = 263:
## pvalue(dip.test) = 0.9111
## Selected k = 1.
##
## Variable 2 has length = 293:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.606 0.489 0.562 0.572
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 185 108
##
## Variable 3 has length = 297:
## pvalue(dip.test) = 0.9912
## Selected k = 1.
##
## Variable 4 has length = 288:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.733 0.571 0.632 0.659
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 228 60
title(main = "standardized Top Gear variables",
cex.main = 1, line = 0.7)
par(mar = c(4, 4, 2, 0.1))
x <- jitter(scars[, 4])
y <- jitter(scars[, 2])
mycol <- rep(NA, nrow(scars))
mycol[scars[, 5] == "Front"] <- "red"
mycol[scars[, 5] == "Rear"] <- "forestgreen"
mycol[scars[, 5] == "4WD"] <- "orange"
plot(x, y, xlim = c(-1.999, 3.999), ylim = c(-2.999, 4.5),
xaxs = "i", yaxs = "i", xlab = "", ylab = "",
pch = 16, cex = 1, col = mycol, las = 1)
title(xlab = "Displacement", line = 2)
title(ylab = "TopSpeed", line = 2)
title(main = "TopSpeed versus Displacement",
cex.main = 1, line = 0.7)
bixplot(scars[, 4], add = TRUE, at = -2.98, horizontal = TRUE,
side = "second", boxwex = 2, bodyOpaque = 0.6)##
## bixplot(scars[, 4], add = TRUE, at = -2.98, horizontal = TRUE,
## bodyOpaque = 0.6, boxwex = 2, side = "second")
##
## Variable 1 has length = 288:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.733 0.571 0.632 0.659
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 228 60
bixplot(scars[, 2], add = TRUE, at = -1.985, horizontal = FALSE,
side = "second", boxwex = 1.4, bodyOpaque = 0.6)##
## bixplot(scars[, 2], add = TRUE, at = -1.985, horizontal = FALSE,
## bodyOpaque = 0.6, boxwex = 1.4, side = "second")
##
## Variable 1 has length = 293:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.606 0.489 0.562 0.572
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 185 108
legend(x = -1.5, y = 4, title = "DriveWheel",
legend = c("Front", "Rear", "4WD"),
fill = c("red", "forestgreen", "orange"), cex = 1)Rear-wheel-drive cars (green) span the full range of displacement and
are concentrated at the higher end of TopSpeed, whereas
front-wheel-drive cars (red) have rather low values of both
TopSpeed and Displacement.
side = "both"The side = "both" option is particularly useful for
two-group comparisons, as it places paired distributions on opposite
sides of a shared axis.
par(las = 0, mfrow = c(1, 2))
par(mar = c(3.5, 2, 2, 1))
bixplot(len ~ supp * dose,
data = ToothGrowth, side = "both",
col = c("orange", "slateblue3"),
main = "", xlab = "Vitamin C dose (mg)",
ylab = "tooth length", bodyOpaque = 0.7,
ylim = c(-1, 44), yaxs = "i", las = 0,
stickCol = "black", stickLwd = 1)##
## bixplot(len ~ supp * dose, data = ToothGrowth, col = c("orange",
## "slateblue3"), bodyOpaque = 0.7, stickCol = "black", stickLwd = 1,
## side = "both", ylim = c(-1, 44), main = "", xlab = "Vitamin C dose (mg)",
## ylab = "tooth length", yaxs = "i", las = 0)
##
## Variable 1 has length = 10:
## Selected k = 1.
##
## Variable 2 has length = 10:
## Selected k = 1.
##
## Variable 3 has length = 10:
## Selected k = 1.
##
## Variable 4 has length = 10:
## Selected k = 1.
##
## Variable 5 has length = 10:
## Selected k = 1.
##
## Variable 6 has length = 10:
## Selected k = 1.
title(main = "guinea pigs' tooth growth by supplement type",
cex.main = 1, line = 0.7)
legend("topleft", title = "supplement",
legend = c("OJ", "AA"),
fill = c("orange", "slateblue3"), cex = 1)
par(mar = c(3.5, 2.8, 2, 0.1))
diris <- data.frame(scale(iris[, 1:4]))
colnames(diris) <- c("Sepal.L", "Sepal.W", "Petal.L", "Petal.W")
bixplot(diris, side = "both", main = "",
bodysize = "width_is_constant",
col = c("orange", "slateblue3"),
modeCol = c("cadetblue3", "hotpink2",
"lawngreen", "gray60", "cyan3"),
stickCol = "red", stickLwd = 1,
bodyOpaque = 0.7, las = 0,
xlab = "standardized measurements")##
## bixplot(diris, col = c("orange", "slateblue3"), bodyOpaque = 0.7,
## bodysize = "width_is_constant", modeCol = c("cadetblue3",
## "hotpink2", "lawngreen", "gray60", "cyan3"), stickCol = "red",
## stickLwd = 1, side = "both", main = "", xlab = "standardized measurements",
## las = 0)
##
## Variable 1 has length = 150:
## pvalue(dip.test) = 0.0789
## Selected k = 1.
##
## Variable 2 has length = 150:
## pvalue(dip.test) = 0.0177
## Selected k = 1.
##
## Variable 3 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.785 0.677 0.670 0.649
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 51 99
##
## Variable 4 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.725 0.727 0.700 0.693
## Selected k = 3.
## Cluster sizes:
## 1 2 3
## [1,] 50 52 48
title(main = "iris data by length and width",
cex.main = 1, line = 0.7)
legend("top", legend = c("length", "width"),
fill = c("orange", "slateblue3"), cex = 1)In the tooth growth data we see that the orange juice (OJ) supplement yields higher tooth growth than ascorbic acid (AA), except at the highest dose, where the medians of OJ and AA align and AA has a larger spread than OJ. In the right panel, the iris data are now shown on two vertical axes instead of four.
Finally, we revisit the iris data with rug colors tied to an external variable, first a continuous one (sepal length) and then a factor (species).
par(las = 1, mfrow = c(1, 1))
par(mar = c(2.2, 2.2, 2, 2.6))
bixplot(diris, main = "",
rugNumeric = iris$Sepal.Length, colorbarW = 0.15,
bodyCol = "grey80", modeCol = "grey80", las = 1)##
## bixplot(diris, bodyCol = "grey80", modeCol = "grey80", rugNumeric = iris$Sepal.Length,
## colorbarW = 0.15, main = "", las = 1)
##
## Variable 1 has length = 150:
## pvalue(dip.test) = 0.0789
## Selected k = 1.
##
## Variable 2 has length = 150:
## pvalue(dip.test) = 0.0177
## Selected k = 1.
##
## Variable 3 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.785 0.677 0.670 0.649
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 51 99
##
## Variable 4 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.725 0.727 0.700 0.693
## Selected k = 3.
## Cluster sizes:
## 1 2 3
## [1,] 50 52 48
Flowers with high sepal length (green) tend to have high petal length and petal width. There is no such relation between sepal length and sepal width.
par(las = 1, mfrow = c(1, 1))
par(mar = c(2.2, 2.2, 2, 0.2))
rugFactorColors <- c("red", "blue", "forestgreen")
bixplot(diris, main = "", rugFactor = iris$Species,
bodyCol = "grey80", modeCol = "grey80",
rugFactorColors = rugFactorColors, las = 1)##
## bixplot(diris, bodyCol = "grey80", modeCol = "grey80", rugFactor = iris$Species,
## rugFactorColors = rugFactorColors, main = "", las = 1)
## The variable rugFactor has the levels:
## [1] "setosa" "versicolor" "virginica"
##
## Variable 1 has length = 150:
## pvalue(dip.test) = 0.0789
## Selected k = 1.
##
## Variable 2 has length = 150:
## pvalue(dip.test) = 0.0177
## Selected k = 1.
##
## Variable 3 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.785 0.677 0.670 0.649
## Selected k = 2.
## Cluster sizes:
## 1 2
## [1,] 51 99
##
## Variable 4 has length = 150:
## pvalue(dip.test) = 0
## mykmax = 5
## Silhouette widths: [1] 0.725 0.727 0.700 0.693
## Selected k = 3.
## Cluster sizes:
## 1 2 3
## [1,] 50 52 48
title(main = "iris data with rug color by species", line = 0.6)
legend("topright",
legend = c("setosa", "versicolor", "virginica"),
fill = rugFactorColors, cex = 1)The three species nicely match the clusters in Petal.W
and Petal.L. The colors inside Sepal.L are
roughly similar, whereas Sepal.W again behaves
differently.
As a final example we use the Titanic dataset from the
classmap package. We compare the standardized age and
log-fare distributions for casualties and survivors using the
side = "both" layout, and then color the rug by cabin
class.
data(data_titanic, package = "classmap")
titanic <- data_titanic[1:100, ]
titanic$logFare <- log(titanic$Fare + 1)
titanic$Pclass <- as.factor(titanic$Pclass)
titanic[, c(5, 14)] <- scale(titanic[, c(5, 14)])
xt <- list(titanic[titanic$y == "casualty", 5],
titanic[titanic$y == "survived", 5],
titanic[titanic$y == "casualty", 14],
titanic[titanic$y == "survived", 14])
names(xt) <- c("standardized Age.C",
"standardized Age.S",
"standardized log(Fare).C",
"standardized log(Fare).S")par(las = 1, mfrow = c(1, 2))
par(mar = c(2.4, 2, 2, 1))
mycol <- c("coral2", "cadetblue3")
bixplot(xt, side = "both", main = "", col = mycol,
stickLwd = 1, stickCol = "purple",
boxW = 0.18, rugW = 0.10, las = 1)##
## bixplot(xt, col = mycol, boxW = 0.18, rugW = 0.1, stickCol = "purple",
## stickLwd = 1, side = "both", main = "", las = 1)
##
## Variable 1 has length = 47:
## pvalue(dip.test) = 0.878
## Selected k = 1.
##
## Variable 2 has length = 31:
## pvalue(dip.test) = 0.8597
## Selected k = 1.
##
## Variable 3 has length = 59:
## pvalue(dip.test) = 0.2383
## Selected k = 1.
##
## Variable 4 has length = 41:
## pvalue(dip.test) = 0.5533
## Selected k = 1.
title(main = "Titanic data by survival", line = 0.6)
legend("top", legend = c("casualty", "survived"),
fill = mycol, cex = 0.9)
par(mar = c(2.4, 3, 2, 0.1))
bixplot(titanic[, c(5, 14)], main = "",
rugFactor = titanic$Pclass, boxW = c(0.22, 0.18),
names = c("standardized Age",
"standardized log(Fare)"), las = 1)##
## bixplot(titanic[, c(5, 14)], names = c("standardized Age", "standardized log(Fare)"),
## boxW = c(0.22, 0.18), rugFactor = titanic$Pclass, main = "",
## las = 1)
## The variable rugFactor has the levels:
## [1] "1" "2" "3"
##
## Variable 1 has length = 78:
## pvalue(dip.test) = 0.8149
## Selected k = 1.
##
## Variable 2 has length = 100:
## pvalue(dip.test) = 0.4382
## Selected k = 1.
title(main = "Titanic data by cabin class", line = 0.6)
legend("top", title = "cabin class",
legend = c("1", "2", "3"),
fill = c("red", "blue", "forestgreen"), cex = 0.9)In the left panel we see that the median of Age and
especially its mode are higher for the survivors, with a smaller effect
of Fare. On the right we see that both Age and
Fare have an effect on cabin class.