Merge pull request #10 from USEPA/develop

CRAN v0.1.1

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: SSN2
 Title: Spatial Modeling on Stream Networks
-Version: 0.1.0
+Version: 0.1.1
 Authors@R: c(
     person(given = "Michael",
            family = "Dumelle",

NAMESPACE

@@ -146,7 +146,7 @@ export(randcov_initial)
 export(randcov_params)
 export(ssn_create_distmat)
 export(ssn_get_data)
-export(ssn_get_netgeometry)
+export(ssn_get_netgeom)
 export(ssn_get_stream_distmat)
 export(ssn_glm)
 export(ssn_import)

NEWS.md

@@ -1,3 +1,23 @@
+# SSN2 0.1.1
+
+## Minor Updates
+
+* Changed network geometry name from `netgeometry` to `netgeom` to avoid exceeding the 10 character limit for column/field names while writing to shapefiles [(#2)](https://github.com/USEPA/SSN2/issues/2).
+* Added an error message when `family` is missing in `ssn_glm()` [(#8)](https://github.com/USEPA/SSN2/issues/8).
+* Added a deprecation warning for `SSN_to_SSN2()`.
+* Minor stability updates.
+* Minor error message updates.
+* Minor documentation updates.
+
+## Bug Fixes
+
+* Fixed a bug in `Torgegram()` that prevented intended computation when `cutoff` was specified.
+* Fixed a bug in `plot.Torgegram()` that occasionally prevented proper spacing of the legend.
+* Fixed a bug that prevented proper printing of the dispersion parameter from `ssn_glm()` model objects (and their summaries) when all covariance parameters were known.
+* Fixed a bug that prevented simulation when `euclid_type` was `"none"`.
+* Fixed a bug that could cause improper prediction behavior when `taildown_type` was `"spherical"`.
+* Fixed a bug that printed response residuals instead of deviance residuals for `ssn_glm()` objects.
+
 # SSN2 0.1.0
 
 * Initial CRAN submission.

R/SSN_to_SSN2.R

@@ -45,6 +45,9 @@
 #' @name SSN_to_SSN2
 #' @export
 SSN_to_SSN2 <- function(object, edge_additive = NULL, site_additive = NULL) {
+
+  .Deprecated(new = "ssn_import")
+
   if (!requireNamespace("sp", quietly = TRUE)) {
     stop("Install the sp package before using SSN_to_SSN2", call. = FALSE)
   } else {
@@ -53,7 +56,7 @@ SSN_to_SSN2 <- function(object, edge_additive = NULL, site_additive = NULL) {
     }
 
     ## ---------------------------------------------------
-    ## Convert edges to sf and add netgeometry column
+    ## Convert edges to sf and add netgeom column
     ## ---------------------------------------------------
 
     sl <- sp::SpatialLines(object@lines, proj4string = object@proj4string)
@@ -62,7 +65,7 @@ SSN_to_SSN2 <- function(object, edge_additive = NULL, site_additive = NULL) {
 
     nl.coords <- object@network.line.coords
 
-    edges$netgeometry <- paste("ENETWORK",
+    edges$netgeom <- paste("ENETWORK",
       paste("(",
         paste(
           nl.coords$NetworkID,
@@ -100,7 +103,7 @@ SSN_to_SSN2 <- function(object, edge_additive = NULL, site_additive = NULL) {
 
 
     ## ------------------------------------------------
-    ## Convert observed sites to sf and add netgeometry
+    ## Convert observed sites to sf and add netgeom
     ## ------------------------------------------------
     ## sites<- st_as_sf(sp::SpatialPointsDataFrame(object@obspoints@SSNPoints[[1]]@point.coords,
     ##       object@obspoints@SSNPoints[[1]]@point.data,
@@ -116,7 +119,7 @@ SSN_to_SSN2 <- function(object, edge_additive = NULL, site_additive = NULL) {
     np.coords <- object@obspoints@SSNPoints[[1]]@network.point.coords
     np.coords <- cbind(np.coords, sites[, c("ratio", "locID")])
     np.coords$pid <- rownames(object@obspoints@SSNPoints[[1]]@network.point.coords)
-    sites$netgeometry <- paste(
+    sites$netgeom <- paste(
       "SNETWORK",
       paste(
         "(",
@@ -161,7 +164,7 @@ SSN_to_SSN2 <- function(object, edge_additive = NULL, site_additive = NULL) {
 
     ## ------------------------------------------------------------- If
     ## prediction sites are present, convert to list of sf data.frames
-    ## and add netgeometry column
+    ## and add netgeom column
     ## -------------------------------------------------------------
 
     if (length(object@predpoints@ID) > 0) {
@@ -184,7 +187,7 @@ SSN_to_SSN2 <- function(object, edge_additive = NULL, site_additive = NULL) {
           tmp.sf[, c("ratio", "locID")]
         )
         np.coords$pid <- rownames(object@predpoints@SSNPoints[[i]]@network.point.coords)
-        tmp.sf$netgeometry <- paste(
+        tmp.sf$netgeom <- paste(
           "SNETWORK",
           paste(
             "(",

R/Torgegram.R

@@ -109,7 +109,7 @@ Torgegram <- function(formula, ssn.object,
 
   # find cutoffs
   if (missing(cutoff) || is.null(cutoff)) {
-    cutoff <- 0.5
+    cutoff <- NULL
   }
 
   if (missing(type) || is.null(type)) {
@@ -136,12 +136,14 @@ Torgegram <- function(formula, ssn.object,
 }
 
 get_esv <- function(dist_vector, resid2_vector, bins, cutoff) {
-  cutoff_val <- max(dist_vector) * cutoff
-  index <- dist_vector > 0 & dist_vector <= cutoff_val
+  if (is.null(cutoff)) {
+    cutoff <- max(dist_vector) * 0.5
+  }
+  index <- dist_vector > 0 & dist_vector <= cutoff
   dist_vector <- dist_vector[index]
   resid2 <- resid2_vector[index]
 
-  dist_classes <- cut(dist_vector, breaks = seq(0, cutoff_val, length.out = bins + 1))
+  dist_classes <- cut(dist_vector, breaks = seq(0, cutoff, length.out = bins + 1))
 
   # compute squared differences within each class
   gamma <- tapply(resid2, dist_classes, function(x) mean(x) / 2)

R/amongSitesDistMat.R

@@ -13,7 +13,7 @@ amongSitesDistMat <- function(ssn, pids, name = "obs", bin.table) {
   if (name != "obs") {
     ind.pids <- ssn$preds[[name]]$ng.pid %in% as.character(pids)
     locID.pid.data <- ssn$preds[[name]]$locID[ind.pids]
-    pid.data <- ssn_get_netgeometry(ssn$preds[[name]][ind.pids, ], c(
+    pid.data <- ssn_get_netgeom(ssn$preds[[name]][ind.pids, ], c(
       "pid", "SegmentID", "locID",
       "DistanceUpstream"
     ))
@@ -22,7 +22,7 @@ amongSitesDistMat <- function(ssn, pids, name = "obs", bin.table) {
   } else {
     ind.pids <- ssn$obs$ng.pid %in% as.character(pids)
     locID.pid.data <- ssn$obs$locID[ind.pids]
-    pid.data <- ssn_get_netgeometry(ssn$obs[ind.pids, ], c(
+    pid.data <- ssn_get_netgeom(ssn$obs[ind.pids, ], c(
       "pid", "SegmentID", "locID",
       "DistanceUpstream"
       ), reformat = TRUE)

R/augment.R

@@ -111,7 +111,7 @@ augment.ssn_lm <- function(x, drop = TRUE, newdata = NULL, se_fit = FALSE,
       augment_data$.se.fit <- preds_data$se.fit
     }
     tibble_out <- tibble::tibble(cbind(data, augment_data, influence(x)))
-    tibble_out$pid <- ssn_get_netgeometry(x$ssn.object$obs, netvars = "pid")$pid
+    tibble_out$pid <- ssn_get_netgeom(x$ssn.object$obs, netvars = "pid")$pid
     coords <- sf::st_coordinates(x$ssn.object$obs)
     tibble_out$.xcoord <- coords[, 1, drop = TRUE]
     tibble_out$.ycoord <- coords[, 2, drop = TRUE]
@@ -158,7 +158,7 @@ augment.ssn_lm <- function(x, drop = TRUE, newdata = NULL, se_fit = FALSE,
       }
       coords <- sf::st_coordinates(newdata)
       tibble_out <- tibble::tibble(cbind(sf::st_drop_geometry(newdata), augment_newdata))
-      augment_newdata$pid <- ssn_get_netgeometry(newdata, netvars = "pid")$pid
+      augment_newdata$pid <- ssn_get_netgeom(newdata, netvars = "pid")$pid
       tibble_out$.xcoord <- coords[, 1, drop = TRUE]
       tibble_out$.ycoord <- coords[, 2, drop = TRUE]
       tibble_out <- sf::st_as_sf(tibble_out, coords = c(".xcoord", ".ycoord"), crs = x$crs)

R/cov_initial_search_glm.R

@@ -4,10 +4,11 @@ cov_initial_search_glm <- function(initial_NA_object, ssn.object, data_object, e
   ns2 <- 1.2 * s2
 
   # create a grid of initial values for the covariance parameters
-  tailup_de <- ns2 * c(rep(9 / 10, 2), rep(1 / 3 * 1 / 10, 2), rep(1 / 3 * 1 / 10, 2), 1 / 3 * 1 / 10, rep(1 / 4, 2))
-  taildown_de <- ns2 * c(rep(1 / 3 * 1 / 10, 2), rep(9 / 10, 2), rep(1 / 3 * 1 / 10, 2), 1 / 3 * 1 / 10, rep(1 / 4, 2))
-  euclid_de <- ns2 * c(rep(1 / 3 * 1 / 10, 2), rep(1 / 3 * 1 / 10, 2), rep(9 / 10, 2), 1 / 3 * 1 / 10, rep(1 / 4, 2))
-  nugget <- ns2 * c(rep(1 / 3 * 1 / 10, 2), rep(1 / 3 * 1 / 10, 2), rep(1 / 3 * 1 / 10, 2), 9 / 10, rep(1 / 4, 2))
+  var_min <- 0.05
+  tailup_de <- pmax(ns2 * c(rep(9 / 10, 2), rep(1 / 3 * 1 / 10, 2), rep(1 / 3 * 1 / 10, 2), 1 / 3 * 1 / 10, rep(1 / 4, 2)), var_min)
+  taildown_de <- pmax(ns2 * c(rep(1 / 3 * 1 / 10, 2), rep(9 / 10, 2), rep(1 / 3 * 1 / 10, 2), 1 / 3 * 1 / 10, rep(1 / 4, 2)), var_min)
+  euclid_de <- pmax(ns2 * c(rep(1 / 3 * 1 / 10, 2), rep(1 / 3 * 1 / 10, 2), rep(9 / 10, 2), 1 / 3 * 1 / 10, rep(1 / 4, 2)), var_min)
+  nugget <- pmax(ns2 * c(rep(1 / 3 * 1 / 10, 2), rep(1 / 3 * 1 / 10, 2), rep(1 / 3 * 1 / 10, 2), 9 / 10, rep(1 / 4, 2)), var_min)
 
   # find the maximum tail and euclidean distances to consider
   tail_max <- data_object$tail_max
@@ -26,6 +27,10 @@ cov_initial_search_glm <- function(initial_NA_object, ssn.object, data_object, e
     euclid_range = euclid_range, rotate = 0, scale = 1, dispersion = 1
   )
 
+  cov_grid2 <- cov_grid
+  cov_grid2$dispersion <- 100
+  cov_grid <- rbind(cov_grid, cov_grid2)
+
   # if there is anisotropy test midpoints
   if (data_object$anisotropy) {
     cov_grid2 <- cov_grid
@@ -79,7 +84,7 @@ cov_initial_search_glm <- function(initial_NA_object, ssn.object, data_object, e
   }
 
   # find the unique rows of the grid
-  cov_grid <- unique(cov_grid_replace(cov_grid, initial_NA_object, data_object))
+  cov_grid <- unique(cov_grid_replace_glm(cov_grid, initial_NA_object, data_object))
   # split into list
   cov_grid_splits <- split(cov_grid, seq_len(NROW(cov_grid)))
   # iterate through list

R/cov_vector.R

@@ -94,11 +94,12 @@ cov_vector.taildown_spherical <- function(params, dist_pred_object, ...) {
   flow_con <- b == 0
   dist_ratio_h <- h / params[["range"]]
   dist_ratio_a <- a / params[["range"]]
+  dist_ratio_b <- b / params[["range"]]
   dist_ratio_h_less1 <- dist_ratio_h <= 1
   dist_ratio_a_less1 <- dist_ratio_a <= 1
 
   h_cor_part <- (1 - 1.5 * dist_ratio_h + 0.5 * dist_ratio_h^3) * dist_ratio_h_less1 * flow_con
-  a_cor_part <- (1 - 1.5 * dist_ratio_a + 0.5 * dist_ratio_a) * (1 - dist_ratio_a)^2 * dist_ratio_a_less1 * (!flow_con)
+  a_cor_part <- (1 - 1.5 * dist_ratio_b + 0.5 * dist_ratio_a) * (1 - dist_ratio_a)^2 * dist_ratio_a_less1 * (!flow_con)
   m * params[["de"]] * (h_cor_part + a_cor_part)
 }
 

R/create_netgeometry.R

@@ -1,12 +1,12 @@
-## Add netgeometry column to sf data.frame sf_data
-create_netgeometry <- function(sf_data, type = NULL) {
+## Add netgeom column to sf data.frame sf_data
+create_netgeom <- function(sf_data, type = NULL) {
   if (type == "point") {
-    sf_data[, "netgeometry"] <- paste0("SNETWORK (", paste(
+    sf_data[, "netgeom"] <- paste0("SNETWORK (", paste(
       sf_data$netID, sf_data$rid, sf_data$upDist,
       sf_data$ratio, sf_data$pid, sf_data$locID
     ), ")", sep = "")
   } else {
-    sf_data[, "netgeometry"] <- paste0("ENETWORK (", paste(
+    sf_data[, "netgeom"] <- paste0("ENETWORK (", paste(
       sf_data$netID,
       sf_data$rid,
       sf_data$upDist
@@ -15,5 +15,5 @@ create_netgeometry <- function(sf_data, type = NULL) {
     sep = ""
     )
   }
-  return(sf_data) ## Return sf data.frame with netgeometry column added
+  return(sf_data) ## Return sf data.frame with netgeom column added
 }

R/get_data_object.R

@@ -184,7 +184,7 @@ get_data_object <- function(formula, ssn.object, additive, anisotropy,
   order <- unlist(split(seq_len(n), local$index), use.names = FALSE)
 
   # store global pid
-  pid <- ssn_get_netgeometry(ssn.object$obs, "pid")$pid
+  pid <- ssn_get_netgeom(ssn.object$obs, "pid")$pid
 
   # restructure ssn
   ssn.object <- restruct_ssn_missing(ssn.object, observed_index, missing_index)

R/get_data_object_glm.R

@@ -212,7 +212,7 @@ get_data_object_glm <- function(formula, ssn.object, family, additive, anisotrop
   order <- unlist(split(seq_len(n), local$index), use.names = FALSE)
 
   # store global pid
-  pid <- ssn_get_netgeometry(ssn.object$obs, "pid")$pid
+  pid <- ssn_get_netgeom(ssn.object$obs, "pid")$pid
 
   # restructure ssn
   ssn.object <- restruct_ssn_missing(ssn.object, observed_index, missing_index)

R/get_dist_object.R

@@ -1,13 +1,13 @@
 # parent function to get the distance object
 get_dist_object <- function(ssn.object, initial_object, additive, anisotropy) {
-  # get netgeometry
-  netgeometry <- ssn_get_netgeometry(ssn.object$obs, reformat = TRUE)
+  # get netgeom
+  netgeom <- ssn_get_netgeom(ssn.object$obs, reformat = TRUE)
 
   # get network index
-  network_index <- netgeometry$NetworkID
+  network_index <- netgeom$NetworkID
 
   # get pid
-  pid <- netgeometry$pid # not needed now but can reorder by it later
+  pid <- netgeom$pid # not needed now but can reorder by it later
 
   # distance order
   dist_order <- order(network_index, pid)

R/get_dist_pred_object.R

@@ -1,28 +1,28 @@
 # parent function to get the distance object
 get_dist_pred_object <- function(object, newdata_name, initial_object) {
-  # get netgeometry
-  netgeometry <- ssn_get_netgeometry(object$ssn.object$obs, reformat = TRUE)
+  # get netgeom
+  netgeom <- ssn_get_netgeom(object$ssn.object$obs, reformat = TRUE)
 
   # get network index
-  network_index <- netgeometry$NetworkID
+  network_index <- netgeom$NetworkID
 
   # get pid
-  pid <- netgeometry$pid
+  pid <- netgeom$pid
 
   # distance order
   dist_order <- order(network_index, pid)
 
   # inverse of distance order
   inv_dist_order <- order(dist_order)
 
-  # get netgeometry
-  netgeometry_pred <- ssn_get_netgeometry(object$ssn.object$preds[[newdata_name]], reformat = TRUE)
+  # get netgeom
+  netgeom_pred <- ssn_get_netgeom(object$ssn.object$preds[[newdata_name]], reformat = TRUE)
 
   # get network pred index
-  network_index_pred <- netgeometry_pred$NetworkID
+  network_index_pred <- netgeom_pred$NetworkID
 
   # get pid
-  pid_pred <- netgeometry_pred$pid
+  pid_pred <- netgeom_pred$pid
 
   # distance order
   dist_order_pred <- order(network_index_pred, pid_pred)

R/get_dist_predbk_object.R

@@ -1,13 +1,13 @@
 # parent function to get the distance object
 get_dist_predbk_object <- function(object, newdata_name, initial_object) {
-  # get netgeometry
-  netgeometry <- ssn_get_netgeometry(object$ssn.object$preds[[newdata_name]], reformat = TRUE)
+  # get netgeom
+  netgeom <- ssn_get_netgeom(object$ssn.object$preds[[newdata_name]], reformat = TRUE)
 
   # get network index
-  network_index <- netgeometry$NetworkID
+  network_index <- netgeom$NetworkID
 
   # get pid
-  pid <- netgeometry$pid
+  pid <- netgeom$pid
 
   # distance order
   dist_order <- order(network_index, pid)

R/laploglik_products.R

@@ -111,6 +111,9 @@ get_w_and_H <- function(data_object, dispersion, SigInv_list, SigInv_X, cov_beta
   if (length(SigInv_list) == 1) {
     while (iter < 50 && wdiffmax > 1e-4) {
       iter <- iter + 1
+      # if (family %in% c("binomial", "beta")) {
+      #   w <- pmax(pmin(w, 8), -8)
+      # }
       # compute the d vector
       d <- get_d(family, w, y, size, dispersion)
       # and then the gradient vector
@@ -134,6 +137,10 @@ get_w_and_H <- function(data_object, dispersion, SigInv_list, SigInv_X, cov_beta
       w <- wnew
     }
 
+    # if (family %in% c("binomial", "beta")) {
+    #   w <- pmax(pmin(w, 8), -8)
+    # }
+
     mHldet <- as.numeric(determinant(-H, logarithm = TRUE)$modulus)
     # mHldet <- 2 * sum(log(diag(mH_upchol)))
     w_and_H_list <- list(w = w, H = NULL, mHldet = mHldet)

R/plot.R

@@ -208,6 +208,7 @@ plot.Torgegram <- function(x, type, separate = FALSE, ...) {
   x <- lapply(type, function(y) {
     x_sub <- x[[y]]
     x_sub[["type"]] <- y
+    x_sub <- x_sub[complete.cases(x_sub), , drop = FALSE]
     x_sub
   })
 
@@ -266,7 +267,7 @@ plot.Torgegram <- function(x, type, separate = FALSE, ...) {
     }))
   } else {
     do.call("plot", args = c(list(
-      x = x$dist, y = x$gamma, cex = x$cex, ylim = c(0, max(x$gamma) * 1.6),
+      x = x$dist, y = x$gamma, cex = x$cex, xlim = c(0, max(x$dist)), ylim = c(0, max(x$gamma) * 1.6),
       col = x$col
     ), dotlist))
     legend(0, max(x$gamma) * 1.6, legend = levels(x$type), col = col_key[levels(x$type)], pch = dotlist$pch)

R/print.R

@@ -334,10 +334,10 @@ print.summary.ssn_glm <- function(x,
   cat("\nCall:\n", paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n", sep = "")
 
   # pasting the residual summary
-  cat("\nResiduals:\n")
+  cat("\nDeviance Residuals:\n")
   resQ <- c(
-    min(x$residuals$response), quantile(x$residuals$response, p = c(0.25, 0.5, 0.75), na.rm = TRUE),
-    max(x$residuals$response)
+    min(x$residuals$deviance), quantile(x$residuals$deviance, p = c(0.25, 0.5, 0.75), na.rm = TRUE),
+    max(x$residuals$deviance)
   )
   names(resQ) <- c("Min", "1Q", "Median", "3Q", "Max")
   print(resQ, digits = digits)