updated package in response to CRAN reviewer comments

.Rbuildignore

@@ -7,7 +7,7 @@
 ^.*\.Rproj$
 
 ## remove default C++ exported functions
-R/RcppExports.R
+R/RcppExports\.R
 
 ## raw data folder
 data-raw/
@@ -34,3 +34,6 @@ README\.html
 build_pipeline/
 builds/
 figure/
+
+# unexported examples
+R/unexported_examples\.R

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: remotePARTS
 Title: Spatiotemporal Autoregression Analyses for Large Data Sets
 Version: 1.0.4
-Authors@R: 
+Authors@R:
     c(person(given = "Clay",
            family = "Morrow",
            role = c("aut", "cre"),
@@ -12,16 +12,16 @@ Authors@R:
            role = c("aut"),
            email = "arives@wisc.edu",
            comment = c(ORCID = "0000-0001-9375-9523")))
-Description: 
-  These tools were created to test map-scale hypotheses about trends in large 
+Description:
+  These tools were created to test map-scale hypotheses about trends in large
   remotely sensed data sets but any data with spatial and temporal variation
   can be analyzed. Tests are conducted using the PARTS method for analyzing spatially
-  autocorrelated time series 
-  (Ives et al., 2021: <doi:10.1016/j.rse.2021.112678>). 
-  The method's unique approach can handle extremely large data sets that other 
-  spatiotemporal models cannot, while still appropriately accounting for 
-  spatial and temporal autocorrelation. This is done by partitioning the data 
-  into smaller chunks, analyzing chunks separately and then combining the 
+  autocorrelated time series
+  (Ives et al., 2021: <doi:10.1016/j.rse.2021.112678>).
+  The method's unique approach can handle extremely large data sets that other
+  spatiotemporal models cannot, while still appropriately accounting for
+  spatial and temporal autocorrelation. This is done by partitioning the data
+  into smaller chunks, analyzing chunks separately and then combining the
   separate analyses into a single, correlated test of the map-scale hypotheses.
 URL: https://github.com/morrowcj/remotePARTS
 BugReports: https://github.com/morrowcj/remotePARTS/issues
@@ -30,16 +30,16 @@ Encoding: UTF-8
 LazyData: TRUE
 RoxygenNote: 7.2.3
 Depends: R (>= 4.0)
-Imports: 
+Imports:
   stats,
-  geosphere (>= 1.5.10), 
-  Rcpp (>= 1.0.5), 
+  geosphere (>= 1.5.10),
+  Rcpp (>= 1.0.5),
   CompQuadForm,
   foreach,
   parallel,
-  iterators, 
+  iterators,
   doParallel
-Suggests: 
+Suggests:
     dplyr (>= 1.0.0),
     data.table,
     knitr,

R/AR_reml.R

@@ -31,8 +31,7 @@
 #' and \eqn{\rho}{rho} is the AR(1) autoregression parameter
 #'
 #' \code{fitAR} estimates the parameter via mathematical optimization
-#' of the restricted log-likelihood function calculated by the workhorse function
-#' \code{remotePARTS:::AR_fun()}.
+#' of the restricted log-likelihood function.
 #'
 #' @references
 #'
@@ -151,32 +150,6 @@ fitAR <- function(formula, data = NULL){
 #' linearly and negatively related to the restricted log likelihood
 #' (i.e., partial log-likelihood).
 #'
-#' @examples
-#'
-#' # simulate dummy data
-#' x = rnorm(31)
-#' time = 1:30
-#' x = x[2:31] + x[1:30] + 0.3*time #AR(1) process + time trend
-#' U = stats::model.matrix(formula(x ~ time))
-#'
-#' # fit an AR
-#' remotePARTS:::AR_fun(par = .2, y = x, X = U, logLik.only = FALSE)
-#'
-#' # get the partial logLik of the AR parameter, given the data.
-#' remotePARTS:::AR_fun(par = .2, y = x, X = U, logLik.only = FALSE)
-#'
-#' # show that minimizing the partial logLik maximizes the true logLik (NOT RUN)
-#' # n = 100
-#' # out.mat = matrix(NA, nrow = n, ncol = 3,
-#' #                  dimnames = list(NULL, c("par", "logLik", "partialLL")))
-#' # out.mat[, "par"] = seq(-10, 10, length.out = n)
-#' # for (i in seq_len(n) ) {
-#' #    p = out.mat[i, "par"]
-#' #    out.mat[i, "logLik"] = remotePARTS:::AR_fun(par = p, y = x, X = U, logLik.only = TRUE)
-#' #    out.mat[i, "partialLL"] = remotePARTS:::AR_fun(par = p, y = x, X = U,
-#' #                                                   logLik.only = FALSE)$logLik
-#' # }
-#' # plot(x = out.mat[, "partialLL"], y = out.mat[, "logLik"])
 AR_fun <- function(par, y, X, logLik.only = TRUE) {
   call = match.call()
 
@@ -319,7 +292,7 @@ AR_fun <- function(par, y, X, logLik.only = TRUE) {
 #'  ## create empty spatiotemporal variables:
 #'  X <- matrix(NA, nrow = nrow(coords), ncol = time.points) # response
 #'  Z <- matrix(NA, nrow = nrow(coords), ncol = time.points) # predictor
-#'  # setup first time point:
+#' # setup first time point:
 #'  Z[, 1] <- .05*coords[,"x"] + .2*coords[,"y"]
 #'  X[, 1] <- .5*Z[, 1] + rnorm(nrow(coords), 0, .05) #x at time t
 #'  ## project through time:
@@ -328,13 +301,13 @@ AR_fun <- function(par, y, X, logLik.only = TRUE) {
 #'    X[, t] <- .2*X[, t-1] + .1*Z[, t] + .05*t + rnorm(nrow(coords), 0 , .25)
 #'  }
 #'
-#' # # visualize dummy data (NOT RUN)
-#' # library(ggplot2);library(dplyr)
-#' # data.frame(coords, X) %>%
-#' #   reshape2::melt(id.vars = c("x", "y")) %>%
-#' #   ggplot(aes(x = x, y = y, fill = value)) +
-#' #   geom_tile() +
-#' #   facet_wrap(~variable)
+#' # visualize dummy data (NOT RUN)
+#' library(ggplot2);library(dplyr)
+#' data.frame(coords, X) %>%
+#'   reshape2::melt(id.vars = c("x", "y")) %>%
+#'   ggplot(aes(x = x, y = y, fill = value)) +
+#'   geom_tile() +
+#'   facet_wrap(~variable)
 #'
 #' # fit AR, showing all output
 #' fitAR_map(X, coords, formula = y ~ t, resids.only = TRUE)

R/check_posdef.R

@@ -5,6 +5,14 @@
 #'
 #' @param M numeric matrix
 #'
+#' @return returns a named logical vector with the following elements:
+#'
+#' \describe{
+#'     \item{sqr}{logical: indicating whether \code{M} is square}
+#'     \item{sym}{logical: indicating whether \code{M} is symmetric}
+#'     \item{posdef}{logical: indicating whether \code{M} is positive-definitive}
+#' }
+#'
 #' @examples
 #'
 #' # distance matrix

R/fitCor.R

@@ -87,12 +87,15 @@
 #' time.points = 30 # time series length
 #' map.width = 8 # square map width
 #' coords = expand.grid(x = 1:map.width, y = 1:map.width) # coordinate matrix
+#'
 #' ## create empty spatiotemporal variables:
 #' X <- matrix(NA, nrow = nrow(coords), ncol = time.points) # response
 #' Z <- matrix(NA, nrow = nrow(coords), ncol = time.points) # predictor
-#' # setup first time point:
+#'
+#' ## setup first time point:
 #' Z[, 1] <- .05*coords[,"x"] + .2*coords[,"y"]
 #' X[, 1] <- .5*Z[, 1] + rnorm(nrow(coords), 0, .05) #x at time t
+#'
 #' ## project through time:
 #' for(t in 2:time.points){
 #'   Z[, t] <- Z[, t-1] + rnorm(map.width^2)
@@ -104,6 +107,7 @@
 #' # using pre-defined covariance function
 #' ## exponential covariance
 #' fitCor(AR.map$residuals, coords, covar_FUN = "covar_exp", start = list(range = .1))
+#'
 #' ## exponential-power covariance
 #' fitCor(AR.map$residuals, coords, covar_FUN = "covar_exppow", start = list(range = .1, shape = .2))
 #'
@@ -115,11 +119,11 @@
 #'
 #' # specify which pixels to use, for reproducibility
 #' fitCor(AR.map$residuals, coords, index = 1:64)$spcor #all
-#' # fitCor(AR.map$residuals, coords, index = 1:20)$spcor #first 20
-#' # fitCor(AR.map$residuals, coords, index = 21:64)$spcor # last 43
+#' fitCor(AR.map$residuals, coords, index = 1:20)$spcor #first 20
+#' fitCor(AR.map$residuals, coords, index = 21:64)$spcor # last 43
 #' # randomly select pixels
-#' # fitCor(AR.map$residuals, coords, fit.n = 20)$spcor #random 20
-#' # fitCor(AR.map$residuals, coords, fit.n = 20)$spcor #random 20 again
+#' fitCor(AR.map$residuals, coords, fit.n = 20)$spcor #random 20
+#' fitCor(AR.map$residuals, coords, fit.n = 20)$spcor # different random 20
 #'
 #' @export
 fitCor <- function(resids, coords, distm_FUN = "distm_scaled", covar_FUN = "covar_exp",
@@ -195,6 +199,8 @@ fitCor <- function(resids, coords, distm_FUN = "distm_scaled", covar_FUN = "cova
 #'
 #' @param x remoteCor object to print
 #' @param ... additional arguments passed to print()
+#'
+#' @return a print-formatted version of key elements of the "remoteCor" object.
 print.remoteCor <- function(x, ...){
   # list(call = call, mod = fit, spcor = spcor, max.distance = max.d, logLik = logLik(fit))
   cat("\nCall:\n")
@@ -215,29 +221,6 @@ print.remoteCor <- function(x, ...){
 #' @param covar.pars vector or list of parameters (other than d) passed to the
 #' covar function
 #'
-#' @examples
-#' #  # distance vector
-#' #  d = seq(0, 1, length.out = 10)
-#' #  # named parameter list
-#' #  test_covar_fun(d = d, covar_FUN = "covar_exppow", covar.pars = list(range = .5))
-#' #  test_covar_fun(d = d, covar_FUN = "covar_exppow", covar.pars = list(range = .5, shape = .5))
-#' #  # unnamed parameter vector
-#' #  test_covar_fun(d = d, covar_FUN = "covar_exppow", covar.pars = .5)
-#' #  test_covar_fun(d = d, covar_FUN = "covar_exppow", covar.pars = c(.5, .5))
-#' #  # different function
-#' #  test_covar_fun(d = d, covar_FUN = "covar_taper", covar.pars = list(theta = .5))
-#' #  # user-defined function, with no extra parameters
-#' #  test_covar_fun(d = d, covar_FUN = function(d){return(d)}, covar.pars = NULL)
-#' #  test_covar_fun(d = d, covar_FUN = function(d){return(d)}, covar.pars = list())
-#' #
-#' #  map.width = 5 # square map width
-#' #  coords = expand.grid(x = 1:map.width, y = 1:map.width) # coordinate matrix
-#' #
-#' #  # calculate distance
-#' #  D = geosphere::distm(coords) # distance matrix
-#' #
-#' #  test_covar_fun(D, covar.pars = list(range = .1*max(D)))
-#'
 test_covar_fun <- function(d, covar_FUN = "covar_exppow", covar.pars = list(range = .5)){
   cov_f <- match.fun(covar_FUN)
   # covar.pars$d = d

R/fitGLS.R

@@ -141,7 +141,6 @@
 #'
 #' @examples
 #'
-#' \donttest{
 #' ## read data
 #' data(ndvi_AK10000)
 #' df = ndvi_AK10000[seq_len(200), ] # first 200 rows
@@ -174,7 +173,7 @@
 #'
 #' ## Log-likelihood (fast)
 #' fitGLS(CLS_coef ~ 0 + land, data = df, V = V, logLik.only = TRUE)
-#' }
+#'
 #' @export
 fitGLS <- function(formula, data, V, nugget = 0, formula0 = NULL, save.xx = FALSE,
                    save.invchol = FALSE, logLik.only = FALSE, no.F = FALSE,

R/fitGLS_opt.R

@@ -100,7 +100,7 @@
 #' \donttest{
 #' ## read data
 #' data(ndvi_AK10000)
-#' df = ndvi_AK10000[seq_len(200), ] # first 500 rows
+#' df = ndvi_AK10000[seq_len(200), ] # first 200 rows
 #'
 #' ## estimate nugget and range (very slow)
 #' fitGLS_opt(formula = CLS_coef ~ 0 + land, data = df,
@@ -259,32 +259,6 @@ fitGLS_opt <- function(formula, data = NULL, coords, distm_FUN = "distm_scaled",
 #' @return \code{fitGLS_opt_FUN} returns the negative log likelihood of a GLS,
 #' given the parameters in \code{op} and \code{fp}
 #'
-#' @examples
-#' \dontrun{
-#' data(ndvi_AK10000)
-#' df = ndvi_AK10000[seq_len(200), ] # first 500 rows
-#' coords = df[, c("lng", "lat")]
-#' remotePARTS:::fitGLS_opt_FUN(op = c(range = .1, nugget = .2),
-#'                              formula = CLS_coef ~ 0 + land, data = df,
-#'                              coords = coords)
-#' remotePARTS:::fitGLS_opt_FUN(op = c(range = .1), fp = c(nugget = 0),
-#'                              formula = CLS_coef ~ 0 + land, data = df,
-#'                              coords = coords)
-#'
-#' logit <- function(p) {log(p / (1 - p))}
-#' inv_logit <- function(l) {1 / (1 + exp(-l))}
-#'
-#' # input logit-transformed range parameters
-#' remotePARTS:::fitGLS_opt_FUN(op = c(range = .1, nugget = logit(.2)),
-#'                              formula = CLS_coef ~ 0 + land, data = df,
-#'                              coords = coords, is.trans = TRUE,
-#'                              backtrans = list(nugget = inv_logit))
-#' # transformed range and nugget
-#' remotePARTS:::fitGLS_opt_FUN(op = c(range = logit(.1), nugget = logit(.2)),
-#'                              formula = CLS_coef ~ 0 + land, data = df,
-#'                              coords = coords, is.trans = TRUE,
-#'                              backtrans = list(nugget = inv_logit, range = inv_logit))
-#' }
 fitGLS_opt_FUN <- function(op, fp, formula, data = NULL, coords, covar_FUN = "covar_exp", distm_FUN = "distm_scaled",
                            is.trans = FALSE, backtrans = list(), ncores = NA){
   ## combine the optimized and fixed parameters into one vector

R/fitGLS_partitionMC.R

@@ -2,32 +2,8 @@
 #'
 #' @rdname partGLS
 #'
-#' @examples
-#' \dontrun{
-#' ## read data
-#' data(ndvi_AK10000)
-#' df = ndvi_AK10000[seq_len(1000), ] # first 1000 rows
+#' @return a "MC_partGLS", which is a precursor to a "partGLS" object
 #'
-#' ## create partition matrix
-#' pm = sample_partitions(nrow(df), npart = 3)
-#'
-#' ## fit GLS with fixed nugget
-#' MCpartGLS = remotePARTS:::MC_GLSpart(formula = CLS_coef ~ 0 + land, partmat = pm,
-#'                                      data = df, nugget = 0, ncores = 2L)
-#'
-#' (partGLS.mc = remotePARTS:::MCGLS_partsummary(MCpartGLS, partsize = nrow(pm)))
-#'
-#' MCpartGLS2 = remotePARTS:::MC_GLSpart(formula = CLS_coef ~ lat, partmat = pm,
-#'                                      data = df, nugget = 0, ncores = 2L)
-#'
-#' (partGLS2.mc = remotePARTS:::MCGLS_partsummary(MCpartGLS2, partsize = nrow(pm)))
-#'
-#' MCpartGLS_int = remotePARTS:::MC_GLSpart(formula = CLS_coef ~ 1, partmat = pm,
-#'                                          data = df, nugget = 0, ncores = 2L)
-#'
-#' (partGLS_int.mc = remotePARTS:::MCGLS_partsummary(MCpartGLS_int, partsize = nrow(pm)))
-#'
-#' }
 MC_GLSpart <- function(formula, partmat, formula0 = NULL, part_FUN = "part_data",
                        distm_FUN = "distm_scaled", covar_FUN = "covar_exp",
                        covar.pars = c(range = .1), nugget = NA, ncross = 6,

R/optimize_nugget.R

@@ -28,21 +28,6 @@
 #'
 #' @seealso \code{?stats::optimize()}
 #'
-#' @examples
-#' \dontrun{
-#' ## read data
-#' data(ndvi_AK10000)
-#' df = ndvi_AK10000[seq_len(200), ] # first 200 rows
-#'
-#' ## format data
-#' X = stats::model.matrix(CLS_coef ~ 0 + land, data = df)
-#'
-#' ## fit covariance matrix
-#' V = covar_exp(distm_scaled(cbind(df$lng, df$lat)), range = .01)
-#'
-#' ## find the ML nugget
-#' remotePARTS:::optimize_nugget(X = X, V = V, y = df$CLS_coef, debug = TRUE)
-#' }
 optimize_nugget <- function(X, y, V, lower = 0.001, upper = 0.999,
                             tol = .Machine$double.eps^.25, debug = FALSE,
                             ncores = NA) {

R/partGLS_methods.R

@@ -7,6 +7,8 @@
 #' @param ... additional arguments passed to print
 #'
 #' @method print partGLS
+#' @return a print-formatted version of key elements of the "partGLS" object.
+#'
 #' @export
 print.partGLS <- function(x, ...){
   cat("\nCoefficients:\n")
@@ -41,8 +43,6 @@ print.partGLS <- function(x, ...){
 #'
 #' @return a p-value for the correlated chisqr test
 #'
-#' @examples
-#' remotePARTS:::part_chisqr(Fmean = 3.6, rSSR = .021, df1 = 2, npart = 5)
 part_chisqr <- function(Fmean, rSSR, df1, npart){
   checks = c(is.na(Fmean) | is.infinite(Fmean),
              is.na(rSSR) | is.infinite(rSSR),
@@ -75,6 +75,7 @@ part_chisqr <- function(Fmean, rSSR, df1, npart){
 #'
 #' @param x object on which to conduct the test
 #' @param ... additional arguments
+#' @return results of the chi-squared test (generic)
 #'
 #' @export
 chisqr <- function(x, ...){