Ben Bolker's suggestions

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: spatPomp
 Type: Package
 Title: Inference for Spatiotemporal Partially Observed Markov Processes
-Version: 0.36.1
-Date: 2024-08-20
+Version: 0.37.0
+Date: 2024-10-25
 Authors@R: c(
   person("Kidus", "Asfaw", email = "kidusasfaw1990@gmail.com", role = c("aut")),
   person("Edward", "Ionides", email = "ionides@umich.edu",role = c("cre","aut")),
@@ -58,6 +58,7 @@ Collate:
     'he10mle.R'
     'ienkf.R'
     'igirf.R'
+    'init.R'
     'iubf.R'
     'loglik.R'
     'lorenz.R'

R/abf.R

@@ -109,7 +109,7 @@ setMethod(
   signature=signature(object="spatPomp"),
   function (object, Nrep, Np, nbhd, params,
     tol = 1.0e-100,
-    ..., verbose=getOption("verbose",FALSE)) {
+    ..., verbose=getOption("spatPomp_verbose",FALSE)) {
 
     ep <- paste0("in ",sQuote("abf"),": ")
 
@@ -188,7 +188,7 @@ setMethod(
   signature=signature(object="abfd_spatPomp"),
   function (object, Nrep, Np, nbhd, tol,
     ...,
-    verbose = getOption("verbose", FALSE)) {
+    verbose = getOption("spatPomp_verbose", FALSE)) {
     if (missing(Np)) Np <- object@Np
     if (missing(Nrep)) Nrep <- object@Nrep
     if (missing(nbhd)) nbhd <- object@nbhd

R/abfir.R

@@ -92,7 +92,7 @@ setMethod(
   signature=signature(object="spatPomp"),
   function (object, Np, Nrep, nbhd,
     Ninter, tol = (1e-100), params, ...,
-    verbose=getOption("verbose",FALSE) ) {
+    verbose=getOption("spatPomp_verbose",FALSE) ) {
     ep <- paste0("in ",sQuote("abfir"),": ")
     ## declare global variable since foreach's u uses non-standard evaluation
     i <- 1

R/bpfilter.R

@@ -127,7 +127,7 @@ setMethod(
 setMethod(
   "bpfilter",
   signature=signature(object="spatPomp"),
-  function (object, Np, block_size, block_list, save_states, filter_traj, ..., verbose=getOption("verbose", FALSE)) {
+  function (object, Np, block_size, block_list, save_states, filter_traj, ..., verbose=getOption("spatPomp_verbose", FALSE)) {
     ep = paste0("in ",sQuote("bpfilter"),": ")
 
     if(missing(save_states)) save_states <- FALSE
@@ -172,7 +172,7 @@ setMethod(
 setMethod(
   "bpfilter",
   signature=signature(object="bpfilterd_spatPomp"),
-  function (object, Np, block_size, block_list, save_states, filter_traj, ..., verbose=getOption("verbose", FALSE)) {
+  function (object, Np, block_size, block_list, save_states, filter_traj, ..., verbose=getOption("spatPomp_verbose", FALSE)) {
     ep = paste0("in ",sQuote("bpfilter"),": ")
 
     if(missing(save_states)) save_states <- FALSE

R/conc.R

@@ -1,4 +1,4 @@
-##' Concatenate
+##' Concatenate spatPomp objects into a listie
 ##'
 ##' @description Internal methods to concatenate objects into useful listie.
 ##' @details Not exported.

R/enkf.R

@@ -95,7 +95,7 @@ setMethod(
   signature=signature(data="spatPomp"),
   function (data,
     Np,
-    ..., verbose = getOption("verbose", FALSE)) {
+    ..., verbose = getOption("spatPomp_verbose", FALSE)) {
     tryCatch(
       enkf.internal(
         data,

R/eunit_measure.R

@@ -1,4 +1,4 @@
-#' eunit_measure
+#' Expectation of the measurement model for one unit
 #'
 #' \code{eunit_measure} evaluates the expectation of a unit's observation given the entire state
 #' @name eunit_measure

R/girf.R

@@ -127,7 +127,7 @@ setMethod(
     kind=c('bootstrap','moment'),
     tol,
     ...,
-    verbose = getOption("verbose", FALSE)) {
+    verbose = getOption("spatPomp_verbose", FALSE)) {
 
     if (missing(tol)) tol <- 1e-100
     if (missing(Ninter)) Ninter <- length(unit_names(object))

R/ibpf.R

@@ -126,7 +126,7 @@ setMethod(
     cooling.type="geometric",
     cooling.fraction.50,
     block_size, block_list,spat_regression,
-    ..., verbose = getOption("verbose", FALSE)
+    ..., verbose = getOption("spatPomp_verbose", FALSE)
   ){
     ep <- paste0("in ",sQuote("ibpf"),": ")
     if (missing(Nbpf)) pStop_(ep, "Nbpf is required")
@@ -189,7 +189,7 @@ setMethod(
     cooling.type="geometric",
     cooling.fraction.50,
     block_size, block_list,spat_regression,
-    ..., verbose = getOption("verbose", FALSE)
+    ..., verbose = getOption("spatPomp_verbose", FALSE)
   ){
     ep <- paste0("in ",sQuote("ibpf"),": ")
     if (!missing(block_list) & !missing(block_size)){
@@ -244,7 +244,7 @@ setMethod(
     cooling.type="geometric",
     cooling.fraction.50,
     block_size, block_list,spat_regression,
-    ..., verbose = getOption("verbose", FALSE)
+    ..., verbose = getOption("spatPomp_verbose", FALSE)
   ) {
     ep <- paste0("in ",sQuote("ibpf"),": ")
     if (!missing(block_list) & !missing(block_size)){

R/ienkf.R

@@ -82,7 +82,7 @@ setMethod(
     Nenkf = 1, rw.sd,
     cooling.type = c("geometric", "hyperbolic"), cooling.fraction.50,
     Np,
-    ..., verbose = getOption("verbose", FALSE)) {
+    ..., verbose = getOption("spatPomp_verbose", FALSE)) {
     tryCatch(
       ienkf.internal(
         data,

R/igirf.R

@@ -98,7 +98,7 @@ setMethod(
   definition=function (data,Ngirf,Np,rw.sd,cooling.type,cooling.fraction.50,
     Ninter,lookahead=1,Nguide,kind=c('bootstrap', 'moment'),
     tol = 1e-100,
-    ..., verbose = getOption("verbose", FALSE)) {
+    ..., verbose = getOption("spatPomp_verbose", FALSE)) {
 
     ep <- paste0("in ", sQuote("igirf") , " : ")
     if(missing(Ninter)) Ninter <- length(unit_names(data))
@@ -129,7 +129,7 @@ setMethod(
   signature=signature(data="igirfd_spatPomp"),
   function (data,Ngirf,Np,rw.sd,cooling.type, cooling.fraction.50, Ninter,
     lookahead,Nguide,kind=c('bootstrap','moment'),tol, ...,
-    verbose = getOption("verbose", FALSE)) {
+    verbose = getOption("spatPomp_verbose", FALSE)) {
     if (missing(Ngirf)) Ngirf <- data@Ngirf
     if (missing(rw.sd)) rw.sd <- data@rw.sd
     if (missing(cooling.type)) cooling.type <- data@cooling.type

R/init.R

@@ -0,0 +1,6 @@
+## set up default verbosity option
+
+.onAttach <- function(libname,pkgname){
+  options(spatPomp_verbose=FALSE)
+}
+

R/iubf.R

@@ -86,7 +86,7 @@ setMethod(
     rw.sd,
     cooling.type = c("geometric","hyperbolic"),
     cooling.fraction.50, tol = (1e-18)^17,
-    verbose = getOption("verbose"),...) {
+    verbose = getOption("spatPomp_verbose"),...) {
 
     ep <- paste0("in ",sQuote("iubf"),": ")
     if(missing(Nubf))

R/listie.R

@@ -1,6 +1,7 @@
-##' listie
+##' listie: List-like objects
 ##'
-##' List-like objects.
+##' A listie allows object-oriented behavior on lists of
+##' \sQuote{spatPomp} objects
 ##'
 ##' @name listie
 ##' @rdname listie

R/loglik.R

@@ -1,6 +1,6 @@
-##' Log likelihood
+##' Log likelihood extractor 
 ##'
-##' Extract the estimated log likelihood.
+##' Extracts the estimated log likelihood from a fitted model
 ##'
 ##' @name logLik
 ##' @rdname loglik

R/lorenz.R

@@ -1,4 +1,4 @@
-#' Lorenz '96 spatPomp simulator
+#' Lorenz '96 spatPomp constructor
 #'
 #' Generate a spatPomp object representing a \code{U}-dimensional stochastic Lorenz '96 process with
 #' \code{N} measurements made at times \eqn{t_n = n * delta_obs}, simulated using an Euler method

R/munit_measure.R

@@ -1,4 +1,4 @@
-#' munit_measure
+#' Matching moments for the unit measurement model
 #'
 #' \code{munit_measure} returns a moment-matched parameter set given an empirically calculated measurement variance and latent states.
 #' This is used in \code{girf()} and \code{igirf()} when they are run with \code{kind='moment'}.

R/plot.R

@@ -1,69 +1,46 @@
-##' Plotting \code{spatPomp} data
+##' Plot methods for \code{spatPomp} objects
+##'
+##' Visualize data in a \code{spatPomp} object or a derived class.
+##' This gives a quick view; the data can be extracted from 
+##' the object to make a customized plot.
 ##'
-##' Visualize \code{spatPomp} data
 ##' @name plot
 ##' @rdname plot
 ##' @include spatPomp_class.R
 ##' @param x a \code{spatPomp} object
-##' @param log should the data be log-transformed before plotting?
+##' @param log should the data be transformed to \code{log10(x+1)} before plotting?
 ##' This helps in contexts where there are spikes that could take away
-##' attention from the dynamics illustrated by the rest of the data.
+##' attention from the dynamics illustrated by the rest of the data. 
 ##' @param ncol the number of columns in the grid plot
 ##' @param type for visualizing an object of class \code{spatPomp}, the user
-##' can obtain a grid of line plots by default (\code{'l'}) or a heat map by
-##' supplying argument \code{'h'}.
+##' can obtain a grid of line plots by default (\code{l}) or a heat map by
+##' supplying argument \code{h}.
 ##' @param plot_unit_names allows suppression of unit names when making
 ##' a heat map for a large number of units
+##' @param params allows selection of a subset of parameters when making
+##' a diagnostic plot for a model with many parameters
 ##' @param ...  for visualizing an object of class \code{spatPomp}, the user
 ##' can add arguments like \code{nrow} to specify the number of rows in the
 ##' grid.
+##' @return a \code{ggplot} plot of class \sQuote{gg} and \sQuote{ggplot}.
+##' @importFrom ggplot2 ggplot geom_line aes facet_wrap
+##' @importFrom graphics par
 NULL
 
 if(getRversion() >= "2.15.1")  utils::globalVariables(c("."))
-##' Plotting output of \code{igirf()}
-##'
-##' Diagnostic plot for \code{igirf()}
-##' @param params the names of the parameters for which the user would like to see a trace plot
-##' @name plot-igirfd_spatPomp
-##' @rdname plot
-##' @importFrom ggplot2 ggplot geom_line aes facet_wrap
-##' @importFrom graphics par
-##' @aliases plot,igirfd_spatPomp-method
-##' @return a \code{ggplot} facet plot of class \sQuote{gg} and \sQuote{ggplot} visualizing
-##' the convergence record of running \code{igirf()} with respect to the likelihood and the parameters of the model.
-setMethod(
-  "plot",
-  signature=signature(x="igirfd_spatPomp"),
-  definition=function (x, params = names(coef(x)), ncol = 3) {
-    plot.df <- data.frame(x@traces[,c("loglik", params)])
-    cn <- colnames(plot.df)
-    plot.df <- cbind(c(seq_len(dim(plot.df)[1])), plot.df)
-    names(plot.df) <- c("iteration", cn)
-    to.gather <- colnames(plot.df)[2:length(colnames(plot.df))]
-    to.plot <- plot.df |> tidyr::gather(key = "param", val = "value", tidyr::all_of(to.gather)) |> tail(-1)
-    ggplot2::ggplot(data = to.plot) +
-      ggplot2::geom_line(mapping = ggplot2::aes(x = .data$iteration, y = .data$value)) +
-      ggplot2::facet_wrap(~param, ncol = ncol, scales = "free")
-  }
-)
 
-##' Plotting \code{spatPomp} data
-##'
-##' Visualize \code{spatPomp} data
 ##' @name plot-spatPomp
 ##' @rdname plot
 ##' @aliases plot,spatPomp-method
-##' @return a \code{ggplot} plot of class \sQuote{gg} and \sQuote{ggplot} visualizing
-##' the time series data over multiple spatial units via a tile-plot.
 ##' @export
 setMethod(
   "plot",
   signature=signature(x="spatPomp"),
-  definition=function (x, type = c('l','h'), log=F, plot_unit_names=T,...) {
+  definition=function (x, type = c('l','h'), log=FALSE, plot_unit_names=TRUE,...) {
     df <- as.data.frame(x)
-    if(log) df[x@unit_obsnames] <- log10(df[x@unit_obsnames]+1)
-    type = match.arg(type)
+    type <- match.arg(type)
     if(type == 'l'){
+      if(log) df[x@unit_obsnames] <- df[x@unit_obsnames]+1
       unit_nm <- rlang::sym(x@unitname)
       df[[unit_nm]] <- factor(df[[unit_nm]], levels = x@unit_names)
       g <- ggplot2::ggplot(data = df,
@@ -72,11 +49,14 @@ setMethod(
           y = !!rlang::sym(x@unit_obsnames)
 	)
       ) +
+      ggplot2::labs(y=paste0(x@unit_obsnames,"+1")) +
+      ggplot2::scale_y_continuous(transform=if(log) "log10" else "identity") +
       ggplot2::geom_line() +
       ggplot2::facet_wrap(unit_nm, ...)
       return(g)
     }
     if(type == 'h'){
+      if(log) df[x@unit_obsnames] <- log10(df[x@unit_obsnames]+1)
       g <- ggplot2::ggplot(data = df,
         mapping = ggplot2::aes(x = !!rlang::sym(x@timename),
           y = factor(!!rlang::sym(x@unitname),
@@ -104,3 +84,22 @@ setMethod(
     }
   }
 )
+
+##' @name plot-igirfd_spatPomp
+##' @rdname plot
+##' @aliases plot,igirfd_spatPomp-method
+setMethod(
+  "plot",
+  signature=signature(x="igirfd_spatPomp"),
+  definition=function (x, params = names(coef(x)), ncol = 3) {
+    plot.df <- data.frame(x@traces[,c("loglik", params)])
+    cn <- colnames(plot.df)
+    plot.df <- cbind(c(seq_len(dim(plot.df)[1])), plot.df)
+    names(plot.df) <- c("iteration", cn)
+    to.gather <- colnames(plot.df)[2:length(colnames(plot.df))]
+    to.plot <- plot.df |> tidyr::gather(key = "param", val = "value", tidyr::all_of(to.gather)) |> tail(-1)
+    ggplot2::ggplot(data = to.plot) +
+      ggplot2::geom_line(mapping = ggplot2::aes(x = .data$iteration, y = .data$value)) +
+      ggplot2::facet_wrap(~param, ncol = ncol, scales = "free")
+  }
+)

R/runit_measure.R

@@ -1,4 +1,4 @@
-#' runit_measure
+#' Random draw from the measurement model for one unit
 #'
 #' \code{runit_measure} simulates a unit's observation given the entire state
 #' @name runit_measure

R/spatPomp.R

@@ -93,7 +93,7 @@ spatPomp <- function (data, units, times, covar, t0, ...,
   rprocess, rmeasure, dprocess, dmeasure, skeleton, rinit, rprior, dprior,
   unit_statenames, unit_accumvars, shared_covarnames, globals, paramnames, params,
   cdir,cfile, shlib.args, PACKAGE,
-  partrans, compile=TRUE, verbose = getOption("verbose",FALSE)) {
+  partrans, compile=TRUE, verbose = getOption("spatPomp_verbose",FALSE)) {
 
   ep <- paste0("in ",sQuote("spatPomp"),": ")
 

R/vec_dmeasure.R

@@ -1,4 +1,4 @@
-##' Vector of measurement densities
+##' Vector of unit measurement densities for each unit
 ##'
 ##' Evaluate the unit measurement model density function for each unit.
 ##' This method is used primarily as part of likelihood evaluation and parameter inference algorithms.

R/vec_rmeasure.R

@@ -1,4 +1,4 @@
-##' Vector of simulated measurements
+##' Vector simulating measurements for each unit using \code{runit_measure}
 ##'
 ##' Simulate from the unit measurement model density function for each unit
 ##' @param object a \code{spatPomp} object

R/vunit_measure.R

@@ -1,4 +1,4 @@
-#' vunit_measure
+#' Conditional variance of the measurement on a single unit
 #'
 #' \code{vunit_measure} evaluates the variance of a unit's observation given the entire state
 #' @name vunit_measure