DistributionalDDA easy howto

Introduction

The package ‘’DistributionalDDA’’ is designed for the statistical analysis of density, distribution and likelihood data. It is useful if, apart of numeric or categorical covariables, your sample has properties that encode a whole distributions or datasets, e.g. particle size distributions, density distributions, or waiting time distributions. ‘’DistributionalDDA’’ uses the Bayes-Space approach to distributional data analysis. A Bayes-Space is a vector space of distributions and likelihoods. This vector space structure allows to define and compute things like means, variances, distances (for clustering), orthogonal directions (for principal component analysis), and linear models on distributions, hence to analyse distributions a lot like you would analyse multivariate datasets. Not all of this is already available in the package, but it ultimately will and the package already provides the computational infrastructure for all necessary computations.

This document provides you a basic introduction into the package on

Loading and Managing Distributional data

Distributional data can come in many different forms. The most common are:

To encapsulate these different kinds of data and work with them in always transparent ways, the package contains several classes organised in a hierarchy.

Example Data

For illustration purposes we will work with a simulated particle dataset provided with the package:

library(DistributionalDDA)
data(DTD)
dim(DTD)
## [1] 20  9
names(DTD)
## [1] "sample"   "region"   "d50"      "shapePDF" "shapeCDF" "sizePDF"  "sizeCDF" 
## [8] "rhoCDF"   "data"
head(DTD)
##   sample region       d50   shapePDF   shapeCDF    sizePDF    sizeCDF
## 1      1      C 1.2246452 f_14c409ec f_e3dbf400 f_04c0fe1c f_04c0fe1c
## 2      2      A 1.1446933 f_4d54eddf f_82186c02 f_d745d47c f_d745d47c
## 3      3      B 1.3731150 f_0cc7b3c1 f_c38b321c f_5d1da56f f_5d1da56f
## 4      4      C 0.8444512 f_0e58fdb3 f_c1147c6e f_a06822dd f_a06822dd
## 5      5      A 0.8782300 f_196642b1 f_d62ac36c f_8b9c108b f_8b9c108b
## 6      6      B 0.9435296 f_1d45abe9 f_d2092a34 f_59d8b882 f_59d8b882
##       rhoCDF       data
## 1 f_47803d5d d_350aae13
## 2 f_57519aa1 d_d2a47b67
## 3 f_72f86b3d d_ca654a5b
## 4 f_6fe499a5 d_3045f7d2
## 5 f_3688a9ee d_9ca63499
## 6 f_987aab12 d_a88ac821

This dataset contains 20 samples of particulate material from three regions (A,B,C), where for each of the there is a sample of individual particles, in the variable data

head(DTD$data[[1]])
##   sample      size     shape      rho
## 1      1 0.3629829 0.8863653 4.746585
## 2      1 1.0936880 0.2089592 2.319610
## 3      1 2.1049808 0.9721672 2.276780
## 4      1 0.5254111 0.5950116 3.257725
## 5      1 0.5260704 0.7247163 4.205198
## 6      1 1.5506693 0.8302345 4.227003

For each of these particles there is a measurement of

For demonstration purposes the dataset also contains functional representations of these distributions as Probability Density Functions (PDF) and cummulative distribution functions (CDF).

head(DTD$shapePDF)
## [[1]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3f0448>
## <environment: 0x5ad22a60fc98>
## 
## [[2]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3efb88>
## <environment: 0x5ad22a6153d8>
## 
## [[3]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3ef2c8>
## <environment: 0x5ad22a638f78>
## 
## [[4]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3eea08>
## <environment: 0x5ad22a63f4b8>
## 
## [[5]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3f1f78>
## <environment: 0x5ad22a643ae8>
## 
## [[6]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3f16b8>
## <environment: 0x5ad22a649f48>
## 
## attr(,"class")
## [1] "ListData" "list"
length(DTD$shapePDF)
## [1] 20
plot(DTD$shapePDF[[1]])

The package allows us to store such complex data in cells of a ‘data.frame’ through the ‘ListData’ functionality. If displayed as entries of a data.frame they are however not displayed in full, but by a series of coded names,

DTD["shapePDF"]
##      shapePDF
## 1  f_14c409ec
## 2  f_4d54eddf
## 3  f_0cc7b3c1
## 4  f_0e58fdb3
## 5  f_196642b1
## 6  f_1d45abe9
## 7  f_d1cbbf34
## 8  f_60637277
## 9  f_9a973a10
## 10 f_0c28ef18
## 11 f_b91738bb
## 12 f_6133b7d6
## 13 f_6f48f753
## 14 f_6c43e7a1
## 15 f_704d654c
## 16 f_57d72076
## 17 f_84223012
## 18 f_ed216a77
## 19 f_b6c363ea
## 20 f_9331b1c1

The names start with a letter, which is the first letter of the class of the object followed by an underscore, e.g. functions get an f_, data.frames a d_, etc. After that follows a hexadecimal code, hopefully unique for the object. Whenever an object is displayed in this way, you get the object by its name from the an environment called “ListDataEnv”

theName = format(DTD$shapePDF)[[1]]
theName
## [1] "f_14c409ec"
get(theName,ListDataEnv)
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3f0448>
## <environment: 0x5ad22a60fc98>

You can also attach ListDataEnv to get it directly by the names displayed, like in

> attach(ListDataEnv)
> f_3e902998

function (x) 
dbeta(x, shape1, shape2, ncp)
<bytecode: 0x636f972dc4b0>
<environment: 0x636f929f7318>

In the next sections we learn how to analyse such distributional data.

BayesSpace and BayesSpaceBasis

The key approach of “DistributionalDDA” is treating distributions as elements of a vector space, a so called Bayes Space encoded in the package as “BayesSpace” objects. A “BayesSpace” is (a typically infinite) dimensional vector space, which contains distributions as its vectors.

A Bayes space \(B(P_0)\) is defined by a reference measure \(P_0\), which forms the origin of the space, i.e. the zero vector. \(B(P_0)\) contains all distributions (and improper priors) \(P\), which have the same zero sets as \(P_0\), i.e. that have essentially the same support as \(P_0\): \(P(A)>0\Leftrightarrow P_0(A)>0\).

The package uses a version of Bayes spaces called Bayes-Hilbert-Spaces, denoted as \(B^2(P_0)\), in which distributions can be operated by means of addition and scaling, as well as by scalar product and Hilbert-norm. This space can be mapped to the subspace of zero-integrable functions of \(L^2\) by means of the clr() mapping. The package “DistributionalDDA”, thus, contains function to create not only a “BayesSpace” but also the underlying “L2Space”, each with its reference measure.

In our first example we consider the shape distribution. From the documentation of the data set above, we know that shape parameters are between 0 and 1 and thus that the support of the distribution is the interval [0,1]. A simple distribution on that interval is the Uniform distribution.

The package already provides a simple Bayes Space, B2Unif01, which has the uniform distribution on the interval [0,1] as a reference measure.

B2Unif01
## BayesSpace with reference measure:ReferenceMeasureUniform(a = 0, b = 1)

We will discuss how to define and use more advanced reference measures for this and other supports later in this vignette. The package can e.g. also deal with discrete, inflated, multivariate, and conditional distributions, if the right BayesSpace is selected.

However, many readily available statistical methods in R, like e.g. k-means clustering or linear models, cannot yet deal with infinite dimensional objects directly. The package thus also asks you to select a basis for a finite dimensional subspace of the Bayes Space that is used to approximate the distributions by a lower dimensional representation for such statistical methods.

B2Unif01Basis5
## BayesSpaceBasis.ReferenceMeasureUniform(basis = basis$referenceMeasure, 
##     bayesSpace = basis, nbasis = 5, order = 4)
plot(B2Unif01Basis5,as="PDF",main="PDF 5-basis-elements", col=2:6)

plot(B2Unif01Basis5,as="logPDF",main="logDensity 5-basis-elements", col=2:6)

B2Unif01Basis7
## BayesSpaceBasis.ReferenceMeasureUniform(basis = basis$referenceMeasure, 
##     bayesSpace = basis, nbasis = 7, order = 4)
plot(B2Unif01Basis7,as="PDF",main="PDF 7-basis-elements", col=2:8)

plot(B2Unif01Basis7,as="logPDF",main="logDensities 7-basis-elements", col=2:8)

Bayes spaces get their name from the fact that their addition is the Bayes-Rule, when interpreting the \(P_0\)-density of one the distributions as a likelihood. Also the other operations are closely related to statistics: multiplying an element of a Bayes space with a scalar is the is the same as statistical weighting and the square norm and scalar product of the Bayes-Space compute Fisher-Information.

One advantage of setting a “BayesSpace” is that we can use the Bayes Space-approach for all kinds of distributions, including those with bounded and infinite support, continuous and discrete and mixed supports, univariate or multivariate. The disadvantage is that we need to first select the right Bayes space that does actually contain the distributions in our data set.

A first step is selecting the right “ReferenceMeasure” object for the “BayesSpace” being created. This can be done by figuring out the support of our distributions. In our first example we first look at the variable shape, which can only have values between 0 and 1.

Parallel analysis of data in lists

In order to analyse datasets in a list, we can apply the same analysis, like e.g. computing a summary dataset to each element and collect the results with the R command sapply():

sapply(DTD$data, \(x) summary(x$shape))
##              [,1]      [,2]      [,3]       [,4]      [,5]       [,6]      [,7]
## Min.    0.2089592 0.3064649 0.1979523 0.04534112 0.1316759 0.04466869 0.2845216
## 1st Qu. 0.7345914 0.8086137 0.7449684 0.60241397 0.7854688 0.53967075 0.8338400
## Median  0.8602592 0.8932408 0.8534113 0.77270798 0.8892876 0.72297043 0.9108957
## Mean    0.8196383 0.8635671 0.8217759 0.72536722 0.8529305 0.68885968 0.8788934
## 3rd Qu. 0.9403832 0.9546061 0.9350704 0.90233503 0.9511142 0.87680304 0.9611216
## Max.    0.9998038 0.9998409 0.9996198 0.99994237 0.9999374 0.99853777 0.9999565
##               [,8]      [,9]     [,10]     [,11]      [,12]     [,13]     [,14]
## Min.    0.03580803 0.1689263 0.4395992 0.1794321 0.09424169 0.3740768 0.2159673
## 1st Qu. 0.69352822 0.7211050 0.8286333 0.7051685 0.67444816 0.8183901 0.7774408
## Median  0.82746728 0.8349482 0.9061280 0.8359707 0.82578524 0.9015911 0.8785918
## Mean    0.79118272 0.7995912 0.8778276 0.7944975 0.77788465 0.8730923 0.8441356
## 3rd Qu. 0.92078326 0.9286724 0.9608074 0.9264019 0.92343431 0.9572970 0.9480600
## Max.    0.99972510 0.9991586 0.9998396 0.9998221 0.99949395 0.9998627 0.9999610
##             [,15]      [,16]     [,17]     [,18]     [,19]     [,20]
## Min.    0.2247027 0.03285693 0.1253844 0.2750655 0.2951025 0.1768640
## 1st Qu. 0.8317325 0.54965215 0.6977271 0.7930465 0.8139313 0.7027430
## Median  0.9086443 0.73313011 0.8317692 0.8920404 0.9022919 0.8426263
## Mean    0.8774800 0.69269103 0.7898380 0.8568793 0.8704642 0.7983515
## 3rd Qu. 0.9571851 0.87816431 0.9252380 0.9515563 0.9583002 0.9286115
## Max.    0.9996288 0.99996312 0.9997938 0.9998649 0.9989510 0.9998713

Note the usage of a lambda expression \(x)f(x), available since R v4.1, to succinctly define an anonymous function for the second argument of sapply().

As applying the same operation to each member of a list is a very common operation in distributional data analysis, the package provides an abbreviation operator %.% for this operation

head(DTD$data) %.% summary(.$shape)
## [[1]]
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2090  0.7346  0.8603  0.8196  0.9404  0.9998 
## 
## [[2]]
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3065  0.8086  0.8932  0.8636  0.9546  0.9998 
## 
## [[3]]
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1980  0.7450  0.8534  0.8218  0.9351  0.9996 
## 
## [[4]]
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.04534 0.60241 0.77271 0.72537 0.90233 0.99994 
## 
## [[5]]
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1317  0.7855  0.8893  0.8529  0.9511  0.9999 
## 
## [[6]]
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.04467 0.53967 0.72297 0.68886 0.87680 0.99854

Here we have reported only the first 6 elements to save space. The operator %.% applies the given expression on the right hand side to each element of the list on the left hand side, in the spirit of the pipe %>% from package “magrittr”. This operator %.% can access fields of elements in lists, as well as evaluate any expression given, where ‘.’ serves as a drop-in replacement for the elements of the list. If users want to compile all this information into a matrix, often the functions simplify2array() and t() can be used to obtain a cleaner output, or else, make use of the alternative operator %..%

# DTD$data %.% summary(.$shape) |> simplify2array() |> t() # works as well
DTD$data %..% summary(.$shape) |> t()
##             Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
##  [1,] 0.20895922 0.7345914 0.8602592 0.8196383 0.9403832 0.9998038
##  [2,] 0.30646493 0.8086137 0.8932408 0.8635671 0.9546061 0.9998409
##  [3,] 0.19795226 0.7449684 0.8534113 0.8217759 0.9350704 0.9996198
##  [4,] 0.04534112 0.6024140 0.7727080 0.7253672 0.9023350 0.9999424
##  [5,] 0.13167587 0.7854688 0.8892876 0.8529305 0.9511142 0.9999374
##  [6,] 0.04466869 0.5396708 0.7229704 0.6888597 0.8768030 0.9985378
##  [7,] 0.28452164 0.8338400 0.9108957 0.8788934 0.9611216 0.9999565
##  [8,] 0.03580803 0.6935282 0.8274673 0.7911827 0.9207833 0.9997251
##  [9,] 0.16892626 0.7211050 0.8349482 0.7995912 0.9286724 0.9991586
## [10,] 0.43959923 0.8286333 0.9061280 0.8778276 0.9608074 0.9998396
## [11,] 0.17943206 0.7051685 0.8359707 0.7944975 0.9264019 0.9998221
## [12,] 0.09424169 0.6744482 0.8257852 0.7778847 0.9234343 0.9994939
## [13,] 0.37407678 0.8183901 0.9015911 0.8730923 0.9572970 0.9998627
## [14,] 0.21596729 0.7774408 0.8785918 0.8441356 0.9480600 0.9999610
## [15,] 0.22470265 0.8317325 0.9086443 0.8774800 0.9571851 0.9996288
## [16,] 0.03285693 0.5496521 0.7331301 0.6926910 0.8781643 0.9999631
## [17,] 0.12538437 0.6977271 0.8317692 0.7898380 0.9252380 0.9997938
## [18,] 0.27506546 0.7930465 0.8920404 0.8568793 0.9515563 0.9998649
## [19,] 0.29510252 0.8139313 0.9022919 0.8704642 0.9583002 0.9989510
## [20,] 0.17686400 0.7027430 0.8426263 0.7983515 0.9286115 0.9998713

this time on the whole data set. This operator ‘%..%’ converts the resulting list into an array if all elements have the same length:

DTD$data %..% mean(.$shape) 
##  [1] 0.8196383 0.8635671 0.8217759 0.7253672 0.8529305 0.6888597 0.8788934
##  [8] 0.7911827 0.7995912 0.8778276 0.7944975 0.7778847 0.8730923 0.8441356
## [15] 0.8774800 0.6926910 0.7898380 0.8568793 0.8704642 0.7983515
DTD$data %..% range(.$shape)
##           [,1]      [,2]      [,3]       [,4]      [,5]       [,6]      [,7]
## [1,] 0.2089592 0.3064649 0.1979523 0.04534112 0.1316759 0.04466869 0.2845216
## [2,] 0.9998038 0.9998409 0.9996198 0.99994237 0.9999374 0.99853777 0.9999565
##            [,8]      [,9]     [,10]     [,11]      [,12]     [,13]     [,14]
## [1,] 0.03580803 0.1689263 0.4395992 0.1794321 0.09424169 0.3740768 0.2159673
## [2,] 0.99972510 0.9991586 0.9998396 0.9998221 0.99949395 0.9998627 0.9999610
##          [,15]      [,16]     [,17]     [,18]     [,19]     [,20]
## [1,] 0.2247027 0.03285693 0.1253844 0.2750655 0.2951025 0.1768640
## [2,] 0.9996288 0.99996312 0.9997938 0.9998649 0.9989510 0.9998713

Selecting a Bayes-Space and a BayesSpaceBasis

The package “DistributionalDDA” provides a selection of ‘ReferenceMeasure’ classes and objects. Two standard reference measures are the uniform distribution on \([0,1]\)

StdReferenceMeasureUnif01
## ReferenceMeasureUniform(a = 0, b = 1)

and the standard normal distribution \(N(0,1)\) on the whole real axis

StdReferenceMeasureN01
## ReferenceMeasureTransformed(cdf = "pnorm", mean = 0, sd = 1)

There are also predefined ‘BayesSpaces’ \(B^2(Unif(0,1))\) and \(B^2(N(0,1))\) for these reference measures,

B2Unif01
## BayesSpace with reference measure:ReferenceMeasureUniform(a = 0, b = 1)
B2N01
## BayesSpace with reference measure:ReferenceMeasureTransformed(cdf = "pnorm", mean = 0, sd = 1)

However, as mentioned above, to select the Bayes Space we need to keep in mind that the reference measure should: (i) have the same support as the distributions you are analysing and (ii) ideally should also have mass, where the mass of the analysed distributions are.

In case of the observed shapes, we know that the support is \([0,1]\), which would suggest that B2Unif01 might be a reasonable choice, but we can see in a global histogram (obtained by unlist-ing the set of selected values), that the mass is concentrated towards the left border of the interval,

# hist(unlist(lapply(DTD$data, function(x) x$shape))) # also works
DTD$data %.% shape |> unlist() |> hist()

It would thus be more sensible to have a reference measure with a similar property, like e.g. a beta distribution with corresponding shape parameters

rbeta(10000,shape1=8,shape2=0.5) |> hist()

We can define our own reference measure and Bayes space accordingly:

myRM = ReferenceMeasureTransformed("pbeta",shape1=8,shape2=0.5)

The function ReferenceMeasureTransformed() allows to transform the reference of the standard BayesSpace B2Unif01 into to any continuous distribution by means of the cdf, quantile function and density.

Alternatively, ReferenceMeasureUniform() allows to construct uniform reference measures on interval supports \([a,b]\) different from \([0,1]\) and there are more “ReferenceMeasure” classes available. These will soon include Reference Measures also for Discrete, Mixed and Multivariate Supports.

Once a “ReferenceMeasure” object is available, the corresponding “BayesSpace” object is created just with a call to BayesSpace()

myB2 = BayesSpace(myRM)
myB2
## BayesSpace with reference measure:ReferenceMeasureTransformed(cdf = "pbeta", shape1 = 8, shape2 = 0.5)

To exploit fast, approximate computations using subspace basis of the \(B^2(P_0)\) Bayes-Hilbert-Space we can simply construct such a basis with command BayesSpaceBasis(),

myBasis = BayesSpaceBasis(myB2, nbasis=7)

This constructs a 7 dimensional subspace basis. The documentation of its parameters together with the documentation for the uniform reference measure, called by ?BayesSpaceBasis.ReferenceMeasureUniform. Users can also implement custom basis, check ? BayesSpaceBasis.function for details.

A “BayesSpaceBasis” object (as created by a call to function BayesSpaceBasis()) contains three fundamental elements:

Thus, in spite of its name, “BayesSpaceBasis” can be used to encapsulate other bases of functions besides B-splines, by using fitting values for basis, A and S.

To understand the basis we have just created we can plot it in different forms, e.g. as distribution densities. In their natural form as zero-integral functions (spanning the image of the clr() mapping)

plot(myBasis, as="is", col=2:8)

as probability densities (with as="PDF"), as log-densities (with as="logPDF") or any of the representations mentioned above.

Users can implement custom basis. See the help to the ? BayesSpaceBasis.function for details.

–>

Finally, some words about accessing support, reference measure, Bayes space and Bayes space basis: these are accessible from all objects of the classes mentioned in the next section by means of the functions getSupport(), getReferenceMeasure(), getBayesSpace() and getBayesSpaceBasis(). Additional introspection functions are

Data Containing Classes

With this preparation we can now include our data into a B2-space. Datasets of such objects are all of class “BData” but can have a specific subclass based on how the data is provided and stored. For data given as samples this is done using the “BDataSamples” class. We first store the samples in a list:

shapeSampleList = DTD$data %.% shape
shapeSampleList %.% head(.)
## [[1]]
## [1] 0.8863653 0.2089592 0.9721672 0.5950116 0.7247163 0.8302345
## 
## [[2]]
## [1] 0.7758731 0.9425645 0.7922046 0.9624506 0.7172939 0.9847528
## 
## [[3]]
## [1] 0.7405442 0.9747334 0.5204016 0.8888155 0.8611618 0.7921460
## 
## [[4]]
## [1] 0.5186528 0.8765541 0.5731402 0.7010106 0.9430818 0.6678860
## 
## [[5]]
## [1] 0.9324004 0.8724167 0.9333947 0.8824227 0.9674738 0.8710433
## 
## [[6]]
## [1] 0.3290565 0.3698881 0.8721967 0.6832863 0.9870057 0.7152036
## 
## [[7]]
## [1] 0.8940945 0.7371354 0.9527969 0.9410313 0.9105519 0.9364430
## 
## [[8]]
## [1] 0.8290853 0.9450872 0.5273792 0.3732960 0.7600679 0.3693583
## 
## [[9]]
## [1] 0.9822911 0.3870021 0.6618772 0.9552959 0.8671555 0.9743838
## 
## [[10]]
## [1] 0.9124911 0.9885316 0.8546842 0.9417093 0.5900723 0.9492331
## 
## [[11]]
## [1] 0.4692731 0.5042157 0.9925732 0.7546001 0.7834937 0.6461458
## 
## [[12]]
## [1] 0.8307735 0.7500557 0.8085610 0.5660862 0.8993581 0.8017951
## 
## [[13]]
## [1] 0.8877072 0.8891281 0.7402295 0.8998500 0.9328819 0.9901007
## 
## [[14]]
## [1] 0.9941455 0.8646323 0.9968254 0.8913990 0.9367338 0.7619233
## 
## [[15]]
## [1] 0.9460251 0.9608802 0.8349781 0.9245737 0.8386814 0.8207921
## 
## [[16]]
## [1] 0.4926794 0.4913018 0.7516545 0.4737378 0.9583709 0.8674214
## 
## [[17]]
## [1] 0.9323455 0.6278002 0.5012123 0.9323985 0.5791718 0.8172869
## 
## [[18]]
## [1] 0.7698503 0.7855638 0.9058426 0.8173555 0.8935076 0.9968967
## 
## [[19]]
## [1] 0.8865360 0.9359341 0.9940524 0.8421927 0.9549555 0.8379152
## 
## [[20]]
## [1] 0.3624885 0.8014512 0.8356658 0.9906906 0.9031482 0.9536974

We now want to understand these values as a dataset of samples describing a distribution in our BayesSpace “myB2”, which we would like to analyse with basis “myBasis”. For that we construct our first BData object, a BDataSamples-Object:

shapeSamples = BDataSamples(DTD$data %.% shape, basis=myBasis)
shapeSamples %..% head(.) 
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
## [1,] 0.8863653 0.7758731 0.7405442 0.5186528 0.9324004 0.3290565 0.8940945
## [2,] 0.2089592 0.9425645 0.9747334 0.8765541 0.8724167 0.3698881 0.7371354
## [3,] 0.9721672 0.7922046 0.5204016 0.5731402 0.9333947 0.8721967 0.9527969
## [4,] 0.5950116 0.9624506 0.8888155 0.7010106 0.8824227 0.6832863 0.9410313
## [5,] 0.7247163 0.7172939 0.8611618 0.9430818 0.9674738 0.9870057 0.9105519
## [6,] 0.8302345 0.9847528 0.7921460 0.6678860 0.8710433 0.7152036 0.9364430
##           [,8]      [,9]     [,10]     [,11]     [,12]     [,13]     [,14]
## [1,] 0.8290853 0.9822911 0.9124911 0.4692731 0.8307735 0.8877072 0.9941455
## [2,] 0.9450872 0.3870021 0.9885316 0.5042157 0.7500557 0.8891281 0.8646323
## [3,] 0.5273792 0.6618772 0.8546842 0.9925732 0.8085610 0.7402295 0.9968254
## [4,] 0.3732960 0.9552959 0.9417093 0.7546001 0.5660862 0.8998500 0.8913990
## [5,] 0.7600679 0.8671555 0.5900723 0.7834937 0.8993581 0.9328819 0.9367338
## [6,] 0.3693583 0.9743838 0.9492331 0.6461458 0.8017951 0.9901007 0.7619233
##          [,15]     [,16]     [,17]     [,18]     [,19]     [,20]
## [1,] 0.9460251 0.4926794 0.9323455 0.7698503 0.8865360 0.3624885
## [2,] 0.9608802 0.4913018 0.6278002 0.7855638 0.9359341 0.8014512
## [3,] 0.8349781 0.7516545 0.5012123 0.9058426 0.9940524 0.8356658
## [4,] 0.9245737 0.4737378 0.9323985 0.8173555 0.8421927 0.9906906
## [5,] 0.8386814 0.9583709 0.5791718 0.8935076 0.9549555 0.9031482
## [6,] 0.8207921 0.8674214 0.8172869 0.9968967 0.8379152 0.9536974

we can use this object for a first plot of the data:

plot(shapeSamples)

All “BData” classes are in plural because they always contain multiple objects, in this case multiple samples. It is sufficient to provide the basis, since the basis knows its Bayes space and its reference measure.

In the same way we can transfer other representations of our distributional data into “BData” objects for further statistical analysis. E.g. the column ‘shapePDF’ contains functions representing the PDF. These can be associated with our BayesSpace in the following form:

shapePDFs = BDataPDFs(DTD$shapePDF, basis=myBasis)
head(shapePDFs)
## [[1]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3f0448>
## <environment: 0x5ad22a60fc98>
## 
## [[2]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3efb88>
## <environment: 0x5ad22a6153d8>
## 
## [[3]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3ef2c8>
## <environment: 0x5ad22a638f78>
## 
## [[4]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3eea08>
## <environment: 0x5ad22a63f4b8>
## 
## [[5]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3f1f78>
## <environment: 0x5ad22a643ae8>
## 
## [[6]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3f16b8>
## <environment: 0x5ad22a649f48>
## 
## attr(,"class")
## [1] "BDataPDFs" "BData"    
## attr(,"basis")
## BayesSpaceBasis.ReferenceMeasureTransformed(basis = basis$referenceMeasure, 
##     bayesSpace = basis, nbasis = 7)
## attr(,"closed")
## [1] FALSE
## attr(,"centered")
## [1] FALSE

This creates a an object of class “BDataPDFs”. The final part of the printed message reports information about the Bayes space structure in use, including reference measure and basis dimension,

BayesSpaceBasis.ReferenceMeasureTransformed(basis = basis$referenceMeasure, 
     bayesSpace = basis, nbasis = 7)

In their “BData” form the columns can again be added to the data.frame, where they retain knowledge of their class, Bayes space, reference measure and basis

DTD["shapePDFs"] = as.data.frame(shapePDFs)
DTD
##    sample region       d50   shapePDF   shapeCDF    sizePDF    sizeCDF
## 1       1      C 1.2246452 f_14c409ec f_e3dbf400 f_04c0fe1c f_04c0fe1c
## 2       2      A 1.1446933 f_4d54eddf f_82186c02 f_d745d47c f_d745d47c
## 3       3      B 1.3731150 f_0cc7b3c1 f_c38b321c f_5d1da56f f_5d1da56f
## 4       4      C 0.8444512 f_0e58fdb3 f_c1147c6e f_a06822dd f_a06822dd
## 5       5      A 0.8782300 f_196642b1 f_d62ac36c f_8b9c108b f_8b9c108b
## 6       6      B 0.9435296 f_1d45abe9 f_d2092a34 f_59d8b882 f_59d8b882
## 7       7      A 1.0222856 f_d1cbbf34 f_1e873ee9 f_0a216287 f_0a216287
## 8       8      C 1.0695661 f_60637277 f_af2ff3aa f_e03a5a17 f_e03a5a17
## 9       9      B 1.2648387 f_9a973a10 f_55dbbbcd f_41144f30 f_41144f30
## 10     10      A 1.0292579 f_0c28ef18 f_c3646ec5 f_a2c4ac4c f_a2c4ac4c
## 11     11      B 1.1020401 f_b91738bb f_765bb966 f_eb2623ee f_eb2623ee
## 12     12      B 1.0638129 f_6133b7d6 f_ae7f360b f_34f0ba19 f_34f0ba19
## 13     13      A 1.2413247 f_6f48f753 f_a004768e f_df72ad48 f_df72ad48
## 14     14      A 0.8427662 f_6c43e7a1 f_a30f667c f_1f79b9a3 f_1f79b9a3
## 15     15      A 1.1766068 f_704d654c f_bf01e491 f_d7b8ed37 f_d7b8ed37
## 16     16      C 0.9597744 f_57d72076 f_989ba1ab f_c4589f21 f_c4589f21
## 17     17      C 1.0532792 f_84223012 f_4b6eb1cf f_f89aab09 f_f89aab09
## 18     18      A 0.9784898 f_ed216a77 f_226debaa f_9f8ec950 f_9f8ec950
## 19     19      A 0.9812755 f_b6c363ea f_798fe237 f_0162ffd9 f_0162ffd9
## 20     20      B 0.9910797 f_9331b1c1 f_5c7d301c f_d593d62d f_d593d62d
##        rhoCDF       data  shapePDFs
## 1  f_47803d5d d_350aae13 f_14c409ec
## 2  f_57519aa1 d_d2a47b67 f_4d54eddf
## 3  f_72f86b3d d_ca654a5b f_0cc7b3c1
## 4  f_6fe499a5 d_3045f7d2 f_0e58fdb3
## 5  f_3688a9ee d_9ca63499 f_196642b1
## 6  f_987aab12 d_a88ac821 f_1d45abe9
## 7  f_63a47087 d_c6834ea2 f_d1cbbf34
## 8  f_4ae1c23d d_7a8cba7c f_60637277
## 9  f_694fee3e d_87ecfc19 f_9a973a10
## 10 f_743aee26 d_943d02b1 f_0c28ef18
## 11 f_e39ca658 d_ff4b94aa f_b91738bb
## 12 f_295e7731 d_cb731d95 f_6133b7d6
## 13 f_babeef4b d_af090f90 f_6f48f753
## 14 f_0f237f8a d_2ca74546 f_6c43e7a1
## 15 f_7b369128 d_fe5cdfea f_704d654c
## 16 f_574bc28a d_6591798a f_57d72076
## 17 f_4c105198 d_e96405aa f_84223012
## 18 f_2ffd689d d_987700d9 f_ed216a77
## 19 f_23268ec1 d_27a68fe3 f_b6c363ea
## 20 f_5c9b4265 d_4da5cdfa f_9331b1c1
DTD$shapePDFs |> head() # to save space
## [[1]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3f0448>
## <environment: 0x5ad22a60fc98>
## 
## [[2]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3efb88>
## <environment: 0x5ad22a6153d8>
## 
## [[3]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3ef2c8>
## <environment: 0x5ad22a638f78>
## 
## [[4]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3eea08>
## <environment: 0x5ad22a63f4b8>
## 
## [[5]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3f1f78>
## <environment: 0x5ad22a643ae8>
## 
## [[6]]
## function (x) 
## stats::dbeta(x, shape1, shape2, ncp)
## <bytecode: 0x5ad22b3f16b8>
## <environment: 0x5ad22a649f48>
## 
## attr(,"class")
## [1] "BDataPDFs" "BData"    
## attr(,"basis")
## BayesSpaceBasis.ReferenceMeasureTransformed(basis = basis$referenceMeasure, 
##     bayesSpace = basis, nbasis = 7)
## attr(,"closed")
## [1] FALSE
## attr(,"centered")
## [1] FALSE
plot(DTD$shapePDFs)

In their “BData” form we can now convert the data between different representations, which might be required to apply certain algorithms. The following representations are available, each implemented as a subclass of “BData”:

Basic analysis and change of representation

Descriptive statistics

In “BData” form we can now start to operate with these distributional objects within the Bayes space structire, for instance to compute summaries. A simple example would be the mean or the median

shapeDistMean = mean(DTD$shapePDFs)
shapeDistMedian = median(DTD$shapePDFs)
plot(DTD$shapePDFs)
plot(shapeDistMean, add=TRUE, lwd=4, col="green")
plot(shapeDistMedian, add=TRUE, lty=2, lwd=3, col="red")

or the metric variance:

mvar(DTD$shapePDFs)

For clarification, these statistics operate on a sample of distributions and produce distributional statistics: the mean() and median() are distributions themselves, while the metric variance, mvar(), is a scalar representing the average distance between each distribution of the sample and the mean distribution. If you wish to obtain statistics of the random variables described by each of these distributions, functions such as caseMeans(), caseSds(), caseVars() or caseQuantiles() computes (numerical approximations of) mean values, standard deviations, variance or quantiles.

Change of representation

More complex statistics and methods than these often require specific representation. Indeed, the selection of the representation often determines the exact algorithm used for certain computations, which can have considerable influence on execution speed and numerical precision. Not all representations will be suitable for all kinds of analysis.

The current representation of a “BData” object can be obtained by using the function class(),

class(DTD$shapePDFs)
## [1] "BDataPDFs" "BData"
class(shapeSamples)
## [1] "BDataUnivariateSamples" "BDataSamples"           "BData"

This shows that the classes described before are organised in a hierarchy, with “BDataPDFs” and “BDataSamples” being subclasses of “BData”, and “BDataUnivariateSamples” being a subclass of “BDataSamples”.

The package provides a(n almost complete) set of converters, implemented following the standard as.NEWCLASS(object) approach,

shapeHists = as.BDataFrequencyHistograms(shapeSamples)
class(shapeHists)
## [1] "BDataFrequencyHistograms" "BDataHistograms"         
## [3] "BData"
plot(shapeHists)

The following diagram documents all currently existing converters

## Loading required package: ggplot2
## 
## Attaching package: 'igraph'
## The following object is masked from 'package:compositions':
## 
##     normalize
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
## Using "stress" as default layout

This includes as well a couple of converters from base R classes (e.g., “list”, “histograms”), which can be used to create “BData” objects from scratch.

In line with current trends in functional data analysis, one of the main ideas of “distributionalDDA” is to define a subspace of the Bayes space and represent the distributions as coordinates in that subspace. The representation “BDataSplines” does exactly that, picking the basis given by getBayesSpaceBasis()

shapeSplineCoefs = as.BDataSplines(shapePDFs)
getBayesSpaceBasis(shapePDFs)
## BayesSpaceBasis.ReferenceMeasureTransformed(basis = basis$referenceMeasure, 
##     bayesSpace = basis, nbasis = 7)

Statistics

Now that the distributions are represented as a discrete set of real coefficients, statistical methods can be used to analyse them. For instance, the following code

plot(shapeSplineCoefs, col=2:6, lty=1:4, xlab="basis coefficients")

produces a parallel plot, where each individual distribution corresponds to a line connecting the values of the coordinates of that distribution. Other useful diagrams could be obtained with boxplot(), stars() or pairs().

Statistical analysis can also be applied to the coordinates. For instance,

l1 = lm(shapeSplineCoefs[,-1]~shapeSplineCoefs[,1])
l1
## 
## Call:
## lm(formula = shapeSplineCoefs[, -1] ~ shapeSplineCoefs[, 1])
## 
## Coefficients:
##                        [,1]     [,2]     [,3]     [,4]     [,5]     [,6]   
## (Intercept)            0.09247  0.20409  0.31401  0.38227  0.28037  0.20736
## shapeSplineCoefs[, 1]  1.20931  1.63344  1.23328  0.84365  0.37388  0.10774

fits a regression line for each coordinate as a function of the first coordinate. Collecting all \(R^2\) coefficients we see that the data set is largely one-dimensional.

summary(l1) |> sapply(\(x)x$r.squared)
## Response Y1 Response Y2 Response Y3 Response Y4 Response Y5 Response Y6 
##   0.9972544   0.9925996   0.9698541   0.9099060   0.7864177   0.3578453

This is supported by a look at a pairs() plot (not included to save space). A full principal component analysis is left for a specific vignette.

Distance based methods such as hierarchical clustering are also readily available via methods for dist(),

hc = shapeSplineCoefs |> dist() |> hclust()
plot(hc)

Similar strategies can be used for many other statistical methods, either: (i) extract a distance and use a metric method (such as hierarchical clustering, multidimensional scaling or Frechet moments, etc); or (ii) compute coefficients with respect to the pre-specified basis, derive covariance matrices within these coefficients and between them and external covariables and compute statistical outcomes from these covariances (as in generalized linear models, principal component analysis, linear or quadratic discriminant analysis, etc).

Computations with distributions

Generally, both the calculation of coordinates as well as the calculation of distances require the same machinery: computing or approximating integrals. The function

expectation(measure, f, ..., method)

can be considered the workhorse of all these calculations, as it calculates \[ TODO: decide-if-we-put-formula-E[f(x)]-here \] the expected value of the function f(...) with respect to the measure measure, with method controlling the way this integral is computed or approximated (the current default method="sample" follows an importance sampling approach). The same idea underlies the calculations in functions such as

CovarianceOperator(measure, X, Y)
scalarProduct(x, y, space)

for all case statistics, i.e. several Case*() functions described in the section Descriptive statistics, or even for the intermediate calculations for clr(), norm(), dist(), cen() or clo().

Besides this, note that arithmetic operators +, -, * and / are overwritten to work with “BData” as perturbation (+, -) and powering (*, /) corresponding to the operations giving a vector space structure to the Bayes space.