--- title: "DistributionalDDA easy howto" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{DistributionalDDA easy howto} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "##", fig.width = 6, fig.height = 4 ) library(DistributionalDDA) ``` # 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 - How to load distributional data and how to store it in variables - How to display distributional data graphically - How to compute summary statistics like means and variances of distributional data. - How to use the tools of the package to apply more advanced statistical tools available in R to distributions. Another vignette of this package deals specifically with principal component analysis for distributional data. - An outlook to working with more complex distributional data like multivariate data, distributions with mixed continuous and discrete supports. # Loading and Managing Distributional data Distributional data can come in many different forms. The most common are: - Statistical samples - Histograms - Discretized cummulative distributions - Density functions (i.e. analytical or formulaic descriptions) 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: ```{r Rintrospection} library(DistributionalDDA) data(DTD) dim(DTD) names(DTD) head(DTD) ``` 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 ```{r Rintrospection2} head(DTD$data[[1]]) ``` For each of these particles there is a measurement of - size: the volume of the particle; - shape: a shape parameter with values between 0 and 1, with 1 describing perfect sphericity; - rho: the mass density of the particle. For demonstration purposes the dataset also contains functional representations of these distributions as Probability Density Functions (PDF) and cummulative distribution functions (CDF). ```{r Rintrospection3} head(DTD$shapePDF) length(DTD$shapePDF) 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, ```{r DFsubsetting} DTD["shapePDF"] ``` 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" ```{r getSubsetter} theName = format(DTD$shapePDF)[[1]] theName get(theName,ListDataEnv) ``` 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) ``` 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. ```{r B2Unif01} B2Unif01 ``` 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. ```{r defaultUnifBases} B2Unif01Basis5 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 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()`: ```{r parallelSummaries1} sapply(DTD$data, \(x) summary(x$shape)) ``` 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 ```{r parallelSummaries2} head(DTD$data) %.% summary(.$shape) ``` 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 `%..%` ```{r parallelSummaries3} # DTD$data %.% summary(.$shape) |> simplify2array() |> t() # works as well DTD$data %..% summary(.$shape) |> t() ``` this time on the whole data set. This operator '%..%' converts the resulting list into an array if all elements have the same length: ```{r parallelSummaries4} DTD$data %..% mean(.$shape) DTD$data %..% range(.$shape) ``` ## 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]$ ```{r StdRefMeasUnif01} StdReferenceMeasureUnif01 ``` and the standard normal distribution $N(0,1)$ on the whole real axis ```{r StdRefMeasN01} StdReferenceMeasureN01 ``` There are also predefined 'BayesSpaces' $B^2(Unif(0,1))$ and $B^2(N(0,1))$ for these reference measures, ```{r defaultB2spaces} B2Unif01 B2N01 ``` 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, ```{r stats2constructB2} # 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 ```{r model2constructB2} rbeta(10000,shape1=8,shape2=0.5) |> hist() ``` We can define our own reference measure and Bayes space accordingly: ```{r myRef2constructB2} 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()` ```{r constructB2} myB2 = BayesSpace(myRM) myB2 ``` 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()`, ```{r constructB2basis} 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: * `basis` contains an R function returning a vector of values of the basis functions in use for the underlying $L^2$-space; the defaults are inherited from `fda::basisfd()` and correspond to B-splines. * `A` is a matrix mapping the evaluations of the functions in `basis` to a set of zero-integral basis functions, i.e. the bases of the $L^2_0$-subspace isomorphic to the current Bayes space; * `S` is the Gram matrix, giving the matrix of scalar products between the bases of $L^2_0$ resulting of operating `A` with evaluations of the `basis` function; this object is necessary because this basis is not necessarily orthonormal, and it can be reached programmatically by means of the function `getMetricTensor()`. 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) ```{r constructB2basisPlot} 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](#bayesspace-and-bayesspacebasis). 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 * `getSample()` giving the sample representing the reference measure (possibly weighted), * `getL2Space()` to access the underlying $L^2$-space below the Bayes space of study. * `getBayesSpaceBasisDimension()` and its alias `getDimension()` (TODO: is this true???), providing the dimension of the basis. ## 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: ```{r BDataSamples0} shapeSampleList = DTD$data %.% shape shapeSampleList %.% head(.) ``` 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: ```{r BDataSamples} shapeSamples = BDataSamples(DTD$data %.% shape, basis=myBasis) shapeSamples %..% head(.) ``` we can use this object for a first plot of the data: ```{r dotPlot} 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: ```{r BDataPDFs} shapePDFs = BDataPDFs(DTD$shapePDF, basis=myBasis) head(shapePDFs) ``` 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 ```{r BDataPDFsPlot} DTD["shapePDFs"] = as.data.frame(shapePDFs) DTD DTD$shapePDFs |> head() # to save space 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": * "BDataSamples" and "BDataUnivariateSamples", to contain lists of sample data, each element of the list corresponding to one set of realisations of a distribution. "BDataUnivariateSamples" is a subclass ensuring the samples to be univariate, which is required for some further plotting and analysis methods. * "BDataHistograms", the general class, and specification "BDataFrequencyHistograms" and "BDataDensityHistograms", giving a matrix of absolute frequencies resp density values on a set of common bins partitioning the support of all distributions of the data set (rows by individual, columns by bin). * "BDataPDFs" and "BDataLogPDFs" are lists, with each element being an R "function" evaluating to the Lebesgue density function resp. log-density function of one distribution. * "BDataDensities" and "BDataLogDensities", conversely, are lists of functions, each evaluating to the Radon-Nykodim density function resp. log-density function of one distribution with respect to the reference measure of the Bayes space under consideration. The only different between these two and the previous pair is indeed the reference measure: while "BDataPDFs" give densities with respect to the Lebesgue measure, "BDataDensities" provide densities with respect to the space reference measure (and the same happens with their log-counterparts). * "BDataSplines", which give a matrix of coefficients of the distributions of the data set with respect to a series of basis functions, typically ZB-splines, but potentially useful for empirical orthogonal functions, orthogonal functions with respect to the reference measure or any other basis of a subspace of the Hilbert space structure of the Bayes space (rows by individual, columns by basis function). * "BDataGeneral", finally, is a facility to construct a distributional data set which individual distributions are not necessarily of the same data type. # 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 ```{r BDataPDFCenters} 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: ```{r BDataPDFsMVar, eval=FALSE} 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()`, ```{r classes} class(DTD$shapePDFs) class(shapeSamples) ``` This shows that the classes described [before](#data-containing-classes) 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, ```{r as.Histograms} shapeHists = as.BDataFrequencyHistograms(shapeSamples) class(shapeHists) plot(shapeHists) ``` The following diagram documents all currently existing converters ```{r convertersGraph, fig.height=6, fig.width=8, echo=FALSE} asfuns = ls(envir = asNamespace("DistributionalDDA"), pattern = "^as\\.") aslist = strsplit(asfuns, split = ".", fixed = TRUE) aslist = aslist[sapply(aslist, \(x) length(x)==3)] asfuns = sapply(aslist, \(x)x[-1]) |> t() asfuns = asfuns[grep("BData", asfuns[,1]),] library(ggraph) library(igraph) g <- graph_from_data_frame(data.frame(from=asfuns[,2], to=asfuns[,1]), directed = TRUE) ggraph(g, layout = "auto") + geom_edge_link(arrow = arrow(length = unit(3, "mm")), end_cap = circle(5, "mm"), colour = "gray40") + geom_node_point(size = 8, colour = "lightblue") + geom_node_text(aes(label = name), size = 4) + theme_void() ``` 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()` ```{r as.Splines} shapeSplineCoefs = as.BDataSplines(shapePDFs) getBayesSpaceBasis(shapePDFs) ``` ## 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 ```{r splinePlots} 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, ```{r lmSplines} l1 = lm(shapeSplineCoefs[,-1]~shapeSplineCoefs[,1]) l1 ``` 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. ```{r lmSplinesR2} summary(l1) |> sapply(\(x)x$r.squared) ``` 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()`, ```{r dendroSplines} 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 ```{r expectation, eval=FALSE} 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 ```{r otherOperators, eval=FALSE} CovarianceOperator(measure, X, Y) scalarProduct(x, y, space) ``` for all case statistics, i.e. several `Case*()` functions described in the section [Descriptive statistics](#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.