Feature/spca biplot

.Rbuildignore

@@ -10,4 +10,5 @@ extra/*
 ^docs$
 ^pkgdown$
 ^CRAN-SUBMISSION$
-revdep
+revdep
+^\.gitattributes$

.gitattributes

@@ -0,0 +1,17 @@
+# Normalise line endings: LF in the repository, native on disk.
+# Without this, files edited on Windows show up as entirely rewritten and
+# diffs become unreadable.
+* text=auto eol=lf
+
+# Files git must not treat as text
+*.png binary
+*.jpg binary
+*.jpeg binary
+*.pdf binary
+*.rda binary
+*.RData binary
+*.rds binary
+*.gpkg binary
+*.shp binary
+*.shx binary
+*.dbf binary

DESCRIPTION

@@ -21,7 +21,6 @@ Description: Precision agriculture spatial data
 License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.3.3
 Imports: 
     data.table,
     e1071,
@@ -43,3 +42,4 @@ Suggests:
 URL: https://ppaccioretti.github.io/paar/, https://github.com/PPaccioretti/paar
 VignetteBuilder: knitr, rmarkdown
 BugReports: https://github.com/PPaccioretti/paar/issues
+Config/roxygen2/version: 8.0.0

NAMESPACE

@@ -1,11 +1,14 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(cbind,paar)
+S3method(print,multispati)
 S3method(print,paar)
 S3method(print,summary.paar)
+S3method(summary,multispati)
 S3method(summary,paar)
 export(compare_zone)
 export(depurate)
 export(fuzzy_k_means)
 export(kmspc)
 export(spatial_t_test)
+export(spca_biplot)

NEWS.md

@@ -1,5 +1,29 @@
 # paar (development version)
 
+* Fixed: in `spca_summary_results`, the `Prop` column of `kmspc()` was a
+  proportion between 0 and 1 while `Acum. Prop.` in the same table was a
+  percentage between 0 and 100. `spca_biplot()` formats `Prop` with a per-cent
+  sign, so the axis labels reported one hundredth of the explained variability,
+  for example `CS1 (0.5%)` where the axis actually explains 50.6%. Both columns
+  are now percentages. The denominator is the total inertia of the analysis
+  rather than the number of retained eigenvalues, which differ when the data
+  matrix is rank deficient. The selection of components from
+  `explainedVariance` is unaffected, as it already operated on the cumulative
+  percentage. The `Prop` and `Acum. Prop.` columns of the classical PCA summary
+  returned by `only_spca_results = FALSE` were also proportions and are now
+  percentages, so that both tables can be read side by side.
+* Fixed: `kmspc()` computed the spatial PCA results but discarded them unless
+  `all_results = TRUE`. With the default arguments `pca_results` came back
+  empty, even though `only_spca_results = TRUE` is meant to return precisely
+  those results. The spatial PCA summary and the loadings are now always
+  returned.
+* `kmspc()` also returns `scores`, the site coordinates on the retained
+  spatial components. These are the values the clustering step operates on and
+  were previously computed but not exposed.
+* New `spca_biplot()`, which draws a biplot of the spatial PCA from a
+  `kmspc()` result, optionally colouring the sites by one of the computed
+  classifications.
+
 # paar 1.0.2
 
 * Help pages were revised for clarity, consistency, and improved user guidance [#21](https://github.com/PPaccioretti/paar/issues/21)

R/kmspc.R

@@ -67,9 +67,21 @@
 #'   \item{cluster}{\code{data.frame} with cluster assignments for each evaluated number of clusters}
 #'   \item{indices}{\code{data.frame} with clustering validity indices}
 #'   \item{summaryResults}{\code{data.frame} with clustering metrics (iterations, SSDW)}
-#'   \item{pca_results}{(optional) PCA and/or spatial PCA summaries depending on arguments}
+#'   \item{pca_results}{a list with the spatial PCA results, always present:
+#'     \code{spca_summary_results} (eigenvalues, spatial variance, explained
+#'     percentage of the total inertia, cumulative percentage and Moran index
+#'     per axis), \code{eigenvectors_used}
+#'     (variable loadings on the retained axes) and \code{scores} (site
+#'     coordinates on the retained axes, the input of the clustering step).
+#'     If \code{only_spca_results = FALSE} it also contains
+#'     \code{pca_results} with the classical PCA summary, and if
+#'     \code{all_results = TRUE}, the complete \code{pca_results_all} and
+#'     \code{spca_results_all} objects.}
 #' }
 #'
+#' @seealso [spca_biplot()] to display `eigenvectors_used` and `scores`
+#'   together in a biplot.
+#'
 #' @export
 #' @example inst/examples/kmspc.R
 #'
@@ -159,7 +171,9 @@ kmspc <- function(
 
   if (!only_spca_results) {
     autov_pca <- pca$eig
-    propvar_pca <- autov_pca / sum(autov_pca)
+    # Percentages, on the same scale as the spatial PCA table below, so that
+    # both summaries can be read side by side.
+    propvar_pca <- 100 * autov_pca / sum(autov_pca)
     propvaracum_pca <- cumsum(propvar_pca)
 
     my_pca <- lapply(pca$l1, spdep::moran.mc, lw, 999)
@@ -168,15 +182,15 @@ kmspc <- function(
 
     nfila_pca <- length(variables)
     eje_pca <- seq_len(nfila_pca)
-    resultado_pca = data.frame(
+    resultado_pca <- data.frame(
       eje_pca,
       autov_pca,
       propvar_pca,
       propvaracum_pca,
       my_pca
     )
     resultado_pca$eje_pca <- as.factor(resultado_pca$eje_pca)
-    names(resultado_pca) = c(
+    names(resultado_pca) <- c(
       "Axis",
       "Eigenvalue",
       "Prop",
@@ -190,8 +204,14 @@ kmspc <- function(
   invisible(utils::capture.output(resms <- summary(ms)))
   var_ms <- resms[, 2, drop = F]
   nfila_ms <- length(ms$eig)
-  propvar_ms <- var_ms / nfila_ms
-  propvaracum_ms <- cumsum(propvar_ms) * 100
+
+  # Both `Prop` and `Acum. Prop.` are expressed as percentages, on the same
+  # scale as the `explainedVariance` threshold. The denominator is the total
+  # inertia of the analysis, which equals the number of variables only when
+  # the data matrix has full rank.
+  total_inertia <- sum(pca$eig)
+  propvar_ms <- 100 * var_ms / total_inertia
+  propvaracum_ms <- cumsum(propvar_ms)
 
   eje_ms <- seq_len(nfila_ms)
   resultado_ms <-
@@ -221,22 +241,26 @@ kmspc <- function(
       unlist(propvaracum_ms)
     ))
 
+  data_clust <- ms$li[num_sPC]
+
+  if (inherits(data_clust, "sf")) {
+    data_clust <- sf::st_drop_geometry(data_clust)
+  }
+
   spca_results <- list(
     spca_summary_results = resultado_ms,
-    eigenvectors_used = ms$c1[, num_sPC, drop = FALSE]
+    eigenvectors_used = ms$c1[, num_sPC, drop = FALSE],
+    scores = data_clust
   )
 
   if (all_results) {
     spca_results <- append(spca_results, list(spca_results_all = ms))
-
-    pca_results <- append(pca_results, spca_results)
   }
 
-  data_clust <- ms$li[num_sPC]
-
-  if (inherits(data_clust, "sf")) {
-    data_clust <- sf::st_drop_geometry(data_clust)
-  }
+  # The spatial PCA results are always returned: they are the output the
+  # `only_spca_results` argument refers to, and `scores` together with
+  # `eigenvectors_used` are what `spca_biplot()` needs.
+  pca_results <- append(pca_results, spca_results)
 
   my_results <- make_clasification(
     data_clust,
@@ -355,7 +379,7 @@ summarize_indices <- function(indices, number_cluster) {
   }
 
   indicesresults <- data.frame(number_cluster, indices, IndN)
-  names(indicesresults) = c(
+  names(indicesresults) <- c(
     "Num. Cluster",
     "Xie Beni",
     # "Fukuyama Sugeno",
@@ -659,6 +683,7 @@ multispati <-
 
 
 #' @noRd
+#' @exportS3Method
 summary.multispati <- function(object, ...) {
   norm.w <- function(X, w) {
     f2 <- function(v) {
@@ -731,6 +756,7 @@ summary.multispati <- function(object, ...) {
 
 
 #' @noRd
+#' @exportS3Method
 print.multispati <- function(x, ...) {
   # cat("Multispati object \n")
   # cat("class: ")