From the simplex to the plane and back

Consider points of the form \((x, y, z)\) where \(x, y, z \geq 0\).

Here are some points

# Generate some points
rbind(
  expand.grid(x=1:10,    y=1:10,    z=1:10)    %>% mutate(set="10x10x10", ID=NA),
  expand.grid(x=c(1,10), y=c(1,10), z=c(1,10)) %>% mutate(set="extremes", ID=LETTERS[row_number()])
) %>% 
  as_tibble() %>%
  mutate(
    # Generate the length from the origin
    length = sqrt(x^2+y^2+z^2),
    
    # Generate their closure
    sum=x+y+z,
    x.clo=x/sum,
    y.clo=y/sum,
    z.clo=z/sum,
    
    # Generate their clr
    x.log=log10(x),
    y.log=log10(y),
    z.log=log10(z),
    gm=(x.log + y.log + z.log)/3,
    x.clr=x.log-gm,
    y.clr=y.log-gm,
    z.clr=z.log-gm,
    
    # Generate a projection of the 3D simplex onto the 2D plane
    # See https://mathworld.wolfram.com/TernaryDiagram.html
    # I've centred this so that (1,1,1) maps to (0,0,0)
    t1 = (2*y + z)/(2*sum) - 0.5,
    t2 = (sqrt(3) * z)/(2*sum) - sqrt(3)/6,
    
    # Generate a projection of the 3D simplex onto ilr coordinates
    # This is the coordinates from Eq 5.1 of Compositional Data Analysis, Theory and Practice
    ilr1=-sqrt(1/2)*log10(x/y),
    ilr2=-sqrt(1/6)*log10(x*y/z^2),
    
    # Do the inverse ilr transform
    x.inv = 1,
    y.inv = 10^(ilr1*sqrt(2)),
    z.inv = 10^(ilr1/sqrt(2) + ilr2*sqrt(3/2)),
    sum.inv= x.inv + y.inv + z.inv,
    x.inv = x.inv/sum.inv,
    y.inv = y.inv/sum.inv,
    z.inv = z.inv/sum.inv
  ) -> X

# Create these handy subsets
Points   <- subset(X, set=="10x10x10")
Extremes <- subset(X, set=="extremes")

select(Points, x.clo, y.clo, z.clo, ilr1, ilr2, x.inv, y.inv, z.inv)  

## # A tibble: 1,000 x 8
##    x.clo  y.clo  z.clo   ilr1   ilr2 x.inv  y.inv  z.inv
##    <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>  <dbl>
##  1 0.333 0.333  0.333   0      0     0.333 0.333  0.333 
##  2 0.5   0.25   0.25   -0.213 -0.123 0.5   0.250  0.25  
##  3 0.6   0.2    0.2    -0.337 -0.195 0.6   0.2    0.2   
##  4 0.667 0.167  0.167  -0.426 -0.246 0.667 0.167  0.167 
##  5 0.714 0.143  0.143  -0.494 -0.285 0.714 0.143  0.143 
##  6 0.75  0.125  0.125  -0.550 -0.318 0.75  0.125  0.125 
##  7 0.778 0.111  0.111  -0.598 -0.345 0.778 0.111  0.111 
##  8 0.8   0.1    0.1    -0.639 -0.369 0.8   0.100  0.1   
##  9 0.818 0.0909 0.0909 -0.675 -0.390 0.818 0.0909 0.0909
## 10 0.833 0.0833 0.0833 -0.707 -0.408 0.833 0.0833 0.0833
## # ... with 990 more rows

Here are the points in 3D

marker.RdBu       <- list(color=~length, colorscale="RdBu", size=2, showscale=FALSE)
marker.extremes   <- list(color="red", size=3, showscale=FALSE, opacity=1)
text.red          <- list(color="red",   size = 20)

plot_ly(
  data=Points, 
  type="scatter3d", 
  name="points",
  mode="markers",
  x=~x, y=~y, z=~z,
  marker=marker.RdBu
) %>%
  add_trace(
    data=Extremes,                             
    name="label", 
    mode="markers+text",  
    x=~x, y=~y, z=~z, text=~ID, 
    marker=marker.extremes,  
    textfont = text.red  
  )

Points are transformed to the simplex \(\mathcal{S}^3\) by dividing them by their sum \(x+y+z\). This operation is called closure in compositional data analysis. Here are the points on the simplex:

plot_ly(
  data=Points, 
  type="scatter3d", 
  name="points",
  mode="markers",
  x=~x.clo, y=~y.clo, z=~z.clo,
  marker=marker.RdBu
) %>%
  add_trace(
    data=Extremes,                             
    name="label", 
    mode="markers+text",  
    x=~x.clo, y=~y.clo, z=~z.clo, text=~ID, 
    marker=marker.extremes,  
    textfont = text.red  
  )

Here is a projection of the simplex onto the 2D plane (see https://mathworld.wolfram.com/TernaryDiagram.html):

ggplot(
  data=Points,  aes(x=t1, y=t2)) + 
  geom_point(aes(color=length)) + 
  scale_color_gradient2(
    low = ("blue"),
    midpoint=median(Points$length),
    high = ("red")) +
  coord_equal() +
  geom_point(data=Extremes, color="red") + 
  geom_text(
    data=Extremes,
    aes(x=1.1*t1, y=1.1*t2, label=ID), color="red"
    ) +
  xlab(expression(ternary[1])) + ylab(expression(ternary[2])) + 
  theme(legend.position="none")

Here is their centred logratio transform

plot_ly(
  data=Points, 
  type="scatter3d", 
  name="points",
  mode="markers",
    x=~x.clr, y=~y.clr, z=~z.clr,
  marker=marker.RdBu
) %>%
  add_trace(
    data=Extremes,                             
    name="label", 
    mode="markers+text",  
      x=~x.clr, y=~y.clr, z=~z.clr, text=~ID, 
    marker=marker.extremes,  
    textfont = text.red  
  )

Here is (an) isometric logratio transform (there are infinitely many) onto the 2D plane:

ggplot(
  data=Points,  aes(x=ilr1, y=ilr2)) + 
  geom_point(aes(color=length)) + 
  scale_color_gradient2(
    low = ("blue"),
    midpoint=median(Points$length),
    high = ("red")) +
  coord_equal() +
  geom_point(data=Extremes, color="red") + 
  geom_text(
    data=Extremes,
    aes(x=1.1*ilr1, y=1.1*ilr2, label=ID), color="red"
    ) +
  xlab(expression(ilr[1])) + ylab(expression(ilr[2])) + 
  theme(legend.position="none")

The isometric logratio transform, and its inverse

There are infinitely many isometric logratio transformations that map the simplex \(\mathcal{S}^D\) to \(\mathbb{R}^{D-1}\). Here’s one from (Pawlowsky-Glahn and Buccianti 2011, 35) which maps a point \(\mathbf{x} = (x_1, x_2, \dots, x_D)\) on the \(D\)-dimensional simplex onto a point \(\mathbf{y} = (y_1, y_2, \dots, y_{D-1})\) on the \(D-1\)-dimensional Euclidean plane: \[ y_i = \frac{1}{\sqrt{i(i+1)}} \ln \left[\frac{x_1 x_2 \cdots x_i}{(x_{i+1})^i} \right], \qquad i = 1, 2, \dots D - 1. \] Let’s see what that looks like for a few values of \(i\): \[ \begin{align} y_1 &= \frac{1}{\sqrt{2}} \ln\frac{x_1}{x_2} \\ &= \frac{1}{\sqrt{2}} (\ln x_1 - \ln x_2)\\ y_2 &= \frac{1}{\sqrt{6}} \ln\frac{x_1 x_2}{x_3^2} \\ &= \frac{1}{\sqrt{2}} (\ln x_1 + \ln x_2 - 2\ln x_3) \end{align} \]

We can use this system of equations to show that \[ \begin{align} \ln x_2 &= \ln x_1 - \sqrt{2} y_1\\ \ln x_3 &= \ln x_1 - \frac{\sqrt{2}}{2} y_1 - \frac{\sqrt{6}}{2} y_2\\ \end{align} \] so that we can write the inverse transform of \((y_1, y_2)\) from the \(D-1\)-dimensional Euclidean plane back to the \(D\)-dimensional simplex as \[ \mathcal{C}(\exp(1, -\sqrt{2}y_i, -\frac{\sqrt{2}}{2} y_1 - \frac{\sqrt{6}}{2} y_2)) \] where \(\mathcal{C}\) is the closure operator.

Example

Here we map an ellipse from an isometric logratio plane in \(\mathbb{R}^{2}\) to the probability simplex in \(\mathcal{S}^D\)

First, we make a function to generate \(n\) points on an ellipse centred on \((x_c, y_c)\) with major axis length \(a\), minor axis length \(b\) and angle to the \(x\)-axis \(\phi\):

# Adapted from https://stackoverflow.com/questions/41820683/how-to-plot-ellipse-given-a-general-equation-in-r
ellipse <- function(xc=0.1, yc=-0.2, a=1, b=0.5, phi=pi/3, n=100){
  t <- seq(0, 2*pi, 0.01)
  tibble(
  x = xc + a*cos(t)*cos(phi) - b*sin(t)*sin(phi),
  y = yc + a*cos(t)*cos(phi) + b*sin(t)*cos(phi)
  )
}

Then we treat the points of that ellipse as though they are ilr coordinates \((\mathrm{ilr}_1, \mathrm{ilr}_2)\), then apply an inverse ilr mapping to \((x_\mathrm{inv}, y_\mathrm{inv}, z_\mathrm{inv})\), and project that onto a 2-D ternary diagram with points \((t_1, t_2)\)

ellipse() %>%
  rename(ilr1=x, ilr2=y) %>%
  mutate(
    x.inv = 1,
    y.inv = 10^(ilr1*sqrt(2)),
    z.inv = 10^(ilr1/sqrt(2) + ilr2*sqrt(3/2)),
    sum.inv= x.inv + y.inv + z.inv,
    x.inv = x.inv/sum.inv,
    y.inv = y.inv/sum.inv,
    z.inv = z.inv/sum.inv,
    t1 = (2*y.inv + z.inv)/2 - 0.5,
    t2 = (sqrt(3) * z.inv)/2 - sqrt(3)/6
  ) -> ilr.ellipse

Here is the view in ilr coordinates, with the points from our cube of counts for reference:

ggplot(
  data=Points,  aes(x=ilr1, y=ilr2)) + 
  geom_point(aes(color=length)) + 
  scale_color_gradient2(
    low = ("blue"),
    midpoint=median(Points$length),
    high = ("red")) +
  coord_equal() +
  geom_point(data=Extremes, color="red") + 
  geom_text(
    data=Extremes,
    aes(x=1.1*ilr1, y=1.1*ilr2, label=ID), color="red"
    ) +
  geom_path(data=ilr.ellipse) +
  xlab(expression(ilr[1])) + ylab(expression(ilr[2])) + 
  theme(legend.position="none")

And here is the view from the simplex

tribble(
  ~x.inv, ~y.inv, ~z.inv, ~ID, ~angle,
  1, 0, 0,     "x = 1", 30,
  1/2, 1/2, 0, "(1/2, 1/2, 0)", 0,
  0, 1, 0,     "y = 1", -30,
  0, 1/2, 1/2, "(0, 1/2, 1/2)", -60,
  0, 0, 1,     "z = 1", 90,
  1/2, 0, 1/2, "(1/2, 0, 1/2)", 60,
  1, 0, 0, NA, 0
) %>%
  mutate(
    t1 = (2*y.inv + z.inv)/2 - 0.5,
    t2 = (sqrt(3) * z.inv)/2 - sqrt(3)/6
  ) -> S3

ggplot(
  data=Points,  aes(x=t1, y=t2)) + 
  # The simplex boundary
  geom_path(data=S3, color="grey") +
  geom_text(
    data=filter(S3, !is.na(ID)), color="grey",size=3, hjust=0.5, vjust=0.5,
    aes(x=1.1*t1, y=1.1*t2, label=ID, angle=angle)
  ) +
  # The reference points
  geom_point(aes(color=length)) + 
  scale_color_gradient2(
    low = ("blue"),
    midpoint=median(Points$length),
    high = ("red")) +
  # The extremes
  geom_point(data=Extremes, color="red") + 
  geom_text(
    data=Extremes,
    aes(x=1.1*t1, y=1.1*t2, label=ID), color="red"
  ) +
  # The transformed ilr ellipse
  geom_path(data=ilr.ellipse) +
  coord_equal() +
  xlab(expression(t[1])) + ylab(expression(t[2])) + 
  theme(legend.position="none")

References

Pawlowsky-Glahn, Vera, and Antonella Buccianti, eds. 2011. Compositional Data Analysis: Theory and Applications. http://onlinelibrary.wiley.com/book/10.1002/9781119976462.