Add log-ratio transforms (#213)

* Add log-ratio transforms

* Update tests

* Fix revert

* Fix typo

Author: eliascarv

Project.toml

@@ -6,6 +6,7 @@ version = "1.16.0"
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
+CoDa = "5900dafe-f573-5c72-b367-76665857777b"
 ColumnSelectors = "9cc86067-7e36-4c61-b350-1ac9833d277f"
 DataScienceTraits = "6cb2f572-2d2b-4ba6-bdb3-e710fa044d6c"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -23,6 +24,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 AbstractTrees = "0.4"
 CategoricalArrays = "0.10"
+CoDa = "1.2"
 ColumnSelectors = "0.1"
 DataScienceTraits = "0.1"
 Distributions = "0.25"

docs/src/transforms.md

@@ -200,6 +200,24 @@ Closure
 Remainder
 ```
 
+## ALR
+
+```@docs
+ALR
+```
+
+## CLR
+
+```@docs
+CLR
+```
+
+## ILR
+
+```@docs
+ILR
+```
+
 ## RowTable
 
 ```@docs

src/TableTransforms.jl

@@ -13,6 +13,7 @@ using LinearAlgebra
 using DataScienceTraits
 using CategoricalArrays
 using Random
+using CoDa
 
 using TransformsBase: Transform, Identity, →
 using DataScienceTraits: SciType, Continuous, coerce
@@ -80,6 +81,9 @@ export
   ProjectionPursuit,
   Closure,
   Remainder,
+  ALR,
+  CLR,
+  ILR,
   RowTable,
   ColTable,
   →,

src/transforms.jl

@@ -290,6 +290,7 @@ include("transforms/eigenanalysis.jl")
 include("transforms/projectionpursuit.jl")
 include("transforms/closure.jl")
 include("transforms/remainder.jl")
+include("transforms/logratio.jl")
 include("transforms/rowtable.jl")
 include("transforms/coltable.jl")
 include("transforms/parallel.jl")

src/transforms/logratio.jl

@@ -0,0 +1,93 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    LogRatio
+
+Parent type of all log-ratio transforms.
+
+See also [`ALR`](@ref), [`CLR`](@ref), [`ILR`](@ref).
+"""
+abstract type LogRatio <: StatelessFeatureTransform end
+
+# log-ratio transform interface
+function refvar end
+function newvars end
+function oldvars end
+function applymatrix end
+function revertmatrix end
+
+isrevertible(::Type{<:LogRatio}) = true
+
+assertions(::LogRatio) = [SciTypeAssertion{Continuous}()]
+
+function applyfeat(transform::LogRatio, feat, prep)
+  cols = Tables.columns(feat)
+  names = Tables.columnnames(cols) |> collect
+
+  # reference variable
+  rvar = refvar(transform, names)
+  @assert rvar ∈ names "invalid reference variable"
+  rind = findfirst(==(rvar), names)
+
+  # permute columns if necessary
+  perm = rind ≠ lastindex(names)
+  pfeat = if perm
+    popat!(names, rind)
+    push!(names, rvar)
+    feat |> Select(names)
+  else
+    feat
+  end
+
+  # apply transform
+  X = Tables.matrix(pfeat)
+  Y = applymatrix(transform, X)
+
+  # new variable names
+  newnames = newvars(transform, names)
+
+  # return same table type
+  𝒯 = (; zip(newnames, eachcol(Y))...)
+  newfeat = 𝒯 |> Tables.materializer(feat)
+
+  newfeat, (rvar, rind, perm)
+end
+
+function revertfeat(transform::LogRatio, newfeat, fcache)
+  cols = Tables.columns(newfeat)
+  names = Tables.columnnames(cols)
+
+  # revert transform
+  Y = Tables.matrix(newfeat)
+  X = revertmatrix(transform, Y)
+
+  # retrieve cache
+  rvar, rind, perm = fcache
+
+  # original variable names
+  onames = oldvars(transform, names, rvar)
+
+  # revert the permutation if necessary
+  if perm
+    n = length(onames)
+    inds = collect(1:(n - 1))
+    insert!(inds, rind, n)
+    onames = onames[inds]
+    X = X[:, inds]
+  end
+
+  # return same table type
+  𝒯 = (; zip(onames, eachcol(X))...)
+  𝒯 |> Tables.materializer(newfeat)
+end
+
+
+# ----------------
+# IMPLEMENTATIONS
+# ----------------
+
+include("logratio/alr.jl")
+include("logratio/clr.jl")
+include("logratio/ilr.jl")

src/transforms/logratio/alr.jl

@@ -0,0 +1,27 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    ALR([refvar])
+
+Additive log-ratio transform.
+
+Optionally, specify the reference variable `refvar` for the ratios.
+Default to the last column of the input table.
+"""
+struct ALR{T<:Union{Symbol,Nothing}} <: LogRatio
+  refvar::T
+end
+
+ALR() = ALR(nothing)
+
+refvar(transform::ALR, names) = isnothing(transform.refvar) ? last(names) : transform.refvar
+
+newvars(::ALR, names) = collect(names)[begin:(end - 1)]
+
+oldvars(::ALR, names, rvar) = [collect(names); rvar]
+
+applymatrix(::ALR, X) = mapslices(alr ∘ Composition, X, dims=2)
+
+revertmatrix(::ALR, Y) = mapslices(components ∘ alrinv, Y, dims=2)

src/transforms/logratio/clr.jl

@@ -0,0 +1,20 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    CLR()
+
+Centered log-ratio transform.
+"""
+struct CLR <: LogRatio end
+
+refvar(::CLR, names) = last(names)
+
+newvars(::CLR, names) = collect(names)
+
+oldvars(::CLR, names, rvar) = collect(names)
+
+applymatrix(::CLR, X) = mapslices(clr ∘ Composition, X, dims=2)
+
+revertmatrix(::CLR, Y) = mapslices(components ∘ clrinv, Y, dims=2)

src/transforms/logratio/ilr.jl

@@ -0,0 +1,27 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    ILR([refvar])
+
+Isometric log-ratio transform.
+
+Optionally, specify the reference variable `refvar` for the ratios.
+Default to the last column of the input table.
+"""
+struct ILR{T<:Union{Symbol,Nothing}} <: LogRatio
+  refvar::T
+end
+
+ILR() = ILR(nothing)
+
+refvar(transform::ILR, names) = isnothing(transform.refvar) ? last(names) : transform.refvar
+
+newvars(::ILR, names) = collect(names)[begin:(end - 1)]
+
+oldvars(::ILR, names, rvar) = [collect(names); rvar]
+
+applymatrix(::ILR, X) = mapslices(ilr ∘ Composition, X, dims=2)
+
+revertmatrix(::ILR, Y) = mapslices(components ∘ ilrinv, Y, dims=2)

test/transforms.jl

@@ -26,6 +26,7 @@ transformfiles = [
   "projectionpursuit.jl",
   "closure.jl",
   "remainder.jl",
+  "logratio.jl",
   "rowtable.jl",
   "coltable.jl",
   "sequential.jl",

test/transforms/logratio.jl

@@ -0,0 +1,41 @@
+@testset "LogRatio" begin
+  @test isrevertible(ALR())
+  @test isrevertible(CLR())
+  @test isrevertible(ILR())
+
+  a = [1.0, 0.0, 1.0]
+  b = [2.0, 2.0, 2.0]
+  c = [3.0, 3.0, 0.0]
+  t = Table(; a, b, c)
+
+  T = ALR()
+  n, c = apply(T, t)
+  @test n == t |> ALR(:c)
+  talr = revert(T, n, c)
+  T = CLR()
+  n, c = apply(T, t)
+  tclr = revert(T, n, c)
+  T = ILR()
+  n, c = apply(T, t)
+  @test n == t |> ILR(:c)
+  tilr = revert(T, n, c)
+  @test Tables.matrix(talr) ≈ Tables.matrix(tclr)
+  @test Tables.matrix(tclr) ≈ Tables.matrix(tilr)
+  @test Tables.matrix(talr) ≈ Tables.matrix(tilr)
+
+  # permute columns
+  a = [1.0, 0.0, 1.0]
+  b = [2.0, 2.0, 2.0]
+  c = [3.0, 3.0, 0.0]
+  t1 = Table(; a, c, b)
+  t2 = Table(; c, a, b)
+
+  T = ALR(:c)
+  n1, c1 = apply(T, t1)
+  n2, c2 = apply(T, t2)
+  @test n1 == n2
+  tₒ = revert(T, n1, c1)
+  @test Tables.schema(tₒ).names == (:a, :c, :b)
+  tₒ = revert(T, n2, c2)
+  @test Tables.schema(tₒ).names == (:c, :a, :b)
+end