Merge branch 'master' of https://github.com/markvanderloo/rtrim

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/)
 
 # rtrim
 Reimplementation of [TRIM](https://www.cbs.nl/en-gb/society/nature-and-environment/indices-and-trends--trim--) for R

pkg/NAMESPACE

@@ -1,12 +1,14 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(check_observations,character)
 S3method(check_observations,data.frame)
 S3method(check_observations,trimcommand)
 S3method(coef,trim)
 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)
@@ -21,7 +23,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)

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 changepoint 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

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

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.
 #'

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]
@@ -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"))

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)

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,...){
@@ -155,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,

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,...)) }

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")) ))
 })