From a88ce03d6507ad9e4aa62a25723a2b9e0c9841ac Mon Sep 17 00:00:00 2001 From: OuGuangyu Date: Tue, 3 Jun 2025 11:07:38 +0800 Subject: [PATCH 1/2] feature: gwpca interface --- NAMESPACE | 2 + R/RcppExports.R | 4 + R/gwpca.R | 165 ++++++++++++++++++++++++++++++++++++ man/gwpca.Rd | 54 ++++++++++++ man/print.Rd | 10 ++- src/Makevars.in | 2 + src/Makevars.win | 2 + src/RcppExports.cpp | 20 +++++ src/gwpca.cpp | 70 +++++++++++++++ tests/testthat/test-gwpca.R | 16 ++++ 10 files changed, 343 insertions(+), 2 deletions(-) create mode 100644 R/gwpca.R create mode 100644 man/gwpca.Rd create mode 100644 src/gwpca.cpp create mode 100644 tests/testthat/test-gwpca.R diff --git a/NAMESPACE b/NAMESPACE index 7e0b4f6..6e6a20d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -12,6 +12,7 @@ S3method(plot,gwrmultiscalem) S3method(plot,modelselcritl) S3method(predict,gwrm) S3method(print,gtdrm) +S3method(print,gwpcam) S3method(print,gwrm) S3method(print,gwrmultiscalem) S3method(residuals,gtdrm) @@ -22,6 +23,7 @@ S3method(step,gtdrm) S3method(step,gwrm) export(gtdr) export(gtdr_config) +export(gwpca) export(gwr_basic) export(gwr_multiscale) export(mgwr_config) diff --git a/R/RcppExports.R b/R/RcppExports.R index 6c1c586..25cc1a2 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -5,6 +5,10 @@ gtdr_fit <- function(x, y, coords, bw, adaptive, kernel, intercept, hatmatrix, p .Call(`_GWmodel3_gtdr_fit`, x, y, coords, bw, adaptive, kernel, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw, optim_bw_criterion, optim_threashold, optim_step, optim_max_iter, select_model, select_model_threshold, variable_names, verbose) } +gwpca_cal <- function(x, coords, bw, adaptive, kernel, longlat, p, theta, keep_components) { + .Call(`_GWmodel3_gwpca_cal`, x, coords, bw, adaptive, kernel, longlat, p, theta, keep_components) +} + gwr_basic_fit <- function(x, y, coords, bw, adaptive, kernel, longlat, p, theta, optim_bw_lower, optim_bw_upper, hatmatrix, intercept, parallel_type, parallel_arg, optim_bw, optim_bw_criterion, select_model, select_model_criterion, select_model_threshold, variable_names, verbose) { .Call(`_GWmodel3_gwr_basic_fit`, x, y, coords, bw, adaptive, kernel, longlat, p, theta, optim_bw_lower, optim_bw_upper, hatmatrix, intercept, parallel_type, parallel_arg, optim_bw, optim_bw_criterion, select_model, select_model_criterion, select_model_threshold, variable_names, verbose) } diff --git a/R/gwpca.R b/R/gwpca.R new file mode 100644 index 0000000..b6d9f25 --- /dev/null +++ b/R/gwpca.R @@ -0,0 +1,165 @@ +#' Calibrate a GWPCA model +#' +#' @param formula Regresison model. +#' @param data A `sf` objects. +#' @param bw Should provide a value to set the size of bandwidth, +#' @param adaptive Whether the bandwidth value is adaptive or not. +#' @param kernel Kernel function used. +#' @param longlat Whether the coordinates +#' @param p Power of the Minkowski distance, +#' default to 2, i.e., Euclidean distance. +#' @param theta Angle in radian to roate the coordinate system, default to 0. +#' @param components How many components want to keep, default to 2. +#' +#' @return A `gwpcam` object. +#' +#' @examples +#' data(LondonHP) +#' +#' # Basic usage +#' gwpca(~FLOORSZ + UNEMPLOY + PROF, LondonHP, 50, TRUE, components = 2) +#' +#' @seealso `browseVignettes("")` +#' +#' @importFrom stats na.action model.frame model.extract model.matrix terms +#' @export +gwpca <- function( + formula, + data, + bw = NA, + adaptive = FALSE, + kernel = c("gaussian", "exp", "bisquare", "tricube", "boxcar"), + longlat = FALSE, + p = 2.0, + theta = 0.0, + components = 2 +) { + ### Check args + kernel = match.arg(kernel) + attr(data, "na.action") <- getOption("na.action") + + ### Extract coords + data <- do.call(na.action(data), args = list(data)) + coords <- as.matrix(sf::st_coordinates(sf::st_centroid(data))) + if (is.null(coords) || nrow(coords) != nrow(data)) + stop("Missing coordinates.") + + ### Extract variables + mc <- match.call(expand.dots = FALSE) + + mf <- model.frame(formula = formula(update(formula, ~ . - 1)), data = sf::st_drop_geometry(data)) + mt <- attr(mf, "terms") + x <- model.matrix(mt, mf) + indep_vars <- colnames(x) + + ### Check whether bandwidth is valid. + if (missing(bw)) { + stop("Please input bandwidth!") + } + if (!(is.numeric(bw) || is.integer(bw))) { + stop("Please input correct bandwidth!") + } + + components <- as.integer(floor(components)) + if (components < 1) { + stop("Components must be an interger greater than 0 !!") + } + if(length(indep_vars) < components) { + stop("Components to keep must be lower than variable counts!") + } + + + + c_result <- tryCatch(gwpca_cal( + x, coords, + bw, adaptive, enum(kernel), longlat, p, theta, + components + ), error = function(e) { + stop("Error:", conditionMessage(e)) + }) + + local_loadings <- c_result$local_loadings + # local_scores <- c_result$local_scores + + local_PV <- c_result$local_PV + + pc_names <- paste0("PC", seq_len(components), "_PV") + colnames(local_PV) <- pc_names + sdf_data <- as.data.frame(local_PV) + sdf_data$geometry <- sf::st_geometry(data) + sdf <- sf::st_sf(sdf_data) + + ### Return result + gwpcam <- list( + SDF = sdf, + local_loadings = local_loadings, + # local_scores = local_scores, + args = list( + x = x, + coords = coords, + bw = bw, + adaptive = adaptive, + kernel = kernel, + longlat = longlat, + p = p, + theta = theta, + keep_components = components + ), + call = mc, + indep_vars = indep_vars + ) + class(gwpcam) <- "gwpcam" + gwpcam +} + +#' Print description of a `gwpcam` object +#' +#' @param x An `gwpcam` object returned by [gwpca()]. +#' @param decimal_fmt The format string passing to [base::sprintf()]. +#' @inheritDotParams print_table_md +#' +#' @method print gwpcam +#' @importFrom stats coef fivenum +#' @rdname print +#' @export +print.gwpcam <- function(x, decimal_fmt = "%.3f", ...) { + if (!inherits(x, "gwpcam")) { + stop("It's not a gwpcam object.") + } + + ### Basic Information + cat("Geographically Weighted Principal Component Analysis", fill = T) + cat("========================================", fill = T) + cat(" Formula:", deparse(x$call$formula), fill = T) + cat(" Data:", deparse(x$call$data), fill = T) + cat(" Kernel:", x$args$kernel, fill = T) + cat("Bandwidth:", x$args$bw, ifelse(x$args$adaptive, "(Nearest Neighbours)", "(Meters)"), fill = T) + cat(" Keep components:", x$args$keep_components, fill = T) + cat("\n", fill = T) + + cat("Summary of Local PV", fill = T) + cat("--------------------------------", fill = T) + local_PV <- sf::st_drop_geometry(x$SDF) + local_PV_fivenum <- t(apply(local_PV, 2, fivenum)) + colnames(local_PV_fivenum) <- c("Min.", "1st Qu.", "Median", "3rd Qu.", "Max.") + rownames(local_PV_fivenum) <- colnames(local_PV) + loadings_1st_str <- rbind( + c("Local loadings", colnames(local_PV_fivenum)), + cbind(rownames(local_PV_fivenum), matrix2char(local_PV_fivenum, decimal_fmt)) + ) + print_table_md(loadings_1st_str, ...) + cat("\n", fill = T) + + cat("Summary of Local Loadings in PC1", fill = T) + cat("--------------------------------", fill = T) + loadings_1st <- x$local_loadings[, , 1] + loadings_1st_fivenum <- t(apply(loadings_1st, 2, fivenum)) + colnames(loadings_1st_fivenum) <- c("Min.", "1st Qu.", "Median", "3rd Qu.", "Max.") + rownames(loadings_1st_fivenum) <- x$indep_vars + loadings_1st_str <- rbind( + c("Local loadings", colnames(loadings_1st_fivenum)), + cbind(rownames(loadings_1st_fivenum), matrix2char(loadings_1st_fivenum, decimal_fmt)) + ) + print_table_md(loadings_1st_str, ...) + cat("\n", fill = T) +} \ No newline at end of file diff --git a/man/gwpca.Rd b/man/gwpca.Rd new file mode 100644 index 0000000..5d79269 --- /dev/null +++ b/man/gwpca.Rd @@ -0,0 +1,54 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gwpca.R +\name{gwpca} +\alias{gwpca} +\title{Calibrate a GWPCA model} +\usage{ +gwpca( + formula, + data, + bw = NA, + adaptive = FALSE, + kernel = c("gaussian", "exp", "bisquare", "tricube", "boxcar"), + longlat = FALSE, + p = 2, + theta = 0, + components = 2 +) +} +\arguments{ +\item{formula}{Regresison model.} + +\item{data}{A \code{sf} objects.} + +\item{bw}{Should provide a value to set the size of bandwidth,} + +\item{adaptive}{Whether the bandwidth value is adaptive or not.} + +\item{kernel}{Kernel function used.} + +\item{longlat}{Whether the coordinates} + +\item{p}{Power of the Minkowski distance, +default to 2, i.e., Euclidean distance.} + +\item{theta}{Angle in radian to roate the coordinate system, default to 0.} + +\item{components}{How many components want to keep, default to 2.} +} +\value{ +A \code{gwpcam} object. +} +\description{ +Calibrate a GWPCA model +} +\examples{ +data(LondonHP) + +# Basic usage +gwpca(~FLOORSZ + UNEMPLOY + PROF, LondonHP, 50, TRUE, components = 2) + +} +\seealso{ +\code{browseVignettes("")} +} diff --git a/man/print.Rd b/man/print.Rd index e7e67ee..ffd70b8 100644 --- a/man/print.Rd +++ b/man/print.Rd @@ -1,13 +1,17 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/gtdr.R, R/gwr_basic.R, R/gwr_multiscale.R +% Please edit documentation in R/gtdr.R, R/gwpca.R, R/gwr_basic.R, +% R/gwr_multiscale.R \name{print.gtdrm} \alias{print.gtdrm} +\alias{print.gwpcam} \alias{print.gwrm} \alias{print.gwrmultiscalem} \title{Print description of a \code{gtdrm} object} \usage{ \method{print}{gtdrm}(x, decimal_fmt = "\%.3f", ...) +\method{print}{gwpcam}(x, decimal_fmt = "\%.3f", ...) + \method{print}{gwrm}(x, decimal_fmt = "\%.3f", ...) \method{print}{gwrmultiscalem}(x, decimal_fmt = "\%.3f", ...) @@ -18,7 +22,7 @@ \item{decimal_fmt}{The format string passing to \code{\link[base:sprintf]{base::sprintf()}}.} \item{...}{ - Arguments passed on to \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}} + Arguments passed on to \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}} \describe{ \item{\code{col.sep}}{Column seperator. Default to \code{""}.} \item{\code{header.sep}}{Header seperator. Default to \code{"-"}.} @@ -33,6 +37,8 @@ Possible values are \code{"plain"}, \code{"md"} or \code{"latex"}. Default to \c \description{ Print description of a \code{gtdrm} object +Print description of a \code{gwpcam} object + Print description of a \code{gwrm} object Print description of a \code{gwrmultiscalem} object diff --git a/src/Makevars.in b/src/Makevars.in index f5897d0..996f6bd 100644 --- a/src/Makevars.in +++ b/src/Makevars.in @@ -45,6 +45,7 @@ OBJECTS_LIBGWMODEL = \ libgwmodel/src/gwmodelpp/GWRBase.o \ libgwmodel/src/gwmodelpp/GWRBasic.o \ libgwmodel/src/gwmodelpp/GWRMultiscale.o \ + libgwmodel/src/gwmodelpp/GWPCA.o \ libgwmodel/src/gwmodelpp/SpatialAlgorithm.o \ libgwmodel/src/gwmodelpp/SpatialMonoscaleAlgorithm.o \ libgwmodel/src/gwmodelpp/SpatialMultiscaleAlgorithm.o \ @@ -61,6 +62,7 @@ OBJECTS_GWMODEL = \ gwr_basic.o \ gwr_multiscale.o \ gtdr.o \ + gwpca.o \ RcppExports.o OBJECTS_CXX = $(OBJECTS_LIBGWMODEL) $(OBJECTS_TELEGRAM) $(OBJECTS_GWMODEL) diff --git a/src/Makevars.win b/src/Makevars.win index 8ac438f..c0a316e 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -47,6 +47,7 @@ OBJECTS_LIBGWMODEL = \ libgwmodel/src/gwmodelpp/GWRBase.o \ libgwmodel/src/gwmodelpp/GWRBasic.o \ libgwmodel/src/gwmodelpp/GWRMultiscale.o \ + libgwmodel/src/gwmodelpp/GWPCA.o \ libgwmodel/src/gwmodelpp/SpatialAlgorithm.o \ libgwmodel/src/gwmodelpp/SpatialMonoscaleAlgorithm.o \ libgwmodel/src/gwmodelpp/SpatialMultiscaleAlgorithm.o \ @@ -63,6 +64,7 @@ OBJECTS_GWMODEL = \ gwr_basic.o \ gwr_multiscale.o \ gtdr.o \ + gwpca.o \ RcppExports.o diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index d0cedbc..60b7886 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -40,6 +40,25 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// gwpca_cal +List gwpca_cal(const arma::mat& x, const arma::mat& coords, double bw, bool adaptive, size_t kernel, bool longlat, double p, double theta, size_t keep_components); +RcppExport SEXP _GWmodel3_gwpca_cal(SEXP xSEXP, SEXP coordsSEXP, SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP, SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP, SEXP keep_componentsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type coords(coordsSEXP); + Rcpp::traits::input_parameter< double >::type bw(bwSEXP); + Rcpp::traits::input_parameter< bool >::type adaptive(adaptiveSEXP); + Rcpp::traits::input_parameter< size_t >::type kernel(kernelSEXP); + Rcpp::traits::input_parameter< bool >::type longlat(longlatSEXP); + Rcpp::traits::input_parameter< double >::type p(pSEXP); + Rcpp::traits::input_parameter< double >::type theta(thetaSEXP); + Rcpp::traits::input_parameter< size_t >::type keep_components(keep_componentsSEXP); + rcpp_result_gen = Rcpp::wrap(gwpca_cal(x, coords, bw, adaptive, kernel, longlat, p, theta, keep_components)); + return rcpp_result_gen; +END_RCPP +} // gwr_basic_fit List gwr_basic_fit(const arma::mat& x, const arma::vec& y, const arma::mat& coords, double bw, bool adaptive, size_t kernel, bool longlat, double p, double theta, double optim_bw_lower, double optim_bw_upper, bool hatmatrix, bool intercept, size_t parallel_type, const IntegerVector& parallel_arg, bool optim_bw, size_t optim_bw_criterion, bool select_model, size_t select_model_criterion, size_t select_model_threshold, const CharacterVector& variable_names, int verbose); RcppExport SEXP _GWmodel3_gwr_basic_fit(SEXP xSEXP, SEXP ySEXP, SEXP coordsSEXP, SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP, SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP, SEXP optim_bw_lowerSEXP, SEXP optim_bw_upperSEXP, SEXP hatmatrixSEXP, SEXP interceptSEXP, SEXP parallel_typeSEXP, SEXP parallel_argSEXP, SEXP optim_bwSEXP, SEXP optim_bw_criterionSEXP, SEXP select_modelSEXP, SEXP select_model_criterionSEXP, SEXP select_model_thresholdSEXP, SEXP variable_namesSEXP, SEXP verboseSEXP) { @@ -134,6 +153,7 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_GWmodel3_gtdr_fit", (DL_FUNC) &_GWmodel3_gtdr_fit, 19}, + {"_GWmodel3_gwpca_cal", (DL_FUNC) &_GWmodel3_gwpca_cal, 9}, {"_GWmodel3_gwr_basic_fit", (DL_FUNC) &_GWmodel3_gwr_basic_fit, 22}, {"_GWmodel3_gwr_basic_predict", (DL_FUNC) &_GWmodel3_gwr_basic_predict, 14}, {"_GWmodel3_gwr_multiscale_fit", (DL_FUNC) &_GWmodel3_gwr_multiscale_fit, 25}, diff --git a/src/gwpca.cpp b/src/gwpca.cpp new file mode 100644 index 0000000..e8a4af3 --- /dev/null +++ b/src/gwpca.cpp @@ -0,0 +1,70 @@ +// [[Rcpp::depends(RcppArmadillo)]] +#include +#include "utils.h" +#include "gwmodel.h" + +using namespace std; +using namespace Rcpp; +using namespace arma; +using namespace gwm; + +// [[Rcpp::export]] +List gwpca_cal( + const arma::mat& x, + const arma::mat& coords, + double bw, + bool adaptive, + size_t kernel, + bool longlat, + double p, + double theta, + size_t keep_components +) { + // Make Spatial Weight + BandwidthWeight bandwidth(bw, adaptive, BandwidthWeight::KernelFunctionType((size_t)kernel)); + Distance* distance = nullptr; + if (longlat) + { + distance = new CRSDistance(true); + } + else + { + if (p == 2.0 && theta == 0.0) + { + distance = new CRSDistance(false); + } + else + { + distance = new MinkwoskiDistance(p, theta); + } + } + SpatialWeight spatial(&bandwidth, distance); + + // Make Algorithm Object + GWPCA algorithm; + algorithm.setCoords(coords); + algorithm.setVariables(x); + algorithm.setSpatialWeight(spatial); + algorithm.setKeepComponents(keep_components); + try + { + algorithm.run(); + } + catch(const std::exception& e) + { + stop(e.what()); + } + + // Return Results + cube loadings = algorithm.loadings(); + + // defined but not calculate in libgwmodel + // cube scores = algorithm.scores(); + + List result_list = List::create( + Named("local_loadings") = loadings, + Named("local_PV") = algorithm.localPV() + // Named("local_scores") = scores + ); + return result_list; +} \ No newline at end of file diff --git a/tests/testthat/test-gwpca.R b/tests/testthat/test-gwpca.R new file mode 100644 index 0000000..17a4881 --- /dev/null +++ b/tests/testthat/test-gwpca.R @@ -0,0 +1,16 @@ +data(LondonHP) +m <- NULL + +test_that("GWPCA: works", { + m <<- expect_no_error(gwpca(~PURCHASE + FLOORSZ + UNEMPLOY + PROF, LondonHP, bw = 50, adaptive = TRUE, components = 3)) +}) + +test_that("GWPCA: error test", { + expect_error(gwpca(~FLOORSZ+UNEMPLOY, LondonHP, bw = 50, adaptive = TRUE, components = 3), + "Components to keep must be lower than variable counts!") +}) + +test_that("GWPCA: error components", { + expect_error(gwpca(~FLOORSZ+UNEMPLOY, LondonHP, bw = 50, adaptive = TRUE, components = 0), + "Components must be an interger greater than 0 !!") +}) \ No newline at end of file From 8906641c8abd4ee83dd376805d66bd322cb6eb99 Mon Sep 17 00:00:00 2001 From: OuGuangyu Date: Sat, 7 Jun 2025 11:35:50 +0800 Subject: [PATCH 2/2] add: glyph.plot for gwpca loadings, add scale option --- NAMESPACE | 2 + R/gwpca.R | 99 ++++++++++++++++++++++++++++++++++++++-- man/glyph.plot.Rd | 11 +++++ man/glyph.plot.gwpcam.Rd | 23 ++++++++++ man/gwpca.Rd | 5 +- 5 files changed, 136 insertions(+), 4 deletions(-) create mode 100644 man/glyph.plot.Rd create mode 100644 man/glyph.plot.gwpcam.Rd diff --git a/NAMESPACE b/NAMESPACE index 6e6a20d..519a71b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,6 +6,7 @@ S3method(coef,gwrmultiscalem) S3method(fitted,gtdrm) S3method(fitted,gwrm) S3method(fitted,gwrmultiscalem) +S3method(glyph.plot,gwpcam) S3method(plot,gtdrm) S3method(plot,gwrm) S3method(plot,gwrmultiscalem) @@ -21,6 +22,7 @@ S3method(residuals,gwrmultiscalem) S3method(step,default) S3method(step,gtdrm) S3method(step,gwrm) +export(glyph.plot) export(gtdr) export(gtdr_config) export(gwpca) diff --git a/R/gwpca.R b/R/gwpca.R index b6d9f25..5694601 100644 --- a/R/gwpca.R +++ b/R/gwpca.R @@ -10,7 +10,8 @@ #' default to 2, i.e., Euclidean distance. #' @param theta Angle in radian to roate the coordinate system, default to 0. #' @param components How many components want to keep, default to 2. -#' +#' @param scale Whether to standardised data. +#' #' @return A `gwpcam` object. #' #' @examples @@ -32,7 +33,8 @@ gwpca <- function( longlat = FALSE, p = 2.0, theta = 0.0, - components = 2 + components = 2, + scale = FALSE ) { ### Check args kernel = match.arg(kernel) @@ -50,6 +52,7 @@ gwpca <- function( mf <- model.frame(formula = formula(update(formula, ~ . - 1)), data = sf::st_drop_geometry(data)) mt <- attr(mf, "terms") x <- model.matrix(mt, mf) + if(scale) x <- scale(x) indep_vars <- colnames(x) ### Check whether bandwidth is valid. @@ -69,7 +72,6 @@ gwpca <- function( } - c_result <- tryCatch(gwpca_cal( x, coords, bw, adaptive, enum(kernel), longlat, p, theta, @@ -162,4 +164,95 @@ print.gwpcam <- function(x, decimal_fmt = "%.3f", ...) { ) print_table_md(loadings_1st_str, ...) cat("\n", fill = T) +} + +#' Glyph plot generic +#' @export +glyph.plot <- function(x, ...) { + UseMethod("glyph.plot") +} + +#' Glyph plot for local loadings in PC1 of GWPCA model. +#' +#' @param x A "gwpcam" object. +#' @param y Ignored. +#' @param columns Column names to plot. +#' If it is missing or non-character value, all coefficient columns are plotted. +#' @param sign whether to plot according to the loadings is beyond 0 or not. +#' @param \dots Additional arguments passed to [sf::plot()]. +#' @method glyph.plot gwpcam +#' @export +glyph.plot.gwpcam <- function(x, y, ..., columns, sign = FALSE) { + if (!inherits(x, "gwpcam")) { + stop("It's not a gwpcam object.") + } + + sdf <- sf::st_as_sf(x$SDF) + coords <- sf::st_coordinates(sdf) + r <- 0.05 * max(diff(range(coords[, 1])), diff(range(coords[, 2]))) + + ld <- x$local_loadings[, , 1] + colnames(ld) <- x$indep_vars + n.col <- ncol(ld) + + colors <- rainbow(n.col) + rgb_colors <- col2rgb(colors) / 255 + + rowmax <- function(z) z[cbind(1:nrow(z), max.col(abs(z)))] + ld <- sweep(ld, 1, sign(rowmax(ld)), "*") + + angles <- (0:(n.col - 1)) * 2 * pi / n.col + J <- 0 + (0 + 1i) + disp <- exp((pi / 2 - angles) * J) * r + loc2 <- coords[, 1] + coords[, 2] * J + ld.max <- max(ld) + ld.scaled <- abs(ld) / ld.max + + plot(coords, asp = 1, type = "n") + points(coords, pch = 14, cex = 0.1, col = "black") + + for (i in 1:nrow(ld)) { + for (j in 1:ncol(ld)) { + l.from <- loc2[i] + l.to <- loc2[i] + disp[j] * ld.scaled[i, j] + + if (sign) { + alpha <- 1 + col <- if (ld[i, j] > 0) { + rgb(1, 0, 0, alpha) # red + } else { + rgb(0, 0, 1, alpha) # blue + } + } else { + col_index <- (j - 1) %% n.col + 1 + col <- rgb(rgb_colors[1, col_index], + rgb_colors[2, col_index], + rgb_colors[3, col_index], + alpha = 1) + } + + lines(Re(c(l.from, l.to)), Im(c(l.from, l.to)), col = col) + } + } + + par(cex = 1.3) + legend_labels <- colnames(ld) + + if (!sign) { + text.w <- max(strwidth(legend_labels)) * 1.2 + legend("bottomleft", + legend = legend_labels, + col = colors, + lty = 1, + lwd = 2, + text.width = text.w, + bg = "white") + } else { + legend("bottomleft", + legend = c("Positive loading", "Negative loading"), + col = c("red", "blue"), + lty = 1, + lwd = 2, + bg = "white") + } } \ No newline at end of file diff --git a/man/glyph.plot.Rd b/man/glyph.plot.Rd new file mode 100644 index 0000000..a5ebac3 --- /dev/null +++ b/man/glyph.plot.Rd @@ -0,0 +1,11 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gwpca.R +\name{glyph.plot} +\alias{glyph.plot} +\title{Glyph plot generic} +\usage{ +glyph.plot(x, ...) +} +\description{ +Glyph plot generic +} diff --git a/man/glyph.plot.gwpcam.Rd b/man/glyph.plot.gwpcam.Rd new file mode 100644 index 0000000..bc264f0 --- /dev/null +++ b/man/glyph.plot.gwpcam.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gwpca.R +\name{glyph.plot.gwpcam} +\alias{glyph.plot.gwpcam} +\title{Glyph plot for local loadings in PC1 of GWPCA model.} +\usage{ +\method{glyph.plot}{gwpcam}(x, y, ..., columns, sign = FALSE) +} +\arguments{ +\item{x}{A "gwpcam" object.} + +\item{y}{Ignored.} + +\item{\dots}{Additional arguments passed to \code{\link[sf:plot]{sf::plot()}}.} + +\item{columns}{Column names to plot. +If it is missing or non-character value, all coefficient columns are plotted.} + +\item{sign}{whether to plot according to the loadings is beyond 0 or not.} +} +\description{ +Glyph plot for local loadings in PC1 of GWPCA model. +} diff --git a/man/gwpca.Rd b/man/gwpca.Rd index 5d79269..488cf6b 100644 --- a/man/gwpca.Rd +++ b/man/gwpca.Rd @@ -13,7 +13,8 @@ gwpca( longlat = FALSE, p = 2, theta = 0, - components = 2 + components = 2, + scale = FALSE ) } \arguments{ @@ -35,6 +36,8 @@ default to 2, i.e., Euclidean distance.} \item{theta}{Angle in radian to roate the coordinate system, default to 0.} \item{components}{How many components want to keep, default to 2.} + +\item{scale}{Whether to standardised data.} } \value{ A \code{gwpcam} object.