Merge joss development to master

.Rbuildignore

@@ -1,9 +1,12 @@
 ^inst/extdata$
-^joss$
+^paper$
 ^data-raw$
 ^LICENSE\.txt$
 ^doc$
 ^Meta$
 ^README\.Rmd$
 ^\.github$
 ^examples$
+vignettes/SpatialGEV-vignette.Rmd.orig
+vignettes/rebuild_vignette.R
+vignettes/vignette-helper-functions.R

.github/workflows/R-CMD-check.yaml

@@ -2,9 +2,9 @@
 # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches: [main, master]
+    branches: [main, master, joss]
   pull_request:
-    branches: [main, master]
+    branches: [main, master, joss]
 
 name: R-CMD-check
 
@@ -22,30 +22,32 @@ jobs:
           - {os: windows-latest, r: 'release'}
           - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
           - {os: ubuntu-latest,   r: 'release'}
-          - {os: ubuntu-latest,   r: 'oldrel-1'}
 
     env:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
       R_KEEP_PKG_SOURCE: yes
 
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
 
-      - uses: r-lib/actions/setup-pandoc@v1
+      - uses: r-lib/actions/setup-pandoc@v2
 
-      - uses: r-lib/actions/setup-r@v1
+      - uses: r-lib/actions/setup-r@v2
         with:
           r-version: ${{ matrix.config.r }}
           http-user-agent: ${{ matrix.config.http-user-agent }}
           use-public-rspm: true
-      
+
       - name: Add inla repo
         run: |
           writeLines('options(repos = c(INLA = "https://inla.r-inla-download.org/R/stable/", getOption(\"repos\")))', file("~/.Rprofile", "ab"))
         shell: Rscript {0}
 
-      - uses: r-lib/actions/setup-r-dependencies@v1
+      - uses: r-lib/actions/setup-r-dependencies@v2
         with:
-          extra-packages: rcmdcheck
+          extra-packages: any::rcmdcheck
+          needs: check
 
-      - uses: r-lib/actions/check-r-package@v1
+      - uses: r-lib/actions/check-r-package@v2
+        with:
+          upload-snapshots: true

.github/workflows/draft-pdf.yml

@@ -0,0 +1,23 @@
+on: [push]
+
+jobs:
+  paper:
+    runs-on: ubuntu-latest
+    name: Paper Draft
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+      - name: Build draft PDF
+        uses: openjournals/openjournals-draft-action@master
+        with:
+          journal: joss
+          # This should be the path to the paper within your repo.
+          paper-path: paper/paper.md
+      - name: Upload
+        uses: actions/upload-artifact@v1
+        with:
+          name: paper
+          # This is the output path where Pandoc will write the compiled
+          # PDF. Note, this should be the same directory as the input
+          # paper.md
+          path: paper/paper.pdf

.github/workflows/test-coverage.yaml

@@ -0,0 +1,35 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/master/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main, master, joss]
+  pull_request:
+    branches: [main, master, joss]
+
+name: test-coverage
+
+jobs:
+  test-coverage:
+    runs-on: ubuntu-latest
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+
+    steps:
+      - uses: actions/checkout@v2
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
+
+      - name: Add inla repo
+        run: |
+          writeLines('options(repos = c(INLA = "https://inla.r-inla-download.org/R/stable/", getOption(\"repos\")))', file("~/.Rprofile", "ab"))
+        shell: Rscript {0}
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: covr
+
+      - name: Test coverage
+        run: covr::codecov(token = "${{ secrets.CODECOV_TOKEN }}")
+        shell: Rscript {0}

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: SpatialGEV
 Title: Fit Spatial Generalized Extreme Value Models
-Version: 1.0.0
+Version: 1.0.1
 Authors@R: c(
 	person(given = "Meixi",
            family = "Chen",
@@ -15,7 +15,7 @@ Authors@R: c(
            role = "ctb",
            email = "rramezan@uwaterloo.ca"))
 Maintainer: Meixi Chen <meixi.chen@uwaterloo.ca>
-Description: Fit latent variable models with the GEV distribution as the data likelihood and the GEV parameters following latent Gaussian processes. The models in this package are built using the template model builder 'TMB' in R, which has the fast ability to integrate out the latent variables using Laplace approximation. This package allows the users to choose in the fit function which GEV parameter(s) is considered as a spatially varying random effect following a Gaussian process, so the users can fit spatial GEV models with different complexities to their dataset without having to write the models in 'TMB' by themselves. This package also offers methods to sample from both fixed and random effects posteriors as well as the posterior predictive distributions at different spatial locations. Methods for fitting this class of models are described in Chen, Ramezan, and Lysy (2021) \href{https://arxiv.org/abs/2110.07051}{arXiv:2110.07051}.
+Description: Fit latent variable models with the GEV distribution as the data likelihood and the GEV parameters following latent Gaussian processes. The models in this package are built using the template model builder 'TMB' in R, which has the fast ability to integrate out the latent variables using Laplace approximation. This package allows the users to choose in the fit function which GEV parameter(s) is considered as a spatially varying random effect following a Gaussian process, so the users can fit spatial GEV models with different complexities to their dataset without having to write the models in 'TMB' by themselves. This package also offers methods to sample from both fixed and random effects posteriors as well as the posterior predictive distributions at different spatial locations. Methods for fitting this class of models are described in Chen, Ramezan, and Lysy (2024) <doi:10.48550/arXiv.2110.07051>.
 License: GPL-3
 Encoding: UTF-8
 LazyData: true
@@ -39,3 +39,5 @@ Suggests:
     rmarkdown
 Additional_repositories: https://inla.r-inla-download.org/R/stable/
 VignetteBuilder: knitr
+URL: https://github.com/meixichen/SpatialGEV
+BugReports: https://github.com/meixichen/SpatialGEV/issues

NAMESPACE

@@ -10,16 +10,18 @@ export(grid_location)
 export(kernel_exp)
 export(kernel_matern)
 export(matern_pc_prior)
-export(r_nll)
 export(sim_cond_normal)
 export(spatialGEV_fit)
 export(spatialGEV_model)
 export(spatialGEV_predict)
 export(spatialGEV_sample)
-import(evd)
-import(mvtnorm)
+importFrom(evd,dgev)
+importFrom(evd,qgev)
+importFrom(evd,rgev)
 importFrom(methods,as)
 importFrom(methods,is)
+importFrom(mvtnorm,dmvnorm)
+importFrom(mvtnorm,rmvnorm)
 importFrom(stats,dist)
 importFrom(stats,nlminb)
 importFrom(stats,quantile)

R/SpatialGEV-package.R

@@ -1,6 +1,6 @@
 #' @keywords internal
-#' @import mvtnorm
-#' @import evd
+#' @importFrom mvtnorm rmvnorm dmvnorm
+#' @importFrom evd dgev rgev qgev
 #' @importFrom stats dist nlminb rnorm runif rexp setNames quantile
 #' @importFrom methods as is
 #' @rawNamespace useDynLib(SpatialGEV, .registration=TRUE);useDynLib(SpatialGEV_TMBExports)

R/old/spatialGEV_fit.R

@@ -262,7 +262,7 @@ spatialGEV_fit <- function(y, locs, random, init_param, reparam_s, kernel="exp",
   mod <- paste("model", random, sep="_")
   mod <- paste(mod, kernel, sep="_")
   n_loc <- length(y)
-  n_obs <- sapply(y, length)
+  n_obs <- vapply(y, length, 1L)
   if (missing(sp_thres)) sp_thres <- -1
   y <- unlist(y)
   loc_ind <- rep(1:n_loc, times=n_obs) # location ind associated with each obs

R/old/spatialGEV_sample.R

@@ -40,15 +40,15 @@ spatialGEV_sample <- function(model, n_draw, observation=FALSE, loc_ind=NULL) {
     meshidxloc <- model$meshidxloc
     if (mod == "a") {
       mean_random <- mean_random[meshidxloc]
-      ind2rm <- setdiff(1:length(random_ind), meshidxloc)
+      ind2rm <- setdiff(seq_along(random_ind), meshidxloc)
       joint_cov <- joint_cov[-ind2rm, -ind2rm]
     } else if (mod == "ab") {
-      ind2rm <- setdiff(1:length(random_ind), c(meshidxloc, length(random_ind)/2 + meshidxloc))
+      ind2rm <- setdiff(seq_along(random_ind), c(meshidxloc, length(random_ind)/2 + meshidxloc))
       mean_random <- mean_random[-ind2rm]
       joint_cov <- joint_cov[-ind2rm, -ind2rm]
     } else {
       # if mod="abs"
-      ind2rm <- setdiff(1:length(random_ind),
+      ind2rm <- setdiff(seq_along(random_ind),
                         c(meshidxloc, # indices of a in the random effects vec
 			  length(random_ind)/3 + meshidxloc, # indices of log_b in the random vec
 			  length(random_ind)/3*2 + meshidxloc)) # indices of s in  the random vec
@@ -65,7 +65,7 @@ spatialGEV_sample <- function(model, n_draw, observation=FALSE, loc_ind=NULL) {
   # Determine the positions of the samples based on location indices to be sampled
   if (is.null(loc_ind)) {
     loc_ind <- 1:n_loc
-    sam_ind <- 1:length(joint_mean)
+    sam_ind <- seq_along(joint_mean)
   } else {
     loc_ind <- sort(loc_ind)
     patterns <- c(paste0("a", loc_ind), paste0("log_b", loc_ind), paste0("s", loc_ind))

R/spatialGEV_fit.R

@@ -9,7 +9,7 @@
 #' @param random Either "a", "ab", or "abs", where `a` indicates the location parameter,
 #' `b` indicates the scale parameter, `s` indicates the shape parameter.  This tells the model
 #' which GEV parameters are considered as random effects.
-#' @param method Either "laplace" or "maxsmooth".
+#' @param method Either "laplace" or "maxsmooth". Default is "laplace". See details.
 #' @param init_param A list of initial parameters. See details.
 #' @param reparam_s A flag indicating whether the shape parameter is "zero", "unconstrained",
 #' constrained to be "negative", or constrained to be "positive". If model "abs" is used,
@@ -85,6 +85,19 @@
 #'
 #' @details
 #' This function adopts Laplace approximation using TMB model to integrate out the random effects.
+#' 
+#' Specifying `method="laplace"` means integrating out the random effects \eqn{u} in the joint likelihood 
+#' via the Laplace approximation: \eqn{p_{\mathrm{LA}}(y \mid \theta) \approx \int p(y, u \mid \theta) \ \mathrm{d}u}.
+#' Then the random effects posterior is constructed via a Normal approximation centered at the Laplace-approximated
+#' marginal likelihood mode with the covariance being the quadrature of it.
+#' If `method="maxsmooth"`, the inference is carried out in two steps. First, the user provide the MLEs
+#' and variance estimates of `a`, `b` and `s` at each location to `data`, which is known as the max step.
+#' The max-step estimates are denoted as \eqn{\hat{u}}, and the likelihood function at each location is approximated
+#' by a Normal distribution at \eqn{\mathcal{N}(\hat{u}, \widehat{Var}(u))}. 
+#' Second, the Laplace approximation is used to integrate out the random effects in the joint likelihood
+#' \eqn{p_{\mathrm{LA}}(\hat{u} \mid \theta) \approx \int p(\hat{u},u \mid \theta) \ \mathrm{d}u}, followed by a Normal
+#' approximation at mode and quadrature of the approximated marginal likelihood \eqn{p_{\mathrm{LA}}(\hat{u} \mid \theta)}. 
+#' This is known as the smooth step.
 #'
 #' The random effects are assumed to follow Gaussian processes with mean 0 and covariance matrix
 #' defined by the chosen kernel function. E.g., using the exponential kernel function:

R/spatialGEV_model.R

@@ -89,13 +89,13 @@ parse_data <- function(data, locs, random,
   if(method == "laplace") {
     y <- data
     if(!is.list(y) ||
-       !all(sapply(y, is.numeric))) {
+       !all(vapply(y, is.numeric, TRUE))) {
       stop("For `method == 'laplace', `data must be a numeric list.")
     } else if(length(y) != n_loc) {
       stop("For `method == 'laplace', must have `length(data) == nrow(locs)`.")
     }
     ## n_loc <- length(y)
-    n_obs <- sapply(y, length)
+    n_obs <- vapply(y, length, 1L)
     y <- unlist(y)
     loc_ind <- rep(1:n_loc, times=n_obs) # location ind associated with each obs
     out <- list(y = y, loc_ind = loc_ind)

R/spatialGEV_predict.R

@@ -118,7 +118,7 @@ spatialGEV_predict <- function(model, locs_new, n_draw, type="response",
       new_a <- a_sim_fun(1) # sample parameter a one time
       if (type == "response") {
         new_y <- t(apply(X = new_a, MARGIN = 1, FUN = function(row) {
-          unlist(Map(rgev, n=1, loc=row, scale=b, shape=s))
+          unlist(Map(evd::rgev, n=1, loc=row, scale=b, shape=s))
         })) # a `1 x n_test` matrix
         pred_y_draws[i, ] <- new_y
       }
@@ -180,7 +180,7 @@ spatialGEV_predict <- function(model, locs_new, n_draw, type="response",
       new_ab <- cbind(new_a, exp(new_logb)) # A `1 x (2*n_test)` matrix constructed by putting the matrix of exp(logb) to the right of the matrix of a
       if (type == "response") {
         new_y <- t(apply(X = new_ab, MARGIN = 1, FUN = function(row) {
-          unlist(Map(rgev, n=1, loc=row[1:n_test],
+          unlist(Map(evd::rgev, n=1, loc=row[1:n_test],
                      scale=row[(n_test+1):length(row)], shape=s))
         })) # a `1 x n_test` matrix
         pred_y_draws[i, ] <- new_y
@@ -265,7 +265,7 @@ spatialGEV_predict <- function(model, locs_new, n_draw, type="response",
       new_abs <- cbind(new_a, exp(new_logb), new_s) # A `1 x (3*n_test)` matrix
       if (type == "response") {
         new_y <- t(apply(X = new_abs, MARGIN = 1, FUN = function(row) {
-          unlist(Map(rgev, n=1, loc=row[1:n_test], scale=row[(n_test+1):(2*n_test)],
+          unlist(Map(evd::rgev, n=1, loc=row[1:n_test], scale=row[(n_test+1):(2*n_test)],
                      shape=row[(2*n_test+1):length(row)]))
         })) # a `1 x n_test` matrix
         pred_y_draws[i, ] <- new_y

R/spatialGEV_sample.R

@@ -32,8 +32,8 @@ spatialGEV_sample <- function(model, n_draw, observation=FALSE, loc_ind=NULL) {
   mean_joint[names(mean_joint) %in% names(mean_fixed)] <- mean_fixed
   # sampling
   prec_joint <- rep$jointPrecision
-  if(!all(sapply(dimnames(prec_joint),
-                 function(x) identical(x, names(mean_joint))))) {
+  if(!all(vapply(dimnames(prec_joint),
+                 function(x) identical(x, names(mean_joint)), TRUE))) {
     stop("Dimension name mismatch between `mean_joint` and `prec_joint`. Please file a bug report.")
   }
   joint_post_draw <- rmvn_prec(n_draw,

R/summary.R

@@ -2,7 +2,7 @@
 #'
 #' @param x Model object of class `spatialGEVfit` returned by `spatialGEV_fit`.
 #' @param ... More arguments for `print`.
-#' @return Information about the fitted model containing number of fixed/random effects, 
+#' @return Information about the fitted model containing number of fixed/random effects,
 #' fitting time, convergence information, etc.
 #' @export
 
@@ -21,7 +21,7 @@ print.spatialGEVfit <- function(x, ...){
   } else {
     cat("The model has not converged and the convergence message output by nlminb is: \n",
         x$fit$message)}
-  
+
   # Kernel info
   cat("The model uses a", x$kernel, "kernel \n")
 
@@ -30,7 +30,7 @@ print.spatialGEVfit <- function(x, ...){
   cat("Number of random effects in the model is", length(x$report$par.random), "\n")
 
   # Hessian info
-  ifelse(x$pdHess_avail, 
+  ifelse(x$pdHess_avail,
 	 mes <- "Hessian matrix is positive definite. Use spatialGEV_sample to obtain posterior samples \n",
          mes <- "Hessian matrix is NOT available or NOT positive definite. spatialGEV_sample and spatialGEV_predict cannot be used \n")
   cat(mes)
@@ -59,19 +59,20 @@ summary.spatialGEVfit <- function(object, ...){
                               do.call(cbind, object$return_levels_sd))
     rl_names <- names(object$return_levels)
     colnames(quantile_summary) <- as.vector(
-      sapply(c("Estimate", "Std.Error"), function(name) paste0(name, rl_names)))
+      vapply(c("Estimate", "Std.Error"), function(name) paste0(name, rl_names),
+	     letters[seq_along(rl_names)]))
     if (object$kernel == "spde") quantile_summary <- quantile_summary[loc_ind,]
     out$return_levels <- quantile_summary
   }
   out
 }
 
 #' Print method for spatialGEVsam
-#' 
+#'
 #' @param x Object of class `spatialGEVsam` returned by `spatialGEV_sample`.
 #' @param ... Additional arguments for `print`.
 #' @return Information about the object including dimension and direction to use `summary` on the object.
-#' @export 
+#' @export
 
 print.spatialGEVsam <- function(x,...){
   # Dimension info
@@ -96,14 +97,14 @@ print.spatialGEVsam <- function(x,...){
 summary.spatialGEVsam <- function(object, q=c(0.025, 0.25, 0.5, 0.75, 0.975), ...){
   # Summary of all parameters
   param_samps <- object[["parameter_draws"]]
-  param_summary <- t(apply(param_samps, 2, quantile, 
+  param_summary <- t(apply(param_samps, 2, quantile,
 	             probs=q))
   param_summary <- cbind(param_summary, apply(param_samps, 2, mean))
   colnames(param_summary)[ncol(param_summary)] <- "mean"
-  # Summary of data resamples 
+  # Summary of data resamples
   if ("y_draws" %in% names(object)){
     y_samps <- object[["y_draws"]]
-    y_summary <- t(apply(y_samps, 2, quantile, 
+    y_summary <- t(apply(y_samps, 2, quantile,
 	           probs=q))
     y_summary <- cbind(y_summary, apply(y_samps, 2, mean))
     colnames(y_summary)[ncol(y_summary)] <- "mean"
@@ -142,7 +143,7 @@ print.spatialGEVpred <- function(x, ...){
 summary.spatialGEVpred <- function(object, q=c(0.025, 0.25, 0.5, 0.75, 0.975), ...){
   # Summary of all parameters
   y_draws <- object[["pred_y_draws"]]
-  out <- t(apply(y_draws, 2, quantile, 
+  out <- t(apply(y_draws, 2, quantile,
 	             probs=q))
   out <- cbind(out, apply(y_draws, 2, mean))
   colnames(out)[ncol(out)] <- "mean"