diff --git a/NAMESPACE b/NAMESPACE index fd701aa..be4ce1d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -38,6 +38,7 @@ S3method(cov_initial_search,cubic) S3method(cov_initial_search,exponential) S3method(cov_initial_search,gaussian) S3method(cov_initial_search,gravity) +S3method(cov_initial_search,ie) S3method(cov_initial_search,jbessel) S3method(cov_initial_search,magnetic) S3method(cov_initial_search,matern) @@ -57,6 +58,7 @@ S3method(cov_initial_search_glm,cubic) S3method(cov_initial_search_glm,exponential) S3method(cov_initial_search_glm,gaussian) S3method(cov_initial_search_glm,gravity) +S3method(cov_initial_search_glm,ie) S3method(cov_initial_search_glm,jbessel) S3method(cov_initial_search_glm,magnetic) S3method(cov_initial_search_glm,matern) @@ -95,6 +97,7 @@ S3method(get_initial_range,cubic) S3method(get_initial_range,exponential) S3method(get_initial_range,gaussian) S3method(get_initial_range,gravity) +S3method(get_initial_range,ie) S3method(get_initial_range,jbessel) S3method(get_initial_range,magnetic) S3method(get_initial_range,matern) @@ -125,6 +128,7 @@ S3method(gloglik_products,cubic) S3method(gloglik_products,exponential) S3method(gloglik_products,gaussian) S3method(gloglik_products,gravity) +S3method(gloglik_products,ie) S3method(gloglik_products,jbessel) S3method(gloglik_products,magnetic) S3method(gloglik_products,matern) @@ -156,6 +160,7 @@ S3method(laploglik_products,cubic) S3method(laploglik_products,exponential) S3method(laploglik_products,gaussian) S3method(laploglik_products,gravity) +S3method(laploglik_products,ie) S3method(laploglik_products,jbessel) S3method(laploglik_products,magnetic) S3method(laploglik_products,matern) @@ -236,6 +241,7 @@ S3method(spcov_matrix,cubic) S3method(spcov_matrix,exponential) S3method(spcov_matrix,gaussian) S3method(spcov_matrix,gravity) +S3method(spcov_matrix,ie) S3method(spcov_matrix,jbessel) S3method(spcov_matrix,magnetic) S3method(spcov_matrix,matern) @@ -257,6 +263,7 @@ S3method(spcov_optim2orig,cubic) S3method(spcov_optim2orig,exponential) S3method(spcov_optim2orig,gaussian) S3method(spcov_optim2orig,gravity) +S3method(spcov_optim2orig,ie) S3method(spcov_optim2orig,jbessel) S3method(spcov_optim2orig,magnetic) S3method(spcov_optim2orig,matern) @@ -276,6 +283,7 @@ S3method(spcov_orig2optim,cubic) S3method(spcov_orig2optim,exponential) S3method(spcov_orig2optim,gaussian) S3method(spcov_orig2optim,gravity) +S3method(spcov_orig2optim,ie) S3method(spcov_orig2optim,jbessel) S3method(spcov_orig2optim,magnetic) S3method(spcov_orig2optim,matern) @@ -294,6 +302,7 @@ S3method(spcov_vector,cubic) S3method(spcov_vector,exponential) S3method(spcov_vector,gaussian) S3method(spcov_vector,gravity) +S3method(spcov_vector,ie) S3method(spcov_vector,jbessel) S3method(spcov_vector,magnetic) S3method(spcov_vector,matern) @@ -312,6 +321,7 @@ S3method(sprnorm,cubic) S3method(sprnorm,exponential) S3method(sprnorm,gaussian) S3method(sprnorm,gravity) +S3method(sprnorm,ie) S3method(sprnorm,jbessel) S3method(sprnorm,magnetic) S3method(sprnorm,matern) @@ -423,6 +433,7 @@ importFrom(stats,model.frame) importFrom(stats,model.matrix) importFrom(stats,model.offset) importFrom(stats,model.response) +importFrom(stats,na.fail) importFrom(stats,na.omit) importFrom(stats,na.pass) importFrom(stats,pchisq) diff --git a/NEWS.md b/NEWS.md index a76a8be..b793531 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,7 +2,20 @@ ## Major Updates -* Added the `range_constrain` argument to `splm()` and `spglm()` to contrain the range parameter to enhance numerical stability. +* Added the `range_constrain` argument to `splm()` and `spglm()` to constrain the range parameter to enhance numerical stability. The default for `range_constrain` is `FALSE`, implying the range is not constrained. +* Updated the `seal` data with additional polygons and a factor variable, `stock`, with two levels (`8` and `10`) that indicates seal stock (i.e., seal type). + +## Minor Updates + +* Changed diagonal tolerance threshold for `spglm()` and `spgautor()` model objects. See [this link](https://usepa.github.io/spmodel/articles/technical.html#sec:computational) for details. +* Added the `"ie"` spatial covariance type to `splm()` and `spglm()` models. For `splm()` models, `"ie"` is an alias for `"none"`. For `spglm()` models, `"none"` now fixes both the `de` and `ie` covariance parameters at zero, while `"ie"` fixes the `de` covariance parameter at zero but allows the `ie` covariance parameter to vary. Thus, `"none"` from `spmodel $\le$ v0.8.0` matches `"ie"` from `spmodel` v0.9.0 and but is different from `"none"` from `spmodel v0.9.0`. +* Added the `na.action` argument to `predict.spmodel()` functions to clarify that missing values in `newdata` return an error. +* Minor documentation updates. + +## Bug Fixes + +* Fixed a bug that caused incorrect degrees of freedom for the likelihood ratio test (`anova(model1, model2)`) when `estmethod` is `"ml"` for both models. +* Fixed a bug that caused an error in `anova(object1, object2)` when the name of `object1` had special characters (e.g., `$`). # spmodel 0.8.0 @@ -149,7 +162,7 @@ * Fixed a bug in `spautor()` that prevented an error from occurring when a partition factor was not categorical or not a factor * Fixed a bug in `covmatrix(object, newdata)` that returned a matrix with improper dimensions when `spcov_type` was `"none"`. * Fixed a bug in `predict()` that caused an error when at least one level of a fixed effect factor was not observed within a local neighborhood (when the `local` method was `"covariance"` or `"distance")`. -* Fixed a bug in `cooks.distance()` that used the Pearson residuals instead of the standarized residuals. +* Fixed a bug in `cooks.distance()` that used the Pearson residuals instead of the standardized residuals. # spmodel 0.3.0 diff --git a/R/anova.R b/R/anova.R index 28e755c..cf1d467 100644 --- a/R/anova.R +++ b/R/anova.R @@ -81,7 +81,6 @@ #' tidy(anova(spmod, lmod)) anova.splm <- function(object, ..., test = TRUE, Terms, L) { - # see if one or two models object2_list <- list(...) @@ -139,13 +138,19 @@ anova.splm <- function(object, ..., test = TRUE, Terms, L) { stop("The fixed effect coefficients must be the same when performing a likeihood ratio test using the reml estimation method. To perform the likelihood ratio tests for different fixed effect and covariance coefficients simultaneously, refit the models using the ml estimation method.", call. = FALSE) } Chi2_stat <- abs(-2 * (logLik(object2) - logLik(object))) - df_diff <- abs(object2$npar - object$npar) + + # df for ml vs reml + df1 <- object$npar + df2 <- object2$npar + if (object$estmethod == "ml") df1 <- df1 + object$p + if (object2$estmethod == "ml") df2 <- df2 + object2$p + df_diff <- abs(df1 - df2) p_value <- pchisq(Chi2_stat, df_diff, lower.tail = FALSE) if (object2$npar < object$npar) { - full_name <- as.character(substitute(object)) + full_name <- deparse(substitute(object)) # replace as.character with deparse reduced_name <- as.character(as.list(substitute(list(...)))[-1]) } else { - reduced_name <- as.character(substitute(object)) + reduced_name <- deparse(substitute(object)) # replace as.character with deparse full_name <- as.character(as.list(substitute(list(...)))[-1]) } if (test) { diff --git a/R/cov_betahat_adjust.R b/R/cov_betahat_adjust.R index 70cc98e..8debc05 100644 --- a/R/cov_betahat_adjust.R +++ b/R/cov_betahat_adjust.R @@ -18,7 +18,7 @@ cov_betahat_adjust <- function(invcov_betahat_list, betahat_list, randcov_params, cov_betahat_noadjust, var_adjust) { P <- length(betahat_list) # reset var_adjust if only one partition - if (P == 1 || inherits(spcov_params, "none")) { + if (P == 1 || inherits(spcov_params, c("none", "ie"))) { var_adjust <- "none" } # reset var_adjust if partitioning used but no local option used diff --git a/R/cov_estimate_gloglik.R b/R/cov_estimate_gloglik.R index 721393a..dccb175 100644 --- a/R/cov_estimate_gloglik.R +++ b/R/cov_estimate_gloglik.R @@ -19,7 +19,7 @@ cov_estimate_gloglik_splm <- function(data_object, formula, spcov_initial, estme if (data_object$anisotropy) { dist_matrix_list <- NULL } else { - if (inherits(spcov_initial, "none") && is.null(data_object$randcov_initial)) { + if (inherits(spcov_initial, c("none", "ie")) && is.null(data_object$randcov_initial)) { dist_matrix_list <- NULL } else { dist_matrix_list <- lapply(data_object$obdata_list, function(x) spdist(x, data_object$xcoord, data_object$ycoord)) diff --git a/R/cov_initial_search.R b/R/cov_initial_search.R index b38775d..35aedd5 100644 --- a/R/cov_initial_search.R +++ b/R/cov_initial_search.R @@ -585,6 +585,9 @@ cov_initial_search.none <- function(spcov_initial_NA, estmethod, data_object, best_params } +#' @export +cov_initial_search.ie <- cov_initial_search.none + #' @export cov_initial_search.matern <- function(spcov_initial_NA, estmethod, data_object, dist_matrix_list, weights, diff --git a/R/cov_initial_search_glm.R b/R/cov_initial_search_glm.R index 93fc645..76bc6b8 100644 --- a/R/cov_initial_search_glm.R +++ b/R/cov_initial_search_glm.R @@ -461,6 +461,9 @@ cov_initial_search_glm.none <- function(spcov_initial_NA, dispersion_initial_NA, best_params } +#' @export +cov_initial_search_glm.ie <- cov_initial_search_glm.none + #' @export cov_initial_search_glm.matern <- function(spcov_initial_NA, dispersion_initial_NA, estmethod, data_object, dist_matrix_list, weights, diff --git a/R/data.R b/R/data.R index 657679c..5ee5143 100644 --- a/R/data.R +++ b/R/data.R @@ -62,9 +62,12 @@ #' #' @description Estimated harbor-seal trends from abundance data in southeast Alaska, USA. #' -#' @format A \code{sf} object with 62 rows and 2 columns: +#' @format A \code{sf} object with 149 rows and 2 columns: #' \itemize{ #' \item log_trend: The log of the estimated harbor-seal trends from abundance data. +#' \item stock: A seal stock factor with two levels: 8 and 10. The factor levels indicate the +#' type of seal stock (i.e., type of seal). Stocks 8 and 10 are two distinct stocks +#' (out of 13 total stocks) in southeast Alaska. #' \item geometry: \code{POLYGON} geometry representing polygons in an Alaska #' Albers projection (EPSG: 3338). #' } @@ -124,7 +127,7 @@ #' \item strat: A factor representing strata (used for sampling). Can take values \code{L} and \code{M}. #' \item count: The count (number) of moose observed. #' \item presence: A binary factor representing whether no moose were observed (value \code{0}) or at least one moose was observed -#' (va ue \code{1}). +#' (value \code{1}). #' \item geometry: \code{POINT} geometry representing coordinates in an Alaska #' Albers projection (EPSG: 3338). Distances between points are in meters. #' } diff --git a/R/get_data_object.R b/R/get_data_object.R index f3bad45..f19f5b8 100644 --- a/R/get_data_object.R +++ b/R/get_data_object.R @@ -31,7 +31,7 @@ get_data_object_splm <- function(formula, data, spcov_initial, xcoord, ycoord, e data_sf <- NULL } - if (!is_sf && missing(xcoord) && !inherits(spcov_initial, "none")) { + if (!is_sf && missing(xcoord) && !inherits(spcov_initial, c("none", "ie"))) { stop("The xcoord argument must be specified.", call. = FALSE) } @@ -52,7 +52,7 @@ get_data_object_splm <- function(formula, data, spcov_initial, xcoord, ycoord, e ycoord_orig_name <- NULL ycoord_orig_val <- NULL # find coordinate dimension and set defaults - if (inherits(spcov_initial, "none") && estmethod %in% c("reml", "ml")) { + if (inherits(spcov_initial, c("none", "ie")) && estmethod %in% c("reml", "ml")) { dim_coords <- 0 if (missing(xcoord)) { xcoord <- ".xcoord" @@ -213,7 +213,7 @@ get_data_object_splm <- function(formula, data, spcov_initial, xcoord, ycoord, e # } # range constrain - max_range_scale <- 5 + max_range_scale <- 4 range_constrain_value <- 2 * max_halfdist * max_range_scale if ("range" %in% names(spcov_initial$is_known)) { if (spcov_initial$is_known[["range"]] || (spcov_initial$initial[["range"]] > range_constrain_value)) { @@ -221,7 +221,7 @@ get_data_object_splm <- function(formula, data, spcov_initial, xcoord, ycoord, e } } - if (inherits(spcov_initial, "none")) { + if (inherits(spcov_initial, c("none", "ie"))) { range_constrain <- FALSE } diff --git a/R/get_data_object_glm.R b/R/get_data_object_glm.R index 223cd0f..bae24a0 100644 --- a/R/get_data_object_glm.R +++ b/R/get_data_object_glm.R @@ -32,7 +32,7 @@ get_data_object_spglm <- function(formula, family, data, spcov_initial, xcoord, data_sf <- NULL } - if (!is_sf && missing(xcoord) && !inherits(spcov_initial, "none")) { + if (!is_sf && missing(xcoord) && !inherits(spcov_initial, c("none", "ie"))) { stop("The xcoord argument must be specified.", call. = FALSE) } @@ -53,7 +53,7 @@ get_data_object_spglm <- function(formula, family, data, spcov_initial, xcoord, ycoord_orig_name <- NULL ycoord_orig_val <- NULL # find coordinate dimension and set defaults - if (inherits(spcov_initial, "none") && estmethod %in% c("reml", "ml")) { + if (inherits(spcov_initial, c("none", "ie")) && estmethod %in% c("reml", "ml")) { dim_coords <- 0 if (missing(xcoord)) { xcoord <- ".xcoord" @@ -222,8 +222,8 @@ get_data_object_spglm <- function(formula, family, data, spcov_initial, xcoord, betahat <- backsolve(R_val, qr.qty(qr_val, y_trans)) resid <- y_trans - X %*% betahat s2 <- sum(resid^2) / (n - p) - # diagtol <- 1e-4 - diagtol <- min(1e-4, 1e-4 * s2) + diagtol <- 1e-4 + # diagtol <- min(1e-4, 1e-4 * s2) @@ -233,7 +233,7 @@ get_data_object_spglm <- function(formula, family, data, spcov_initial, xcoord, max_halfdist <- sqrt((max(x_range) - min(x_range))^2 + (max(y_range) - min(y_range))^2) / 2 # range constrain - max_range_scale <- 5 + max_range_scale <- 4 range_constrain_value <- 2 * max_halfdist * max_range_scale if ("range" %in% names(spcov_initial$is_known)) { if (spcov_initial$is_known[["range"]] || (spcov_initial$initial[["range"]] > range_constrain_value)) { @@ -241,7 +241,7 @@ get_data_object_spglm <- function(formula, family, data, spcov_initial, xcoord, } } - if (inherits(spcov_initial, "none")) { + if (inherits(spcov_initial, c("none", "ie"))) { range_constrain <- FALSE } diff --git a/R/get_initial_range.R b/R/get_initial_range.R index de38f0e..c2f877c 100644 --- a/R/get_initial_range.R +++ b/R/get_initial_range.R @@ -33,6 +33,10 @@ get_initial_range.circular <- function(spcov_type, max_halfdist, ...) { get_initial_range.none <- function(spcov_type, max_halfdist, ...) { Inf } + +#' @export +get_initial_range.ie <- get_initial_range.none + #' @export get_initial_range.cubic <- function(spcov_type, max_halfdist, ...) { max_halfdist diff --git a/R/get_spcov_params.R b/R/get_spcov_params.R index 0cd1239..29c0084 100644 --- a/R/get_spcov_params.R +++ b/R/get_spcov_params.R @@ -7,7 +7,7 @@ #' #' @noRd get_spcov_params <- function(spcov_type, spcov_orig_val) { - if (spcov_type %in% c("exponential", "spherical", "gaussian", "triangular", "circular", "none", "cubic", "pentaspherical", "cosine", "wave", "jbessel", "gravity", "rquad", "magnetic")) { + if (spcov_type %in% c("exponential", "spherical", "gaussian", "triangular", "circular", "none", "ie", "cubic", "pentaspherical", "cosine", "wave", "jbessel", "gravity", "rquad", "magnetic")) { spcov_params_val <- spcov_params( spcov_type = spcov_type, de = spcov_orig_val[["de"]], diff --git a/R/gloglik_products.R b/R/gloglik_products.R index 9b61b20..ab0abb1 100644 --- a/R/gloglik_products.R +++ b/R/gloglik_products.R @@ -85,6 +85,8 @@ gloglik_products.circular <- gloglik_products.exponential #' @export gloglik_products.none <- gloglik_products.exponential #' @export +gloglik_products.ie <- gloglik_products.none +#' @export gloglik_products.cubic <- gloglik_products.exponential #' @export gloglik_products.pentaspherical <- gloglik_products.exponential diff --git a/R/laploglik_products.R b/R/laploglik_products.R index 7e7ebaa..ba23289 100644 --- a/R/laploglik_products.R +++ b/R/laploglik_products.R @@ -121,6 +121,8 @@ laploglik_products.circular <- laploglik_products.exponential #' @export laploglik_products.none <- laploglik_products.exponential #' @export +laploglik_products.ie <- laploglik_products.none +#' @export laploglik_products.cubic <- laploglik_products.exponential #' @export laploglik_products.pentaspherical <- laploglik_products.exponential diff --git a/R/loocv.R b/R/loocv.R index 52bb312..2b3136f 100644 --- a/R/loocv.R +++ b/R/loocv.R @@ -64,7 +64,7 @@ loocv.splm <- function(object, cv_predict = FALSE, se.fit = FALSE, local, ...) { } # iid if relevant otherwise pass - if (inherits(coef(object, type = "spcov"), "none") && is.null(object$random)) { + if (inherits(coef(object, type = "spcov"), c("none", "ie")) && is.null(object$random)) { return(loocv_iid(object, cv_predict, se.fit, local)) } diff --git a/R/predict.R b/R/predict.R index 54b9d28..20ad711 100644 --- a/R/predict.R +++ b/R/predict.R @@ -21,7 +21,13 @@ #' (ignored if \code{scale} is not specified). The default is \code{Inf}. #' @param interval Type of interval calculation. The default is \code{"none"}. #' Other options are \code{"confidence"} (for confidence intervals) and -#' \code{"prediction"} (for prediction intervals). +#' \code{"prediction"} (for prediction intervals). When \code{interval} +#' is \code{"none"} or \code{"prediction"}, predictions are returned (and when +#' requested, their corresponding uncertainties). When \code{interval} +#' is \code{"confidence"}, mean estimates are returned (and when +#' requested, their corresponding uncertainties). This \code{"none"} behavior +#' differs from that of \code{lm()}, as \code{lm()} returns confidence +#' uncertainties (in \code{.$se.fit}). #' @param level Tolerance/confidence level. The default is \code{0.95}. #' @param type The prediction type, either on the response scale, link scale (only for #' \code{spglm()} or \code{spgautor()} model objects), or terms scale. @@ -65,6 +71,8 @@ #' \code{list(size = 100, method = "covariance", parallel = FALSE)}. #' @param terms If \code{type} is \code{"terms"}, the type of terms to be returned, #' specified via either numeric position or name. The default is all terms are included. +#' @param na.action Missing (\code{NA}) values in \code{newdata} will return an error and should +#' be removed before proceeding. #' @param ... Other arguments. Only used for models fit using \code{splmRF()} #' or \code{spautorRF()} where \code{...} indicates other #' arguments to \code{ranger::predict.ranger()}. @@ -114,7 +122,7 @@ #' predict(spmod, sulfate_preds, interval = "prediction") #' augment(spmod, newdata = sulfate_preds, interval = "prediction") predict.splm <- function(object, newdata, se.fit = FALSE, scale = NULL, df = Inf, interval = c("none", "confidence", "prediction"), - level = 0.95, type = c("response", "terms"), local, terms = NULL, ...) { + level = 0.95, type = c("response", "terms"), local, terms = NULL, na.action = na.fail, ...) { # match interval argument so the three display interval <- match.arg(interval) @@ -232,6 +240,12 @@ predict.splm <- function(object, newdata, se.fit = FALSE, scale = NULL, df = Inf attr(newdata_model, "assign") <- attr_assign[keep_cols] attr(newdata_model, "contrasts") <- attr_contrasts + # finding rows w/out NA + ob_predictors <- complete.cases(newdata_model) + if (any(!ob_predictors)) { + stop("Cannot have NA values in predictors.", call. = FALSE) + } + # call terms if needed if (type == "terms") { return(predict_terms(object, newdata_model, se.fit, scale, df, interval, level, add_newdata_rows, terms, ...)) @@ -380,7 +394,7 @@ predict.splm <- function(object, newdata, se.fit = FALSE, scale = NULL, df = Inf # ) cov_matrix_val <- covmatrix(object) # handling closed form of none covariance - if (inherits(spcov_params_val, "none") && is.null(randcov_params_val)) { + if (inherits(spcov_params_val, c("none", "ie")) && is.null(randcov_params_val)) { cov_lowchol <- cov_matrix_val diag(cov_lowchol) <- sqrt(diag(cov_lowchol)) # already diagonal don't need transpose } else { @@ -538,7 +552,7 @@ predict.splm <- function(object, newdata, se.fit = FALSE, scale = NULL, df = Inf #' @order 2 #' @export predict.spautor <- function(object, newdata, se.fit = FALSE, scale = NULL, df = Inf, interval = c("none", "confidence", "prediction"), - level = 0.95, type = c("response", "terms"), local, terms = NULL, ...) { + level = 0.95, type = c("response", "terms"), local, terms = NULL, na.action = na.fail, ...) { # match interval argument so the three display interval <- match.arg(interval) @@ -603,6 +617,15 @@ predict.spautor <- function(object, newdata, se.fit = FALSE, scale = NULL, df = attr(newdata_model, "assign") <- attr_assign[keep_cols] attr(newdata_model, "contrasts") <- attr_contrasts + # finding rows w/out NA + # this isn't really needed, because the error should come on model building + # but someone could accidentally write over their newdata object after fitting + # so it is good to keep for completeness + ob_predictors <- complete.cases(newdata_model) + if (any(!ob_predictors)) { + stop("Cannot have NA values in predictors.", call. = FALSE) + } + # call terms if needed if (type == "terms") { # may want to add add_newdata_rows if we allow spautor prediction for the observed model matrix @@ -906,8 +929,10 @@ get_pred_spautor_parallel <- function(cluster_list, cov_matrix_lowchol, betahat, #' @method predict splm_list #' @order 3 #' @export -predict.splm_list <- function(object, newdata, se.fit = FALSE, interval = c("none", "confidence", "prediction"), - level = 0.95, local, ...) { +predict.splm_list <- function(object, newdata, se.fit = FALSE, scale = NULL, + df = Inf, interval = c("none", "confidence", "prediction"), + level = 0.95, type = c("response", "terms"), local, + terms = NULL, na.action = na.fail, ...) { # match interval argument so the three display interval <- match.arg(interval) @@ -918,11 +943,34 @@ predict.splm_list <- function(object, newdata, se.fit = FALSE, interval = c("non if (missing(newdata)) { preds <- lapply(object, function(x) { - predict(x, se.fit = se.fit, interval = interval, level = level, local = local, ...) + predict( + x, + se.fit = se.fit, + scale = scale, + df = df, + interval = interval, + level = level, + type = type, + local = local, + terms = terms, # don't need na.action as it is fixed in predict + ... + ) }) } else { preds <- lapply(object, function(x) { - predict(x, newdata = newdata, se.fit = se.fit, interval = interval, level = level, local = local, ...) + predict( + x, + newdata = newdata, + se.fit = se.fit, + scale = scale, + df = df, + interval = interval, + level = level, + type = type, + local = local, + terms = terms, # don't need na.action as it is fixed in predict + ... + ) }) } names(preds) <- names(object) diff --git a/R/predict_glm.R b/R/predict_glm.R index b3f6edf..57ec0a4 100644 --- a/R/predict_glm.R +++ b/R/predict_glm.R @@ -23,7 +23,7 @@ #' augment(spgmod, newdata = moose_preds, interval = "prediction") #' } predict.spglm <- function(object, newdata, type = c("link", "response", "terms"), se.fit = FALSE, interval = c("none", "confidence", "prediction"), - level = 0.95, dispersion = NULL, terms = NULL, local, var_correct = TRUE, newdata_size, ...) { + level = 0.95, dispersion = NULL, terms = NULL, local, var_correct = TRUE, newdata_size, na.action = na.fail, ...) { @@ -154,6 +154,12 @@ predict.spglm <- function(object, newdata, type = c("link", "response", "terms") attr(newdata_model, "assign") <- attr_assign[keep_cols] attr(newdata_model, "contrasts") <- attr_contrasts + # finding rows w/out NA + ob_predictors <- complete.cases(newdata_model) + if (any(!ob_predictors)) { + stop("Cannot have NA values in predictors.", call. = FALSE) + } + # call terms if needed if (type == "terms") { # glm supports standard errors for terms objects but not intervals (no interval argument) @@ -616,7 +622,7 @@ get_pred_spglm <- function(newdata_list, se.fit, interval, formula, obdata, xcoo #' @export predict.spgautor <- function(object, newdata, type = c("link", "response", "terms"), se.fit = FALSE, interval = c("none", "confidence", "prediction"), - level = 0.95, dispersion = NULL, terms = NULL, local, var_correct = TRUE, newdata_size, ...) { + level = 0.95, dispersion = NULL, terms = NULL, local, var_correct = TRUE, newdata_size, na.action = na.fail, ...) { # match type argument so the two display type <- match.arg(type) @@ -697,6 +703,15 @@ predict.spgautor <- function(object, newdata, type = c("link", "response", "term attr(newdata_model, "assign") <- attr_assign[keep_cols] attr(newdata_model, "contrasts") <- attr_contrasts + # finding rows w/out NA + # this isn't really needed, because the error should come on model building + # but someone could accidentally write over their newdata object after fitting + # so it is good to keep for completeness + ob_predictors <- complete.cases(newdata_model) + if (any(!ob_predictors)) { + stop("Cannot have NA values in predictors.", call. = FALSE) + } + # call terms if needed if (type == "terms") { # scale df not used for glms @@ -923,8 +938,10 @@ get_pred_spgautor_parallel <- function(cluster_list, cov_matrix_lowchol, betahat #' @method predict spglm_list #' @order 11 #' @export -predict.spglm_list <- function(object, newdata, type = c("link", "response"), se.fit = FALSE, interval = c("none", "confidence", "prediction"), - newdata_size, level = 0.95, local, var_correct = TRUE, ...) { +predict.spglm_list <- function(object, newdata, type = c("link", "response", "terms"), se.fit = FALSE, + interval = c("none", "confidence", "prediction"), level = 0.95, + dispersion = NULL, terms = NULL, local, var_correct = TRUE, newdata_size, + na.action = na.fail, ...) { type <- match.arg(type) # match interval argument so the three display interval <- match.arg(interval) @@ -941,11 +958,36 @@ predict.spglm_list <- function(object, newdata, type = c("link", "response"), se if (missing(newdata)) { preds <- lapply(object, function(x) { - predict(x, type = type, se.fit = se.fit, interval = interval, newdata_size = newdata_size, level = level, local = local, var_correct = var_correct, ...) + predict( + x, + type = type, + se.fit = se.fit, + interval = interval, + level = level, + dispersion = dispersion, + terms = terms, + local = local, + var_correct = var_correct, + newdata_size = newdata_size, # don't need na.action as it is fixed in predict + ... + ) }) } else { preds <- lapply(object, function(x) { - predict(x, newdata = newdata, type = type, se.fit = se.fit, interval = interval, newdata_size = newdata_size, level = level, local = local, var_correct = var_correct, ...) + predict( + x, + newdata = newdata, + type = type, + se.fit = se.fit, + interval = interval, + level = level, + dispersion = dispersion, + terms = terms, + local = local, + var_correct = var_correct, + newdata_size = newdata_size, # don't need na.action as it is fixed in predict + ... + ) }) } names(preds) <- names(object) diff --git a/R/print.R b/R/print.R index 9047711..7bec939 100644 --- a/R/print.R +++ b/R/print.R @@ -45,7 +45,7 @@ print.splm <- function(x, digits = max(3L, getOption("digits") - 3L), if (!x$anisotropy) { spcoef <- spcoef[-which(names(spcoef) %in% c("rotate", "scale"))] } - if (inherits(coef(x, type = "spcov"), "none")) { + if (inherits(coef(x, type = "spcov"), c("none", "ie"))) { spcoef <- spcoef["ie"] } @@ -165,7 +165,7 @@ print.summary.splm <- function(x, if (!x$anisotropy) { spcoef <- spcoef[-which(names(spcoef) %in% c("rotate", "scale"))] } - if (inherits(x$coefficients$spcov, "none")) { + if (inherits(x$coefficients$spcov, c("none", "ie"))) { spcoef <- spcoef["ie"] } diff --git a/R/print_glm.R b/R/print_glm.R index fa919d0..47017f9 100644 --- a/R/print_glm.R +++ b/R/print_glm.R @@ -24,7 +24,7 @@ print.spglm <- function(x, digits = max(3L, getOption("digits") - 3L), if (!x$anisotropy) { spcoef <- spcoef[-which(names(spcoef) %in% c("rotate", "scale"))] } - if (inherits(coef(x, type = "spcov"), "none")) { + if (inherits(coef(x, type = "spcov"), c("none", "ie"))) { spcoef <- spcoef["ie"] } @@ -160,7 +160,7 @@ print.summary.spglm <- function(x, if (!x$anisotropy) { spcoef <- spcoef[-which(names(spcoef) %in% c("rotate", "scale"))] } - if (inherits(x$coefficients$spcov, "none")) { + if (inherits(x$coefficients$spcov, c("none", "ie"))) { spcoef <- spcoef["ie"] } diff --git a/R/spautor.R b/R/spautor.R index 4619ede..6251e00 100644 --- a/R/spautor.R +++ b/R/spautor.R @@ -47,7 +47,7 @@ #' specifying the partition factor. The partition factor assumes observations #' from different levels of the partition factor are uncorrelated. #' @param W Weight matrix specifying the neighboring structure used. -#' Not required if \code{data} is an \code{sf} object wtih \code{POLYGON} geometry, +#' Not required if \code{data} is an \code{sf} object with \code{POLYGON} geometry, #' as \code{W} is calculated internally using queen contiguity. If calculated internally, #' \code{W} is computed using \code{sf::st_intersects()}. Also not required if \code{data} #' is an \code{sf} object with \code{POINT} geometry as long as \code{cutoff} is specified. diff --git a/R/spautor_checks.R b/R/spautor_checks.R index 24cfdc1..aa43ea3 100644 --- a/R/spautor_checks.R +++ b/R/spautor_checks.R @@ -12,7 +12,7 @@ spautor_checks <- function(spcov_type, W_given, data, estmethod) { "exponential", "spherical", "gaussian", "triangular", "circular", "cubic", "pentaspherical", "cosine", "wave", "jbessel", "gravity", "rquad", "magnetic", - "matern", "cauchy", "pexponential", "none" + "matern", "cauchy", "pexponential", "none", "ie" )) { stop("Invalid spatial covariance type for spautor(). To fit models for point-referenced data, use splm().", call. = FALSE) } diff --git a/R/spcov_initial.R b/R/spcov_initial.R index 13ee39a..86cf524 100644 --- a/R/spcov_initial.R +++ b/R/spcov_initial.R @@ -71,7 +71,7 @@ #' #' All spatial covariance functions are valid in one spatial dimension. All #' spatial covariance functions except \code{triangular} and \code{cosine} are -#' valid in two dimensions. +#' valid in two dimensions. An alias for \code{none} is \code{ie}. #' #' When the spatial covariance function is \code{car} or \code{sar}, \code{extra} #' represents the variance parameter for the observations in \code{W} without @@ -100,7 +100,7 @@ spcov_initial <- function(spcov_type, de, ie, range, extra, rotate, scale, known stop("spcov_type must be specified", call. = FALSE) } else if (!spcov_type %in% c( "exponential", "spherical", "gaussian", "triangular", "circular", - "none", "cubic", "pentaspherical", "cosine", "wave", "matern", "car", "sar", "jbessel", + "none", "ie", "cubic", "pentaspherical", "cosine", "wave", "matern", "car", "sar", "jbessel", "gravity", "rquad", "magnetic", "cauchy", "pexponential" )) { stop(paste(spcov_type, "is not a valid spatial covariance function", sep = " "), call. = FALSE) diff --git a/R/spcov_initial_NA.R b/R/spcov_initial_NA.R index 89d9b5e..abe4739 100644 --- a/R/spcov_initial_NA.R +++ b/R/spcov_initial_NA.R @@ -36,7 +36,7 @@ spcov_initial_NA <- function(spcov_initial, anisotropy = FALSE, is_W_connected = spcov_val_default <- c(de = NA, ie = 0, range = NA, extra = NA) spcov_known_default <- c(de = FALSE, ie = TRUE, range = FALSE, extra = FALSE) } - } else if (inherits(spcov_initial, "none")) { + } else if (inherits(spcov_initial, c("none", "ie"))) { spcov_names <- c("de", "ie", "range", "rotate", "scale") spcov_val_default <- c(de = 0, ie = NA, range = Inf, rotate = 0, scale = 1) spcov_known_default <- c(de = TRUE, ie = FALSE, range = TRUE, rotate = TRUE, scale = TRUE) @@ -50,7 +50,7 @@ spcov_initial_NA <- function(spcov_initial, anisotropy = FALSE, is_W_connected = # put in is_known not in initial spcov_initial$is_known[spcov_out] <- spcov_known_default[spcov_out] # reorder names - if (inherits(spcov_initial, "none")) { + if (inherits(spcov_initial, c("none", "ie"))) { # reset if none covariance spcov_initial$initial[c("de", "range", "rotate", "scale")] <- spcov_val_default[c("de", "range", "rotate", "scale")] # put in is_known not in initial diff --git a/R/spcov_initial_NA_glm.R b/R/spcov_initial_NA_glm.R index 9272a7c..b3cb051 100644 --- a/R/spcov_initial_NA_glm.R +++ b/R/spcov_initial_NA_glm.R @@ -2,6 +2,11 @@ spcov_initial_NA_glm <- function(family, spcov_initial, anisotropy = FALSE, is_W spcov_initial_NA_val <- spcov_initial_NA(spcov_initial, anisotropy, is_W_connected) spcov_initial_NA_glm_val <- spcov_initial_NA_val + if (inherits(spcov_initial_NA_glm_val, "none")) { + spcov_initial_NA_glm_val$initial["ie"] <- 0 + spcov_initial_NA_glm_val$is_known["ie"] <- TRUE + } + # fix ie can be confounded with dispersion (problems with ie error) # if (family %in% c("nbinomial", "Gamma", "inverse.gaussian", "beta")) { # if (is.na(spcov_initial_NA_glm_val$initial[["ie"]])) { diff --git a/R/spcov_matrix.R b/R/spcov_matrix.R index d365487..4e7dbd5 100644 --- a/R/spcov_matrix.R +++ b/R/spcov_matrix.R @@ -70,6 +70,10 @@ spcov_matrix.none <- function(spcov_params, dist_matrix, diagtol = 0, ...) { spcov_matrix_val } +# spcov_matrix ie +#' @export +spcov_matrix.ie <- spcov_matrix.none + # spcov_matrix cubic #' @export spcov_matrix.cubic <- function(spcov_params, dist_matrix, diagtol = 0, ...) { diff --git a/R/spcov_optim2orig.R b/R/spcov_optim2orig.R index 5e7d516..3ade2b5 100644 --- a/R/spcov_optim2orig.R +++ b/R/spcov_optim2orig.R @@ -48,6 +48,8 @@ spcov_optim2orig.circular <- spcov_optim2orig.exponential #' @export spcov_optim2orig.none <- spcov_optim2orig.exponential #' @export +spcov_optim2orig.ie <- spcov_optim2orig.none +#' @export spcov_optim2orig.cubic <- spcov_optim2orig.exponential #' @export spcov_optim2orig.pentaspherical <- spcov_optim2orig.exponential diff --git a/R/spcov_orig2optim.R b/R/spcov_orig2optim.R index 2f693a0..492cddf 100644 --- a/R/spcov_orig2optim.R +++ b/R/spcov_orig2optim.R @@ -92,6 +92,8 @@ spcov_orig2optim.circular <- spcov_orig2optim.exponential #' @export spcov_orig2optim.none <- spcov_orig2optim.exponential #' @export +spcov_orig2optim.ie <- spcov_orig2optim.none +#' @export spcov_orig2optim.cubic <- spcov_orig2optim.exponential #' @export spcov_orig2optim.pentaspherical <- spcov_orig2optim.exponential diff --git a/R/spcov_params.R b/R/spcov_params.R index bc635f7..da857ba 100644 --- a/R/spcov_params.R +++ b/R/spcov_params.R @@ -9,7 +9,7 @@ #' \code{"pentaspherical"}, \code{"cosine"}, \code{"wave"}, #' \code{"jbessel"}, \code{"gravity"}, \code{"rquad"}, #' \code{"magnetic"}, \code{"matern"}, \code{"cauchy"}, \code{"pexponential"}, -#' \code{"car"}, \code{"sar"}, and \code{"none"}. +#' \code{"car"}, \code{"sar"}, \code{"none"}, and \code{"ie"} (an alias for \code{"none"}). #' @param de The spatially dependent (correlated) random error variance. Commonly referred to as #' a partial sill. #' @param ie The spatially independent (uncorrelated) random error variance. Commonly referred to as @@ -44,7 +44,7 @@ spcov_params <- function(spcov_type, de, ie, range, extra, rotate = 0, scale = 1 stop("spcov_type must be specified", call. = FALSE) } else if (!spcov_type %in% c( "exponential", "spherical", "gaussian", "triangular", "circular", - "none", "cubic", "pentaspherical", "cosine", "wave", "matern", "car", "sar", "jbessel", + "none", "ie", "cubic", "pentaspherical", "cosine", "wave", "matern", "car", "sar", "jbessel", "gravity", "rquad", "magnetic", "cauchy", "pexponential" )) { stop(paste(spcov_type), "is not a valid spatial covariance function.") @@ -54,7 +54,7 @@ spcov_params <- function(spcov_type, de, ie, range, extra, rotate = 0, scale = 1 stop("spcov_type must be specified.", call. = FALSE) } - if (spcov_type == "none") { + if (spcov_type %in% c("none", "ie")) { de <- 0 range <- Inf } @@ -69,11 +69,11 @@ spcov_params <- function(spcov_type, de, ie, range, extra, rotate = 0, scale = 1 } # some parameter specification checks - if (spcov_type != "none" && any(missing(de), missing(ie), missing(range))) { + if ((!spcov_type %in% c("none", "ie")) && any(missing(de), missing(ie), missing(range))) { stop("de, ie, and range must be specified.", call. = FALSE) } - if (spcov_type == "none" && missing(ie)) { + if (spcov_type %in% c("none", "ie") && missing(ie)) { stop("ie must be specified.", call. = FALSE) } @@ -117,7 +117,7 @@ spcov_params <- function(spcov_type, de, ie, range, extra, rotate = 0, scale = 1 stop("extra must be positive and no larger than 2.", call. = FALSE) } - if (spcov_type %in% c("exponential", "spherical", "gaussian", "triangular", "circular", "none", "cubic", "pentaspherical", "cosine", "wave", "jbessel", "gravity", "rquad", "magnetic")) { + if (spcov_type %in% c("exponential", "spherical", "gaussian", "triangular", "circular", "none", "ie", "cubic", "pentaspherical", "cosine", "wave", "jbessel", "gravity", "rquad", "magnetic")) { extra <- NULL } if (spcov_type %in% c("car", "sar")) { diff --git a/R/spcov_vector.R b/R/spcov_vector.R index eb60c81..9eb941a 100644 --- a/R/spcov_vector.R +++ b/R/spcov_vector.R @@ -59,6 +59,10 @@ spcov_vector.none <- function(spcov_params, dist_vector) { spcov_vector_val } +# spcov_vector ie +#' @export +spcov_vector.ie <- spcov_vector.none + # spcov_vector cubic #' @export spcov_vector.cubic <- function(spcov_params, dist_vector) { diff --git a/R/spgautor_checks.R b/R/spgautor_checks.R index 0873dfd..d33acf3 100644 --- a/R/spgautor_checks.R +++ b/R/spgautor_checks.R @@ -12,7 +12,7 @@ spgautor_checks <- function(family, spcov_type, W_given, data, estmethod) { "exponential", "spherical", "gaussian", "triangular", "circular", "cubic", "pentaspherical", "cosine", "wave", "jbessel", "gravity", "rquad", "magnetic", - "matern", "cauchy", "pexponential", "none" + "matern", "cauchy", "pexponential", "none", "ie" )) { stop("Invalid spatial covariance type for spautor(). To fit models for point-referenced data, use splm().", call. = FALSE) } diff --git a/R/spglm.R b/R/spglm.R index 6bf9098..c7ebf0e 100644 --- a/R/spglm.R +++ b/R/spglm.R @@ -28,7 +28,7 @@ #' \code{"pentaspherical"}, \code{"cosine"}, \code{"wave"}, #' \code{"jbessel"}, \code{"gravity"}, \code{"rquad"}, #' \code{"magnetic"}, \code{"matern"}, \code{"cauchy"}, \code{"pexponential"}, -#' and \code{"none"}. Parameterizations of each spatial covariance type are +#' \code{"ie"}, and \code{"none"}. Parameterizations of each spatial covariance type are #' available in Details. Multiple spatial covariance types can be provided as #' a character vector, and then \code{spglm()} is called iteratively for each #' element and a list is returned for each model fit. The default for @@ -123,7 +123,7 @@ #' \code{list(size = 100, method = "kmeans", var_adjust = "theoretical", parallel = FALSE)}. #' @param range_constrain An optional logical that indicates whether the range #' should be constrained to enhance numerical stability. If \code{range_constrain = TRUE}, -#' the maximum possible range value is 5 times the maximum distance in the domain. +#' the maximum possible range value is 4 times the maximum distance in the domain. #' If \code{range_constrain = FALSE}, then maximum possible range is unbounded. #' The default is \code{FALSE}. #' Note that if \code{range_constrain = TRUE} and the value of \code{range} in \code{spcov_initial} @@ -225,8 +225,19 @@ #' \item cauchy: \eqn{(1 + \eta^2)^{-extra}}, \eqn{extra > 0} #' \item pexponential: \eqn{exp(h^{extra}/range)}, \eqn{0 < extra \le 2} #' \item none: \eqn{0} +#' \item ie: \eqn{0} (see below for differences compared to none) #' } #' +#' There is a slight difference between \code{none} and \code{ie} (the \code{spcov_type}) from above. \code{none} +#' fixes the \code{"ie"} covariance parameter at zero, which should yield a model that has very similar +#' parameter estimates as one from [stats::glm]. Note that \code{spglm()} uses +#' a different likelihood, so direct likelihood-based comparisons (e.g., \code{AIC}) +#' should not be made to compare an \code{spglm()} model to a \code{glm()} one +#' (though deviance can be used). \code{ie} allows the \code{"ie"} covariance +#' parameter to vary. If the \code{"ie"} covariance parameter is estimated near +#' zero, then \code{none} and \code{ie} yield similar models. +#' +#' #' All spatial covariance functions are valid in one spatial dimension. All #' spatial covariance functions except \code{triangular} and \code{cosine} are #' valid in two dimensions. diff --git a/R/splm.R b/R/splm.R index d513284..0133ee4 100644 --- a/R/splm.R +++ b/R/splm.R @@ -120,7 +120,7 @@ #' \code{list(size = 100, method = "kmeans", var_adjust = "theoretical", parallel = FALSE)}. #' @param range_constrain An optional logical that indicates whether the range #' should be constrained to enhance numerical stability. If \code{range_constrain = TRUE}, -#' the maximum possible range value is 5 times the maximum distance in the domain. +#' the maximum possible range value is 4 times the maximum distance in the domain. #' If \code{range_constrain = FALSE}, then maximum possible range is unbounded. #' The default is \code{FALSE}. #' Note that if \code{range_constrain = TRUE} and the value of \code{range} in \code{spcov_initial} @@ -164,7 +164,7 @@ #' #' All spatial covariance functions are valid in one spatial dimension. All #' spatial covariance functions except \code{triangular} and \code{cosine} are -#' valid in two dimensions. +#' valid in two dimensions. An alias for \code{none} is \code{ie}. #' #' \code{estmethod} Details: The various estimation methods are #' \itemize{ @@ -304,7 +304,7 @@ splm <- function(formula, data, spcov_type, xcoord, ycoord, spcov_initial, } # set local explicitly to FALSE if iid - if (inherits(spcov_initial, "none") && is.null(random)) { + if (inherits(spcov_initial, c("none", "ie")) && is.null(random)) { local <- FALSE } @@ -355,7 +355,7 @@ splm <- function(formula, data, spcov_type, xcoord, ycoord, spcov_initial, - if (inherits(cov_est_object$spcov_params_val, "none") && is.null(random)) { + if (inherits(cov_est_object$spcov_params_val, c("none", "ie")) && is.null(random)) { model_stats <- get_model_stats_splm_iid(cov_est_object, data_object, estmethod) } else { model_stats <- get_model_stats_splm(cov_est_object, data_object, estmethod) diff --git a/R/spmodel-package.R b/R/spmodel-package.R index 1ab19e7..b39e4be 100644 --- a/R/spmodel-package.R +++ b/R/spmodel-package.R @@ -9,7 +9,7 @@ #' @importFrom stats AIC anova BIC coef coefficients complete.cases confint cooks.distance cor #' dbeta dbinom dgamma dnbinom dpois #' delete.response deviance dist fitted fitted.values formula .getXlevels hatvalues -#' influence kmeans lm logLik model.frame model.matrix model.offset model.response na.omit na.pass pchisq +#' influence kmeans lm logLik model.frame model.matrix model.offset model.response na.fail na.omit na.pass pchisq #' pnorm predict printCoefmat pt qnorm qqnorm qqline qt quantile #' rbeta rbinom rgamma rnbinom rpois resid #' residuals reformulate rnorm rstandard terms var vcov diff --git a/R/sprnorm.R b/R/sprnorm.R index 942cf25..d04f60b 100644 --- a/R/sprnorm.R +++ b/R/sprnorm.R @@ -266,6 +266,12 @@ sprnorm.none <- function(spcov_params, mean = 0, samples = 1, data, randcov_para sprnorm_val } + +#' @rdname sprnorm +#' @method sprnorm ie +#' @export +sprnorm.ie <- sprnorm.none + #' @rdname sprnorm #' @method sprnorm car #' @export diff --git a/R/tidy.R b/R/tidy.R index 58377c9..4019eaa 100644 --- a/R/tidy.R +++ b/R/tidy.R @@ -69,7 +69,7 @@ tidy.splm <- function(x, conf.int = FALSE, which_scale <- which(result$term == "scale") result <- result[-c(which_rotate, which_scale), , drop = FALSE] } - if (inherits(spcoef, "none")) { + if (inherits(spcoef, c("none", "ie"))) { which_ie <- which(result$term == "ie") result <- result[which_ie, , drop = FALSE] } diff --git a/R/tidy_glm.R b/R/tidy_glm.R index 6b4e0d6..312019d 100644 --- a/R/tidy_glm.R +++ b/R/tidy_glm.R @@ -41,7 +41,7 @@ tidy.spglm <- function(x, conf.int = FALSE, which_scale <- which(result$term == "scale") result <- result[-c(which_rotate, which_scale), , drop = FALSE] } - if (inherits(spcoef, "none")) { + if (inherits(spcoef, c("none", "ie"))) { which_ie <- which(result$term == "ie") result <- result[which_ie, , drop = FALSE] } diff --git a/data/seal.rda b/data/seal.rda index ded3457..5d0cbf1 100644 Binary files a/data/seal.rda and b/data/seal.rda differ diff --git a/docs/articles/SPGLMs.html b/docs/articles/SPGLMs.html index f1304c3..d515f05 100644 --- a/docs/articles/SPGLMs.html +++ b/docs/articles/SPGLMs.html @@ -789,7 +789,7 @@

Simulation -

An application to binary (presence/abscence) data +

An application to binary (presence/absence) data

In the next two sections we to show how to use spmodel to analyze data collected on moose in the Togiak region of Alaska, USA. @@ -883,7 +883,7 @@

An application to binary (presence/abscence) dat 
 AIC(bin_mod, bin_spmod, bin_spmod_anis)
#>                df      AIC
-#> bin_mod         1 719.1628
+#> bin_mod         0 717.1627
 #> bin_spmod       3 686.9031
 #> bin_spmod_anis  5 677.9403

The spatial models (, ) have a much lower AIC than the non-spatial @@ -996,7 +996,7 @@

An application to binary (presence/abscence) dat #> <chr> <int> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> #> 1 bin_spmod_anis 218 4 5 668. 678. 678. 695. -334. 161. #> 2 bin_spmod 218 4 3 681. 687. 687. 697. -340. 161. -#> 3 bin_mod 218 4 1 717. 719. 719. 723. -359. 294. +#> 3 bin_mod 218 4 0 717. 717. 717. 717. -359. 294. #> # ℹ 1 more variable: pseudo.r.squared <dbl>

The function augments the data with model diagnostics:

@@ -1206,7 +1206,7 @@ 
An application to count data
 
Estimated mean moose counts for the data (obtained via ) and
 predicted mean moose counts for the data (obtained via or ) share
-similiar patterns and are overlain by running

+similar patterns and are overlain by running

 
 aug_mod <- augment(count_mod_nb)
 aug_pred <- augment(count_mod_nb, newdata = moose_preds, type.predict = "response")
diff --git a/docs/articles/guide.html b/docs/articles/guide.html
index bffbb2f..ccb9462 100644
--- a/docs/articles/guide.html
+++ b/docs/articles/guide.html
@@ -753,25 +753,25 @@ 
An Areal Data Exampleseal by running

 
 seal

-
#> Simple feature collection with 62 features and 1 field
+
#> Simple feature collection with 149 features and 2 fields
 #> Geometry type: POLYGON
 #> Dimension:     XY
-#> Bounding box:  xmin: 913618.8 ymin: 1007542 xmax: 1116002 ymax: 1145054
+#> Bounding box:  xmin: 913618.8 ymin: 855730.2 xmax: 1221859 ymax: 1145054
 #> Projected CRS: NAD83 / Alaska Albers
-#> # A tibble: 62 × 2
-#>    log_trend                                                            geometry
-#>        <dbl>                                                       <POLYGON [m]>
-#>  1  NA       ((1035002 1054710, 1035002 1054542, 1035002 1053542, 1035002 10525…
-#>  2  -0.282   ((1037002 1039492, 1037006 1039490, 1037017 1039492, 1037035 10394…
-#>  3  -0.00121 ((1070158 1030216, 1070185 1030207, 1070187 1030207, 1070211 10302…
-#>  4   0.0354  ((1054906 1034826, 1054931 1034821, 1054936 1034822, 1055001 10348…
-#>  5  -0.0160  ((1025142 1056940, 1025184 1056889, 1025222 1056836, 1025256 10567…
-#>  6   0.0872  ((1026035 1044623, 1026037 1044605, 1026072 1044610, 1026083 10446…
-#>  7  -0.266   ((1100345 1060709, 1100287 1060706, 1100228 1060706, 1100170 10607…
-#>  8   0.0743  ((1030247 1029637, 1030248 1029637, 1030265 1029642, 1030328 10296…
-#>  9  NA       ((1043093 1020553, 1043097 1020550, 1043101 1020550, 1043166 10205…
-#> 10  -0.00961 ((1116002 1024542, 1116002 1023542, 1116002 1022542, 1116002 10215…
-#> # ℹ 52 more rows

+#> # A tibble: 149 × 3
+#>    log_trend stock                                                      geometry
+#>  *     <dbl> <fct>                                                 <POLYGON [m]>
+#>  1  NA       8     ((1035002 1054710, 1035002 1054542, 1035002 1053542, 1035002…
+#>  2  -0.282   8     ((1037002 1039492, 1037006 1039490, 1037017 1039492, 1037035…
+#>  3  -0.00121 8     ((1070158 1030216, 1070185 1030207, 1070187 1030207, 1070211…
+#>  4   0.0354  8     ((1054906 1034826, 1054931 1034821, 1054936 1034822, 1055001…
+#>  5  -0.0160  8     ((1025142 1056940, 1025184 1056889, 1025222 1056836, 1025256…
+#>  6   0.0872  8     ((1026035 1044623, 1026037 1044605, 1026072 1044610, 1026083…
+#>  7  -0.266   8     ((1100345 1060709, 1100287 1060706, 1100228 1060706, 1100170…
+#>  8   0.0743  8     ((1030247 1029637, 1030248 1029637, 1030265 1029642, 1030328…
+#>  9  NA       8     ((1043093 1020553, 1043097 1020550, 1043101 1020550, 1043166…
+#> 10  -0.00961 8     ((1116002 1024542, 1116002 1023542, 1116002 1022542, 1116002…
+#> # ℹ 139 more rows

 
We can learn more about the data by running
 help("seal", "spmodel").

 
We can visualize the distribution of log seal trends in the
@@ -817,51 +817,50 @@ 
An Areal Data Example#> spautor(formula = log_trend ~ 1, data = seal, spcov_type = "car")
 #> 
 #> Residuals:
-#>      Min       1Q   Median       3Q      Max 
-#> -0.34455 -0.10417  0.04410  0.07338  0.20475 
+#>       Min        1Q    Median        3Q       Max 
+#> -0.402385 -0.073439  0.004368  0.073465  0.848290 
 #> 
 #> Coefficients (fixed):
-#>             Estimate Std. Error z value Pr(>|z|)   
-#> (Intercept) -0.07090    0.02497  -2.839  0.00452 **
-#> ---
-#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#>             Estimate Std. Error z value Pr(>|z|)
+#> (Intercept) -0.01306    0.02004  -0.652    0.514
 #> 
 #> Coefficients (car spatial covariance):
 #>      de   range   extra 
-#> 0.03252 0.42037 0.02177

+#> 0.04936 0.56167 0.01750
 tidy(sealmod)
#> # A tibble: 1 × 5
 #>   term        estimate std.error statistic p.value
 #>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
-#> 1 (Intercept)  -0.0709    0.0250     -2.84 0.00452
+#> 1 (Intercept) -0.0131 0.0200 -0.652 0.514 
 glance(sealmod)
#> # A tibble: 1 × 10
 #>       n     p  npar value   AIC  AICc   BIC logLik deviance pseudo.r.squared
 #>   <int> <dbl> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>    <dbl>            <dbl>
-#> 1    34     1     3 -36.9 -30.9 -30.1 -26.3   18.4     33.1                0
+#> 1 94 1 3 -74.5 -68.5 -68.3 -60.9 37.3 92.9 0 
 augment(sealmod)
-
#> Simple feature collection with 34 features and 6 fields
+
#> Simple feature collection with 94 features and 6 fields
 #> Geometry type: POLYGON
 #> Dimension:     XY
-#> Bounding box:  xmin: 980001.5 ymin: 1010815 xmax: 1116002 ymax: 1145054
+#> Bounding box:  xmin: 980001.5 ymin: 855730.2 xmax: 1221859 ymax: 1145054
 #> Projected CRS: NAD83 / Alaska Albers
-#> # A tibble: 34 × 7
-#>    log_trend .fitted  .resid   .hat .cooksd .std.resid                  geometry
-#>  *     <dbl>   <dbl>   <dbl>  <dbl>   <dbl>      <dbl>             <POLYGON [m]>
-#>  1  -0.282   -0.0709 -0.211  0.0161 0.0209      -1.13  ((1037002 1039492, 10370…
-#>  2  -0.00121 -0.0709  0.0697 0.0473 0.0256       0.718 ((1070158 1030216, 10701…
-#>  3   0.0354  -0.0709  0.106  0.0290 0.0183       0.782 ((1054906 1034826, 10549…
-#>  4  -0.0160  -0.0709  0.0549 0.0228 0.00157      0.260 ((1025142 1056940, 10251…
-#>  5   0.0872  -0.0709  0.158  0.0276 0.0383       1.16  ((1026035 1044623, 10260…
-#>  6  -0.266   -0.0709 -0.195  0.0287 0.0530      -1.34  ((1100345 1060709, 11002…
-#>  7   0.0743  -0.0709  0.145  0.0491 0.0901       1.32  ((1030247 1029637, 10302…
-#>  8  -0.00961 -0.0709  0.0613 0.0122 0.00242      0.442 ((1116002 1024542, 11160…
-#>  9  -0.182   -0.0709 -0.111  0.0225 0.0224      -0.986 ((1079864 1025088, 10798…
-#> 10   0.00351 -0.0709  0.0744 0.0316 0.0100       0.555 ((1110363 1037056, 11103…
-#> # ℹ 24 more rows

+#> # A tibble: 94 × 7
+#>    log_trend .fitted   .resid    .hat   .cooksd .std.resid
+#>  *     <dbl>   <dbl>    <dbl>   <dbl>     <dbl>      <dbl>
+#>  1  -0.282   -0.0131 -0.269   0.00579 0.00705      -1.10  
+#>  2  -0.00121 -0.0131  0.0118  0.0156  0.000527      0.182 
+#>  3   0.0354  -0.0131  0.0485  0.0103  0.00100       0.310 
+#>  4  -0.0160  -0.0131 -0.00295 0.00749 0.0000640    -0.0921
+#>  5   0.0872  -0.0131  0.100   0.00939 0.00361       0.617 
+#>  6  -0.266   -0.0131 -0.253   0.0229  0.0877       -1.93  
+#>  7   0.0743  -0.0131  0.0873  0.0171  0.00620       0.597 
+#>  8  -0.00961 -0.0131  0.00345 0.00909 0.000254      0.166 
+#>  9  -0.182   -0.0131 -0.169   0.00785 0.00916      -1.08  
+#> 10   0.00351 -0.0131  0.0166  0.0111  0.000144      0.113 
+#> # ℹ 84 more rows
+#> # ℹ 1 more variable: geometry <POLYGON [m]>

Note that for spautor() models, the ie spatial covariance parameter is assumed zero by default (and omitted from the summary() output). This default behavior can be @@ -1063,25 +1062,25 @@

Predictionaugment() to make the predictions:

 augment(sealmod, newdata = sealmod$newdata)
-
#> Simple feature collection with 28 features and 2 fields
+
#> Simple feature collection with 55 features and 3 fields
 #> Geometry type: POLYGON
 #> Dimension:     XY
-#> Bounding box:  xmin: 913618.8 ymin: 1007542 xmax: 1115097 ymax: 1132682
+#> Bounding box:  xmin: 913618.8 ymin: 866541.5 xmax: 1214641 ymax: 1132682
 #> Projected CRS: NAD83 / Alaska Albers
-#> # A tibble: 28 × 3
-#>    log_trend  .fitted                                                   geometry
-#>  *     <dbl>    <dbl>                                              <POLYGON [m]>
-#>  1        NA -0.115   ((1035002 1054710, 1035002 1054542, 1035002 1053542, 1035…
-#>  2        NA -0.00908 ((1043093 1020553, 1043097 1020550, 1043101 1020550, 1043…
-#>  3        NA -0.0602  ((1099737 1054310, 1099752 1054262, 1099788 1054278, 1099…
-#>  4        NA -0.0359  ((1099002 1036542, 1099134 1036462, 1099139 1036431, 1099…
-#>  5        NA -0.0723  ((1076902 1053189, 1076912 1053179, 1076931 1053179, 1076…
-#>  6        NA -0.0548  ((1070501 1046969, 1070317 1046598, 1070308 1046542, 1070…
-#>  7        NA -0.0976  ((1072995 1054942, 1072996 1054910, 1072997 1054878, 1072…
-#>  8        NA -0.0714  ((960001.5 1127667, 960110.8 1127542, 960144.1 1127495, 9…
-#>  9        NA -0.0825  ((1031308 1079817, 1031293 1079754, 1031289 1079741, 1031…
-#> 10        NA -0.0592  ((998923.7 1053647, 998922.5 1053609, 998950 1053631, 999…
-#> # ℹ 18 more rows

+#> # A tibble: 55 × 4
+#>    log_trend stock  .fitted                                             geometry
+#>  *     <dbl> <fct>    <dbl>                                        <POLYGON [m]>
+#>  1        NA 8     -0.106   ((1035002 1054710, 1035002 1054542, 1035002 1053542…
+#>  2        NA 8      0.0359  ((1043093 1020553, 1043097 1020550, 1043101 1020550…
+#>  3        NA 8     -0.0303  ((1099737 1054310, 1099752 1054262, 1099788 1054278…
+#>  4        NA 8      0.00167 ((1099002 1036542, 1099134 1036462, 1099139 1036431…
+#>  5        NA 8     -0.0360  ((1076902 1053189, 1076912 1053179, 1076931 1053179…
+#>  6        NA 8     -0.0133  ((1070501 1046969, 1070317 1046598, 1070308 1046542…
+#>  7        NA 8     -0.0704  ((1072995 1054942, 1072996 1054910, 1072997 1054878…
+#>  8        NA 8     -0.0159  ((960001.5 1127667, 960110.8 1127542, 960144.1 1127…
+#>  9        NA 8     -0.0558  ((1031308 1079817, 1031293 1079754, 1031289 1079741…
+#> 10        NA 8     -0.0296  ((998923.7 1053647, 998922.5 1053609, 998950 105363…
+#> # ℹ 45 more rows

Advanced Features diff --git a/docs/articles/guide_files/figure-html/log_trend-1.png b/docs/articles/guide_files/figure-html/log_trend-1.png index a0875f8..1148459 100644 Binary files a/docs/articles/guide_files/figure-html/log_trend-1.png and b/docs/articles/guide_files/figure-html/log_trend-1.png differ diff --git a/docs/articles/introduction.html b/docs/articles/introduction.html index fd400a7..6b48dfe 100644 --- a/docs/articles/introduction.html +++ b/docs/articles/introduction.html @@ -862,7 +862,7 @@

Model Summaries#> model n p npar value AIC AICc BIC logLik deviance #> <chr> <int> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> #> 1 binmod 218 2 3 692. 698. 698. 708. -346. 190. -#> 2 glmod 218 2 1 715. 717. 717. 721. -358. 302. +#> 2 glmod 218 2 0 715. 715. 715. 715. -358. 302. #> # ℹ 1 more variable: pseudo.r.squared <dbl>

The lower AIC and AICc for the spatial binomial regression model indicates it is a much better fit to the data.

diff --git a/docs/articles/technical.html b/docs/articles/technical.html index 61263d7..bc27fac 100644 --- a/docs/articles/technical.html +++ b/docs/articles/technical.html @@ -1248,7 +1248,7 @@

The forms of R for each spatial covariance type available in splm(). All spatial covariance functions are valid in two dimensions except the triangular and cosine functions, which are only valid in one -dimension. +dimension. An alias for none is ie. @@ -1345,6 +1345,10 @@

\(\texttt{"none"}\) \(0\) + +\(\texttt{"ie"}\) +\(0\) +

@@ -2756,7 +2760,9 @@

spglm() Spatial Covariance Functions

Covariance functions for spglm() are defined the same as -covariance functions for splm().

+covariance functions for splm() except that for +"none", the ie parameter is also set to zero +(analogous to glm() models).

Model-fitting @@ -3471,11 +3477,13 @@

A Note on Computational Stability, where \(\sigma^2_{ie, up}\) denotes an “updated” version of \(\sigma^2_{ie}\). For spatial generalized linear models, \(\sigma^2_{ie, up} = \text{max}(\sigma^2_{ie}, -\sigma^2_{de}/10^4, d)\), where \(d = -\text{max}(1/10^4, s^2/10^4)\), where \(s^2\) is the sample variance of a linear -regression of \(\ln(\mathbf{y} + 1)\) -on \(\mathbf{X}\). This value of \(\sigma^2_{ie, up}\) is also added to the -diagonal of \(\mathbf{X}^\top +\sigma^2_{de}/10^4, d)\), where d = \(1/10^4\). Previously (spmodel +v\(\le 0.8.0\)), \(d = \text{min}(1/10^4, s^2/10^4)\), where +\(s^2\) is the sample variance of a +linear regression of \(\ln(\mathbf{y} + +1)\) on \(\mathbf{X}\). This +value of \(\sigma^2_{ie, up}\) is also +added to the diagonal of \(\mathbf{X}^\top \boldsymbol{\Sigma}^{-1} \mathbf{X} + \mathbf{X}^\top \boldsymbol{\Sigma}^{-1} (\mathbf{D} - \boldsymbol{\Sigma}^{-1})^{-1} \boldsymbol{\Sigma}^{-1} \mathbf{X}\), used via the diff --git a/docs/news/index.html b/docs/news/index.html index 53764ad..9ff01dc 100644 --- a/docs/news/index.html +++ b/docs/news/index.html @@ -70,10 +70,23 @@

Changelog

-

spmodel 0.9.0

+

spmodel 0.9.02024-11-06

Major Updates

-
  • Added the range_constrain argument to splm() and spglm() to contrain the range parameter to enhance numerical stability.
    • +
    • Added the range_constrain argument to splm() and spglm() to constrain the range parameter to enhance numerical stability. The default for range_constrain is FALSE, implying the range is not constrained.
      • +
    • Updated the seal data with additional polygons and a factor variable, stock, with two levels (8 and 10) that indicates seal stock (i.e., seal type).
      • +
+
+

Minor Updates

+
  • Changed diagonal tolerance threshold for spglm() and spgautor() model objects. See this link for details.
    • +
  • Added the "ie" spatial covariance type to splm() and spglm() models. For splm() models, "ie" is an alias for "none". For spglm() models, "none" now fixes both the de and ie covariance parameters at zero, while "ie" fixes the de covariance parameter at zero but allows the ie covariance parameter to vary. Thus, "none" from spmodel $\le$ v0.8.0 matches "ie" from spmodel v0.9.0 and but is different from "none" from spmodel v0.9.0.
    • +
  • Added the na.action argument to predict.spmodel() functions to clarify that missing values in newdata return an error.
    • +
  • Minor documentation updates.
    • +
+
+

Bug Fixes

+
  • Fixed a bug that caused incorrect degrees of freedom for the likelihood ratio test (anova(model1, model2)) when estmethod is "ml" for both models.
    • +
  • Fixed a bug that caused an error in anova(object1, object2) when the name of object1 had special characters (e.g., $).
@@ -235,7 +248,7 @@

Bug fixesspautor() that prevented an error from occurring when a partition factor was not categorical or not a factor
  • Fixed a bug in covmatrix(object, newdata) that returned a matrix with improper dimensions when spcov_type was "none".
  • Fixed a bug in predict() that caused an error when at least one level of a fixed effect factor was not observed within a local neighborhood (when the local method was "covariance" or "distance").
    • -
  • Fixed a bug in cooks.distance() that used the Pearson residuals instead of the standarized residuals.
    • +
  • Fixed a bug in cooks.distance() that used the Pearson residuals instead of the standardized residuals.

    • diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 6753687..b045613 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -7,7 +7,7 @@ articles: introduction: introduction.html SPGLMs: SPGLMs.html technical: technical.html -last_built: 2024-08-29T19:42Z +last_built: 2024-11-07T00:40Z urls: reference: https://usepa.github.io/spmodel/reference article: https://usepa.github.io/spmodel/articles diff --git a/docs/reference/augment.spmodel.html b/docs/reference/augment.spmodel.html index 5ce1640..23fde45 100644 --- a/docs/reference/augment.spmodel.html +++ b/docs/reference/augment.spmodel.html @@ -319,45 +319,46 @@

    Examples

    # missingness in original data spmod_seal <- spautor(log_trend ~ 1, data = seal, spcov_type = "car") augment(spmod_seal) -#> Simple feature collection with 34 features and 6 fields +#> Simple feature collection with 94 features and 6 fields #> Geometry type: POLYGON #> Dimension: XY -#> Bounding box: xmin: 980001.5 ymin: 1010815 xmax: 1116002 ymax: 1145054 +#> Bounding box: xmin: 980001.5 ymin: 855730.2 xmax: 1221859 ymax: 1145054 #> Projected CRS: NAD83 / Alaska Albers -#> # A tibble: 34 × 7 -#> log_trend .fitted .resid .hat .cooksd .std.resid geometry -#> * <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <POLYGON [m]> -#> 1 -0.282 -0.0709 -0.211 0.0161 0.0209 -1.13 ((1037002 1039492, 10370… -#> 2 -0.00121 -0.0709 0.0697 0.0473 0.0256 0.718 ((1070158 1030216, 10701… -#> 3 0.0354 -0.0709 0.106 0.0290 0.0183 0.782 ((1054906 1034826, 10549… -#> 4 -0.0160 -0.0709 0.0549 0.0228 0.00157 0.260 ((1025142 1056940, 10251… -#> 5 0.0872 -0.0709 0.158 0.0276 0.0383 1.16 ((1026035 1044623, 10260… -#> 6 -0.266 -0.0709 -0.195 0.0287 0.0530 -1.34 ((1100345 1060709, 11002… -#> 7 0.0743 -0.0709 0.145 0.0491 0.0901 1.32 ((1030247 1029637, 10302… -#> 8 -0.00961 -0.0709 0.0613 0.0122 0.00242 0.442 ((1116002 1024542, 11160… -#> 9 -0.182 -0.0709 -0.111 0.0225 0.0224 -0.986 ((1079864 1025088, 10798… -#> 10 0.00351 -0.0709 0.0744 0.0316 0.0100 0.555 ((1110363 1037056, 11103… -#> # ℹ 24 more rows +#> # A tibble: 94 × 7 +#> log_trend .fitted .resid .hat .cooksd .std.resid +#> * <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> +#> 1 -0.282 -0.0131 -0.269 0.00579 0.00705 -1.10 +#> 2 -0.00121 -0.0131 0.0118 0.0156 0.000527 0.182 +#> 3 0.0354 -0.0131 0.0485 0.0103 0.00100 0.310 +#> 4 -0.0160 -0.0131 -0.00295 0.00749 0.0000640 -0.0921 +#> 5 0.0872 -0.0131 0.100 0.00939 0.00361 0.617 +#> 6 -0.266 -0.0131 -0.253 0.0229 0.0877 -1.93 +#> 7 0.0743 -0.0131 0.0873 0.0171 0.00620 0.597 +#> 8 -0.00961 -0.0131 0.00345 0.00909 0.000254 0.166 +#> 9 -0.182 -0.0131 -0.169 0.00785 0.00916 -1.08 +#> 10 0.00351 -0.0131 0.0166 0.0111 0.000144 0.113 +#> # ℹ 84 more rows +#> # ℹ 1 more variable: geometry <POLYGON [m]> augment(spmod_seal, newdata = spmod_seal$newdata) -#> Simple feature collection with 28 features and 2 fields +#> Simple feature collection with 55 features and 3 fields #> Geometry type: POLYGON #> Dimension: XY -#> Bounding box: xmin: 913618.8 ymin: 1007542 xmax: 1115097 ymax: 1132682 +#> Bounding box: xmin: 913618.8 ymin: 866541.5 xmax: 1214641 ymax: 1132682 #> Projected CRS: NAD83 / Alaska Albers -#> # A tibble: 28 × 3 -#> log_trend .fitted geometry -#> * <dbl> <dbl> <POLYGON [m]> -#> 1 NA -0.115 ((1035002 1054710, 1035002 1054542, 1035002 1053542, 1035… -#> 2 NA -0.00908 ((1043093 1020553, 1043097 1020550, 1043101 1020550, 1043… -#> 3 NA -0.0602 ((1099737 1054310, 1099752 1054262, 1099788 1054278, 1099… -#> 4 NA -0.0359 ((1099002 1036542, 1099134 1036462, 1099139 1036431, 1099… -#> 5 NA -0.0723 ((1076902 1053189, 1076912 1053179, 1076931 1053179, 1076… -#> 6 NA -0.0548 ((1070501 1046969, 1070317 1046598, 1070308 1046542, 1070… -#> 7 NA -0.0976 ((1072995 1054942, 1072996 1054910, 1072997 1054878, 1072… -#> 8 NA -0.0714 ((960001.5 1127667, 960110.8 1127542, 960144.1 1127495, 9… -#> 9 NA -0.0825 ((1031308 1079817, 1031293 1079754, 1031289 1079741, 1031… -#> 10 NA -0.0592 ((998923.7 1053647, 998922.5 1053609, 998950 1053631, 999… -#> # ℹ 18 more rows +#> # A tibble: 55 × 4 +#> log_trend stock .fitted geometry +#> * <dbl> <fct> <dbl> <POLYGON [m]> +#> 1 NA 8 -0.106 ((1035002 1054710, 1035002 1054542, 1035002 1053542… +#> 2 NA 8 0.0359 ((1043093 1020553, 1043097 1020550, 1043101 1020550… +#> 3 NA 8 -0.0303 ((1099737 1054310, 1099752 1054262, 1099788 1054278… +#> 4 NA 8 0.00167 ((1099002 1036542, 1099134 1036462, 1099139 1036431… +#> 5 NA 8 -0.0360 ((1076902 1053189, 1076912 1053179, 1076931 1053179… +#> 6 NA 8 -0.0133 ((1070501 1046969, 1070317 1046598, 1070308 1046542… +#> 7 NA 8 -0.0704 ((1072995 1054942, 1072996 1054910, 1072997 1054878… +#> 8 NA 8 -0.0159 ((960001.5 1127667, 960110.8 1127542, 960144.1 1127… +#> 9 NA 8 -0.0558 ((1031308 1079817, 1031293 1079754, 1031289 1079741… +#> 10 NA 8 -0.0296 ((998923.7 1053647, 998922.5 1053609, 998950 105363… +#> # ℹ 45 more rows
    diff --git a/docs/reference/formula.spmodel.html b/docs/reference/formula.spmodel.html index 25d5e61..fa035be 100644 --- a/docs/reference/formula.spmodel.html +++ b/docs/reference/formula.spmodel.html @@ -113,7 +113,7 @@

    Examples

    ) formula(spmod) #> z ~ water + tarp -#> <environment: 0x000002d7a071edf8> +#> <environment: 0x000001fb1bc49688> diff --git a/docs/reference/moose.html b/docs/reference/moose.html index 0747498..831fbd4 100644 --- a/docs/reference/moose.html +++ b/docs/reference/moose.html @@ -84,7 +84,7 @@

    Format

  • strat: A factor representing strata (used for sampling). Can take values L and M.

  • count: The count (number) of moose observed.

  • presence: A binary factor representing whether no moose were observed (value 0) or at least one moose was observed -(va ue 1).

    • +(value 1).

  • geometry: POINT geometry representing coordinates in an Alaska Albers projection (EPSG: 3338). Distances between points are in meters.

    • diff --git a/docs/reference/predict.spmodel.html b/docs/reference/predict.spmodel.html index be304af..4070d8b 100644 --- a/docs/reference/predict.spmodel.html +++ b/docs/reference/predict.spmodel.html @@ -87,6 +87,7 @@

    Model predictions (Kriging)

    type = c("response", "terms"), local, terms = NULL, + na.action = na.fail, ... ) @@ -102,6 +103,7 @@

    Model predictions (Kriging)

    type = c("response", "terms"), local, terms = NULL, + na.action = na.fail, ... ) @@ -110,9 +112,14 @@

    Model predictions (Kriging)

    object, newdata, se.fit = FALSE, + scale = NULL, + df = Inf, interval = c("none", "confidence", "prediction"), level = 0.95, + type = c("response", "terms"), local, + terms = NULL, + na.action = na.fail, ... ) @@ -121,9 +128,14 @@

    Model predictions (Kriging)

    object, newdata, se.fit = FALSE, + scale = NULL, + df = Inf, interval = c("none", "confidence", "prediction"), level = 0.95, + type = c("response", "terms"), local, + terms = NULL, + na.action = na.fail, ... ) @@ -152,6 +164,7 @@

    Model predictions (Kriging)

    local, var_correct = TRUE, newdata_size, + na.action = na.fail, ... ) @@ -168,6 +181,7 @@

    Model predictions (Kriging)

    local, var_correct = TRUE, newdata_size, + na.action = na.fail, ... ) @@ -175,13 +189,16 @@

    Model predictions (Kriging)

    predict( object, newdata, - type = c("link", "response"), + type = c("link", "response", "terms"), se.fit = FALSE, interval = c("none", "confidence", "prediction"), - newdata_size, level = 0.95, + dispersion = NULL, + terms = NULL, local, var_correct = TRUE, + newdata_size, + na.action = na.fail, ... ) @@ -189,13 +206,16 @@

    Model predictions (Kriging)

    predict( object, newdata, - type = c("link", "response"), + type = c("link", "response", "terms"), se.fit = FALSE, interval = c("none", "confidence", "prediction"), - newdata_size, level = 0.95, + dispersion = NULL, + terms = NULL, local, var_correct = TRUE, + newdata_size, + na.action = na.fail, ... ) @@ -237,7 +257,13 @@

    Arguments

    interval

    Type of interval calculation. The default is "none". Other options are "confidence" (for confidence intervals) and -"prediction" (for prediction intervals).

    +"prediction" (for prediction intervals). When interval +is "none" or "prediction", predictions are returned (and when +requested, their corresponding uncertainties). When interval +is "confidence", mean estimates are returned (and when +requested, their corresponding uncertainties). This "none" behavior +differs from that of lm(), as lm() returns confidence +uncertainties (in .$se.fit).

    level
    @@ -292,6 +318,11 @@

    Arguments

    specified via either numeric position or name. The default is all terms are included.

    +
    na.action
    +

    Missing (NA) values in newdata will return an error and should +be removed before proceeding.

    + +
    ...

    Other arguments. Only used for models fit using splmRF() or spautorRF() where ... indicates other diff --git a/docs/reference/seal.html b/docs/reference/seal.html index fa98c07..c7142e8 100644 --- a/docs/reference/seal.html +++ b/docs/reference/seal.html @@ -80,7 +80,10 @@

    Estimated harbor-seal trends from abundance data in southeast Alaska, USA

    Format

    -

    A sf object with 62 rows and 2 columns:

    All spatial covariance functions are valid in one spatial dimension. All spatial covariance functions except triangular and cosine are -valid in two dimensions.

    +valid in two dimensions. An alias for none is ie.

    When the spatial covariance function is car or sar, extra represents the variance parameter for the observations in W without at least one neighbor (other than itself) -- these are called unconnected diff --git a/docs/reference/spcov_params.html b/docs/reference/spcov_params.html index 5b793aa..165d4af 100644 --- a/docs/reference/spcov_params.html +++ b/docs/reference/spcov_params.html @@ -89,7 +89,7 @@

    Arguments

    "pentaspherical", "cosine", "wave", "jbessel", "gravity", "rquad", "magnetic", "matern", "cauchy", "pexponential", -"car", "sar", and "none".

    +"car", "sar", "none", and "ie" (an alias for "none").

    de
    diff --git a/docs/reference/spgautor.html b/docs/reference/spgautor.html index b0eef06..003a3c4 100644 --- a/docs/reference/spgautor.html +++ b/docs/reference/spgautor.html @@ -364,22 +364,22 @@

    Examples

    #> spcov_type = "car") #> #> Deviance Residuals: -#> Min 1Q Median 3Q Max -#> -4.484 -2.461 -1.030 0.386 2.823 +#> Min 1Q Median 3Q Max +#> -4.5034 -2.1822 -1.2659 -0.3373 6.2996 #> #> Coefficients (fixed): #> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -3.7975 0.3396 -11.18 <2e-16 *** +#> (Intercept) -3.6675 0.1937 -18.94 <2e-16 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Coefficients (car spatial covariance): -#> de range extra -#> 0.001738 0.995833 0.002374 +#> de range extra +#> 0.00657 0.99331 0.00657 #> #> Coefficients (Dispersion for Gamma family): #> dispersion -#> 0.3051 +#> 0.3374 diff --git a/docs/reference/spglm.html b/docs/reference/spglm.html index f2427c7..140fe22 100644 --- a/docs/reference/spglm.html +++ b/docs/reference/spglm.html @@ -136,7 +136,7 @@

    Arguments

    "pentaspherical", "cosine", "wave", "jbessel", "gravity", "rquad", "magnetic", "matern", "cauchy", "pexponential", -and "none". Parameterizations of each spatial covariance type are +"ie", and "none". Parameterizations of each spatial covariance type are available in Details. Multiple spatial covariance types can be provided as a character vector, and then spglm() is called iteratively for each element and a list is returned for each model fit. The default for @@ -261,7 +261,7 @@

    Arguments

    range_constrain

    An optional logical that indicates whether the range should be constrained to enhance numerical stability. If range_constrain = TRUE, -the maximum possible range value is 5 times the maximum distance in the domain. +the maximum possible range value is 4 times the maximum distance in the domain. If range_constrain = FALSE, then maximum possible range is unbounded. The default is FALSE. Note that if range_constrain = TRUE and the value of range in spcov_initial @@ -374,7 +374,16 @@

    Details

  • cauchy: \((1 + \eta^2)^{-extra}\), \(extra > 0\)

  • pexponential: \(exp(h^{extra}/range)\), \(0 < extra \le 2\)

  • none: \(0\)

    • -

    All spatial covariance functions are valid in one spatial dimension. All +

  • ie: \(0\) (see below for differences compared to none)

    • +

    There is a slight difference between none and ie (the spcov_type) from above. none +fixes the "ie" covariance parameter at zero, which should yield a model that has very similar +parameter estimates as one from stats::glm. Note that spglm() uses +a different likelihood, so direct likelihood-based comparisons (e.g., AIC) +should not be made to compare an spglm() model to a glm() one +(though deviance can be used). ie allows the "ie" covariance +parameter to vary. If the "ie" covariance parameter is estimated near +zero, then none and ie yield similar models.

    +

    All spatial covariance functions are valid in one spatial dimension. All spatial covariance functions except triangular and cosine are valid in two dimensions.

    estmethod Details: The various estimation methods are

    All spatial covariance functions are valid in one spatial dimension. All spatial covariance functions except triangular and cosine are -valid in two dimensions.

    +valid in two dimensions. An alias for none is ie.

    estmethod Details: The various estimation methods are