This vignette gives you a quick introduction to data.tree applications. We took care to keep the examples simple enough so non-specialists can follow them. The price for this is, obviously, that the examples are often simple compared to real-life applications.
If you are using data.tree for things not listed here, and if you believe this is of general interest, then please do drop us a note, so we can include your application in a future version of this vignette.
This example is inspired by the examples of the treemap package.
You’ll learn how to
Aggregate
and Cumulate
Prune
methodThe original example visualizes the world population as a tree map.
library(treemap)
data(GNI2014)
treemap(GNI2014,
index=c("continent", "iso3"),
vSize="population",
vColor="GNI",
type="value")
As there are many countries, the chart gets clustered with many very small boxes. In this example, we will limit the number of countries and sum the remaining population in a catch-all country called “Other”.
We use data.tree to do this aggregation.
First, let’s convert the population data into a data.tree structure:
library(data.tree)
GNI2014$continent <- as.character(GNI2014$continent)
GNI2014$pathString <- paste("world", GNI2014$continent, GNI2014$country, sep = "/")
tree <- as.Node(GNI2014[,])
print(tree, pruneMethod = "dist", limit = 20)
## levelName
## 1 world
## 2 ¦--North America
## 3 ¦ ¦--Bermuda
## 4 ¦ ¦--United States
## 5 ¦ °--... 22 nodes w/ 0 sub
## 6 ¦--Europe
## 7 ¦ ¦--Norway
## 8 ¦ ¦--Switzerland
## 9 ¦ °--... 39 nodes w/ 0 sub
## 10 ¦--Asia
## 11 ¦ ¦--Qatar
## 12 ¦ ¦--Macao SAR, China
## 13 ¦ °--... 45 nodes w/ 0 sub
## 14 ¦--Oceania
## 15 ¦ ¦--Australia
## 16 ¦ ¦--New Zealand
## 17 ¦ °--... 11 nodes w/ 0 sub
## 18 ¦--South America
## 19 ¦ ¦--Uruguay
## 20 ¦ ¦--Chile
## 21 ¦ °--... 10 nodes w/ 0 sub
## 22 ¦--Seven seas (open ocean)
## 23 ¦ ¦--Seychelles
## 24 ¦ ¦--Mauritius
## 25 ¦ °--... 1 nodes w/ 0 sub
## 26 °--Africa
## 27 °--... 48 nodes w/ 0 sub
We can also navigate the tree to find the population of a specific
country. Luckily, RStudio is quite helpful with its code completion (use
CTRL + SPACE
):
## [1] 7604467
Or, we can look at a sub-tree:
northAm <- tree$`North America`
Sort(northAm, "GNI", decreasing = TRUE)
print(northAm, "iso3", "population", "GNI", limit = 12)
## levelName iso3 population GNI
## 1 North America NA NA
## 2 ¦--Bermuda BMU 67837 106140
## 3 ¦--United States USA 313973000 55200
## 4 ¦--Canada CAN 33487208 51630
## 5 ¦--Bahamas, The BHS 309156 20980
## 6 ¦--Trinidad and Tobago TTO 1310000 20070
## 7 ¦--Puerto Rico PRI 3971020 19310
## 8 ¦--Barbados BRB 284589 15310
## 9 ¦--St. Kitts and Nevis KNA 40131 14920
## 10 ¦--Antigua and Barbuda ATG 85632 13300
## 11 ¦--Panama PAN 3360474 11130
## 12 °--... 14 nodes w/ 0 sub NA NA
Or, we can find out what is the country with the largest GNI:
maxGNI <- Aggregate(tree, "GNI", max)
#same thing, in a more traditional way:
maxGNI <- max(sapply(tree$leaves, function(x) x$GNI))
tree$Get("name", filterFun = function(x) x$isLeaf && x$GNI == maxGNI)
## Bermuda
## "Bermuda"
We aggregate the population. For non-leaves, this will recursively
iterate through children, and cache the result in the
population
field.
tree$Do(function(x) {
x$population <- Aggregate(node = x,
attribute = "population",
aggFun = sum)
},
traversal = "post-order")
Next, we sort each node by population:
Finally, we cumulate among siblings, and store the running sum in an
attribute called cumPop
:
The tree now looks like this:
## levelName population cumPop
## 1 world 6683146875 6683146875
## 2 ¦--Asia 4033277009 4033277009
## 3 ¦ ¦--China 1338612970 1338612970
## 4 ¦ ¦--India 1166079220 2504692190
## 5 ¦ °--... 45 nodes w/ 0 sub NA NA
## 6 ¦--Africa 962382035 4995659044
## 7 ¦ ¦--Nigeria 149229090 149229090
## 8 ¦ ¦--Ethiopia 85237338 234466428
## 9 ¦ °--... 46 nodes w/ 0 sub NA NA
## 10 ¦--Europe 728669949 5724328993
## 11 ¦ ¦--Russian Federation 140041247 140041247
## 12 ¦ ¦--Germany 82329758 222371005
## 13 ¦ °--... 39 nodes w/ 0 sub NA NA
## 14 ¦--North America 528748158 6253077151
## 15 ¦ ¦--United States 313973000 313973000
## 16 ¦ ¦--Mexico 111211789 425184789
## 17 ¦ °--... 22 nodes w/ 0 sub NA NA
## 18 ¦--South America 394352338 6647429489
## 19 ¦ ¦--Brazil 198739269 198739269
## 20 ¦ ¦--Colombia 45644023 244383292
## 21 ¦ °--... 10 nodes w/ 0 sub NA NA
## 22 ¦--Oceania 33949312 6681378801
## 23 ¦ ¦--Australia 21262641 21262641
## 24 ¦ ¦--Papua New Guinea 6057263 27319904
## 25 ¦ °--... 11 nodes w/ 0 sub NA NA
## 26 °--Seven seas (open ocean) 1768074 6683146875
## 27 °--... 3 nodes w/ 0 sub NA NA
The previous steps were done to define our threshold: big countries should be displayed, while small ones should be grouped together. This lets us define a pruning function that will allow a maximum of 7 countries per continent, and that will prune all countries making up less than 90% of a continent’s population.
We would like to store the original number of countries for further use:
We are now ready to prune. This is done by defining a pruning function, returning ‘FALSE’ for all countries that should be combined:
myPruneFun <- function(x, cutoff = 0.9, maxCountries = 7) {
if (isNotLeaf(x)) return (TRUE)
if (x$position > maxCountries) return (FALSE)
return (x$cumPop < (x$parent$population * cutoff))
}
We clone the tree, because we might want to play around with different parameters:
treeClone <- Clone(tree, pruneFun = myPruneFun)
print(treeClone$Oceania, "population", pruneMethod = "simple", limit = 20)
## levelName population
## 1 Oceania 33949312
## 2 ¦--Australia 21262641
## 3 °--Papua New Guinea 6057263
Finally, we need to sum countries that we pruned away into a new “Other” node:
treeClone$Do(function(x) {
missing <- x$population - sum(sapply(x$children, function(x) x$population))
other <- x$AddChild("Other")
other$iso3 <- paste0("OTH(", x$origCount, ")")
other$country <- "Other"
other$continent <- x$name
other$GNI <- 0
other$population <- missing
},
filterFun = function(x) x$level == 2
)
print(treeClone$Oceania, "population", pruneMethod = "simple", limit = 20)
## levelName population
## 1 Oceania 33949312
## 2 ¦--Australia 21262641
## 3 ¦--Papua New Guinea 6057263
## 4 °--Other 6629408
In order to plot the treemap, we need to convert the data.tree structure back to a data.frame:
Obviously, we should also aggregate the GNI as a weighted average. Namely, we should do this for the OTH catch-all countries that we add to the tree.
In this example, we show how to display an investment portfolio as a hierarchic breakdown into asset classes. You’ll see:
Aggregate
fileName <- system.file("extdata", "portfolio.csv", package="data.tree")
pfodf <- read.csv(fileName, stringsAsFactors = FALSE)
head(pfodf)
## ISIN Name Ccy Type Duration
## 1 LI0015327682 LGT Money Market Fund (CHF) - B CHF Fund NA
## 2 LI0214880598 CS (Lie) Money Market Fund EUR EB EUR Fund NA
## 3 LI0214880689 CS (Lie) Money Market Fund USD EB USD Fund NA
## 4 LU0243957825 Invesco Euro Corporate Bond A EUR Acc EUR Fund 5.10
## 5 LU0408877412 JPM Euro Gov Sh. Duration Bd A (acc)-EUR EUR Fund 2.45
## 6 LU0376989207 Aberdeen Global Sel Emerg Mkt Bd A2 HEUR EUR Fund 6.80
## Weight AssetCategory AssetClass SubAssetClass
## 1 0.030 Cash CHF
## 2 0.060 Cash EUR
## 3 0.020 Cash USD
## 4 0.120 Fixed Income EUR Sov. and Corp. Bonds
## 5 0.065 Fixed Income EUR Sov. and Corp. Bonds
## 6 0.030 Fixed Income EUR Em. Mkts Bonds
Let us convert the data.frame to a data.tree structure. Here, we use
again the path string method. For other options, see
?as.Node.data.frame
To calculate the weight per asset class, we use the
Aggregate
method:
t <- Traverse(pfo, traversal = "post-order")
Do(t, function(x) x$Weight <- Aggregate(node = x, attribute = "Weight", aggFun = sum))
We now calculate the WeightOfParent
,
Duration is a bit more complicated, as this is a concept that applies only to the fixed income asset class. Note that, in the second statement, we are reusing the traversal from above.
We can add default formatters to our data.tree structure. Here, we add them to the root, but we might as well add them to any Node in the tree.
SetFormat(pfo, "WeightOfParent", function(x) FormatPercent(x, digits = 1))
SetFormat(pfo, "Weight", FormatPercent)
FormatDuration <- function(x) {
if (x != 0) res <- FormatFixedDecimal(x, digits = 1)
else res <- ""
return (res)
}
SetFormat(pfo, "Duration", FormatDuration)
These formatter functions will be used when printing a data.tree structure.
## levelName Weight WeightOfParent Duration
## 1 portfolio 100.00 % 0.8
## 2 ¦--Cash 11.00 % 11.0 %
## 3 ¦ ¦--CHF 3.00 % 27.3 %
## 4 ¦ ¦--EUR 6.00 % 54.5 %
## 5 ¦ °--USD 2.00 % 18.2 %
## 6 ¦--Fixed Income 28.50 % 28.5 % 3.0
## 7 ¦ ¦--EUR 26.00 % 91.2 % 3.1
## 8 ¦ ¦ ¦--Sov. and Corp. Bonds 18.50 % 71.2 % 2.4
## 9 ¦ ¦ ¦--Em. Mkts Bonds 3.00 % 11.5 % 6.8
## 10 ¦ ¦ °--High Yield Bonds 4.50 % 17.3 % 3.4
## 11 ¦ °--USD 2.50 % 8.8 % 1.6
## 12 ¦ °--High Yield Bonds 2.50 % 100.0 % 1.6
## 13 ¦--Equities 40.00 % 40.0 %
## 14 ¦ ¦--Switzerland 6.00 % 15.0 %
## 15 ¦ ¦--Euroland 14.50 % 36.2 %
## 16 ¦ ¦--US 8.10 % 20.2 %
## 17 ¦ ¦--UK 0.90 % 2.2 %
## 18 ¦ ¦--Japan 3.00 % 7.5 %
## 19 ¦ ¦--Australia 2.00 % 5.0 %
## 20 ¦ °--Emerging Markets 5.50 % 13.7 %
## 21 °--Alternative Investments 20.50 % 20.5 %
## 22 ¦--Real Estate 5.50 % 26.8 %
## 23 ¦ °--Eurozone 5.50 % 100.0 %
## 24 ¦--Hedge Funds 10.50 % 51.2 %
## 25 °--Commodities 4.50 % 22.0 %
This example shows you the following:
Thanks a lot for all the helpful comments made by Holger von Jouanne-Diedrich.
Classification trees are very popular these days. If you have never come across them, you might be interested in classification trees. These models let you classify observations (e.g. things, outcomes) according to the observations’ qualities, called features. Essentially, all of these models consist of creating a tree, where each node acts as a router. You insert your mushroom instance at the root of the tree, and then, depending on the mushroom’s features (size, points, color, etc.), you follow along a different path, until a leaf node spits out your mushroom’s class, i.e. whether it’s edible or not.
There are two different steps involved in using such a model: training (i.e. constructing the tree), and predicting (i.e. using the tree to predict whether a given mushroom is poisonous). This example provides code to do both, using one of the very early algorithms to classify data according to discrete features: ID3. It lends itself well for this example, but of course today there are much more elaborate and refined algorithms available.
During the prediction step, each node routes our mushroom according to a feature. But how do we chose the feature? Should we first separate our set according to color or size? That is where classification models differ.
In ID3, we pick, at each node, the feature with the highest Information Gain. In a nutshell, this is the feature which splits the sample in the possibly purest subsets. For example, in the case of mushrooms, dots might be a more sensible feature than organic.
The entropy is a measure of the purity of a dataset.
Mathematically, the information gain IG is defined as:
$$ IG(T,a) = H(T)-\sum_{v\in vals(a)}\frac{|\{\textbf{x}\in T|x_a=v\}|}{|T|} \cdot H(\{\textbf{x}\in T|x_a=v\}) $$
In words, the information gain measures the difference between the entropy before the split, and the weighted sum of the entropies after the split.
So, let’s rewrite that in R:
We are all set for the ID3 training algorithm.
We start with the entire training data, and with a root. Then:
For the following implementation, we assume that the classifying features are in columns 1 to n-1, whereas the class (the edibility) is in the last column.
TrainID3 <- function(node, data) {
node$obsCount <- nrow(data)
#if the data-set is pure (e.g. all toxic), then
if (IsPure(data)) {
#construct a leaf having the name of the pure feature (e.g. 'toxic')
child <- node$AddChild(unique(data[,ncol(data)]))
node$feature <- tail(names(data), 1)
child$obsCount <- nrow(data)
child$feature <- ''
} else {
#calculate the information gain
ig <- sapply(colnames(data)[-ncol(data)],
function(x) InformationGain(
table(data[,x], data[,ncol(data)])
)
)
#chose the feature with the highest information gain (e.g. 'color')
#if more than one feature have the same information gain, then take
#the first one
feature <- names(which.max(ig))
node$feature <- feature
#take the subset of the data-set having that feature value
childObs <- split(data[ ,names(data) != feature, drop = FALSE],
data[ ,feature],
drop = TRUE)
for(i in 1:length(childObs)) {
#construct a child having the name of that feature value (e.g. 'red')
child <- node$AddChild(names(childObs)[i])
#call the algorithm recursively on the child and the subset
TrainID3(child, childObs[[i]])
}
}
}
Our training data looks like this:
## color size points edibility
## 1 red small yes toxic
## 2 brown small no edible
## 3 brown large yes edible
## 4 green small no edible
## 5 red large no edible
Indeed, a bit small. But you get the idea.
We are ready to train our decision tree by running the function:
## levelName feature obsCount
## 1 mushroom color 5
## 2 ¦--brown edibility 2
## 3 ¦ °--edible 2
## 4 ¦--green edibility 1
## 5 ¦ °--edible 1
## 6 °--red size 2
## 7 ¦--large edibility 1
## 8 ¦ °--edible 1
## 9 °--small edibility 1
## 10 °--toxic 1
We need a predict function, which will route data through our tree and make a prediction based on the leave where it ends up:
This demo calculates and plots a simple decision tree. It demonstrates the following:
YAML is similar to JSON, but targeted towards humans (as opposed to computers). It’s consise and easy to read. YAML can be a neat format to store your data.tree structures, as you can use it across different software and systems, you can edit it with any text editor, and you can even send it as an email.
This is how our YAML file looks:
fileName <- system.file("extdata", "jennylind.yaml", package="data.tree")
cat(readChar(fileName, file.info(fileName)$size))
## name: Jenny Lind
## type: decision
## Sign with Movie Company:
## type: chance
## Small Box Office:
## type: terminal
## p: 0.3
## payoff: 200000
## Medium Box Office:
## type: terminal
## p: 0.6
## payoff: 1000000
## Large Box Office:
## type: terminal
## p: 0.1
## payoff: 3000000
## Sign with TV Network:
## type: chance
## Small Box Office:
## type: terminal
## p: 0.3
## payoff: 900000
## Medium Box Office:
## type: terminal
## p: 0.6
## payoff: 900000
## Large Box Office:
## type: terminal
## p: 0.1
## payoff: 900000
Let’s convert the YAML into a data.tree structure. First, we load it
with the yaml package into a list of lists. Then we use
as.Node
to convert the list into a data.tree structure:
library(data.tree)
library(yaml)
lol <- yaml.load_file(fileName)
jl <- as.Node(lol)
print(jl, "type", "payoff", "p")
## levelName type payoff p
## 1 Jenny Lind decision NA NA
## 2 ¦--Sign with Movie Company chance NA NA
## 3 ¦ ¦--Small Box Office terminal 200000 0.3
## 4 ¦ ¦--Medium Box Office terminal 1000000 0.6
## 5 ¦ °--Large Box Office terminal 3000000 0.1
## 6 °--Sign with TV Network chance NA NA
## 7 ¦--Small Box Office terminal 900000 0.3
## 8 ¦--Medium Box Office terminal 900000 0.6
## 9 °--Large Box Office terminal 900000 0.1
Next, we define our payoff function, and apply it to the tree. Note that we use post-order traversal, meaning that we calculate the tree from leaf to root:
payoff <- function(node) {
if (node$type == 'chance') node$payoff <- sum(sapply(node$children, function(child) child$payoff * child$p))
else if (node$type == 'decision') node$payoff <- max(sapply(node$children, function(child) child$payoff))
}
jl$Do(payoff, traversal = "post-order", filterFun = isNotLeaf)
The decision function is the next step. Note that we filter on decision nodes:
The data tree plotting facility uses GraphViz / DiagrammeR. You can provide a function as a style:
GetNodeLabel <- function(node) switch(node$type,
terminal = paste0( '$ ', format(node$payoff, scientific = FALSE, big.mark = ",")),
paste0('ER\n', '$ ', format(node$payoff, scientific = FALSE, big.mark = ",")))
GetEdgeLabel <- function(node) {
if (!node$isRoot && node$parent$type == 'chance') {
label = paste0(node$name, " (", node$p, ")")
} else {
label = node$name
}
return (label)
}
GetNodeShape <- function(node) switch(node$type, decision = "box", chance = "circle", terminal = "none")
SetEdgeStyle(jl, fontname = 'helvetica', label = GetEdgeLabel)
SetNodeStyle(jl, fontname = 'helvetica', label = GetNodeLabel, shape = GetNodeShape)
Note that the fontname
is inherited as is by all
children, whereas e.g. the label
argument is a function,
it’s called on each inheriting child node.
Another alternative is to set the style per node:
jl$Do(function(x) SetEdgeStyle(x, color = "red", inherit = FALSE),
filterFun = function(x) !x$isRoot && x$parent$type == "decision" && x$parent$decision == x$name)
Finally, we direct our plot from left-to-right, and use the plot function to display:
In this example, we will replicate Mike Bostock’s bubble example. See
here for details: http://bl.ocks.org/mbostock/4063269.
We use Joe Cheng’s bubbles package. All of
this is inspired by Timelyportfolio, the king
of htmlwidgets.
You’ll learn how to convert a complex JSON into a data.frame, and how to use this to plot hierarchic visualizations.
The data represents the Flare class hierarchy, which is a code library for creating visualizations. The JSON is long, deeply nested, and complicated.
fileName <- system.file("extdata", "flare.json", package="data.tree")
flareJSON <- readChar(fileName, file.info(fileName)$size)
cat(substr(flareJSON, 1, 300))
## {
## "name": "flare",
## "children": [
## {
## "name": "analytics",
## "children": [
## {
## "name": "cluster",
## "children": [
## {"name": "AgglomerativeCluster", "size": 3938},
## {"name": "CommunityStructure", "size": 3812},
## {"name": "HierarchicalCluster", "size": 6714},
## {"name
So, let’s convert it into a data.tree structure:
library(jsonlite)
flareLoL <- fromJSON(file(fileName),
simplifyDataFrame = FALSE
)
flareTree <- as.Node(flareLoL, mode = "explicit", check = "no-warn")
flareTree$attributesAll
## [1] "size"
## levelName size
## 1 flare NA
## 2 ¦--analytics NA
## 3 ¦ ¦--cluster NA
## 4 ¦ ¦ ¦--AgglomerativeCluster 3938
## 5 ¦ ¦ ¦--CommunityStructure 3812
## 6 ¦ ¦ ¦--HierarchicalCluster 6714
## 7 ¦ ¦ °--MergeEdge 743
## 8 ¦ ¦--graph NA
## 9 ¦ ¦ ¦--BetweennessCentrality 3534
## 10 ¦ ¦ ¦--LinkDistance 5731
## 11 ¦ ¦ ¦--MaxFlowMinCut 7840
## 12 ¦ ¦ ¦--ShortestPaths 5914
## 13 ¦ ¦ °--SpanningTree 3416
## 14 ¦ °--optimization NA
## 15 ¦ °--AspectRatioBanker 7074
## 16 ¦--animate NA
## 17 ¦ ¦--Easing 17010
## 18 ¦ ¦--FunctionSequence 5842
## 19 ¦ ¦--interpolate NA
## 20 ¦ ¦ ¦--ArrayInterpolator 1983
## 21 ¦ ¦ ¦--ColorInterpolator 2047
## 22 ¦ ¦ ¦--DateInterpolator 1375
## 23 ¦ ¦ ¦--Interpolator 8746
## 24 ¦ ¦ ¦--MatrixInterpolator 2202
## 25 ¦ ¦ ¦--NumberInterpolator 1382
## 26 ¦ ¦ ¦--ObjectInterpolator 1629
## 27 ¦ ¦ ¦--PointInterpolator 1675
## 28 ¦ ¦ °--RectangleInterpolator 2042
## 29 ¦ ¦--ISchedulable 1041
## 30 ¦ °--... 8 nodes w/ 0 sub NA
## 31 °--... 8 nodes w/ 215 sub NA
Finally, we can convert it into a data.frame. The
ToDataFrameTable
only converts leafs, but inherits
attributes from ancestors:
flare_df <- ToDataFrameTable(flareTree,
className = function(x) x$parent$name,
packageName = "name",
"size")
head(flare_df)
## className packageName size
## 1 cluster AgglomerativeCluster 3938
## 2 cluster CommunityStructure 3812
## 3 cluster HierarchicalCluster 6714
## 4 cluster MergeEdge 743
## 5 graph BetweennessCentrality 3534
## 6 graph LinkDistance 5731
This does not look spectacular. But take a look at this stack overflow question to see how people struggle to do this type of operation.
Here, it was particularly simple, because the underlying JSON
structure is regular. If it were not (e.g. some nodes contain different
attributes than others), the conversion from JSON to data.tree would
still work. And then, as a second step, we could modify the data.tree
structure before converting it into a data.frame. For example, we could
use Prune
and Remove
to remove unwanted nodes,
use Set
to remove or add default values, etc.
What follows has nothing to do with data.tree anymore. We simply provide the bubble chart printing for your enjoyment. In order to run it yourself, you need to install the bubbles package from github:
devtools::install_github("jcheng5/bubbles@6724e43f5e")
library(scales)
library(bubbles)
library(RColorBrewer)
bubbles(
flare_df$size,
substr(flare_df$packageName, 1, 2),
tooltip = flare_df$packageName,
color = col_factor(
brewer.pal(9,"Set1"),
factor(flare_df$className)
)(flare_df$className),
height = 800,
width = 800
)
In this example, we print the files that exist in the folder structure of the file system. As a special goodie, we’ll show code that lets you build your own R File Explorer, an interactive tree / list widget that lets you expand folders and browse through your file system.
First, let’s read the files in a directory tree into R. In this
example, the root path “..” is the parent of the vignettes
folder, i.e. the data.tree package folder itself:
path <- ".."
files <- list.files(path = path,
recursive = TRUE,
include.dirs = FALSE)
df <- data.frame(
filename = sapply(files,
function(fl) paste0("data.tree","/",fl)
),
file.info(paste(path, files, sep = "/")),
stringsAsFactors = FALSE
)
print(head(df)[c(1,2,3,4)], row.names = FALSE)
## filename size isdir mode
## data.tree/CRAN-SUBMISSION 91 FALSE 644
## data.tree/DESCRIPTION 2158 FALSE 644
## data.tree/NAMESPACE 1467 FALSE 644
## data.tree/NEWS 12112 FALSE 644
## data.tree/R/data.tree-package.R 4990 FALSE 644
## data.tree/R/data_doc.R 952 FALSE 644
We now convert this into a data.tree:
fileStructure <- as.Node(df, pathName = "filename")
fileStructure$leafCount / (fileStructure$totalCount - fileStructure$leafCount)
## [1] 9.384615
## levelName mode size
## 1 data.tree NA NA
## 2 ¦--CRAN-SUBMISSION 420 91
## 3 ¦--DESCRIPTION 420 2158
## 4 ¦--NAMESPACE 420 1467
## 5 ¦--NEWS 420 12112
## 6 ¦--R NA NA
## 7 ¦ ¦--data.tree-package.R 420 4990
## 8 ¦ ¦--data_doc.R 420 952
## 9 ¦ ¦--node.R 420 40117
## 10 ¦ ¦--node_actives.R 420 2184
## 11 ¦ ¦--node_conversion.R 420 1801
## 12 ¦ ¦--node_conversion_ape.R 420 4127
## 13 ¦ ¦--node_conversion_dataframe.R 420 15451
## 14 ¦ ¦--node_conversion_dendrogram.R 420 3748
## 15 ¦ ¦--node_conversion_igraph.R 420 1554
## 16 ¦ ¦--node_conversion_list.R 420 10479
## 17 ¦ ¦--node_conversion_party.R 420 5261
## 18 ¦ ¦--node_conversion_rpart.R 420 1801
## 19 ¦ ¦--node_methods.R 420 16378
## 20 ¦ ¦--node_methods_sideeffect.R 420 4443
## 21 ¦ ¦--node_methods_traversal.R 420 11599
## 22 ¦ ¦--node_plot.R 420 10568
## 23 ¦ ¦--register-s3.R 420 3729
## 24 ¦ ¦--release.R 420 368
## 25 ¦ °--... 2 nodes w/ 0 sub NA NA
## 26 °--... 12 nodes w/ 99 sub NA NA
Finally, we can display the files by timelyportfolio’s listviewer. As it’s not on CRAN, we only display a screenshot of the widget in in this vignette. This is not half as fun as the interactive widget, of course. So please try it out for yourself to see it in action.
#This requires listviewer, which is available only on github
devtools::install_github("timelyportfolio/listviewer")
library(listviewer)
l <- ToListSimple(fileStructure)
jsonedit(l)
(Run the code yourself to see the widget in action)
This is a simplistic example from the area of genetics. Similar models are found in many attributes, namely wherever you have multi-generation models and probabilities.
The code generates 100 simulations of a 3 generation population. Individuals can inherit or develop a certain feature (e.g. colour blindness). The probability to develop the feature is based on sex. We then plot the probability distribution of the feature in the last generation.
You’ll learn how to build a data.tree structure according to probabilistic rules, and how to use the structure to infer a probability distribution.
First, we generate a family tree of a population exhibiting a certain feature (e.g. colour blindness).
#' @param children the number of children each population member has
#' @param probSex the probability of the sex of a descendant
#' @param probInherit the probability the feature is inherited, depending on the sex of the descendant
#' @param probDevelop the probability the feature is developed (e.g. a gene defect), depending on the sex
#' of the descendant
#' @param generations the number of generations our simulated population should have
#' @param parent for recursion
GenerateChildrenTree <- function(children = 2,
probSex = c(male = 0.52, female = 0.48),
probInherit = c(male = 0.8, female = 0.5),
probDevelop = c(male = 0.05, female = 0.01),
generations = 3,
parent = NULL) {
if (is.null(parent)) {
parent <- Node$new("1")
parent$sex <- 1
parent$feature <- TRUE
parent$develop <- FALSE
}
#sex of descendants
#1 = male
#2 = female
sex <- sample.int(n = 2, size = children, replace = TRUE, prob = probSex)
for (i in 1:children) child <- parent$AddChild(i)
Set(parent$children, sex = sex)
#inherit
if (parent$feature == TRUE) {
for (i in 1:2) {
subPop <- Traverse(parent, filterFun = function(x) x$sex == i)
inherit <- sample.int(n = 2,
size = length(subPop),
replace = TRUE,
prob = c(1 - probInherit[i], probInherit[i]))
Set(subPop, feature = as.logical(inherit - 1))
}
} else {
Set(parent$children, feature = FALSE)
}
#develop
Set(parent$children, develop = FALSE)
for (i in 1:2) {
subPop <- Traverse(parent, filterFun = function(x) x$sex == i && !x$feature)
develop <- sample.int(n = 2,
size = length(subPop),
replace = TRUE,
prob = c(1 - probDevelop[i], probDevelop[i]))
Set(subPop, feature = as.logical((develop - 1)), develop = as.logical((develop - 1)))
}
#recursion to next generation
if (generations > 0) for (i in 1:children) GenerateChildrenTree(children,
probSex,
probInherit,
probDevelop,
generations - 1,
parent$children[[i]])
return (parent)
}
Just for demonstration purpose, this is what a tree looks like:
## levelName sex feature develop
## 1 1 1 TRUE FALSE
## 2 ¦--1 2 FALSE FALSE
## 3 ¦ ¦--1 1 TRUE TRUE
## 4 ¦ ¦ ¦--1 1 TRUE FALSE
## 5 ¦ ¦ ¦ ¦--1 2 FALSE FALSE
## 6 ¦ ¦ ¦ °--2 1 FALSE FALSE
## 7 ¦ ¦ °--2 2 FALSE FALSE
## 8 ¦ ¦ ¦--1 2 TRUE FALSE
## 9 ¦ ¦ °--2 2 FALSE FALSE
## 10 ¦ °--2 1 FALSE FALSE
## 11 ¦ ¦--1 2 FALSE FALSE
## 12 ¦ ¦ ¦--1 2 FALSE FALSE
## 13 ¦ ¦ °--2 2 FALSE FALSE
## 14 ¦ °--2 2 FALSE FALSE
## 15 ¦ ¦--1 2 FALSE FALSE
## 16 ¦ °--2 2 FALSE FALSE
## 17 °--2 2 FALSE FALSE
## 18 ¦--1 2 FALSE FALSE
## 19 ¦ ¦--1 1 FALSE FALSE
## 20 ¦ ¦ °--... 2 nodes w/ 0 sub NA NA NA
## 21 ¦ °--... 1 nodes w/ 4 sub NA NA NA
## 22 °--... 1 nodes w/ 11 sub NA NA NA
How big is our population after three generations?
## [1] 31
For a given tree, how many have the feature?
## [1] 4
How many males have developed the feature without inheritance?
## [1] 1
What is the occurrence of the feature in the last generation?
FreqLastGen <- function(tree) {
l <- tree$leaves
sum(sapply(l, function(x) x$feature))/length(l)
}
FreqLastGen(tree)
## [1] 0.0625
Generate 100 sample trees and get the frequency of the feature in the last generation
## user system elapsed
## 1.385 0.000 1.386
Plot a histogram of the frequency of the defect in the last generation:
For larger populations, you might consider parallelisation, of course. See below for some hints.
It is straight forward to parallelise the simulation. If, as in this example, you do not need to pass around a data.tree structure from one process (fork) to another, it is also rather efficient.
library(foreach)
library(doParallel)
registerDoParallel(makeCluster(3))
#On Linux, there are other alternatives, e.g.: library(doMC); registerDoMC(3)
system.time(x <- foreach (i = 1:100, .packages = "data.tree") %dopar% FreqLastGen(GenerateChildrenTree()))
stopImplicitCluster()
## user system elapsed
## 0.07 0.02 1.40
For the more complicated case where you want to parallelise operations on a single tree, see below.
In this example, we do a brute force solution of Tic-Tac-Toe, the well-known 3*3 game.
You’ll learn how data.tree can be used to build a tree of game history, and how the resulting data.tree structure can be used to analyze the game.
In addition, this example shows you how parallelisation can speed up data.tree.
We want to set up the problem in a way such that each
Node
is a move of a player, and each path describes the
entire history of a game.
We number the attributes from 1 to 9. Additionally, for easy readability, we label the Nodes in an Excel-like manner, such that field 9, say, is ‘c3’:
## Var1 Var2
## 1 a 1
## 2 b 1
## 3 c 1
## 4 a 2
## 5 b 2
## 6 c 2
## 7 a 3
## 8 b 3
## 9 c 3
To speed up things a bit, we consider rotation, so that, say, the first move in a3 and a1 are considered equal, because they could be achieved with a 90 degree rotation of the board. This leaves us with only a3, b3, and b2 for the first move of player 1:
ttt <- Node$new("ttt")
#consider rotation, so first move is explicit
ttt$AddChild("a3")
ttt$a3$f <- 7
ttt$AddChild("b3")
ttt$b3$f <- 8
ttt$AddChild("b2")
ttt$b2$f <- 5
ttt$Set(player = 1, filterFun = isLeaf)
Now we recurse through the tree, and add possible moves to the
leaves, growing it eventually to hold all possible games. To do this, we
define a method which, based on a Node's
path, adds
possible moves as children.
AddPossibleMoves <- function(node) {
t <- Traverse(node, traversal = "ancestor", filterFun = isNotRoot)
available <- rownames(attributes)[!rownames(attributes) %in% Get(t, "f")]
for (f in available) {
child <- node$AddChild(paste0(attributes[f, 1], attributes[f, 2]))
child$f <- as.numeric(f)
child$player <- ifelse(node$player == 1, 2, 1)
hasWon <- HasWon(child)
if (!hasWon && child$level <= 10) AddPossibleMoves(child)
if (hasWon) {
child$result <- child$player
print(paste("Player ", child$player, "wins!"))
} else if(child$level == 10) {
child$result <- 0
print("Tie!")
}
}
return (node)
}
Note that we store additional info along the way. For example, in the
line child$player <- ifelse(node$player == 1, 2, 1)
, the
player is deferred from the parent Node
, and set as an
attribute in the Node
.
Our algorithm stops whenever either player has won, or when all 9 attributes are taken. Whether a player has won is determined by this function:
HasWon <- function(node) {
t <- Traverse(node, traversal = "ancestor", filterFun = function(x) !x$isRoot && x$player == node$player)
mine <- Get(t, "f")
mineV <- rep(0, 9)
mineV[mine] <- 1
mineM <- matrix(mineV, 3, 3, byrow = TRUE)
result <- any(rowSums(mineM) == 3) ||
any(colSums(mineM) == 3) ||
sum(diag(mineM)) == 3 ||
sum(diag(t(mineM))) == 3
return (result)
}
The following code plays all possible games. Depending on your computer, this might take a few minutes:
## user system elapsed
## 345.645 3.245 346.445
What is the total number of games?
## [1] 89796
How many nodes (moves) does our tree have?
## [1] 203716
What is the average length of a game?
## [1] 8.400775
What is the average branching factor?
## [1] 1.788229
How many games were won by each player?
winnerOne <- Traverse(ttt, filterFun = function(x) x$isLeaf && x$result == 1)
winnerTwo <- Traverse(ttt, filterFun = function(x) x$isLeaf && x$result == 2)
ties <- Traverse(ttt, filterFun = function(x) x$isLeaf && x$result == 0)
c(winnerOne = length(winnerOne), winnerTwo = length(winnerTwo), ties = length(ties))
## winnerOne winnerTwo ties
## 39588 21408 28800
We can, for example, look at any Node, using the
PrintBoard
function. This function prints the game
history:
PrintBoard <- function(node) {
mineV <- rep(0, 9)
t <- Traverse(node, traversal = "ancestor", filterFun = function(x) !x$isRoot && x$player == 1)
field <- Get(t, "f")
value <- Get(t, function(x) paste0("X", x$level - 1))
mineV[field] <- value
t <- Traverse(node, traversal = "ancestor", filterFun = function(x) !x$isRoot && x$player == 2)
field <- Get(t, "f")
value <- Get(t, function(x) paste0("O", x$level - 1))
mineV[field] <- value
mineM <- matrix(mineV, 3, 3, byrow = TRUE)
rownames(mineM) <- letters[1:3]
colnames(mineM) <- as.character(1:3)
mineM
}
The first number denotes the move (1 to 9). The second number is the player:
## 1 2 3
## a "O2" "X3" "O4"
## b "X5" "O6" "X7"
## c "X1" "O8" "X9"
Exercise: Do the same for Chess!
Here, the parallelisation is more challenging as with the Gene Defect example above. The reason is that we have only one tree, albeit a big one. So we need a strategy to do what we call intra-tree parallelisation.
In a perfect world, data.tree and intra-tree parallelisation would tell a love story: Many operations are recursive, and can be called equally well on a subtree or on an entire tree. Therefore, it is very natural to delegate the calculation of multiple sub-trees to different processes.
For example, tic-tac-toe seems almost trivial to parallelise:
Remember that, on level 2, we created manually 3 Nodes
. The
creation of the sub-trees on these Nodes
will be completely
independent on the other sub-trees. Then, each sub-tree can be created
in its own process.
So, in theory, we could use any parallelisation mechanism available
in R. Unfortunately, you need to take into account a few things. As a
matter of fact, to pass the sub-trees from a fork process back to the
main process, R needs to serialize the Nodes
of the
sub-tree, and this results in huge objects. As a result, collecting the
sub-trees would take ages.
So, instead, we can
AnalyseTicTacToe <- function(subtree) {
# 1. create sub-tree
AddPossibleMoves(subtree)
# 2. run the analysis
winnerOne <- Traverse(subtree, filterFun = function(x) x$isLeaf && x$result == 1)
winnerTwo <- Traverse(subtree, filterFun = function(x) x$isLeaf && x$result == 2)
ties <- Traverse(subtree, filterFun = function(x) x$isLeaf && x$result == 0)
res <- c(winnerOne = length(winnerOne),
winnerTwo = length(winnerTwo),
ties = length(ties))
# 3. return the result
return(res)
}
library(foreach)
library(doParallel)
registerDoParallel(makeCluster(3))
#On Linux, there are other alternatives, e.g.: library(doMC); registerDoMC(3)
system.time(
x <- foreach (child = ttt$children,
.packages = "data.tree") %dopar% AnalyseTicTacToe(child)
)
## user system elapsed
## 0.05 0.04 116.86
## winnerOne winnerTwo ties
## 39588 21408 28800