compositions' v2.0: R classes for compositional analysis

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

The "compositions" package

The package "compositions" was released in 2005 to enable R users the analysis of compositional data in R. The main goals were twofold: to provide functions for the logratio analysis of compositions, and the facilitation of the comparison of results and interpretation among several analysis paradigms, in a transparent and consistent way. These paradigms correspond to different statistical scales for the data, thus for different geometries of the underlying sample space (the simplex). Each geometry was encoded in an S3 class. The following table summarizes these classes. Details on how to choose a class, and how to work with them can be found in UsingCompositions.pdf. Readers new to "compositions" are recommended to start there, as this vignette rather focuses on what has changed and how to do certain advanced analysis with the new functionality.

tbl = matrix(c("no", "no", "no", "no", "yes", "no",  "no", "no", "yes", "yes", "yes", "--",   
               "no", "yes","no", "yes","yes", "no",
               "real compositional", "relative compositional", "real amounts", "relative amounts", "count multinomial", "real multivariate",
               "rcomp", "acomp", "rplus", "aplus", "ccomp", "rmult"), ncol=5)
colnames(tbl) = c("counts?", "total?", "relative?", "geometry", "class")
knitr::kable(tbl)

In this table, the first three columns correspond to three binary questions that guide you in choosing the right geometry for your data

Compositional classes, old and new

Once you have chosen your geometry (or a couple of them that you want to compare), you can load "compositions" with

library(compositions)

and load your data with any data reading function (read.csv(), read_delim(), read_excel(), etc). Here we use for illustration purposes the data set "Hydrochem"

data("Hydrochem")
dim(Hydrochem)
colnames(Hydrochem)

having 5 descriptive variables and 14 chemical components (from H to TOC). These last form a composition, which we preliminarily consider of an aplus geometry

xp = aplus(Hydrochem[, c("Ca","K", "Na", "Mg")])
summary(xp)
plot(xp, col=Hydrochem$River, cex=0.5)

The object is of course an "aplus" object

is.aplus(xp)    # check with S3 paradigm
is(xp, "aplus") # check with S4 paradigm

but it is also an amounts object

is(xp, "amounts") # only S4 paradigm

The "amounts" class is an abstract superclass, that cannot be created directly. Also objects of class "rplus" and "ccomp" are amounts

is(rplus(xp), "amounts")
is(ccomp(xp), "amounts")

If we consider the total sum of each row to be irrelevant, the appropriate class is an "acomp"

xc = acomp(Hydrochem[, c("Ca","K", "Na", "Mg")])
summary(xc)
plot(xc, col=Hydrochem$River, cex=0.5)

This can be thus checked

is.acomp(xc)    # check with S3 paradigm
is(xc, "acomp") # check with S4 paradigm

There is also an abstract superclass for compositional objects, i.e. those where the total is irrelevant or an artifact

is(xc, "compositional") 

that contains objects of class "acomp", "rcomp" and "ccomp"

is(rcomp(xc), "compositional")
is(ccomp(xc), "compositional")

Additionally, all comopostional classes have natural ways to be automatically converted to "data.frame" and to "structure" classes. This works like this

as(xc, "data.frame")
as(xc, "structure")

and it allows S4 classes with slots expecting a "data.frame" or one of "vector", "matrix" or "array" to accept a compositional object. This will be discussed in the last section of this vignette.

Subsetting and transformations

Compositional data sets react now to the dollar notation. Since compositions v2 you can extract one particular variable with

summary(xp$Ca)

As discussed in UsingCompositions.pdf, the key idea of compositional analysis is to transform the data prior to the analysis. All transformations are now accessible using the dollar notation

summary(xp$clr)

This includes each of alr, ilr, clr, log, ilt, iit, apt, ipt, cpt, their generic versions cdt and idt, as well as all their inverses (e.g. cptInv or alrInv). Also raw, clo, unclass and pwlr resp. pwlrInv (pairwise logratio and its inverse) work with this trick

summary(xp)
summary(cdt(xp))
summary(xp$cdt)
summary(xp$cdt$cdtInv)

Another (potentially conflictive) novelty is the [-subsetting, in particular for closed classes ("acomp" and "rcomp"). In compositions-v1, subseting invariably destroyed the class, returning the selected rows or columns in a matrix, or even in a vector if only one row or one column was selected. From compositions-v2 on, the nature of the subsetting depends on class of the parent object and the number of selected elements:

xp0 = xp[1:5,]
xp0
xc0 = xc[1:5,]
xc0
xc[10,]
xp0[,"Ca"]
xp0[,c("Ca","K")]
xc0[,c("Ca","K")]
xc0$Ca
xc0[,"Ca"]
xp0[,"Ca"]
xc0[,c("Ca","K"), drop=TRUE]
xp0[,"Ca", drop=TRUE]

You can check which is the current status of sticky classes with getStickyClassOption().

Compositions as columns

Matrix-like compositions can be embedded in other data containers, typically "data.frame"s or tibbles. For this, we just need to define an extra column

Hydrochem$comp <- xc
head(Hydrochem)
dim(Hydrochem)

This composition can be then recovered with the dollar notation, and also used in further calls to statistical methods, e.g.

head(Hydrochem$comp)
codalm = lm(alr(comp)~River, data=Hydrochem) 
summary(codalm)

The transformation can even be called with the dollar notation, and transformed data keep the information necessary to construct their back-transformation.

codalm = lm(comp$idt~River, data=Hydrochem) 
head(predict(codalm)$idtInv)

To recover the column names in the backtransformation you might need to make use of the orig argument of the explicit transformation (that is, without the dollar notation)

head(idtInv(predict(codalm), orig=xc))

Note as well the existence of the new function backtransform, which will take an "rmult" object and return its expression in the appropriate original components. This can even be allowed to be controlled by a third object through an extra argument as

head(backtransform(idt(xp)))
head(backtransform(predict(codalm), as=codalm$residuals))

In the case of tibbles, the column names of the embedded object are even preserved, so that

tbHydro = tibble::as_tibble(Hydrochem)
tbHydro$comp = xc
tbHydro$comp$Ca

would return the same as xc$Ca. This is not included here to redude the dependencies of this package. Note that this trick does not work with "data.frame" objects, because [<- takes care of removing the column names of the composition before attaching it.

This property of self-back-transformability is the single most important addition to "compositions" v2, and will experiment improvements in the future. In particular, we strive to make the recovery of column names automatic in the near future.

Embedding in S4 objects

The final addition to compositions v2 is the ability of compositional data to fake to be "data.frame"s or "structure"s themselves. This allows them to be embedded in S3 elements and S4 classes in slots expecting one of these objects, keeping their additional attributes (hence their compositional nature, backtransformability and so on). This is an example of how this could work, for instance in a spatial data frame of package "sp". As before, this is not actually included here to avoid including one extra package. Readers are invited to check how it works. Packages "gstat" and "sp" are necessary

data("jura", package="gstat")
coords = jura.pred[, 1:2]
compo = jura.pred[, 7:13]
xc = cdt(acomp(compo))
spcompo = SpatialPointsDataFrame(coords=coords, data=xc)
summary(spcompo@data)
class(spcompo@data)

For an S4 class to admit compositional classes in their slots for "data.frame", "vector", "matrix" or "array" objects, though, there is a condition: validity checks must be made after as(x,"data.frame"), as(x,"matrix") or equivalent is applied.



Try the compositions package in your browser

Any scripts or data that you put into this service are public.

compositions documentation built on June 22, 2024, 12:15 p.m.