From 801b7e7211fe3795cb6d6a3ca2dea0fe6b14e195 Mon Sep 17 00:00:00 2001 From: markvanderloo Date: Fri, 11 Nov 2016 09:08:46 +0100 Subject: [PATCH 1/9] documentation --- pkg/R/data_checks.R | 2 +- pkg/R/rtrim-pkg.R | 2 +- pkg/R/utils.R | 9 ++++++--- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/pkg/R/data_checks.R b/pkg/R/data_checks.R index a079c76..0aad71c 100644 --- a/pkg/R/data_checks.R +++ b/pkg/R/data_checks.R @@ -35,7 +35,7 @@ check_observations <- function(x,...){ #' first column indicates the time point, the second column indicates for what covariate value insufficient #' counts are found.} #' \item{For Model 2, \code{$errors} is a list with a single element \code{changepoints}. It -#' points out what changepoint lead to a time slice with zero observations.} +#' points out what changepoints lead to a time slice with zero observations.} #' } #' #' diff --git a/pkg/R/rtrim-pkg.R b/pkg/R/rtrim-pkg.R index 2386585..3c6b02d 100644 --- a/pkg/R/rtrim-pkg.R +++ b/pkg/R/rtrim-pkg.R @@ -1,6 +1,6 @@ #' Trend and Indices for Monitoring Data #' -#' The TRIIM model is used to estimate species populations based on frequent +#' The TRIM model is used to estimate species populations based on frequent #' (anual) counts at fixed sites. The current package is a complete #' re-implementation of the \code{Delphi} based #' \href{https://www.cbs.nl/en-gb/society/nature-and-environment/indices-and-trends--trim--}{TRIM} diff --git a/pkg/R/utils.R b/pkg/R/utils.R index 463c5a7..1d24203 100644 --- a/pkg/R/utils.R +++ b/pkg/R/utils.R @@ -4,18 +4,21 @@ #' @section Details: #' #' Control how much output \code{\link{trim}} writes to the screen while -#' fitting the model. +#' fitting the model. By default, \code{trim} only returns the output +#' and does not write any progress to the screen. After calling +#' \code{set_trim_verbose(TRUE)}, \code{trim} will write information +#' about running iterations and convergence to the screen during optmization. +#' #' #' #' @param verbose \code{[logical]} toggle verbosity. \code{TRUE} means: be #' verbose, \code{FALSE} means be quiet (this is the default). -#' +#' @family modelspec #' @export set_trim_verbose <- function(verbose=FALSE){ stopifnot(isTRUE(verbose)|!isTRUE(verbose)) options(trim_verbose=verbose) } -set_trim_verbose(TRUE) # Convenience function for console output during runs rprintf <- function(fmt,...) { if(getOption("trim_verbose")) cat(sprintf(fmt,...)) } From d54f25c7a80b7f95720985c8da285189c6f18c01 Mon Sep 17 00:00:00 2001 From: markvanderloo Date: Fri, 11 Nov 2016 10:54:19 +0100 Subject: [PATCH 2/9] check_observations now works for model 2 w/covariates --- pkg/NAMESPACE | 2 +- pkg/R/data_checks.R | 44 ++++++++++++++++++++------------ pkg/R/rtrim-pkg.R | 8 +++--- pkg/tests/testthat/test_checks.R | 11 +++++++- 4 files changed, 44 insertions(+), 21 deletions(-) diff --git a/pkg/NAMESPACE b/pkg/NAMESPACE index 6519540..4f96e32 100644 --- a/pkg/NAMESPACE +++ b/pkg/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +S3method(check_observations,character) S3method(check_observations,data.frame) S3method(check_observations,trimcommand) S3method(coef,trim) @@ -21,7 +22,6 @@ S3method(trim,formula) S3method(trim,trimcommand) S3method(vcov,trim) S3method(wald,trim) -export(check_observarions.character) export(check_observations) export(count_summary) export(gof) diff --git a/pkg/R/data_checks.R b/pkg/R/data_checks.R index 0aad71c..8ddf3c6 100644 --- a/pkg/R/data_checks.R +++ b/pkg/R/data_checks.R @@ -16,7 +16,7 @@ check_observations <- function(x,...){ } #' @param model \code{[numeric]} Model 1, 2 or 3? -#' @param covar \code{[character|numeric]} column index of covariates in \code{x} +#' @param covars \code{[character|numeric]} column index of covariates in \code{x} #' @param time.id \code{[character|numeric]} column index of time points in \code{x} #' @param count.id \code{[character|numeric]} column index of the counts in \code{x} #' @param changepoints \code{[numeric]} Changepoints (model 2 only) @@ -32,36 +32,48 @@ check_observations <- function(x,...){ #' points with insufficient counts}. #' \item{For model 3 with covariates, \code{$errors} is a named list with an element for each covariate #' for which insufficients counts are encountered. Each element is a two-column \code{data.frame}. The -#' first column indicates the time point, the second column indicates for what covariate value insufficient +#' first column indicates the time point, the second column indicates for which covariate value insufficient #' counts are found.} -#' \item{For Model 2, \code{$errors} is a list with a single element \code{changepoints}. It -#' points out what changepoints lead to a time slice with zero observations.} +#' \item{For Model 2, without covariates \code{$errors} is a list with a single +#' element \code{changepoints}. It points out what changepoints lead to a time +#' slice with zero observations.} +#' \item{For Model 2, with covariates \code{$errors} is a named list with an +#' element for each covariate for which inssufficients counts are encountered. +#' Each element is a two-column \code{data.frame}, The first colum indicates the +#' changepoint, the second column indicates for which covariate value +#' insufficient counts are found.} #' } #' #' #' #' @export #' @rdname check_observations -check_observations.data.frame <- function(x, model, covar = list() +check_observations.data.frame <- function(x, model, covars = list() , changepoints = numeric(0), time.id="time",count.id="count", eps=1e-8, ...){ stopifnot(model %in% 1:3) out <- list() - if (model==3 && length(covar) == 0 ){ + if (model==3 && length(covars) == 0 ){ time_totals <- tapply(X = x[,count.id], INDEX = x[,time.id], FUN = sum, na.rm=TRUE) ii <- time_totals <= eps out$sufficient <- !any(ii) out$errors <- setNames(list(x[ii,time.id]),time.id) - } else if (model == 3 && length(covar>0)) { - out$errors <- get_cov_count_errlist(x[,count.id],x[,time.id],covars=x[covar],timename=time.id) + } else if (model == 3 && length(covars>0)) { + out$errors <- get_cov_count_errlist(x[,count.id],x[,time.id],covars=x[covars],timename=time.id) out$sufficient <- length(out$errors) == 0 - } else if (model == 2 && length(covar) == 0) { + } else if ( model == 2 ) { pieces <- pieces_from_changepoints(x[,time.id],changepoints) - time_totals <- tapply(X=x[,count.id],INDEX=pieces, FUN = sum, na.rm=TRUE) - ii <- time_totals <= eps - out$sufficient <- !any(ii) - out$errors <- list(changepoint = as.numeric(names(time_totals))[ii]) + ok <- pieces > 0 # allow zero counts for changepoint 0 + if ( length(covars) == 0){ + time_totals <- tapply(X=x[ok,count.id],INDEX=pieces[ok], FUN = sum, na.rm=TRUE) + ii <- time_totals <= eps + out$sufficient <- !any(ii) + out$errors <- list(changepoint = as.numeric(names(time_totals))[ii]) + } else { + out$errors <- get_cov_count_errlist(x[ok,count.id], pieces[ok], x[ok,covars,drop=FALSE], timename="changepoint") + out$sufficient <- length(out$errors) == 0 + } } out @@ -72,13 +84,13 @@ check_observations.data.frame <- function(x, model, covar = list() #' @rdname check_observations check_observations.trimcommand <- function(x,...){ dat <- read_tdf(x$file) - check_observations.data.frame(x=dat,model=x$model, covar=x$labels[x$covariates] + check_observations.data.frame(x=dat,model=x$model, covars=x$labels[x$covariates] , changepoints = x$changepoints) } #' @export #' @rdname check_observations -check_observarions.character <- function(x,...){ +check_observations.character <- function(x,...){ tc <- read_tcf(x) check_observations.trimcommand(tc,...) } @@ -169,7 +181,7 @@ get_cov_count_errlist <- function(count, time, covars, timename="time"){ # CP0 = levels(df$time)[1] # err <- df[df$Freq==0 & df$time!=CP0, 1:2] err <- df[df$Freq==0, 1:2] - + row.names(err) <- NULL if (nrow(err) > 0){ names(err)[2] <- covname ERR[[covname]] <- err diff --git a/pkg/R/rtrim-pkg.R b/pkg/R/rtrim-pkg.R index 19b8fa7..02e7d43 100644 --- a/pkg/R/rtrim-pkg.R +++ b/pkg/R/rtrim-pkg.R @@ -1,8 +1,10 @@ #' Trend and Indices for Monitoring Data #' -#' The TRIIM model is used to estimate species populations based on frequent -#' (anual) counts at fixed sites. The current package is a complete -#' re-implementation of the \code{Delphi} based +#' The TRIIM model is used to estimate species populations based on frequent +#' (annual) counts at a varying collection of sites. The model is able to take +#' account of missing data by imputing prior to estimation of population totals. +#' The current package is a complete re-implementation of the \code{Delphi} +#' based #' \href{https://www.cbs.nl/en-gb/society/nature-and-environment/indices-and-trends--trim--}{TRIM} #' software by J. Pannekoek. #' diff --git a/pkg/tests/testthat/test_checks.R b/pkg/tests/testthat/test_checks.R index 10ae445..c9b5968 100644 --- a/pkg/tests/testthat/test_checks.R +++ b/pkg/tests/testthat/test_checks.R @@ -24,7 +24,7 @@ test_that("check_observations",{ # model 3 with covariates d <- data.frame(count = rep(0:1,2), time = rep(1:2,each=2), cov = rep(1:2,2)) - out <- check_observations(d, model=3,covar="cov") + out <- check_observations(d, model=3,covars="cov") expect_false(out$sufficient) expect_equal(out$errors,list(cov=data.frame(time=factor(1:2),cov=factor(c(1,1),levels=1:2))) ) @@ -36,6 +36,15 @@ test_that("check_observations",{ # model 2, with covariates + d <- data.frame( + time=1:4 + , X = c("A","A","A","B") + , count = c(1,1,1,0) + ) + + out <- check_observations(d,model=2, covars = "X",changepoints=1) + expect_false(out$sufficient) + expect_equal(out$errors$X, data.frame(changepoint=factor(1),X=factor("B",levels=c("A","B")) )) }) From 4beb1828af84d6bc0376d9da1c4e41ac8fca7a13 Mon Sep 17 00:00:00 2001 From: markvanderloo Date: Fri, 11 Nov 2016 10:58:11 +0100 Subject: [PATCH 3/9] fix in trim.formula signature --- pkg/R/trim.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pkg/R/trim.R b/pkg/R/trim.R index 862d748..ae54ae1 100644 --- a/pkg/R/trim.R +++ b/pkg/R/trim.R @@ -248,7 +248,7 @@ trim.data.frame <- function(x, count.id = "count", site.id="site", time.id="time #' @rdname trim #' @param data \code{[data.frame]} Data containing at least counts, times, and sites. #' @export -trim.formula <- function(x, data, model=c(1,2,3), weights=numeric(0) +trim.formula <- function(x, data, model=2, weights=numeric(0) , serialcor=FALSE, overdisp=FALSE, changepoints=integer(0), stepwise=FALSE , autodelete=FALSE, ...){ stopifnot(inherits(data,"data.frame")) From 35e02e2b9e78162ae173228beb3ab925d20411c4 Mon Sep 17 00:00:00 2001 From: markvanderloo Date: Fri, 11 Nov 2016 11:11:57 +0100 Subject: [PATCH 4/9] added print method for trim objects --- pkg/R/trim_post.R | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/pkg/R/trim_post.R b/pkg/R/trim_post.R index b6edeae..8392a37 100644 --- a/pkg/R/trim_post.R +++ b/pkg/R/trim_post.R @@ -1,11 +1,29 @@ # ########################################### TRIM postprocessing functions #### +# ================================================================= Print ====== + + +#' print a 'trim' object +#' +#' @param x a \code{\link{trim}} object +#' @param ... currently unused +#' +#' @export +#' @keywords internal +print.trim <- function(x,...){ + cat("Call:\n") + print(x$call) + cat("\n",x$convergence,"\n") + cat("\nCoefficients:\n") + print(coef.trim(x)) +} + # ================================================================= Summary ==== # ----------------------------------------------------------------- extract ---- -#' Print summary information for a TRIM job +#' Summary information for a TRIM job #' #' Print a summary of a \code{\link{trim}} object. #' @@ -41,6 +59,7 @@ summary.trim <- function(object,...) { ),class="trim.summary") } + #' @export #' @keywords internal print.trim.summary <- function(x,...){ From 629a79b0a3d6b2fceace38afddd65583d8838c3d Mon Sep 17 00:00:00 2001 From: markvanderloo Date: Fri, 11 Nov 2016 11:12:24 +0100 Subject: [PATCH 5/9] trim now does not return invisibly anymore --- pkg/NAMESPACE | 1 + pkg/R/trim.R | 26 +++++++++++++------------- 2 files changed, 14 insertions(+), 13 deletions(-) diff --git a/pkg/NAMESPACE b/pkg/NAMESPACE index 4f96e32..17e31d3 100644 --- a/pkg/NAMESPACE +++ b/pkg/NAMESPACE @@ -8,6 +8,7 @@ S3method(gof,trim) S3method(plot,trim.index) S3method(plot,trim.overall) S3method(print,count.summary) +S3method(print,trim) S3method(print,trim.gof) S3method(print,trim.overall) S3method(print,trim.summary) diff --git a/pkg/R/trim.R b/pkg/R/trim.R index ae54ae1..28ba3ce 100644 --- a/pkg/R/trim.R +++ b/pkg/R/trim.R @@ -183,18 +183,18 @@ trim.trimcommand <- function(x,...){ if (isTRUE(x$covin)) covin <- read_icv(x) else covin <- list() - out <- trim_estimate(count=dat$count - , time.id = dat$time - , site.id = dat$site - , covars = dat[covars] - , model = x$model - , serialcor = x$serialcor - , overdisp = x$overdisp - , changepoints = x$changepoints - , stepwise = x$stepwise - , weights = wgt - , covin = covin - , ...) + trim_estimate(count=dat$count + , time.id = dat$time + , site.id = dat$site + , covars = dat[covars] + , model = x$model + , serialcor = x$serialcor + , overdisp = x$overdisp + , changepoints = x$changepoints + , stepwise = x$stepwise + , weights = wgt + , covin = covin + , ...) } @@ -229,7 +229,7 @@ trim.data.frame <- function(x, count.id = "count", site.id="site", time.id="time stopifnot(all(weights>0), length(weights) == 0 || length(weights) == nrow(x)) # estimate the model and return - m <- trim_estimate( + trim_estimate( count = x[,count.id] , time.id = x[,time.id] , site.id = x[,site.id] From 206e2bf924f8fe86d8daeb0dc8cb62fd0bc90c20 Mon Sep 17 00:00:00 2001 From: markvanderloo Date: Fri, 11 Nov 2016 11:15:06 +0100 Subject: [PATCH 6/9] added badges --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index 6a3f542..50f1818 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,7 @@ [![Build Status](https://travis-ci.org/markvanderloo/rtrim.svg)](https://travis-ci.org/markvanderloo/rtrim) [![Coverage Status](https://coveralls.io/repos/github/markvanderloo/rtrim/badge.svg?branch=master)](https://coveralls.io/github/markvanderloo/rtrim?branch=master) +[![CRAN](http://www.r-pkg.org/badges/version/rtrim)](http://cran.r-project.org/web/packages/rtrim) +[![Downloads](http://cranlogs.r-pkg.org/badges/rtrim)](http://cran.r-project.org/package=rtrim/)[![Research software impact](http://depsy.org/api/package/cran/rtrim/badge.svg)](http://depsy.org/package/r/rtrim) # rtrim Reimplementation of [TRIM](https://www.cbs.nl/en-gb/society/nature-and-environment/indices-and-trends--trim--) for R From 8a278b9492bbb95a61b00e0053a34d113b2d3360 Mon Sep 17 00:00:00 2001 From: Mark van der Loo Date: Fri, 11 Nov 2016 11:16:00 +0100 Subject: [PATCH 7/9] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 50f1818..a4c9dc0 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ [![Build Status](https://travis-ci.org/markvanderloo/rtrim.svg)](https://travis-ci.org/markvanderloo/rtrim) [![Coverage Status](https://coveralls.io/repos/github/markvanderloo/rtrim/badge.svg?branch=master)](https://coveralls.io/github/markvanderloo/rtrim?branch=master) [![CRAN](http://www.r-pkg.org/badges/version/rtrim)](http://cran.r-project.org/web/packages/rtrim) -[![Downloads](http://cranlogs.r-pkg.org/badges/rtrim)](http://cran.r-project.org/package=rtrim/)[![Research software impact](http://depsy.org/api/package/cran/rtrim/badge.svg)](http://depsy.org/package/r/rtrim) +[![Downloads](http://cranlogs.r-pkg.org/badges/rtrim)](http://cran.r-project.org/package=rtrim/) # rtrim Reimplementation of [TRIM](https://www.cbs.nl/en-gb/society/nature-and-environment/indices-and-trends--trim--) for R From 591f4c1bdfcae5afa95c615447ade653024c905a Mon Sep 17 00:00:00 2001 From: markvanderloo Date: Fri, 11 Nov 2016 11:27:22 +0100 Subject: [PATCH 8/9] first arg of index is now 'x' (not z) --- pkg/R/trim_index.R | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/pkg/R/trim_index.R b/pkg/R/trim_index.R index 5bb490a..c8a84a5 100644 --- a/pkg/R/trim_index.R +++ b/pkg/R/trim_index.R @@ -57,7 +57,7 @@ #' Extract time-indices from trim output #' -#' @param z an object of class \code{\link{trim}} +#' @param x an object of class \code{\link{trim}} #' @param which Selector to distinguish between time indices based on the imputed data (default), #' the fitted model, or both. #' @param covars Switch to compute indices for covariate categories as well. @@ -88,12 +88,12 @@ #' z <- trim(count ~ time + site + Habitat, data=skylark, model=2) #' index(z, covars=TRUE) #' -index <- function(z, which=c("imputed","fitted","both"), covars=FALSE, base=1) { - stopifnot(inherits(z,"trim")) +index <- function(x, which=c("imputed","fitted","both"), covars=FALSE, base=1) { + stopifnot(inherits(x,"trim")) # Match base to actual time points - if (base %in% z$time.id) { - base = which(base == z$time.id) + if (base %in% x$time.id) { + base = which(base == x$time.id) stopifnot(length(base)==1) } @@ -102,22 +102,22 @@ index <- function(z, which=c("imputed","fitted","both"), covars=FALSE, base=1) { which <- match.arg(which) if (which=="fitted") { # Call workhorse function to do the actual computation - mod <- .index(z$tt_mod, z$var_tt_mod, base) + mod <- .index(x$tt_mod, x$var_tt_mod, base) # Store results in a data frame - out <- data.frame(time = z$time.id, + out <- data.frame(time = x$time.id, fitted = mod$tau, se_fit = sqrt(mod$var_tau)) } else if (which=="imputed") { # Idem, using the imputed time totals instead - imp <- .index(z$tt_imp, z$var_tt_imp, base) - out = data.frame(time = z$time.id, + imp <- .index(x$tt_imp, x$var_tt_imp, base) + out = data.frame(time = x$time.id, imputed = imp$tau, se_imp = sqrt(imp$var_tau)) } else if (which=="both") { # Idem, using both modelled and imputed time totals. - mod <- .index(z$tt_mod, z$var_tt_mod, base) - imp <- .index(z$tt_imp, z$var_tt_imp, base) - out = data.frame(time = z$time.id, + mod <- .index(x$tt_mod, x$var_tt_mod, base) + imp <- .index(x$tt_imp, x$var_tt_imp, base) + out = data.frame(time = x$time.id, fitted = mod$tau, se_fit = sqrt(mod$var_tau), imputed = imp$tau, @@ -128,7 +128,7 @@ index <- function(z, which=c("imputed","fitted","both"), covars=FALSE, base=1) { if (covars) { out <- cbind(data.frame(covariate="Overall", category=0), out) - tt <- z$covar_tt + tt <- x$covar_tt index <- list() ncovar <- length(tt) for (i in seq_len(ncovar)) { @@ -138,7 +138,7 @@ index <- function(z, which=c("imputed","fitted","both"), covars=FALSE, base=1) { index[[name]] <- vector("list", nclass) for (j in seq_len(nclass)) { ttij <- tti[[j]] - df = data.frame(covariate=ttij$covariate, category=ttij$class, time=z$time.id) + df = data.frame(covariate=ttij$covariate, category=ttij$class, time=x$time.id) # Compute model index+variance if (which %in% c("fitted","both")) { idx <- .index(ttij$mod, ttij$var_mod, base) From 73cb39129f662a8eba07e4f2fc871c74054d3e34 Mon Sep 17 00:00:00 2001 From: markvanderloo Date: Fri, 11 Nov 2016 11:27:44 +0100 Subject: [PATCH 9/9] doc fix in count_summary --- pkg/R/read_tdf.R | 6 +++--- pkg/R/trim_post.R | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/pkg/R/read_tdf.R b/pkg/R/read_tdf.R index 9dddc58..b0f033d 100644 --- a/pkg/R/read_tdf.R +++ b/pkg/R/read_tdf.R @@ -97,9 +97,9 @@ snifreport <- function(file, colclasses){ #' #' @param x A \code{data.frame} with anual counts per site. #' @param eps Numbers smaller then \code{eps} are treated a zero. -#' @param count.id name of the column containing the counts -#' @param time.id name of the column containing the time codes -#' @param site.id name of the column containing the site id's +#' @param count.id \code{[character|numeric]} index of the column containing the counts +#' @param time.id \code{[character|numeric]} index of the column containing the time codes +#' @param site.id \code{[character|numeric]} index of the column containing the site id's #' #' @return A \code{list} of class \code{count.summary} containing individual names. #' @export diff --git a/pkg/R/trim_post.R b/pkg/R/trim_post.R index 8392a37..411e7eb 100644 --- a/pkg/R/trim_post.R +++ b/pkg/R/trim_post.R @@ -174,12 +174,12 @@ overdispersion <- function(x){ #' slope defined by \eqn{\alpha^* + \beta^*d_j} is returned. #' #' Finally, note that both the overall slope and the deviations can be written -#' in multiplicative format. +#' in multiplicative form as well. #' #' #' @param object TRIM output structure (i.e., output of a call to \code{trim}) #' @param representation \code{[character]} Choose the coefficient -#' representation \code{"trend"} and \code{"deviations"} are for model 3 only. +#' representation. Options \code{"trend"} and \code{"deviations"} are for model 3 only. #' @param ... currently unused #' #' @return A \code{data.frame} containing coefficients and their standard errors,