Data exploration is a vital step to gain insights into data sets. Raw data needs to be cleaned, merged, summarized and assessed. This process is resource intensive, accounting for 80% of time spent on data analysis, by one estimate (Hellerstein et al. 2017) . Furthermore, decisions made in this stage can impact scientific rigor and reproducibility. Recently, an appreciation has emerged for systematic and transparent protocols about data inspection steps to be performed prior to formal data analysis (e.g. Huebner et al. (2016) ). Such protocols are designed to provide structure at this key stage while preventing statistical fishing for results.
Tools for data exploration, like tables and figures, have been historically important for science. For instance, in the late 1800s Florence Nightingale used rose plots to discover patterns in data that matched her clinical intuition-that most soldiers in the Crimean War were dying from hygiene-related infections rather than on the battlefield--and subsequently used this to influence the British Parliament (Nelson and Rafferty 2012) . This and other methods were a catalyst for the early-1900's revolution of statistical inference in many scientific fields.
Data exploration tools are more important today than ever. Data is more ubiquitous with a higher volume, velocity and variety than any time in history (Katal, Wazid, and Goudar 2013) . Further, these data are more accessible to analysis due to cheaper and more powerful computation (Waldrop 2016) . Consequently, data literacy and intuitive data exploration tools are required for exploring and communicating findings.
In this paper we introduce variable trees as a tool for exploring subsets of data, and their implementation in the vtree R package. The objectives of this paper are i) to compare variable trees to several established data exploration tools, ii) to review the functionality of the vtree package, and iii) to demonstrate the utility of variable trees in two case studies.
Subsets play an important role in almost any data analysis. Consider the variables relating to the 2207 passengers and crew members of the Titanic, represented in the the data set titanic_data from the EIX R package (https://modeloriented.github.io/EIX/ articles/titanic_data.html). Among other variables, the data set includes each person's home country, which we have grouped into regions, and age, which we have divided into children (under age 13) and adults. In Figure 1 each person's home region is shown, and within each region the number and percentage of children and adults are shown. Missing values are shown as NA. We call this a variable tree. The vtree package provides a general solution for drawing variable trees and describing nested subsets.
Even in simple situations like Figure 1 , it can be a chore to keep track of nested subsets and calculate the corresponding percentages. The denominator used to calculate percentages may also depend on whether the variables have any missing values, as discussed later. Finally, as the number of variables increases, the magnitude of the task balloons, because the number of nested subsets grows exponentially.
A variable tree consists of nodes connected by arrows. At the top of Figure 1 , the root node of the tree contains all 2207 people on the Titanic. The rest of the nodes are arranged in successive layers, where each layer corresponds to a specific variable. This highlights one difference between variable trees and some other kinds of trees: each layer of a variable tree corresponds to just one variable. This is distinct from decision trees, where a layer may include splits based on different variables.
The nodes immediately below the root node in Figure 1 represent values of Region and are referred to as the children 1 of the root node. Inside each of the nodes, the number of people is displayed and-except for in a missing value node-the corresponding percentage is also shown. An example of a missing value node appears in Figure 1 , where Region was missing (NA) for 81 people. Note that, by default, vtree displays "valid" percentages, i.e. the denominator used to calculate the percentage is the total number of non-missing values, in this case 2126. By default, vtree displays the full missing-value structure of the specified variables.
The final layer of the tree corresponds to values of Age. Each of these leaf nodes represents children and adults nested within a subset defined by a value of Region. Let's use the notation Region → Age to represent Age nested within Region. A missing-value node, like any node, can have children. For example, of the 81 people for whom Region is missing, 10 were children and 71 were adults.
A variety of tools have been developed to display joint distributions of discrete variables, the most basic being the contingency table, often enhanced with row, column, or table percentages. For example, Table 1 presents the same information as Figure 1 . Note that through the use of column percentages, the reader is encouraged to focus on age group nested within region. While the contingency table above is more compact than the variable tree in Figure 1 , we find the variable tree to be more intuitive. Furthermore, domain experts often respond well to such visual representations. Now suppose we'd like to examine Region → Age → Survived (i.e. survival within age within region of origin). Multi-way cross classifications (three or more variables) are typically displayed using several two-way tables, referred to as layers or slices. Table 2 shows two-way tables of survival within age group for each of the regions of origin. This is followed by a variable tree showing the same information ( Figure 2 ). Table 2 shows the same information.
Note that by default, vtree shows percentages in each node except for the root. For example, of the 2207 people on board the Titanic, 300 (14%) were from North America, of whom 287 (96%) were adults, of whom 160 (56%) survived. In its simplest form, a contingency table only shows crosstabulated frequencies, corresponding to the frequencies shown in the leaf nodes of a variable tree. Additionally, a variety of marginal and conditional percentages are often shown.
As the number of variables increases, contingency tables for multi-way classifications become increasingly difficult to interpret. In such situations, large variable trees can also become unwieldy, however this can be mitigated by pruning away branches of lesser interest.
Contingency tables are not always more compact than variable trees. When most cells of a large contingency table are empty (in which case the table is said to be sparse), the corresponding variable tree may be much more compact since empty nodes are not shown. In the Titanic data set, there are two missing values of Age, and both are for individuals from Europe. This appears as a single node in Figure 2 , but in Table 2 in addition to the cell showing these 2 missing values, there are 9 cells containing zero.
Like contingency tables, variable trees show numerical values (frequencies and percentages) rather than using graphical elements such as area to encode such quantities. In contrast to contingency tables, which use a tabular layout to represent subsets, variable trees use the graphical arrangement of nodes and arrows to represent the nesting structure.
Several visualization methods have been proposed for discrete multivariate data. Barplots are often used to visually represent the number of observations of each value of a variable. They can also be produced for subsets, defined by values of another variable. A more compact representation is the stacked barplot, however these are harder to read since they there is no common baseline, except for the bottom category in the stack.
An elegant extension of the stacked barplot is the mosaic plot (Hartigan and Kleiner 1981) .
In a mosaic plot, the area of each rectangle represents the number of observations in the corresponding subset of the data. Mosaic plots are available in base R through the mosaicplot function, or via the ggmosaic package or the vcd package. Mosaic plots can provide an intuitive visual representation of the number of observations in subsets of the data, however they tend to become overwhelming when there are more than three variables. Figure 3 is a mosaic plot for Region → Age → Survived for the people onboard the Titanic, as in Table 2 and Figure Table 2 and Figure 2 show the same information.
Visualizations like Figure 3 have advantages and disadvantages compared to text and tabular summaries. On the one hand, they represent quantitative and qualitative information in a way that is quickly decoded by our visual perceptual systems. On the other, visualizations can be unfamiliar and even perplexing compared to the familiarity of numerical and tabular representations. On a practical level, text and tabular information are easier to format and manipulate with current software. Variable trees have characteristics (and hence advantages as well as disadvantages) of both tabular representations and visualizations,
A special type of discrete multivariate data is when all of the variables are binary, in which case they can be interpreted as representing set membership. Venn diagrams use overlapping closed curves such that all intersections between sets are represented by overlapping areas. Euler diagrams are like Venn diagrams but empty intersections need not be shown. Venn and Euler diagrams have long been used to represent the intersection of sets. For datasets, software is available to calculate the number of observations in each of the intersections, for example in R, the VennDiagram and venneuler packages. A further elaboration of these diagrams is to make the areas of the sets and their intersections approximately proportional to the number of observations in each subset. The package eulerr provides this functionality. For example, in Figure 4 a dataset of Wilkinson (2012) is represented using an approximately area-proportional Euler diagram. As the number of sets grows, Venn and Euler diagrams can become unwieldy. An innovative way to represent the intersections of a large number of sets is UpSet (Lex et al. 2014) . The R package UpSetR (Conway, Lex, and Gehlenborg 2017) was used to produce Figure 5 for the dataset of Wilkinson (2012) . UpSet uses a grid layout to represent the intersections (see the dots at the bottom of Figure 5 ), together with bar graphs on each side to represent the size of sets and intersections. Variable trees can also represent the intersection of sets, however unlike UpSet and areaproportional Euler diagrams, they do not use graphical elements to encode quantity. Like non-proportional Venn Diagrams, variable trees graphically depict the relationships between subsets of the data, but represent quantities numerically ( Figure 6 ). Unlike Venn, Euler, and UpSet diagrams, variable trees require a prespecified ordering. For example, Figure 6 uses the ordering A → B → C → D → E → F. vtree can also display a pattern tree, which depicts every intersection. Each row in Figure 7 corresponds to the combination of values represented by a terminal node in Figure 6 . Since the intermediate nodes in Figure 6 are not represented, this is a loss of information. The pattern tree is much easier to read, however. Pattern trees have some of the same structure as an UpSet plot, except that sizes of subsets are not represented graphically as in the bar graphs on the sides of an UpSet plot. 
This section provides an overview of the features of the vtree package. Additional resources are available in the package vignette, a cheatsheet, and video tutorials on YouTube.
Suppose the Titanic data are in a data frame called td. To display a variable tree for a single variable, say Class, use the following command:
R> vtree(td,"Class") The variable Class is specified as a character string. To produce a variable tree for Class → Age, the character string is specified as "Class Age": R> vtree(td,"Class Age",horiz=FALSE) Figure 9 : A two-layer vertical variable tree.
By default, vtree produces horizontal trees. The tree in Figure 9 is vertical because of the specification horiz=FALSE.
When a variable tree gets too big, or you are only interested in certain parts of the tree, it may be useful to remove some nodes along with their descendants. This is known as pruning. For convenience, there are several different ways to prune a tree, described below.
Suppose you don't wish to show the "Europe" node or the "Other" node (which represents people from other parts of the world such as India, the Middle East, etc.). Specifying prune=list(Region=c("Europe","Other")) removes those nodes, and all of their descendants:
R> vtree(td,"Region Age", prune = list(Region = c("Europe", "Other")), + horiz = FALSE) In general, the argument of the prune parameter is a list with an element named for each variable you wish to prune. In the example above, the list has a single element, named Region. In turn, that element is a vector c("Europe","Other") indicating the values of Region to prune.
Note that once a variable tree has been pruned, it is no longer complete. This can sometimes be confusing since not all observations are represented at certain layers of the tree. For example in the tree above, only 1737 observations are shown in the Region nodes and their children.
Sometimes it is more convenient to specify which nodes should be retained rather than which ones should be discarded. The keep parameter is used for this purpose, and can thus be considered the complement of the prune parameter. Suppose we wish to create a variable tree for Region → Class → gender → Age. This tree has four layers, and without any pruning it would be quite large. But suppose we are only interested in certain branches of the tree, say the "Europe" node of Region, the "3rd" node of Class, and the "male" node of gender. Using the keep parameter a compact tree can be produced: R> vtree(td, "Region Class gender Age", + keep = list(Region = "Europe", Class = "3rd", gender = "male")) In the Region layer of Figure 11 , note that along with the "Europe" node, the missing value node has also been retained. In general, whenever valid percentages are used (which is the default), missing value nodes are retained when keep is used. This is because valid percentages are difficult to interpret without knowing the denominator, which requires knowing the Note that the missing value node for Region is no longer present, since the percentage for the "Europe" node can be interpreted without knowing how many missing values are present. Also, note that missing value node for Age includes a percentage, and the percentages for the other nodes of Age are slightly different. (With only two missing values, the difference is slight, but as the proportion of missing data increases, the percentages become substantially different.)
An alternative is to prune below the specified nodes (i.e. to prune their descendants), so that the counts always add up. In the present example, this means that the other nodes will be shown, but not their descendants. The prunebelow parameter is used to do this: R> vtree(td,"Region Age", + prunebelow = list(Region =c("UK and Ireland", "North America", "Other"))) The complement of the prunebelow parameter is the follow parameter. Instead of specifying which nodes should be pruned below, this allows you to specify which nodes should be followed (that is, not pruned below).
As a variable trees grow, it can become difficult to see the forest for the tree. For example, the following variable tree is hard to read.
R> vtree(td, "Class Region", horiz = FALSE) Figure 15 : Using the prunesmaller parameter.
Similar to the keep parameter, when valid percentages are used (vp=TRUE, which is the default), nodes represent missing values will not be pruned. As noted previously, this is because percentages are confusing when missing values are not shown. On the other hand, when vp=FALSE, missing nodes can be pruned.
Readability of a variable tree can be improved by customizing the variable and node names using the labelvar and labelnode parameters. By default, vtree labels variables and nodes exactly as they appear in the data frame. But it is often useful to change these labels.
For example, the embarked variable indicates the port where a passenger or crew member went on board the Titanic. Suppose we wish this variable to appear as Port in the variable tree. The labelvar parameter is used to do this. Figure 16 : Using the labelvar parameter.
By default, vtree labels nodes (except for the root node) using the values of the variable in question. (If the variable is a factor, the levels of the factor are used). Sometimes it is convenient to instead specify custom labels for nodes. The labelnode argument can be used to relabel the values. For example, to relabel the classes as "First Class", "Second Class", and "Third Class":
R> vtree(td,"Class", horiz = FALSE, labelnode = list(Class = c( + "First Class" = "1st", "Second Class" = "2nd", "Third Class" = "3rd"))) 
For convenience, in the call to the vtree function, you can specify variable names (separated by whitespace) in a single character string. (If, however, any of the variable names have internal spaces, the variable names must be specified as a vector of character strings.) Additionally, several modifiers can be used, as detailed below.
If an individual variable name is preceded by is.na:, that variable will be replaced by a missing value indicator in the variable tree. This facilitates exploration of missing data, for example: Figure 18 : Using the is.na: variable specification.
A variety of other specifications are available. For example <, =, and > can be used to dichotomize numeric variables. While this is a powerful tool for data exploration, a word of caution is needed. To ensure scientific rigor, it is essential that this functionality not be used to explore a variety of dichotomizations of a predictor variable in relation to the outcome variable. There is a large literature on the misuse of dichotomization and its detrimental effect on statistical inference (Altman 1994) . It is therefore recommended that any dichotomization using vtree be conducted according to a pre-specified protocol (Huebner et al. 2016 Figure 19 : Using the > variable specification.
It is often useful to display information about other variables (apart from those that define the tree) in the nodes of a variable tree. This is particularly useful for numeric variables, which generally would not be used to build the tree since they have too many distinct values. The summary parameter allows you to show information (for example, the mean of a numeric variable) within each subset of the data frame.
Suppose you are interested in summary information concerning the number of siblings/spouses aboard the Titanic (the sibsp variable) for all of the observations in the data frame (i.e. in the root node). In that case you don't need to specify any variables for the tree itself: Suppose you wish to see the mean of this variable within levels of Region and Class. (To reduce the size of this tree we'll hide the crew and the "Other" region.)
R> vtree(td, "Region Class", summary = "sibsp \nmean %mean%", horiz = FALSE, + prune = list(Region = "Other", Class = "Crew")) 
Each node in a variable tree provides the frequency of a particular combination of values of the variables. The leaf nodes represent the observed combinations of values of all of the variables. For example, in a variable tree gender nested within Class, the leaf nodes correspond to Male and Female. These combinations, or patterns, can be treated as an additional variable. And if this new pattern variable is used as the first variable in a tree, then the branches of the tree will be simplified: each branch will represent a unique pattern, with no sub-branches. A pattern tree can be easily produced by specifying pattern=TRUE. 
Study-flow diagrams provide a visual representation of how participants (or study units) meet or do not meet a sequence of inclusion criteria. These diagrams provide critical information to the reader of published study. Medical research in particular has embraced these data visualizations as part of recommended reporting guidelines. Randomized clinical trials use CONSORT diagrams to show the flow of participants through a single study (Schulz, Altman, and Moher 2010) . Systematic reviews use PRISMA flow diagrams to depict study screening (Page et al. 2020) , (Stovold et al. 2014) . While presenting study-flow diagrams is widely considered to be best practice, preparing these diagrams has traditionally been a slow, resource intensive, manual process, which has to be repeated when small changes are made to the data.
vtree uses an R data frame to make a data-driven study flow diagram. This automates the production of study-flow diagrams. As more data arrives, data cleaning changes the existing data and the analysis plan is modified after initial assessment of the data (Huebner et al. 2016) , the study-flow diagram is easily kept up to date. Not only does this increase efficiency, it minimizes the risk of introducing human error.
Consider, for example, the Remdesivir trial of Spinner et al. (2020) , in which 612 patients with confirmed severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection and moderate COVID-19 pneumonia were screened for inclusion. Although, in this case, the full data set is not publicly available, the variables required for the flow diagram can be reconstructed from Figure 1 of the published paper. The build.data.frame function built into the vtree package makes it easy to construct a data frame indicating which participants were screened, included (and of these, who was eligible, and who consented), the group participants were randomized to, and who started the intervention. (Additional details have been omitted for the sake of brevity.) R> rem <-build.data.frame( + c( "included","elig","consent","randgrp","started"), +
list ( 
Next, let's define node labels:
R> nodelabels <-list( + included=c("Randomized"="1","Excluded"="0"), + randgrp=c( + "Randomized to receive 10 d of remdesivir"="1", + "Randomized to receive 5 d of remdesivir"="2", + "Randomized to continue standard care"="3"), + started=c( + "Did not start remdesivir"="0", + "Started remdesivir"="1"))
Having set up these objects, the code to produce a CONSORT-style diagram is fairly straightforward. In particular, the follow parameter makes it easy to specify which branches of the tree should be retained.
R> vtree(rem,"included randgrp started", + labelnode=nodelabels, + follow=list(included="1",randgrp=c("1","2")), + summary=c( + "elig \nDid not meet eligibility criteria %sum%%var=includedd%%node=1%", + "consent=0 \n(Withdrew consent %sum%%var=included%%node=0%)", + "elig=0 \n(Ineligible %sum%%var=included%%node=0%)"), + cdigits=0,showvarnames=FALSE,title="patients screened", + horiz=FALSE,fillcolor="lightsteelblue1",showpct=FALSE) 
Following a 2005 racial profiling complaint to the Ontario Human Rights Commission, the Ottawa Police Service agreed to collect race data in traffic stops, known as the Traffic Stop Race Data Collection Project (TSRDCP). The TSRDCP required police officers to record their perception of the driver's race, for traffic stops over a two-year period from June 27, 2013 to June 26, 2015. A data set representing these traffic stops was made public (https: //www.ottawapolice.ca/en/news-and-community/race-archive.aspx).
Important questions concern whether some racialized or ethnic groups are stopped at a rate disproportionate to overall makeup of the population. This requires external data, not presented here. See the report by researchers at York University, dated October 2016, for a comprehensive analysis: https://www.ottawapolice.ca/en/about-us/resources/ .TSRDCP_York_Research_Report.pdf
In the York University report, some records from the raw data were removed due to errors. Additionally, since some drivers were stopped more than once, only a single report per driver was included. It was not possible to replicate this last step because driver identifiers were not included in the publicly available data set. Figure 26 : Variable tree for race=white → age → gender. Each node also shows the percentage of traffic stops with "final/no action" outcome.
Variable trees are an intuitive way to represent discrete multivariate data. The vtree package in R provides an implementation of variable trees along with a number of convenient extensions. There are a variety of other methods for displaying discrete multivariate data, and depending on the context, one of these methods be preferable. However, the simple structure of variable trees provides not only ease of interpretation but also considerable generality. We have found that variable trees facilitate iterative data exploration when a statistician is working together with a domain expert.
A key characteristic of variable trees is that the order of variables is important. Sometimes the ordering of variables is natural (e.g. school board → school → teacher), in other cases it is dictated by the research question, and in still other cases the choice of ordering is up to the analyst. Depending on the situation, this may be a strength or a weakness.
While vtree can be used to explore data, it can also be used to generate study-flow diagrams.
In recent years there has been growing concern about the "reproducibility crisis" in science (Baker 2016) . The design of vtree was influenced by the tidyverse philosophy (Wickham et al. 2019) , with its emphasis on reproducible workflows. In order to produce study-flow diagrams using vtree, all of the variables and the corresponding set of inclusion/exclusion steps must be in a single data frame, which encourages a reproducible workflow. A key barrier to the wider adoption of these diagrams has been the difficulty required to produce them. vtree facilitates reproducible research by making it easy to produce accurate study-flow diagrams.
To conclude, variable trees are an intuitive new data exploration tool for visualizing nested subsets. Applications of variable trees include revealing patterns in data, understanding missingness and producing study-flow diagrams for reproducible research.