projection (prediction) and vector sum (interpolation) geoms

[corybrunson]

A geom layer geom_u_projection(from = i, to = j) should render a (by default dashed) line from the ordinates of the ith row of $U$ to the linear subspace containing those of the jth row of $V$, maybe with an optional cute right angle symbol. Both from and to could contain multiple indices. The projection should adopt the axes (primary or secondary) used by $U$.

[jtr13]

First of all, thank you for this fantastic package (as well as for ggalluvial!). A dropped perpendicular feature would be a great addition: combined with the calibrated axes the perpendiculars really help explain the meaning of a biplot. I've been creating them with calibrate package in base R and have taken a stab at using calibrate calculations in ggplot2 but of course it would be nice to have a smoother interface.

library(calibrate)
#> Loading required package: MASS
    df <- data.frame(x = c(5, 8, 9), y = c(4, 13, 6))
    pca <- prcomp(df)
    scores <- pca$x
    loadings <- pca$rotation
    axisrange <- max(df[,"x"]) - min(df[,"x"])
    ticklab <- round(seq(min(df[,"x"])-.5*axisrange,
                         max(df[,"x"])+.5*axisrange,
                         by = axisrange/2))
    plot(scores, pch = 16, asp = 1, ylim = c(-5, 5))
    text(x = scores[,1]-.5, y = scores[,2], labels = 1:nrow(scores))
    c <- calibrate(g = loadings[,"PC1"],
                   y = df[,"x"] - mean(df[,"x"]),
                   tm = ticklab - mean(df[,"x"]),
                   Fr = scores,
                   tmlab = ticklab,
                   tl = .2,
                   lm = TRUE,
                   axislab="x",
                   where = 1,
                   labpos = 4,
                   dp = TRUE)
#> ---------- Calibration Results for  x  ---------------------
#> Length of 1 unit of the original variable =  1  
#> Angle                                     =  76.3 degrees
#> Optimal calibration factor                =  1  
#> Used calibration factor                   =  1  
#> Goodness-of-fit                           =  1  
#> Goodness-of-scale                         =  1  
#> ------------------------------------------------------------
    arrows(0, 0, loadings[,1], loadings[,2], length = .1, col = "red",
           lwd = 1.5)

^{Created on 2021-08-25 by the reprex package (v2.0.1)}

library(calibrate)
#> Loading required package: MASS
    library(ggplot2)
    df <- data.frame(x = c(5, 8, 9), y = c(4, 13, 6))
    pca <- prcomp(df)
    scores <- pca$x
    loadings <- pca$rotation
    axisrange <- max(df[,"x"]) - min(df[,"x"])
    ticklab <- round(seq(min(df[,"x"])-.5*axisrange,
                         max(df[,"x"])+.5*axisrange,
                         by = axisrange/2))
    c <- calibrate(g = loadings[,"PC1"],
                   y = df[,"x"] - mean(df[,"x"]),
                   tm = ticklab - mean(df[,"x"]),
                   Fr = scores,
                   tmlab = ticklab,
                   tl = .3,
                   graphics = FALSE)
#> ---------- Calibration Results for    ----------------------
#> Length of 1 unit of the original variable =  1  
#> Angle                                     =  76.3 degrees
#> Optimal calibration factor                =  1  
#> Used calibration factor                   =  1  
#> Goodness-of-fit                           =  1  
#> Goodness-of-scale                         =  1  
#> ------------------------------------------------------------
    
    
    dfpoints <- data.frame(scores)
    dfpoints$xsdrop <- c$Fpr[,1]
    dfpoints$ysdrop <- c$Fpr[,2]
    dfaxis <- data.frame(x = c$M[1, 1], y = c$M[1, 2],
                         xend = c$M[nrow(c$M),1], yend = c$M[nrow(c$M), 2])
    dfticks <- data.frame(c$M, c$Mn, ticklab)
    colnames(dfticks) <- c("x", "y", "xend", "yend", "label")
    dfarrows <- data.frame(xend = loadings[,1], yend = loadings[,2])
    
    
    ggplot(dfpoints, aes(x = PC1, y = PC2)) +
        geom_point() +
        geom_segment(data = dfticks, aes(x = x, y = y, xend = xend, yend = yend), col = "blue") +
        geom_text(data = dfticks, aes(x = xend, y = yend, label = label), nudge_x = .3) +
        geom_segment(data = dfaxis, aes(x = x, y = y, xend = xend, yend = yend), col = "blue") +
        geom_segment(data = dfpoints, aes(x = PC1, y = PC2, xend = xsdrop, yend = ysdrop),
                     lty = "dashed", col = "cornflowerblue") +
        geom_segment(data = dfarrows, aes(x = 0, y = 0, xend = xend, yend = yend),
                     arrow = arrow(length = unit(.03, "npc")), color = "red", lwd = 1.5) + coord_equal()

^{Created on 2021-08-25 by the reprex package (v2.0.1)}

The closest I can get with ordr:

library(ggplot2)
    library(ggbiplot)
#> Loading required package: plyr
#> Loading required package: scales
#> Loading required package: grid
    library(ordr)
#> 
#> Attaching package: 'ordr'
#> The following object is masked from 'package:ggbiplot':
#> 
#>     ggbiplot
    df <- data.frame(x = c(5, 8, 9), y = c(4, 13, 6))
    pca_ordr <- ordinate(df, cols = 1:2, model = ~ prcomp(., scale. = TRUE))
    ggbiplot(pca_ordr) +
        xlim(c(-2, 2)) +
        ylim(c(-2, 2)) +
        geom_rows_point() +
        geom_rows_text(aes(label = 1:5), nudge_x = .2) +
        geom_cols_vector(color = "red", lwd = 1.5) + 
        geom_cols_axis() +
        geom_cols_axis_ticks(aes(center = .center, scale = .scale)) +
        geom_cols_axis_text(aes(center = .center, scale = .scale)) +
        geom_cols_axis_label(aes(label = .name))
#> Warning: Ignoring unknown aesthetics: center, scale

#> Warning: Ignoring unknown aesthetics: center, scale
#> Warning: Ignoring unknown aesthetics: label

^{Created on 2021-08-25 by the reprex package (v2.0.1)}

I should also note that it would be helpful to be able to plot one axis at a time since the graph gets messy with multiple axes, though that might be challenge in the ggplot2 framework.

[corybrunson]

@jtr13 thank you for these details!

I'm going to reassess what features should go into a first release, and you've got me thinking that the projection geom is one that should be put back on that list. (Meanwhile, the predict = TRUE option for ggbiplot() is perhaps too ambitious, and i'm planning to remove it for the first release.) I'll come back to work on it ASAP.

Ditto the ability to specify what row or column elements should be plotted by each layer, as you mention at the end. An experimental solution is the subset parameter, which is implemented in GeomIsolines$setup_data and recycled by a few other layers, not including GeomAxis. I'll add an issue to expand this to all plot layers.

[corybrunson]

Correction and clarification: Thanks to Gower &al's books, i should have remembered that projection onto an axis is a prediction step that calls for a prediction biplot, whereas vector addition along the axes is an interpolation step that makes use of the (currently implemented) interpolation biplot. (Interpolation axis unit vectors are located analogously to case markers, whereas prediction axes must be rescaled in the linear setting and completely redrawn in the non-linear setting.) This issue should really be about rendering both sets of steps, each being appropriate for one type of biplot. Vector addition should come first, for the first CRAN release.

[corybrunson]

In the conventional use of PCA with rows of cases and columns of variables, the vector sum geom could look like geom_cols_addition(rows = new_data). The layer would inherit the column (standard) coordinates from ggbiplot(), and the rows (or perhaps new_data) parameter would accept rows of data with the same variables as the original table fed to, e.g., princomp(), or a list of variable values, e.g. list(mpg = 28, cyl = 4L, ...). If some values are missing or NA but the original PCA is centered, then the interpolated points would correspond to mean imputations.

It shouldn't be much extra work to include an option type = c("centroid", "sequence") to control whether the vectors along the axes are added in sequence or their centroid is scaled by their number.

[corybrunson]

Experimental work is ongoing in the vecsum branch.

[corybrunson]

Update: Based on further reading, and on the counterintuitiveness and complexity of the experimental implementation, i think both of these annotations should work as follows:

1. Make generics and S3 class methods to calculate layer data (like the output of Stat*$compute_*()) for predictive layers, interpolation vector sums, prediction projections, and other geometric elements derived from matrix algebra.
2. Make a pronoun (.O()?) for the underlying ordination model, analogous to .G() in tidygraph, by which to access these methods from within a ggproto build.
3. Make new statistical transformations StatRows* and StatCols* that use this pronoun to calculate layer data.

The reason for a change in approach is that the 'tbl_ord' class is designed for display (print, summary, annotation, plot), and ggbiplot() itself can only see the fortified tibble, whereas the underlying ordination classes, which are designed for modeling (interpolate, predict), are the proper setting for these tasks. For example, predictive MDS biplots use nonlinear axes while predictive LDA biplots use Voronoi cells, and the S3 methods can be tailored to these needs.

I think this deserves to wait until feedback is collected from a first CRAN release.

A compensatory benefit will be the new generics and methods themselves, which might be used in base R biplots or for other purposes. To my knowledge, there are no standalone implementations of them; they are only built inaccessibly into visualization tools like those of Gower &al.

[corybrunson]

The calibrate branch contains a partial solution that requires integer specification of axes to offset and manual expansion of the plotting window if necessary. See the last example in ?biplot after installing from that branch.

I think doing this properly—specifically, determining the offset and the extent of each axis from the range of plotted elements from the other matrix factor rather than by trial and error—will come with the registration/pronoun solution or some other trick to access the model from the plot build.

[jtr13]

I will check it out!

[corybrunson]

@jtr13 i'm drafting what i hope will be a legit solution in the offset branch. Lots of work remains to be done, but the core transformations seem to work! The functionality is illustrated in the same example.

[corybrunson]

As of 13cad83, the offset branch is feature-complete and mostly debugged. Documentation, examples, and tests remain to be written, but it's ready for stress-testing!

Reminder to self: This work toward delimited and offset calibrated axes was only a prelude to projection and vector sum graphical elements, so this issue will remain open until those are implemented.

[corybrunson]

This issue might be resolved alongside #64. Solutions to both are being drafted in offset and its spinoff branch offset-referent, using a ggplot_add() method to inherit the global aesthetic mapping as suggested in this comment.

A standalone (non-'tbl_ord') convenience parent StatReferent and one child StatProjection are drafted, which expect data passed to referent to have the same position columns (i.e. those mapped to x and y) as the inherited data. The tasks remaining are

1. to default referent to the cols data in stat_rows_projection() shortcuts and vice-versa;
2. to allow referent to accept a subset-like specification; and
3. to adapt build-layers.r to properly construct the stat_*_projection() shortcuts.

Once these tasks are done, the infrastructure will be ready to support StatInterpolation (rather than GeomAddition as in the vecsum branch) and a rewrite of StatRule.

[corybrunson]

Tentatively resolved in the offset branch at 6af4a2d.

Note: Interpolation should not require the referent trick, so the vecsum branch might be merged with offset to finally close this issue.

[corybrunson]

GeomAddition needs some reconciliation and probably a renaming, but the merge is underway in the offset-vecsum branch.

[corybrunson]

@jtr13 i believe this issue is resolved with today's merge. As you have the bandwidth, i'd be glad to know if the following new layers and shortcuts suit your needs—or, if not, what remains to be done.

• geom_interpolation()*
• stat_projection()
• stat_rule() and geom_rule()

Thanks a lot for your continued input!

*Should perhaps be reverted to geom_addition() for consistency with stat_projection(), since "interpolation" is better paired with "prediction".

[jtr13]

Sounds good! I can't wait to finish end-of-semester grading and other admin stuff so I can take a look at this!

[jtr13]

The new additions are great! I've experimented with stat_projection() with very positive results:

library(Lock5withR)
library(ordr)
pca_ordr <- ordinate(Cereal, cols = 3:10,
                     model = ~ prcomp(., scale. = TRUE))

ggbiplot(pca_ordr) +
  geom_rows_point(size = 1, color = "cornflowerblue") +
  geom_rows_text_repel(aes(label = Name), size = 2, color = "cornflowerblue") +
  geom_cols_vector(color = "red", vector_labels = FALSE) +
  geom_cols_axis(aes(label = name), label_size = 2.5, label_color = "red", num = 10) +
  stat_rows_projection(subset = 23, ref_subset = c(2, 4, 6), color = "blue", linetype = "dashed") +
  xlim(c(-5, 5)) +
  theme_minimal()

Questions:

• Is geom_interpolation() meant for predicting with new data?

• referent isn't explained in ?stat_projection() -- is there any reason this might be changed?

• Curious why it's not geom_rows_projection() -- I know there's a larger discussion about stats vs geoms...

A few really minor comments:

• It took me a while to figure out that the row annotations are stored as Name, it's also a little confusing to use name for columns and Name for rows

• I like how geom_cols_vector() puts the vector labels on the positive side; it would be nice if geom_cols_axis() did the same.

The more I explore the package the more I realize how much there is that I haven't tried! It is an impressive project.