diff --git a/NAMESPACE b/NAMESPACE index 396e2d2..08b1435 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,13 +1,16 @@ # Generated by roxygen2: do not edit by hand S3method(coef,gtdrm) +S3method(coef,gwlcrm) S3method(coef,gwrm) S3method(coef,gwrmultiscalem) S3method(fitted,gtdrm) +S3method(fitted,gwlcrm) S3method(fitted,gwrm) S3method(fitted,gwrmultiscalem) S3method(plot,gtdrm) S3method(plot,gwcorrm) +S3method(plot,gwlcrm) S3method(plot,gwrm) S3method(plot,gwrmultiscalem) S3method(plot,modelselcritl) @@ -15,9 +18,11 @@ S3method(predict,gwrm) S3method(print,gtdrm) S3method(print,gwavgm) S3method(print,gwcorrm) +S3method(print,gwlcrm) S3method(print,gwrm) S3method(print,gwrmultiscalem) S3method(residuals,gtdrm) +S3method(residuals,gwlcrm) S3method(residuals,gwrm) S3method(residuals,gwrmultiscalem) S3method(step,default) @@ -29,6 +34,7 @@ export(gwaverage) export(gwcorr_config) export(gwcorrelation) export(gwr_basic) +export(gwr_lcr) export(gwr_multiscale) export(mgwr_config) export(print_table_md) diff --git a/R/RcppExports.R b/R/RcppExports.R index 3f3babd..43f8738 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -21,6 +21,10 @@ gwr_basic_predict <- function(pcoords, x, y, coords, bw, adaptive, kernel, longl .Call(`_GWmodel3_gwr_basic_predict`, pcoords, x, y, coords, bw, adaptive, kernel, longlat, p, theta, intercept, parallel_type, parallel_arg, verbose) } +gwr_lcr_fit <- function(x, y, coords, bw, adaptive, kernel, longlat, p, theta, lambda, lambda_adjust, cn_thresh, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw) { + .Call(`_GWmodel3_gwr_lcr_fit`, x, y, coords, bw, adaptive, kernel, longlat, p, theta, lambda, lambda_adjust, cn_thresh, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw) +} + gwr_multiscale_fit <- function(x, y, coords, bw, adaptive, kernel, longlat, p, theta, optim_bw, optim_bw_criterion, threashold, initial_type, centered, optim_bw_lower, optim_bw_upper, criterion, hatmatrix, intercept, retry_times, max_iterations, parallel_type, parallel_arg, variable_names, verbose) { .Call(`_GWmodel3_gwr_multiscale_fit`, x, y, coords, bw, adaptive, kernel, longlat, p, theta, optim_bw, optim_bw_criterion, threashold, initial_type, centered, optim_bw_lower, optim_bw_upper, criterion, hatmatrix, intercept, retry_times, max_iterations, parallel_type, parallel_arg, variable_names, verbose) } diff --git a/R/gwr_localcollinearity.R b/R/gwr_localcollinearity.R new file mode 100644 index 0000000..1780317 --- /dev/null +++ b/R/gwr_localcollinearity.R @@ -0,0 +1,324 @@ +#' Calibrate a GWR local collinearity model +#' +#' @param formula Regresison model. +#' @param data A `sf` objects. +#' @param bw Either a value to set the size of bandwidth, +#' or `CV` to set the criterion for bandwidth auto-optimization process. +#' Note that if `NA` or other non-numeric value is setted, +#' this parameter will be reset to `Inf`. +#' @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 lambda Option for a globally-defined (constant) ridge parameter. +#' Default is lambda=0, which gives a basic GWR fit +#' @param lambda_adjust A locally-varying ridge parameter.Default FALSE, refers to: +#' -i a basic GWR without a local ridge adjustment +#' (i.e. lambda=0, everywhere); +#' -ii a penalised GWR with a global ridge adjustment +#' (i.e. lambda is user-specified as some constant, other than 0 everywhere); +#' if TRUE, use cn.tresh to set the maximum condition number. +#' Here for locations with a condition number (for its local design matrix) +#' above this user-specified threshold, a local ridge parameter is found +#' @param cn_thresh maximum value for condition number, commonly set between 20 and 30 +#' @param hatmatrix If TRUE, great circle will be caculated. +#' @param parallel_method Parallel method, multithreading (`omp`) is available +#' @param parallel_arg Parallel method argument. +#' @param verbose Whether to print additional information. +#' +#' @return A `gwlcrm` object. +#' +#' @examples +#' data(LondonHP) +#' +#' # Basic usage +#' gwr_lcr(PURCHASE ~ FLOORSZ + UNEMPLOY, LondonHP, 64, TRUE) +#' +#' # Bandwidth Optimization +#' m <- gwr_lcr(PURCHASE ~ FLOORSZ + UNEMPLOY + PROF, LondonHP, 'CV', TRUE) +#' m +#' +#' @importFrom stats na.action model.frame model.extract model.matrix terms +#' @export +gwr_lcr <- function( + formula, + data, + bw = NA, + adaptive = FALSE, + kernel = c("gaussian", "exp", "bisquare", "tricube", "boxcar"), + longlat = FALSE, + p = 2.0, + theta = 0.0, + lambda = 0, + lambda_adjust = FALSE, + cn_thresh = 30, + hatmatrix = TRUE, + parallel_method = c("no", "omp"), + parallel_arg = c(0), + verbose = FALSE +) { + ### 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(formula), data = sf::st_drop_geometry(data)) + mt <- attr(mf, "terms") + y <- model.extract(mf, "response") + x <- model.matrix(mt, mf) + dep_var <- as.character(attr(terms(mf), "variables")[[2]]) + has_intercept <- attr(terms(mf), "intercept") == 1 + indep_vars <- colnames(x) + indep_vars[which(indep_vars == "(Intercept)")] <- "Intercept" + colnames(x) <- indep_vars + if (has_intercept && indep_vars[1] != "Intercept") { + stop("Please put Intercept to the first column.") + } + + ### Check whether bandwidth is valid. + if (missing(bw)) { + optim_bw <- TRUE + bw <- Inf + } else if (is.numeric(bw) || is.integer(bw)) { + optim_bw <- FALSE + } else { + optim_bw <- TRUE + if(is.character(bw)) + if(!match.arg(bw, c("CV"))){ + stop("Parameter bw must be 'CV' or a number") + } + bw <- Inf + } + + ### Check lambda and cnthresh + if (lambda<0 || lambda>1){ + stop("Error: lambda must in [0,1]") + } + if (cn_thresh>30 || cn_thresh<20){ + warning("cn_thresh is recommended in [20,30]") + } + + ### Call solver + c_result <- tryCatch(gwr_lcr_fit( + x, y, coords, + bw, adaptive, enum(kernel), longlat, p, theta, + lambda, lambda_adjust, cn_thresh, + has_intercept, hatmatrix, + enum_list(parallel_method, parallel_types), parallel_arg, + optim_bw + ), error = function (e) { + stop("Error:", conditionMessage(e)) + }) + if (optim_bw) + bw <- c_result$bandwidth + betas <- c_result$betas + yhat <- c_result$yhat + diagnostic <- c_result$diagnostic + resi <- y - yhat + local_cn <- c_result$localCN + local_lambda <- c_result$localLambda + # n_dp <- nrow(coords) + # rss_gw <- sum(resi * resi) + # sigma <- rss_gw / (n_dp - 2 * shat_trace[1] + shat_trace[2]) + # betas_se <- sqrt(sigma * betas_se) + # betas_tv <- betas / betas_se + + ### Create result Layer + colnames(betas) <- indep_vars + # colnames(betas_se) <- paste(indep_vars, "SE", sep = ".") + # colnames(betas_tv) <- paste(indep_vars, "TV", sep = ".") + sdf_data <- as.data.frame(cbind( + betas, + "yhat" = yhat, + "residual" = resi, + "localCN" = local_cn, + "localLambda" = local_lambda + # betas_se, + # betas_tv + )) + sdf_data$geometry <- sf::st_geometry(data) + sdf <- sf::st_sf(sdf_data) + + ### Return result + gwlcrm <- list( + SDF = sdf, + diagnostic = diagnostic, + args = list( + x = x, + y = y, + coords = coords, + bw = bw, + adaptive = adaptive, + kernel = kernel, + longlat = longlat, + p = p, + theta = theta, + lambda = lambda, + lambda_adjust = lambda_adjust, + cn_thresh = cn_thresh, + hatmatrix = hatmatrix, + has_intercept = has_intercept, + parallel_method = parallel_method, + parallel_arg = parallel_arg, + optim_bw = optim_bw + ), + call = mc, + indep_vars = indep_vars, + dep_var = dep_var + ) + class(gwlcrm) <- "gwlcrm" + gwlcrm +} + + + +#' Print description of a `gwlcrm` object +#' +#' @param x An `hgwlcrm` object returned by [gwr_lcr()]. +#' @param decimal_fmt The format string passing to [base::sprintf()]. +#' @inheritDotParams print_table_md +#' +#' @method print gwlcrm +#' @importFrom stats coef fivenum +#' @rdname print +#' @export +print.gwlcrm <- function(x, decimal_fmt = "%.3f", ...) { + if (!inherits(x, "gwlcrm")) { + stop("It's not a gwlcrm object.") + } + + ### Basic Information + cat("Results of Ridge Geographically Weighted Regression", 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)"), + ifelse(x$args$optim_bw, paste0( + "(Optimized accroding to CV)" + ), ""), fill = T) + cat("Lambda(ridge parameter for gwr ridge model):", x$args$lambda, + ifelse(x$args$lambda_adjust, " (Adjusted)", ""), + fill = T) + if(x$args$lambda_adjust) + cat(" cnThresh:", x$args$cn_thresh, fill = T) + + cat("\n", fill = T) + + cat("Summary of Coefficient Estimates", fill = T) + cat("--------------------------------", fill = T) + betas <- coef(x) + beta_fivenum <- t(apply(betas, 2, fivenum)) + colnames(beta_fivenum) <- c("Min.", "1st Qu.", "Median", "3rd Qu.", "Max.") + rownames(beta_fivenum) <- colnames(betas) + beta_str <- rbind( + c("Coefficient", colnames(beta_fivenum)), + cbind(rownames(beta_fivenum), matrix2char(beta_fivenum, decimal_fmt)) + ) + print_table_md(beta_str, ...) + cat("\n", fill = T) + + cat("Diagnostic Information", fill = T) + cat("----------------------", fill = T) + cat(" RSS:", x$diagnostic$RSS, fill = T) + cat(" ENP:", x$diagnostic$ENP, fill = T) + cat(" EDF:", x$diagnostic$EDF, fill = T) + cat(" R2:", x$diagnostic$RSquare, fill = T) + cat("R2adj:", x$diagnostic$RSquareAdjust, fill = T) + cat(" AIC:", x$diagnostic$AIC, fill = T) + cat(" AICc:", x$diagnostic$AICc, fill = T) + cat("\n", fill = T) +} + +#' @describeIn gwr_lcr Plot the result of GWR local collinearity model. +#' +#' @param x A "gwlcrm" object. +#' @param y Ignored. +#' @param columns Column names to plot. +#' If it is missing or non-character value, all coefficient columns are plottd. +#' @param \dots Additional arguments passing to [sf::plot()]. +#' @method plot gwlcrm +#' +#' @examples +#' plot(m) +#' +#' @export +plot.gwlcrm <- function(x, y, ..., columns) { + if (!inherits(x, "gwlcrm")) { + stop("It's not a gwlcrm object.") + } + + sdf <- sf::st_as_sf(x$SDF) + sdf_colnames <- names(sf::st_drop_geometry(x$SDF)) + if (!missing(columns) && is.character(columns)) { + valid_columns <- intersect(columns, sdf_colnames) + if (length(valid_columns) > 0) { + sdf <- sdf[valid_columns] + } + } else { ### Select coefficient columns. + sdf <- sdf[x$indep_vars] + } + plot(sdf, ...) +} + +#' @describeIn gwr_lcr Get coefficients of a GWR local collinearity model. +#' +#' @param object A "gwlcrm" object. +#' @param \dots Additional arguments passing to [coef()]. +#' +#' @examples +#' coef(m) +#' +#' @method coef gwlcrm +#' @export +coef.gwlcrm <- function(object, ...) { + if (!inherits(object, "gwlcrm")) { + stop("It's not a gwlcrm object.") + } + sf::st_drop_geometry(object$SDF[object$indep_vars]) +} + +#' @describeIn gwr_lcr Get fitted values of a GWR local collinearity model. +#' +#' @param object A "gwlcrm" object. +#' @param \dots Additional arguments passing to [fitted()]. +#' +#' @examples +#' fitted(m) +#' +#' @method fitted gwlcrm +#' @export +fitted.gwlcrm <- function(object, ...) { + if (!inherits(object, "gwlcrm")) { + stop("It's not a gwlcrm object.") + } + object$SDF[["yhat"]] +} + +#' @describeIn gwr_lcr Get residuals of a GWR local collinearity model. +#' +#' @param object A "gwlcrm" object. +#' @param \dots Additional arguments passing to [residuals()]. +#' +#' @examples +#' residuals(m) +#' +#' @method residuals gwlcrm +#' @export +residuals.gwlcrm <- function(object, ...) { + if (!inherits(object, "gwlcrm")) { + stop("It's not a gwlcrm object.") + } + object$SDF[["residual"]] +} diff --git a/man/gwr_lcr.Rd b/man/gwr_lcr.Rd new file mode 100644 index 0000000..77e7087 --- /dev/null +++ b/man/gwr_lcr.Rd @@ -0,0 +1,126 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gwr_localcollinearity.R +\name{gwr_lcr} +\alias{gwr_lcr} +\alias{plot.gwlcrm} +\alias{coef.gwlcrm} +\alias{fitted.gwlcrm} +\alias{residuals.gwlcrm} +\title{Calibrate a GWR local collinearity model} +\usage{ +gwr_lcr( + formula, + data, + bw = NA, + adaptive = FALSE, + kernel = c("gaussian", "exp", "bisquare", "tricube", "boxcar"), + longlat = FALSE, + p = 2, + theta = 0, + lambda = 0, + lambda_adjust = FALSE, + cn_thresh = 30, + hatmatrix = TRUE, + parallel_method = c("no", "omp"), + parallel_arg = c(0), + verbose = FALSE +) + +\method{plot}{gwlcrm}(x, y, ..., columns) + +\method{coef}{gwlcrm}(object, ...) + +\method{fitted}{gwlcrm}(object, ...) + +\method{residuals}{gwlcrm}(object, ...) +} +\arguments{ +\item{formula}{Regresison model.} + +\item{data}{A \code{sf} objects.} + +\item{bw}{Either a value to set the size of bandwidth, +or \code{CV} to set the criterion for bandwidth auto-optimization process. +Note that if \code{NA} or other non-numeric value is setted, +this parameter will be reset to \code{Inf}.} + +\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{lambda}{Option for a globally-defined (constant) ridge parameter. +Default is lambda=0, which gives a basic GWR fit} + +\item{lambda_adjust}{A locally-varying ridge parameter.Default FALSE, refers to: +-i a basic GWR without a local ridge adjustment +(i.e. lambda=0, everywhere); +-ii a penalised GWR with a global ridge adjustment +(i.e. lambda is user-specified as some constant, other than 0 everywhere); +if TRUE, use cn.tresh to set the maximum condition number. +Here for locations with a condition number (for its local design matrix) +above this user-specified threshold, a local ridge parameter is found} + +\item{cn_thresh}{maximum value for condition number, commonly set between 20 and 30} + +\item{hatmatrix}{If TRUE, great circle will be caculated.} + +\item{parallel_method}{Parallel method, multithreading (\code{omp}) is available} + +\item{parallel_arg}{Parallel method argument.} + +\item{verbose}{Whether to print additional information.} + +\item{x}{A "gwlcrm" object.} + +\item{y}{Ignored.} + +\item{\dots}{Additional arguments passing to \code{\link[=residuals]{residuals()}}.} + +\item{columns}{Column names to plot. +If it is missing or non-character value, all coefficient columns are plottd.} + +\item{object}{A "gwlcrm" object.} +} +\value{ +A \code{gwlcrm} object. +} +\description{ +Calibrate a GWR local collinearity model +} +\section{Functions}{ +\itemize{ +\item \code{plot(gwlcrm)}: Plot the result of GWR local collinearity model. + +\item \code{coef(gwlcrm)}: Get coefficients of a GWR local collinearity model. + +\item \code{fitted(gwlcrm)}: Get fitted values of a GWR local collinearity model. + +\item \code{residuals(gwlcrm)}: Get residuals of a GWR local collinearity model. + +}} +\examples{ +data(LondonHP) + +# Basic usage +gwr_lcr(PURCHASE ~ FLOORSZ + UNEMPLOY, LondonHP, 64, TRUE) + +# Bandwidth Optimization +m <- gwr_lcr(PURCHASE ~ FLOORSZ + UNEMPLOY + PROF, LondonHP, 'CV', TRUE) +m + +plot(m) + +coef(m) + +fitted(m) + +residuals(m) + +} diff --git a/man/print.Rd b/man/print.Rd index ea7a1e2..b48b76b 100644 --- a/man/print.Rd +++ b/man/print.Rd @@ -1,11 +1,12 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/gtdr.R, R/gwaverage.R, R/gwcorrelation.R, -% R/gwr_basic.R, R/gwr_multiscale.R +% R/gwr_basic.R, R/gwr_localcollinearity.R, R/gwr_multiscale.R \name{print.gtdrm} \alias{print.gtdrm} \alias{print.gwavgm} \alias{print.gwcorrm} \alias{print.gwrm} +\alias{print.gwlcrm} \alias{print.gwrmultiscalem} \title{Print description of a \code{gtdrm} object} \usage{ @@ -17,6 +18,8 @@ \method{print}{gwrm}(x, decimal_fmt = "\%.3f", ...) +\method{print}{gwlcrm}(x, decimal_fmt = "\%.3f", ...) + \method{print}{gwrmultiscalem}(x, decimal_fmt = "\%.3f", ...) } \arguments{ @@ -25,7 +28,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}}, \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}}, \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{"-"}.} @@ -46,5 +49,7 @@ Print description of a \code{gwcorrm} object Print description of a \code{gwrm} object +Print description of a \code{gwlcrm} object + Print description of a \code{gwrmultiscalem} object } diff --git a/src/Makevars.in b/src/Makevars.in index f939301..abce31c 100644 --- a/src/Makevars.in +++ b/src/Makevars.in @@ -47,6 +47,7 @@ OBJECTS_LIBGWMODEL = \ libgwmodel/src/gwmodelpp/GWRMultiscale.o \ libgwmodel/src/gwmodelpp/GWAverage.o \ libgwmodel/src/gwmodelpp/GWCorrelation.o \ + libgwmodel/src/gwmodelpp/GWRLocalCollinearity.o \ libgwmodel/src/gwmodelpp/SpatialAlgorithm.o \ libgwmodel/src/gwmodelpp/SpatialMonoscaleAlgorithm.o \ libgwmodel/src/gwmodelpp/SpatialMultiscaleAlgorithm.o \ @@ -63,6 +64,7 @@ OBJECTS_GWMODEL = \ utils.o \ gwr_basic.o \ gwr_multiscale.o \ + gwr_localcollinearity.o \ gtdr.o \ gwaverage.o \ gwcorrelation.o \ diff --git a/src/Makevars.win b/src/Makevars.win index a3a9120..9545035 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -49,6 +49,7 @@ OBJECTS_LIBGWMODEL = \ libgwmodel/src/gwmodelpp/GWRMultiscale.o \ libgwmodel/src/gwmodelpp/GWAverage.o \ libgwmodel/src/gwmodelpp/GWCorrelation.o \ + libgwmodel/src/gwmodelpp/GWRLocalCollinearity.o \ libgwmodel/src/gwmodelpp/SpatialAlgorithm.o \ libgwmodel/src/gwmodelpp/SpatialMonoscaleAlgorithm.o \ libgwmodel/src/gwmodelpp/SpatialMultiscaleAlgorithm.o \ @@ -65,6 +66,7 @@ OBJECTS_GWMODEL = \ utils.o \ gwr_basic.o \ gwr_multiscale.o \ + gwr_localcollinearity.o \ gtdr.o \ gwaverage.o \ gwcorrelation.o \ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index d4b7979..d784092 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -142,6 +142,33 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// gwr_lcr_fit +List gwr_lcr_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 lambda, bool lambda_adjust, double cn_thresh, bool intercept, bool hatmatrix, size_t parallel_type, const IntegerVector& parallel_arg, bool optim_bw); +RcppExport SEXP _GWmodel3_gwr_lcr_fit(SEXP xSEXP, SEXP ySEXP, SEXP coordsSEXP, SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP, SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP, SEXP lambdaSEXP, SEXP lambda_adjustSEXP, SEXP cn_threshSEXP, SEXP interceptSEXP, SEXP hatmatrixSEXP, SEXP parallel_typeSEXP, SEXP parallel_argSEXP, SEXP optim_bwSEXP) { +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::vec& >::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< double >::type lambda(lambdaSEXP); + Rcpp::traits::input_parameter< bool >::type lambda_adjust(lambda_adjustSEXP); + Rcpp::traits::input_parameter< double >::type cn_thresh(cn_threshSEXP); + Rcpp::traits::input_parameter< bool >::type intercept(interceptSEXP); + Rcpp::traits::input_parameter< bool >::type hatmatrix(hatmatrixSEXP); + Rcpp::traits::input_parameter< size_t >::type parallel_type(parallel_typeSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type parallel_arg(parallel_argSEXP); + Rcpp::traits::input_parameter< bool >::type optim_bw(optim_bwSEXP); + rcpp_result_gen = Rcpp::wrap(gwr_lcr_fit(x, y, coords, bw, adaptive, kernel, longlat, p, theta, lambda, lambda_adjust, cn_thresh, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw)); + return rcpp_result_gen; +END_RCPP +} // gwr_multiscale_fit List gwr_multiscale_fit(const arma::mat& x, const arma::vec& y, const arma::mat& coords, const NumericVector& bw, const LogicalVector& adaptive, const IntegerVector& kernel, const LogicalVector& longlat, const NumericVector& p, const NumericVector& theta, const LogicalVector& optim_bw, const IntegerVector& optim_bw_criterion, const NumericVector& threashold, const IntegerVector& initial_type, const LogicalVector& centered, double optim_bw_lower, double optim_bw_upper, size_t criterion, bool hatmatrix, bool intercept, size_t retry_times, size_t max_iterations, size_t parallel_type, const IntegerVector& parallel_arg, const CharacterVector& variable_names, int verbose); RcppExport SEXP _GWmodel3_gwr_multiscale_fit(SEXP xSEXP, SEXP ySEXP, SEXP coordsSEXP, SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP, SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP, SEXP optim_bwSEXP, SEXP optim_bw_criterionSEXP, SEXP threasholdSEXP, SEXP initial_typeSEXP, SEXP centeredSEXP, SEXP optim_bw_lowerSEXP, SEXP optim_bw_upperSEXP, SEXP criterionSEXP, SEXP hatmatrixSEXP, SEXP interceptSEXP, SEXP retry_timesSEXP, SEXP max_iterationsSEXP, SEXP parallel_typeSEXP, SEXP parallel_argSEXP, SEXP variable_namesSEXP, SEXP verboseSEXP) { @@ -184,6 +211,7 @@ static const R_CallMethodDef CallEntries[] = { {"_GWmodel3_gw_correlation_cal", (DL_FUNC) &_GWmodel3_gw_correlation_cal, 15}, {"_GWmodel3_gwr_basic_fit", (DL_FUNC) &_GWmodel3_gwr_basic_fit, 22}, {"_GWmodel3_gwr_basic_predict", (DL_FUNC) &_GWmodel3_gwr_basic_predict, 14}, + {"_GWmodel3_gwr_lcr_fit", (DL_FUNC) &_GWmodel3_gwr_lcr_fit, 17}, {"_GWmodel3_gwr_multiscale_fit", (DL_FUNC) &_GWmodel3_gwr_multiscale_fit, 25}, {NULL, NULL, 0} }; diff --git a/src/gwr_localcollinearity.cpp b/src/gwr_localcollinearity.cpp new file mode 100644 index 0000000..06dff98 --- /dev/null +++ b/src/gwr_localcollinearity.cpp @@ -0,0 +1,109 @@ +// [[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 gwr_lcr_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 lambda, + bool lambda_adjust, + double cn_thresh, + bool intercept, + bool hatmatrix, + size_t parallel_type, + const IntegerVector& parallel_arg, + bool optim_bw +) { + vector vpar_args = as< vector >(IntegerVector(parallel_arg)); + + // 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); + + GWRLocalCollinearity algorithm; + algorithm.setCoords(coords); + algorithm.setDependentVariable(y); + algorithm.setIndependentVariables(x); + algorithm.setSpatialWeight(spatial); + algorithm.setHasIntercept(intercept); + algorithm.setHasHatMatrix(hatmatrix); + + algorithm.setLambdaAdjust(lambda_adjust); + algorithm.setLambda(lambda); + algorithm.setCnThresh(cn_thresh); + + switch (ParallelType(size_t(parallel_type))) + { + case ParallelType::SerialOnly: + algorithm.setParallelType(ParallelType::SerialOnly); + break; +#ifdef _OPENMP + case ParallelType::OpenMP: + algorithm.setParallelType(ParallelType::OpenMP); + algorithm.setOmpThreadNum(vpar_args[0]); + break; +#endif + default: + algorithm.setParallelType(ParallelType::SerialOnly); + break; + } + + if (optim_bw) + { + algorithm.setIsAutoselectBandwidth(true); + algorithm.setBandwidthSelectionCriterion(GWRLocalCollinearity::BandwidthSelectionCriterionType::CV); + } + + + + algorithm.fit(); + + // Return Results + mat betas = algorithm.betas(); + vec yhat = sum(betas % x, 1); + List result_list = List::create( + Named("betas") = betas, + Named("diagnostic") = mywrap(algorithm.diagnostic()), + Named("yhat") = yhat, + Named("localCN") = algorithm.localCN(), + Named("localLambda") = algorithm.localLambda() + ); + + if (optim_bw){ + double bw_value = algorithm.spatialWeight().weight()->bandwidth(); + result_list["bandwidth"] = wrap(bw_value); + } + + return result_list; +} diff --git a/tests/testthat/test-gwr_localcollinearity.R b/tests/testthat/test-gwr_localcollinearity.R new file mode 100644 index 0000000..b78a598 --- /dev/null +++ b/tests/testthat/test-gwr_localcollinearity.R @@ -0,0 +1,38 @@ +data(LondonHP) +m <- NULL + +test_that("GWR LocalCollinearity: works", { + m1 <<- expect_no_error( + gwr_lcr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, 64, TRUE) + ) +}) + +test_that("GWR LocalCollinearity: bw selection", { + expect_no_error( + gwr_lcr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP) + ) +}) + +test_that("GWR LocalCollinearity: omp", { + expect_no_error( + gwr_lcr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, parallel_method = "omp", parallel_arg = 4) + ) +}) + +test_that("GWR LocalCollinearity: lambda & cnthresh", { + m2 <<- expect_no_error( + gwr_lcr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, 64, TRUE, lambda = 0.3, cn_thresh = 30) + ) +}) + +test_that("GWR LocalCollinearity: lambda & cnthresh", { + m3 <<- expect_no_error( + gwr_lcr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, 64, TRUE, lambda_adjust = TRUE,cn_thresh = 30) + ) +}) + +test_that("GWR LocalCollinearity: lambda error", { + expect_error({ + gwr_lcr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, 64, TRUE, lambda = 10, cn_thresh = 30) + }, "Error: lambda must in \\[0,1\\]") +})