Skip to content

Commit

Permalink
Changed order 'site' and 'time' arguments is all TRIM calls
Browse files Browse the repository at this point in the history
  • Loading branch information
pwbogaart committed Nov 21, 2016
1 parent 7a374d1 commit 95ed9ab
Show file tree
Hide file tree
Showing 13 changed files with 93 additions and 105 deletions.
24 changes: 12 additions & 12 deletions pkg/R/read_tdf.R
Original file line number Diff line number Diff line change
Expand Up @@ -90,37 +90,37 @@ snifreport <- function(file, colclasses){

#' Compute a summary of counts
#'
#'
#'
#' Summarize counts over a trim input dataset. Sites without counts are removed
#' before any counting takes place (since these will not be used when calling
#' \code{\link{trim}}). For the remaining records, the total number of
#' zero-counts, positive counts, total number of observed counts and the total
#' number of missings are reported.
#'
#'
#' @param x A \code{data.frame} with annual counts per site.
#' @param eps \code{[numeric]} Numbers smaller then \code{eps} are treated a zero.
#' @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
#'
#' @param time.id \code{[character|numeric]} index of the column containing the time codes
#' @param count.id \code{[character|numeric]} index of the column containing the counts
#'
#' @return A \code{list} of class \code{count.summary} containing individual names.
#' @export
#' @examples
#' data(skylark)
#' count_summary(skylark)
#'
#'
#' s <- count_summary(skylark)
#' s$zero_counts # obtain number of zero counts
count_summary <- function(x, count.id="count",time.id="time",site.id="site", eps=1e-8){
count_summary <- function(x, count.id="count", site.id="site", time.id="time", eps=1e-8){

site_count <- tapply(X = x[,count.id], INDEX = x[site.id], FUN=sum, na.rm=TRUE)
ii <- abs(site_count) < eps
sites_wout_counts <- character(0)
if (any(ii)){
sites_wout_counts <- names(site_count[ii])
x <- x[!x[,site.id] %in% sites_wout_counts,,drop=FALSE]
}

cnt <- x[,count.id]
L <- list(
sites = length(unique(x[,site.id]))
Expand All @@ -135,10 +135,10 @@ count_summary <- function(x, count.id="count",time.id="time",site.id="site", eps
}

#' print a count summary
#'
#'
#' @param x An R object
#' @param ... unused
#'
#'
#' @export
#' @keywords internal
print.count.summary <- function(x,...){
Expand All @@ -152,7 +152,7 @@ print.count.summary <- function(x,...){
printf("Total number of observed counts %8d\n",x$total_observed)
printf("Number of missing counts %8d\n",x$missing_counts)
printf("Total number of counts %8d\n",x$total_counts)

}


Expand Down
26 changes: 13 additions & 13 deletions pkg/R/trim.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
#'
#'
#' Compute TRIM model parameters.
#'
#'
#'
#' @section Models:
#'
Expand All @@ -14,8 +14,8 @@
#' In \bold{Model 1}, the imputed values are modeled as
#'
#' \eqn{\ln\mu_{ij} = \alpha_i,}
#'
#' where \eqn{\alpha_i} is the site effect. This model implies that the counts
#'
#' where \eqn{\alpha_i} is the site effect. This model implies that the counts
#' vary accross sites, not over time. The model-based \link[=totals]{time totals} are equal to
#' each time point and the model-based \link[=index]{indices} are all equal to one.
#'
Expand Down Expand Up @@ -138,7 +138,7 @@
#' @param x a \code{\link{trimcommand}}, a \code{data.frame}, or a \code{formula}
#' If \code{x} is a \code{formula}, the dependent variable (left-hand-side)
#' is treated as the 'counts' variable. The first and second independent variable
#' are treated as the 'time' and 'site' variable, \bold{in that specific order}. All
#' are treated as the 'site' and 'time' variable, \bold{in that specific order}. All
#' other variables are treated as covariates.
#' @param ... Currently unused
#'
Expand All @@ -151,7 +151,7 @@
#'
#' @examples
#' data(skylark)
#' m <- trim(count ~ time + site, data=skylark, model=2)
#' m <- trim(count ~ site + time, data=skylark, model=2)
#' summary(m)
#' coefficients(m)
#'
Expand All @@ -161,13 +161,13 @@
#' # match weights to sites
#' weights <- w[skylark$site]
#' # run model
#' m <- trim(count ~ time + site, data=skylark, model=3, weights=weights)
#' m <- trim(count ~ site + time, data=skylark, model=3, weights=weights)
#'
#'
#' # An example using change points, a covariate, and overdispersion
#' # 1 is added as cp automatically
#' cp <- c(2,6)
#' m <- trim(count ~ time + site + Habitat, data=skylark, model=2, changepoints=cp, overdisp=TRUE)
#' m <- trim(count ~ site + time + Habitat, data=skylark, model=2, changepoints=cp, overdisp=TRUE)
#' coefficients(m)
#' # check significance of changes in slope
#' wald(m)
Expand All @@ -191,8 +191,8 @@ trim.trimcommand <- function(x,...){
else covin <- list()

trim_estimate(count=dat$count
, time.id = dat$time
, site.id = dat$site
, time.id = dat$time
, covars = dat[covars]
, model = x$model
, serialcor = x$serialcor
Expand All @@ -207,8 +207,8 @@ trim.trimcommand <- function(x,...){


#' @param count.id \code{[character]} name of the column holding species counts
#' @param time.id \code{[character]} name of the column holding the time of counting
#' @param site.id \code{[character]} name of the column holding the site id
#' @param time.id \code{[character]} name of the column holding the time of counting
#' @param covars \code{[character]} name(s) of column(s) holding covariates
#' @param model \code{[numeric]} TRIM model type 1, 2, or 3.
#' @param weights \code{[numeric]} Optional vector of site weights. The length of
Expand Down Expand Up @@ -238,8 +238,8 @@ trim.data.frame <- function(x, count.id = "count", site.id="site", time.id="time
# estimate the model and return
trim_estimate(
count = x[,count.id]
, time.id = x[,time.id]
, site.id = x[,site.id]
, time.id = x[,time.id]
, covars = x[covars]
, model = model
, serialcor=serialcor
Expand All @@ -253,15 +253,15 @@ 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.
#' @param data \code{[data.frame]} Data containing at least counts, sites, and times
#' @export
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"))
L <- parse_formula(x, names(data))
trim.data.frame(x=data
, count.id=L$count, time.id=L$time, site.id=L$site, covars = L$cov
, count.id=L$count, site.id=L$site, time.id=L$time, covars = L$cov
, model=model, weights=weights
, serialcor=serialcor, overdisp=overdisp, changepoints=changepoints
, stepwise=stepwise, autodelete=autodelete, ...)
Expand All @@ -280,7 +280,7 @@ parse_formula <- function(x, vars){
if (!all(valid_vars)){
stop(sprintf("Variables %s not found in data", pr(all_vars[!valid_vars])))
}
list(count = lhs, time = rhs[1], site=rhs[2], cov=rhs[-(1:2)])
list(count = lhs, site=rhs[1], time = rhs[2], cov=rhs[-(1:2)])
}


Expand Down
10 changes: 5 additions & 5 deletions pkg/R/trim_gof.R
Original file line number Diff line number Diff line change
Expand Up @@ -21,18 +21,18 @@
#' and \code{AIC}, for Chi-squared, Likelihoof Ratio and Akaike informatiuon content,
#' respectively.
#' @export
#'
#'
#' @family analyses
#' @examples
#' @examples
#' data(skylark)
#' z <- trim(count ~ time + site, data=skylark, model=2)
#' z <- trim(count ~ site + time, data=skylark, model=2)
#' # prettyprint GOF information
#' gof(z)
#'
#'
#' # get individual elements, e.g. p-value
#' L <- gof(z)
#' LR_p <- L$LR$p # get p-value for likelihood ratio
#'
#'
gof <- function(x) UseMethod("gof")

# Here is a simple wrapper function for TRIM output lists.
Expand Down
6 changes: 3 additions & 3 deletions pkg/R/trim_index.R
Original file line number Diff line number Diff line change
Expand Up @@ -81,14 +81,14 @@
#' @examples
#'
#' data(skylark)
#' z <- trim(count ~ time + site, data=skylark, model=2)
#' z <- trim(count ~ site + time, data=skylark, model=2)
#' index(z)
#' # mimic classic TRIM:
#' index(z, "both")
#' # Extract standard errors for the imputed data
#' SE <- index(z,"imputed")$se_mod
#' # Include covariates
#' z <- trim(count ~ time + site + Habitat, data=skylark, model=2)
#' z <- trim(count ~ site + time + Habitat, data=skylark, model=2)
#' ind <- index(z, covars=TRUE)
#' plot(ind)
index <- function(x, which=c("imputed","fitted","both"), covars=FALSE, base=1) {
Expand Down Expand Up @@ -178,7 +178,7 @@ index <- function(x, which=c("imputed","fitted","both"), covars=FALSE, base=1) {
#' @examples
#'
#' data(skylark)
#' z <- trim(count ~ time + site + Habitat, data=skylark, model=2)
#' z <- trim(count ~ site + time + Habitat, data=skylark, model=2)
#' idx <- index(z, covars=TRUE)
#' plot(idx, covar="Habitat", main="Skylark")
#'
Expand Down
6 changes: 3 additions & 3 deletions pkg/R/trim_overall.R
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@
#'
#' # obtain the overall slope accross all change points.
#' data(skylark)
#' z <- trim(count ~ time + site, data=skylark, model=2)
#' z <- trim(count ~ site + time, data=skylark, model=2)
#' overall(z)
#' plot(overall(z))
#'
Expand All @@ -39,7 +39,7 @@
#' L$coef
#'
#' # Obtain the slope from changepoint to changepoint
#' z <- trim(count ~ time + site, data=skylark, model=2,changepoints=c(1,4,6))
#' z <- trim(count ~ site + time, data=skylark, model=2,changepoints=c(1,4,6))
#' # slope from time point 1 to 5
#' overall(z,changepoints=c(1,5,7))
overall <- function(x, which=c("imputed","fitted"), changepoints=numeric(0)) {
Expand Down Expand Up @@ -255,7 +255,7 @@ print.trim.overall <- function(x,...) {
#'
#' @examples
#' data(skylark)
#' m <- trim(count ~ time + site, data=skylark, model=2)
#' m <- trim(count ~ site + time, data=skylark, model=2)
#' plot(overall(m))
#'
#' @export
Expand Down
12 changes: 6 additions & 6 deletions pkg/R/trim_post.R
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ print.trim <- function(x,...){
#' @examples
#'
#' data(skylark)
#' z <- trim(count ~ time + site,data=skylark,model=2,overdisp=TRUE)
#' z <- trim(count ~ site + time, data=skylark, model=2, overdisp=TRUE)
#'
#' summary(z)
summary.trim <- function(object,...) {
Expand Down Expand Up @@ -190,7 +190,7 @@ overdispersion <- function(x){
#' @family analyses
#' @examples
#' data(skylark)
#' z <- trim(count ~ time + site, data=skylark, model=2,overdisp=TRUE)
#' z <- trim(count ~ site + time, data=skylark, model=2, overdisp=TRUE)
#' coefficients(z)
coef.trim <- function(object,
representation=c("standard","trend","deviations"),...) {
Expand Down Expand Up @@ -244,7 +244,7 @@ coef.trim <- function(object,
#' @family analyses
#' @examples
#' data(skylark)
#' z <- trim(count ~ time + site, data=skylark, model=2,changepoints=c(3,5))
#' z <- trim(count ~ site + time, data=skylark, model=2, changepoints=c(3,5))
#' totals(z)
#'
#' totals(z, "both") # mimics classic TRIM
Expand Down Expand Up @@ -297,7 +297,7 @@ export.trim.totals <- function(x, species, stratum) {
#' @family analyses
#' @examples
#' data(skylark)
#' z <- trim(count ~ time + site, data=skylark, model=2);
#' z <- trim(count ~ site + time, data=skylark, model=2);
#' totals(z)
#' vcv1 <- vcov(z) # Use imputed data
#' vcv2 <- vcov(z,"model") # Use modelled data
Expand Down Expand Up @@ -334,7 +334,7 @@ vcov.trim <- function(object, which = c("imputed","model"), ... ) {
#' @family analyses
#' @examples
#' data(skylark)
#' z <- trim(count ~ time + site, data=skylark, model=2);
#' z <- trim(count ~ site + time, data=skylark, model=2);
#' out <- results(z)
results <- function(z) {
stopifnot(inherits(z,"trim"))
Expand Down Expand Up @@ -380,7 +380,7 @@ plot.trim.results <- function(z, ...) {
#' @examples
#'
#' data(skylark)
#' z <- trim(count ~ time + site,data=skylark,model=2,overdisp=TRUE)
#' z <- trim(count ~ site + time, data=skylark, model=2, overdisp=TRUE)
#' now_what(z)
now_what <- function(z) {
stopifnot(inherits(z,"trim"))
Expand Down
12 changes: 6 additions & 6 deletions pkg/R/trim_refine.R
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@
#' TRIM stepwise refinement
#'
#' @param count a numerical vector of count data.
#' @param time.id a numerical vector time points for each count data point.
#' @param site.id a numerical vector time points for each count data point.
#' @param time.id a numerical vector time points for each count data point.
#' @param covars an optional list of covariates
#' @param model a model type selector
#' @param serialcor a flag indication of autocorrelation has to be taken into account.
Expand All @@ -17,7 +17,7 @@
#' Usually this information is retrieved by a set of postprocessing functions
#'
#' @keywords internal
trim_refine <- function(count, time.id, site.id, covars=list(),
trim_refine <- function(count, site.id, time.id, covars=list(),
model=2, serialcor=FALSE, overdisp=FALSE,
changepoints=integer(0),weights=numeric(0))
{
Expand All @@ -27,13 +27,13 @@ trim_refine <- function(count, time.id, site.id, covars=list(),

# Always start with an estimation using all proposed changepoints
cur_cp <- org_cp
z <- trim_workhorse(count, time.id, site.id, covars, model, serialcor, overdisp, cur_cp, weights)
z <- trim_workhorse(count, site.id, time.id, covars, model, serialcor, overdisp, cur_cp, weights)

# # Hack: remove all except the first changepoints
# n <- length(org_cp)
# active[2:n] <- FALSE
# cur_cp <- org_cp[active]
# z <- trim_workhorse(count, time.id, site.id, covars, model, serialcor, overdisp, cur_cp)
# z <- trim_workhorse(count, site.id, time.id, covars, model, serialcor, overdisp, cur_cp)

for (iter in 1:100) {
# Phase 1: can one of the changepoints be removed?
Expand All @@ -56,7 +56,7 @@ trim_refine <- function(count, time.id, site.id, covars=list(),
# If a changepoint has been removed, we'll need to re-estimate the model
if (removed) {
cur_cp = org_cp[active]
z <- trim_workhorse(count, time.id, site.id, covars, model, serialcor, overdisp, cur_cp, weights)
z <- trim_workhorse(count, site.id, time.id, covars, model, serialcor, overdisp, cur_cp, weights)
}

# Phase 2: try to re-insert previously removed changepoints
Expand Down Expand Up @@ -103,7 +103,7 @@ trim_refine <- function(count, time.id, site.id, covars=list(),
# If a changepoint has been re-inserted, we'll need to re-estimate the model
if (added) {
cur_cp = org_cp[active]
z <- trim_workhorse(count, time.id, site.id, covars, model, serialcor, overdisp, cur_cp, weights)
z <- trim_workhorse(count, site.id, time.id, covars, model, serialcor, overdisp, cur_cp, weights)
}

# Finished refinement?
Expand Down
4 changes: 2 additions & 2 deletions pkg/R/trim_wald.R
Original file line number Diff line number Diff line change
Expand Up @@ -53,10 +53,10 @@
#'
#' @examples
#' data(skylark)
#' z2 <- trim(count ~ time + site, data=skylark, model=2)
#' z2 <- trim(count ~ site + time, data=skylark, model=2)
#' # print info on significance of slope parameters
#' print(z2)
#' z3 <- trim(count ~ time + site, data=skylark, model=3)
#' z3 <- trim(count ~ site + time, data=skylark, model=3)
#' # print info on significance of deviations from linear trend
#' wald(z3)
wald <- function(x) UseMethod("wald")
Expand Down
Loading

0 comments on commit 95ed9ab

Please sign in to comment.