Feature: GWDA

NAMESPACE

@@ -22,6 +22,7 @@ S3method(print,ggwrm)
 S3method(print,gtdrm)
 S3method(print,gwavgm)
 S3method(print,gwcorrm)
+S3method(print,gwdam)
 S3method(print,gwlcrm)
 S3method(print,gwpcam)
 S3method(print,gwrm)
@@ -39,6 +40,7 @@ export(gtdr_config)
 export(gwaverage)
 export(gwcorr_config)
 export(gwcorrelation)
+export(gwda)
 export(gwpca)
 export(gwr_basic)
 export(gwr_generalized)

R/RcppExports.R

@@ -13,6 +13,10 @@ gw_correlation_cal <- function(x1, x2, coords, bw, adaptive, kernel, longlat, p,
     .Call(`_GWmodel3_gw_correlation_cal`, x1, x2, coords, bw, adaptive, kernel, longlat, p, theta, initial_type, optim_bw_criterion, parallel_type, parallel_arg, variable_names, verbose)
 }
 
+gwda_cal <- function(x, y, coords, bw, adaptive, kernel, longlat, p, theta, method, parallel_type, parallel_arg) {
+    .Call(`_GWmodel3_gwda_cal`, x, y, coords, bw, adaptive, kernel, longlat, p, theta, method, parallel_type, parallel_arg)
+}
+
 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)
 }

R/gwda.R

@@ -0,0 +1,176 @@
+#' Geographically Weighted Summary Statistics
+#'
+#' @param formula Regresison model.
+#' @param data A `sf` objects.
+#' @param bw Bandwidth value
+#' @param adaptive Whether the bandwidth value is adaptive or not.
+#' @param quantile Whether to calculate local quantiles.
+#' @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 parallel_method Parallel method.
+#' @param parallel_arg Parallel method argument.
+#'
+#' @return A `gwdam` object.
+#'
+#' @export
+gwda <- function(
+    formula,
+    data,
+    bw,
+    adaptive = FALSE,
+    kernel = c("gaussian", "exp", "bisquare", "tricube", "boxcar"),
+    longlat = FALSE,
+    p = 2.0,
+    theta = 0.0,
+    method = c("wqda", "wlda"),
+    parallel_method = c("no", "omp"),
+    parallel_arg = c(0))
+{
+    ### Check args
+    kernel <- match.arg(kernel)
+    parallel_method <- match.arg(parallel_method)
+    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")
+    y <- model.extract(mf, "response")
+    y <- as.character(y)
+    x <- model.matrix(mt, mf)
+    dep_var <- as.character(attr(terms(mf), "variables")[[2]])
+    indep_vars <- colnames(x)
+
+    ### Check whether bandwidth is valid.
+    if (missing(bw)) {
+        stop("Bandwidth is missing.")
+    } else if (is.numeric(bw) || is.integer(bw)) {
+        if (bw == 0) stop("Bandwidth must be non-zero.")
+    } else {
+        stop("Bandwidth must be a number.")
+    }
+
+    method <- match.arg(method)
+    wqda_method <- ifelse(method == "wqda", TRUE, FALSE)
+
+    c_results <- gwda_cal(
+        x, y, coords, bw,
+        adaptive, enum(kernel),
+        longlat, p, theta,
+        wqda_method,
+        enum_list(parallel_method, parallel_types), parallel_arg
+    )
+
+    sdf_data <- data.frame()
+    res <- c_results$res
+    group <- as.data.frame(c_results$group)
+    probs <- c_results$probs
+    pmax <- c_results$pmax
+    entropy <- c_results$entropy
+
+    colnames(res) <- indep_vars
+    colnames(group) <- "groups"
+    colnames(probs) <- paste0(indep_vars, "_probs")
+    colnames(pmax) <- "pmax"
+    colnames(entropy) <- "entropy"
+    sdf_data <- as.data.frame(cbind(res, group, probs, pmax, entropy))
+
+    sdf_data$geometry <- sf::st_geometry(data)
+    sdf <- sf::st_sf(sdf_data)
+
+    ### Return result
+    gwdam <- list(
+        SDF = sdf,
+        args = list(
+            x = x,
+            y = y,
+            correctRate = c_results$correctRate,
+            coords = coords,
+            bw = bw,
+            adaptive = adaptive,
+            kernel = kernel,
+            longlat = longlat,
+            p = p,
+            theta = theta,
+            method = method,
+            parallel_method = parallel_method,
+            parallel_arg = parallel_arg
+        ),
+        call = mc,
+        indep_vars = indep_vars,
+        dep_var = dep_var
+    )
+    class(gwdam) <- "gwdam"
+    gwdam
+}
+
+#' Print description of a `gwdam` object
+#'
+#' @param x An `gwdam` object returned by [gwr_basic()].
+#' @param decimal_fmt The format string passing to [base::sprintf()].
+#' @inheritDotParams print_table_md
+#'
+#' @method print gwdam
+#' @rdname print
+#' @export
+print.gwdam <- function(x, ..., decimal_fmt) {
+    if (!inherits(x, "gwdam")) {
+        stop("It's not a 'gwdam' object.")
+    }
+
+    cat("   ***********************************************************************\n")
+    cat("   *                         Package   GWmodel3                          *\n")
+    cat("   ***********************************************************************\n")
+    cat("   *            Results of Geographically Discriminant Analysis          *\n")
+    cat("   ***********************************************************************\n")
+    cat("\n   *********************Model Calibration Information*********************\n")
+
+    cat("   Formula:", deparse(x$call$formula), fill = T)
+    cat("   Data:", deparse(x$call$data), fill = T)
+    cat("   Method:", x$args$method, fill = T)
+    cat("   Number of summary points:", nrow(x$args$x), fill = T)
+    cat("   Kernel function:", x$args$kernel, fill = T)
+    cat("   Bandwidth:", x$args$bw,
+        ifelse(x$args$adaptive, "(Nearest Neighbours)", "(Meters)"),
+        fill = T
+    )
+    distance_type <- "Euclidean"
+    if (x$args$longlat) {
+        distance_type <- "Geodetic"
+    } else if (x$args$p == 2) {
+        distance_type <- "Euclidean"
+    } else if (x$args$p == 1) {
+        distance_type <- "Manhattan"
+    } else if (is.infinite(x$args$p)) {
+        distance_type <- "Chebyshev"
+    } else {
+        distance_type <- "Generalized Minkowski"
+    }
+    distance_rotated <- (x$args$theta != 0 && x$args$p != 2 && !x$args$longlat)
+    cat("   Distance:", distance_type, ifelse(distance_rotated, " (rotated)", ""), fill = T)
+    res <- st_drop_geometry(x$SDF)
+    group_counts <- as.data.frame(table(res$groups))
+    cat("   Discrimination Result:", fill = T)
+    colnames(group_counts) <- c("Group", "Count")
+    output <- capture.output(print(group_counts, row.names = FALSE))
+    output <- paste0("  ", output)
+    cat(output, sep = "\n")
+
+    cat("\n   ***********************Local Summary Statistics************************", fill = T)
+    numeric_res <- res[sapply(res, is.numeric)]
+    df <- as.data.frame(t(sapply(numeric_res, summary)))
+    output <- capture.output(print(df))
+    output <- paste0("   ", output)
+    cat(output, sep = "\n")
+    cat("   ***********************************************************************\n")
+}

src/Makevars.in

@@ -50,6 +50,7 @@ OBJECTS_LIBGWMODEL = \
 	libgwmodel/src/gwmodelpp/GWRGeneralized.o \
 	libgwmodel/src/gwmodelpp/GWRMultiscale.o \
 	libgwmodel/src/gwmodelpp/GWAverage.o \
+	libgwmodel/src/gwmodelpp/GWDA.o \
 	libgwmodel/src/gwmodelpp/GWCorrelation.o \
 	libgwmodel/src/gwmodelpp/GWRLocalCollinearity.o \
 	libgwmodel/src/gwmodelpp/GWPCA.o \
@@ -73,6 +74,7 @@ OBJECTS_GWMODEL = \
 	gwr_localcollinearity.o \
 	gtdr.o \
 	gwaverage.o \
+	gwda.o \
 	gwcorrelation.o \
 	gwpca.o \
 	RcppExports.o

src/Makevars.win

@@ -52,6 +52,7 @@ OBJECTS_LIBGWMODEL = \
 	libgwmodel/src/gwmodelpp/GWRGeneralized.o \
 	libgwmodel/src/gwmodelpp/GWRMultiscale.o \
 	libgwmodel/src/gwmodelpp/GWAverage.o \
+	libgwmodel/src/gwmodelpp/GWDA.o \
 	libgwmodel/src/gwmodelpp/GWCorrelation.o \
 	libgwmodel/src/gwmodelpp/GWRLocalCollinearity.o \
 	libgwmodel/src/gwmodelpp/GWPCA.o \
@@ -75,6 +76,7 @@ OBJECTS_GWMODEL = \
 	gwr_localcollinearity.o \
 	gtdr.o \
 	gwaverage.o \
+	gwda.o \
 	gwcorrelation.o \
 	gwpca.o \
 	RcppExports.o

src/RcppExports.cpp

@@ -86,6 +86,28 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// gwda_cal
+List gwda_cal(const arma::mat& x, std::vector<std::string>& y, const arma::mat& coords, double bw, bool adaptive, size_t kernel, bool longlat, double p, double theta, bool method, size_t parallel_type, const IntegerVector& parallel_arg);
+RcppExport SEXP _GWmodel3_gwda_cal(SEXP xSEXP, SEXP ySEXP, SEXP coordsSEXP, SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP, SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP, SEXP methodSEXP, SEXP parallel_typeSEXP, SEXP parallel_argSEXP) {
+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< std::vector<std::string>& >::type y(ySEXP);
+    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< bool >::type method(methodSEXP);
+    Rcpp::traits::input_parameter< size_t >::type parallel_type(parallel_typeSEXP);
+    Rcpp::traits::input_parameter< const IntegerVector& >::type parallel_arg(parallel_argSEXP);
+    rcpp_result_gen = Rcpp::wrap(gwda_cal(x, y, coords, bw, adaptive, kernel, longlat, p, theta, method, parallel_type, parallel_arg));
+    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) {
@@ -279,6 +301,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_GWmodel3_gtdr_fit", (DL_FUNC) &_GWmodel3_gtdr_fit, 19},
     {"_GWmodel3_gwaverage_fit", (DL_FUNC) &_GWmodel3_gwaverage_fit, 11},
     {"_GWmodel3_gw_correlation_cal", (DL_FUNC) &_GWmodel3_gw_correlation_cal, 15},
+    {"_GWmodel3_gwda_cal", (DL_FUNC) &_GWmodel3_gwda_cal, 12},
     {"_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},