diff --git a/DESCRIPTION b/DESCRIPTION index 5f6eebb..34260bf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,9 +1,8 @@ Package: activity Type: Package Title: Animal Activity Statistics -Version: 1.1 -Date: 2016-09-23 -Author: Marcus Rowcliffe +Version: 1.2 +Author: Marcus Rowcliffe Maintainer: Marcus Rowcliffe Description: Provides functions to fit kernel density functions to animal activity time data; plot activity distributions; @@ -11,11 +10,12 @@ Description: Provides functions to fit kernel density functions activity metrics through bootstrapping; and evaluate variation in linear variables with time (or other circular variables). License: GPL-3 +Encoding: UTF-8 Depends: methods -Imports: circular, overlap, pbapply, -LazyData: TRUE -RoxygenNote: 5.0.1 +Imports: circular, pbapply, +LazyData: true +RoxygenNote: 6.1.1 NeedsCompilation: no -Packaged: 2016-09-25 21:15:10 UTC; Rowcliffe.M +Packaged: 2019-01-16 15:47:24 UTC; Rowcliffe.M Repository: CRAN -Date/Publication: 2016-09-26 08:03:36 +Date/Publication: 2019-01-16 16:20:04 UTC diff --git a/MD5 b/MD5 index 616b8f1..ce5ced5 100644 --- a/MD5 +++ b/MD5 @@ -1,23 +1,23 @@ -8c94e5b21d96c2641dd8b120693d48a9 *DESCRIPTION -d6d32bbdd0ef2430432d632b36648906 *NAMESPACE -cfd442f3b42dabb6c62a98c978e38070 *R/activity_code_1_1.r +0110609d6cda094b6dbc52e2d2463ac2 *DESCRIPTION +77b7cb20b88d12d3f93f90a6bd2eddf1 *NAMESPACE +f7e7aa1c906be1242f8f16e3f6711cb2 *R/activity_code.r 601e4a0baf81c965b52eb2e420cccb89 *data/BCIspeed.rda 49d79fe17bc56762619f55a6e1c816b4 *data/BCItime.rda -9b1bab5216976a1a5b8ce652938217df *man/BCIspeed.Rd -f991d5b73e53ccf0725455fd849a1ae8 *man/BCItime.Rd -42d9583064c7b58e3416c6e6d890842d *man/activity.Rd -4a199aeff8a1ecbd05a694543e164771 *man/actmod-class.Rd -ae24fa53bdc5484026e4d61b8ded5669 *man/bwcalc.Rd -772c5b40ca3dd72fa5d53e5d933136e8 *man/compareAct.Rd -9479e2d17d267007e50a176b4deb64cd *man/compareCkern.Rd -90e9a4389ed1a52b83fe1f1331dbd9bc *man/compareTimes.Rd -a498a15607e2e962141a57f7edb4e910 *man/dvmkern.Rd -3580acc5366e36f9d98e125b253d5db0 *man/fitact.Rd -862b3cfdc11583fab590002d872c7dc9 *man/fitlincirc.Rd -b47af65fd7f08455aa45c730eff17f23 *man/lincircKern.Rd -3c09eed335eea2b1c458bcf9ca3967c1 *man/lincircmod-class.Rd -c5760a4d8dc6a69069646262946abe3e *man/plot-actmod-methods.Rd -e2bace67835000dee06dbd65b967a93f *man/plot-lincircmod-methods.Rd -31b70269a919604eb5b8f2570ff193f6 *man/rvmkern.Rd -ed1c4c423549fa4817b7ec0c132358dc *man/show-actmod-method.Rd -f28910d3dbe3188bbc511e447c61e50f *man/yfromx.Rd +15290f976976a3087ebcc33909050685 *man/BCIspeed.Rd +092ff3feba8134e3732ae2d6928709bc *man/BCItime.Rd +f9a8656625b0003da4874ebe9027fe1c *man/activity.Rd +22217069e29f038820dca2cb3b2b48c8 *man/actmod-class.Rd +f03c4513788ddc920d9419bbd6cd5468 *man/bwcalc.Rd +04a3e4dd30d55517635bac3e825c3991 *man/compareAct.Rd +d1f467fc386482ef74022b4525392818 *man/compareCkern.Rd +a177b1b4fd6d565c1a8a27193f7f7161 *man/compareTimes.Rd +83d195d4d282c71c6fd68bf0710f5f51 *man/density2.Rd +657639fee7a58c873750cbcefff482e7 *man/dvmkern.Rd +2e25bb96d338adefa051a662211c60b0 *man/fitact.Rd +0d3239eb12023df32e74f7f003d06acf *man/fitlincirc.Rd +65e154f383e0abdd1c6ea4fc0a49d121 *man/lincircKern.Rd +cbe944f12f84f1713acc949e380f35b5 *man/lincircmod-class.Rd +aa03318797be35f8d805142b15323ac6 *man/ovl5.Rd +c0f8e7372852a652e16f29c1c70ee0eb *man/plot.actmod.Rd +a4776155b8dd4e218f75373feefd5881 *man/plot.lincircmod.Rd +e66ffa98a14a17af257b5c3e49a56497 *man/redf.Rd diff --git a/NAMESPACE b/NAMESPACE index 034ac3f..8a06734 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,17 +1,16 @@ import(methods) -importFrom("graphics", "hist", "lines") -importFrom("stats", "ecdf", "optimise", "pchisq", "quantile", "runif", "sd") +S3method(plot,actmod) +S3method(plot,lincircmod) export(bwcalc) export(compareAct) export(compareCkern) export(compareTimes) +export(density2) export(dvmkern) export(fitact) export(fitlincirc) export(lincircKern) -export(rvmkern) -export(yfromx) +export(ovl5) +export(redf) exportClasses(actmod) exportClasses(lincircmod) -exportMethods(plot) -exportMethods(show) diff --git a/R/activity_code.r b/R/activity_code.r new file mode 100644 index 0000000..b7fc7aa --- /dev/null +++ b/R/activity_code.r @@ -0,0 +1,835 @@ +#' Animal activity statistics +#' +#' Provides functions to estimate and compare activity parameters from sensor data. +#' +#' @details Sensors that record active animals (eg camera traps) build up a record of +#' the distribution of activity over the course of the day. Records are more frequent +#' when animals are more active, and less frequent or absent when animals are inactive. +#' The area under the distribution of records thus contains information on the overall +#' level of activity in a sampled population. This package provides tools for plotting +#' activity distributions, quantifying the overall level of activity with error, and +#' statistically comparing distributions through bootstrapping. +#' +#' The core function is \code{fitact}, which creates an \code{actmod} object containing +#' the circular kernel PDF, and the activity level estimate derived from this. The +#' generic plot function for \code{actmod} objects plots the distribution. Functions +#' starting with \code{compare} make statistical comparisons between distributions or +#' activity estimates. Note that all time or other circular data should be in radians +#' (in the range 0 to 2*pi). +#' +#' @references Rowcliffe, M., Kays, R., Kranstauber, B., Carbone, C., Jansen, P.A. (2014) Quantifying animal activity level using camera trap data. Methods in Ecology and Evolution. +#' @docType package +#' @name activity +NULL + +#' Animal record time of day data +#' +#' Barro Colorado Island 2008 data: times of day at which animal records occured +#' (\code{time}), together with species (\code{species}). +#' +#' @format A dataframe with 17820 observations and 2 variables. +#' @source http://dx.doi.org/10.6084/m9.figshare.1160536 +#' @name BCItime +#' @docType data +NULL + +#' Animal speed data +#' +#' Barro Colorado Island 2008 data: speeds of animal passages past camera traps +#' (\code{speed}), together with species (\code{species}) and time of day (\code{time}) +#' for each record. +#' +#' @format A dataframe with 2204 observations and 3 variables. +#' @source http://dx.doi.org/10.6084/m9.figshare.1160536 +#' @name BCIspeed +#' @docType data +NULL + +#' Activity model class. +#' +#' An S4 class describing activity models fitted to time of observation data. +#' +#' @slot data Object of class \code{"numeric"}, the input data. +#' @slot wt Object of class \code{"numeric"}, weights applied to the data. +#' @slot bw Object of class \code{"numeric"}, kernel bandwidth. +#' @slot adj Object of class \code{"numeric"}, kernel bandwidth adjustment multiplier. +#' @slot pdf Object of class \code{"matrix"} describing fitted probability density function: +#' Column 1: A regular sequence of radian times at which PDF evaluated; range is [0, 2*pi] if unbounded, and sequence steps are range difference divided by 512. +#' Column 2: Corresponding circular kernel PDF values. +#' Additionally if errors bootstrapped: +#' Column 3: PDF standard error. +#' Column 4: PDF lower 95\% confidence limit. Column 5: PDF upper 95\% confidence limit. +#' @slot act Object of class \code{"numeric"} giving activity level estimate and, if errors boostrapped, standard error and 95 percent confidence limits. +#' @export +setClass("actmod", + representation(data="numeric", wt="numeric", bw="numeric", adj="numeric", + pdf="matrix", act="numeric")) + +#' An S4 class describing linear-circular relationships. +#' +#' @slot data Object of class \code{"data.frame"}, the input data, with columns +#' \code{lindat} (linear data) and \code{circdat} (circular data). +#' @slot fit Object of class \code{"data.frame"}, summary of the model fit, with columns: +#' \code{x}: A regular ascending sequence from 0 to 2*pi at which other columns evaluated; +#' \code{fit}: The linear fitted values; +#' \code{p}: The two tailed probability of observing the fitted values under a random (null) circular distribution; +#' \code{nullLCL}: The lower 95\% confidence limit of the null distribution; +#' \code{nullUCL}: The upper 95\% confidence limit of the null distribution. +#' @export +setClass("lincircmod", representation(data="data.frame", fit="data.frame")) + + +#' Index of overlap between circular distributions. +#' +#' Calculates Dhat5 overlap index (see reference) between two kernel distributions . +#' +#' Uses linear interpolation to impute values from kernel distributions. +#' +#' @param fit1,fit2 Fitted activity models of class actmod created using function fitact. +#' @return Scalar overlap index (specifically Dhat5). +#' @references Ridout, M.S. & Linkie, M. (2009) Estimating overlap of daily activity patterns from camera trap data. Journal of Agricultural Biological and Environmental Statistics, 14, 322-337. +#' @examples +#' data(BCItime) +#' oceAct <- fitact(subset(BCItime, species=="ocelot")$time*2*pi) +#' broAct <- fitact(subset(BCItime, species=="brocket")$time*2*pi) +#' ovl5(oceAct, broAct) +#' @export +ovl5 <- function(fit1, fit2){ + f <- stats::approxfun(fit1@pdf[,1], fit1@pdf[,2]) + g <- stats::approxfun(fit2@pdf[,1], fit2@pdf[,2]) + fx <- f(fit1@data) + gx <- g(fit1@data) + fy <- f(fit2@data) + gy <- g(fit2@data) + xr <- gx/fx + yr <- fy/gy + (mean(ifelse(xr>1, 1, xr)) + mean(ifelse(yr>1, 1, yr))) / 2 +} + +#' Calculate circular kernel bandwidth. +#' +#' Uses an optimisation procedure to calculate the circular kernel bandwidth giving the best fit to the data. +#' +#' Mainly for internal use. +#' +#' @param dat Numeric data vector of radian times. +#' @param K Integer number of values of kappa over which to maximise (see references for details). +#' @return Single numeric bandwidth value. +#' @references Ridout, M.S. & Linkie, M. (2009) Estimating overlap of daily activity patterns from camera trap data. Journal of Agricultural Biological and Environmental Statistics, 14, 322-337. +#' @export +bwcalc <- function(dat,K=3) +{ if(!all(dat>=0 & dat<=2*pi)) warning("some dat values are <0 or >2*pi, expecting radian data") + if(max(dat)<1) warning("max(dat) < 1, expecting radian data") + + minfunc <- function(kap,k,dat) + { trigmom <- circular::trigonometric.moment(circular::circular(dat),k,center=T)$rho + (besselI(kap,k)/besselI(kap,0) - trigmom)^2 + } + kapk.calc <- function(k,dat) + stats::optimise(minfunc,c(0,100),k,dat)$minimum + kap <- max(sapply(1:K, kapk.calc, dat)) + ((3*length(dat)*kap^2*besselI(2*kap,2)) / (4*pi^0.5*besselI(kap,0)^2))^(2/5) +} + +#' Circular kernel probability density function. +#' +#' Optionally weighted Von Mises kernel probability densities. +#' +#' If \code{bw} not provided it is calculated internally using \code{bw.calc}. The \code{adj} argument is used to adjust \code{bw} to facilitate exploration of fit flexibility. +#' +#' @param x Numeric vector of radian times at which to evaluate the PDF. +#' @param dat Numeric vector of radian time data to which the PDF is fitted. +#' @param wt A numeric vector of weights for each \code{dat} value. +#' @param bw Numeric value for kernel bandwidth. +#' @param adj Numeric kernel bandwidth multiplier. +#' @return Numeric vector of probability densities evaluated at \code{x}. +#' @seealso \code{\link{density.circular}, \link{bwcalc}} +#' @examples +#' #Example with made up input +#' tt <- runif(100,0,2*pi) +#' xx <- seq(0,2*pi, pi/256) +#' pdf <- dvmkern(xx, tt) +#' plot(xx, pdf, type="l") +#' @export +dvmkern <- function(x,dat,wt=NULL,bw=NULL,adj=1){ + if(!all(dat>=0 & dat<=2*pi)) warning("some dat values are <0 or >2*pi, expecting radian data") + if(max(dat)<1) warning("max(dat) < 1, expecting radian data") + if(!all(x>=0 & x<=2*pi)) warning("some x values are <0 or >2*pi, expecting radian values") + if(!is.null(wt) & length(wt)!=length(dat)) stop("dat and wt have different lengths") + + if(is.null(bw)) bw <- bwcalc(dat) + if(is.null(wt)) wt <- rep(1,length(dat)) + dx <- expand.grid(dat,x) + dif <- abs(dx[,2]-dx[,1]) + i <- dif>pi + dif[i] <- 2*pi-dif[i] + prob <- circular::dvonmises(circular::circular(dif),circular::circular(0),bw*adj) + apply(matrix(prob*wt, nrow=length(dat)),2,sum)/sum(wt) +} + +#' Random numbers from empirical distribution function. +#' +#' Random numbers drawn from an empirical distribution defined by paired values and probabilities. +#' +#' @details The distribution function is defined by \code{fit}, which must be a dataframe containing (at least) columns named: +#' x: a regular sequence of values from which to draw; +#' y: corresponding pdf values. +#' @param n Integer number of random numbers to return. +#' @param fit Data frame defining the emprical distribution (see details). +#' @return A numeric vector. +#' @examples +#' data(BCItime) +#' tm <- 2*pi*subset(BCItime, species=="paca")$time +#' mod <- fitact(tm) +#' rn <- redf(1000, as.data.frame(mod@pdf)) +#' @export +redf <- function(n, fit){ + if(sum(c("x","y") %in% names(fit)) != 2) stop("fit must be a dataframe with (at least) columns named x and y") + if(diff(range(diff(fit$x)))>0.0001) stop("x doesn't seem to be a regular sequence") + + df <- (fit$y[-1]+fit$y[-nrow(fit)])/2 + cdf <- c(0,cumsum(df)/sum(df)) + rn <- stats::runif(n) + stats::approx(cdf, fit$x, rn)$y +} + +#' Modified kernel density function +#' +#' Modifies \code{stats::density} by: +#' Adding SE and 95\% confidence intervals for the density to the output; and +#' Truncating calculation (not just reporting) of density values on from and/or to. +#' +#' @details Truncation copes with cases where where no data are available outside truncation points. +#' Truncation is achieved by fitting the density to the data augmented by reflecting it +#' across each bound using the optimal bandwidth for the unaugmented data, and returning +#' the resulting densities for the region between the bounds. +#' +#' @param x numeric data vector +#' @param reps bootstrap iterations for SE/interval calculation; set to NULL to suppress +#' @param ... Additional arguments passed to \code{stas::density} +#' @return A list with the same components as \code{stats::density} output plus: +#' \code{se}: standard error of the density +#' \code{lcl}, \code{ucl}: lower and upper 95\% confidence intervals of the density + +#' @examples +#' data(BCItime) +#' tm <- subset(BCItime, species=="ocelot")$time +#' dens <- density2(tm, from=0.25, to=0.75) +#' plot(dens$x, dens$y, type="l") +#' @export +density2 <- function(x, reps=999, ...){ + prm <- list(...) + prmnm <- names(prm) + + if(sum(c("from","to") %in% prmnm)==2) + if(prm$from>prm$to) stop("When double-bounded, from must be less than to") + + warn <- FALSE + dmult <- 1 + xx <- x + wt <- prm$weights + if("from" %in% prmnm){ + if(any(xprm$from] + if(!is.null(wt)) prm$weights <- prm$weights[x>prm$from] + } + dmult <- dmult+1 + xx <- c(xx, 2*prm$from-x) + if(!is.null(wt)) wt <- c(wt, prm$weights) + } + if("to" %in% prmnm){ + if(any(x>prm$to)){ + warn <- TRUE + x <- x[xpi/2 & tm2*pi) stop("bounds must be radian (between 0 and 2*pi)") + bdiff <- diff(bounds) + if(diff(bounds)<0){ + dat <- dat-ifelse(dat>pi, 2*pi, 0) + bounds[1] <- bounds[1]-2*pi + bdiff <- diff(bounds) + } + oob <- datbounds[2] + if(sum(oob>0)){ + warning("Some x values outside bounds were removed") + dat <- dat[!oob] + if(!is.null(wt)){ + wt <- wt[!oob] + wt <- wt/sum(wt) + } + } + dens <- density2(dat, from=bounds[1], to=bounds[2], weights=wt, adjust=adj, + bw=if(is.null(bw)) "nrd0" else bw, reps=NULL) + bw <- dens$bw + x <- dens$x + pdf <- dens$y + act <- 1/(bdiff*max(pdf)) + } + + if(sample=="none") + sepdf <- lclpdf <- uclpdf <- seact <- lclact <- uclact <- numeric(0) else{ + if(sample=="model") + samp <- matrix(redf(reps*length(dat), data.frame(x=x,y=pdf)), ncol=reps) else + samp <- matrix(sample(dat, reps*length(dat), replace=TRUE, prob=wt), ncol=reps) + + if(is.null(bounds)){ + if(show) + pdfs <- pbapply::pbapply(samp, 2, function(dat) dvmkern(x,dat,wt,adj=adj)) else + pdfs <- apply(samp, 2, function(dat) dvmkern(x,dat,wt,adj=adj)) + } else + if(show) + pdfs <- pbapply::pbapply(samp, 2, function(dat) density2(dat, from=bounds[1], to=bounds[2], + weights=wt, adjust=adj, bw=if(is.null(bw)) "nrd0" else bw, + reps=NULL)$y) else + pdfs <- apply(samp, 2, function(dat) density2(dat, from=bounds[1], to=bounds[2], weights=wt, + adjust=adj, bw=if(is.null(bw)) "nrd0" else bw, reps=NULL)$y) + + sepdf <- apply(pdfs,1,stats::sd) + lclpdf <- apply(pdfs,1,stats::quantile,probs=0.025) + uclpdf <- apply(pdfs,1,stats::quantile,probs=0.975) + if(is.null(bounds)) acts <- 1/(2*pi*apply(pdfs,2,max)) else + acts <- 1/(bdiff*apply(pdfs,2,max)) + seact <- stats::sd(acts) + lclact <- stats::quantile(acts,0.025) + uclact <- stats::quantile(acts,0.975) + } + + if(is.null(wt)) wt <- 1 + if(min(x)<0){ + dat <- dat+ifelse(dat<0, 2*pi, 0) + x <- x+ifelse(x<0, 2*pi, 0) + } + pdftab <- cbind(x=x, y=pdf, se=sepdf, lcl=lclpdf, ucl=uclpdf)[order(x), ] + + methods::new("actmod", data=dat, wt=wt, bw=bw, adj=adj, pdf=pdftab, + act=c(act=act, se=seact, lcl=lclact, ucl=uclact)) +} + +#' Compare circular distributions. +#' +#' Randomisation test for the probability that two sets of circular observations come from the same distribution. +#' +#' Calculates overlap index Dhat5 (see references) for the two fitted distributions, then generates a null distribution of overlap indices using data sampled randomly with replacement from the combined data. +#' This randomised distribution is then used to define an empirical probability distribution against which the probability that the observed overlap arose by chance is judged. +#' When one or both fitted models use weighted distributions, sampling probabilities are taken from the weights. If both models are weighted, the weights must therefore be on the same scale. +#' +#' @param fit1,fit2 Fitted activity models of class actmod created using function fitact. +#' @param reps Number of bootstrap iterations. +#' @return A named 4-element vector: obs = observed overlap index; null = mean null overlap index; seNull = standard error of the null distribution; pNull = probability observed index arose by chance. +#' @references Ridout, M.S. & Linkie, M. (2009) Estimating overlap of daily activity patterns from camera trap data. Journal of Agricultural Biological and Environmental Statistics, 14, 322-337. +#' @examples +#' #Example with bootstrap reps limited to reduce run time +#' data(BCItime) +#' tPaca <- 2*pi*BCItime$time[BCItime$species=="paca"] +#' tRat <- 2*pi*BCItime$time[BCItime$species=="rat"] +#' fPaca <- fitact(tPaca) +#' fRat <- fitact(tRat) +#' compareCkern(fPaca,fRat,reps=10) +#' @export +compareCkern <- function(fit1, fit2, reps=999){ + if(class(fit1)!="actmod" | class(fit1)!="actmod") + stop("Input must be fitted activity models (class actmod)") + bnd <- range(fit1@pdf[,1]) + if(!all(bnd==range(fit2@pdf[,1]))) + stop("Distribution bounds are not identical") + + if(diff(bnd)==2*pi) bnd <- NULL + olp <- ovl5(fit1,fit2) + y1 <- fit1@data + y2 <- fit2@data + w1 <- fit1@wt + w2 <- fit2@wt + if(length(w1)==1) w1 <- rep(1, length(y1)) + if(length(w2)==1) w2 <- rep(1, length(y2)) + y <- c(y1,y2) + w <- c(w1,w2) + samp <- matrix(sample(1:length(y), reps*length(y), replace=T, prob=w), nrow=reps) + f <- function(s){ + m1 <- fitact(y[s[1:length(y1)]], sample="n", bw=fit1@bw, adj=fit1@adj, bounds=bnd) + m2 <- fitact(y[s[(length(y1)+1):length(y)]], sample="n", bw=fit1@bw, adj=fit1@adj, bounds=bnd) + ovl5(m1,m2) + } + res <- pbapply::pbapply(samp, 1, f) + fun <- stats::ecdf(res) + c(obs=olp, null=mean(res), seNull=stats::sd(res), pNull=fun(olp)) +} + + +#' Compare activity level estimates +#' +#' Wald test for the statistical difference between two or more activitiy level estimates. +#' +#' Uses a Wald test to ask whether the difference between estimates a1 and a2 is +#' significantly different from 0: statistic W = (a1-a2)^2 / (SE1^2+SE2^2) tested +#' on chi-sq distribution with 1 degree of freedom. +#' +#' @param fits A list of fitted \code{actmod} objects +#' @return A matrix with 4 columns: 1. differences between estimates; 2. SEs of the differences; 3. Wald statistics; 4. p-values (H0 is no difference between estimates). Matrix rows give all possible pairwise comparisons, numbered in the order in which they entered in the list \code{fits}. +#' @examples +#' #Test whether paca have a sigificantly different activity level from rat. +#' #Bootstrap reps limited to speed up example. +#' data(BCItime) +#' tPaca <- 2*pi*BCItime$time[BCItime$species=="ocelot"] +#' tRat <- 2*pi*BCItime$time[BCItime$species=="rat"] +#' fPaca <- fitact(tPaca, sample="data", reps=10) +#' fRat <- fitact(tRat, sample="data", reps=10) +#' fPaca@act +#' fRat@act +#' compareAct(list(fPaca,fRat)) +#' @export +compareAct <- function(fits) +{ if(class(fits)!="list" | !all(unlist(lapply(fits,class))=="actmod")) + stop("fits must be a list of actmod objects") + if(min(unlist(lapply(fits, function(x) length(x@act))))==1) + stop("all input model fits must be boostrapped") + + len <- length(fits) + i <- rep(1:(len-1), (len-1):1) + j <- unlist(sapply(2:len, function(i) i:len)) + acts <- unlist(lapply(fits, function(fit) fit@act[1])) + seacts <- unlist(lapply(fits, function(fit) fit@act[2])) + dif <- acts[i]-acts[j] + vardif <- seacts[i]^2 + seacts[j]^2 + W <- dif^2/vardif + prob <- 1-stats::pchisq(W,1) + res <- cbind(Difference=dif, SE=sqrt(vardif), W=W, p=prob) + dimnames(res)[[1]] <- paste(i,j,sep="v") + res +} + +#' Compare activity between times of day +#' +#' Uses a Wald test to statistically compare activity levels at given radian times of day for a fitted activity distribution. +#' +#' Bootrapping the activity model yields standard error estimates for the PDF. This function uses these SEs to compute a Wald statistic for the difference between PDF values (by inference activity levels) at given times of day: statistic W = (a1-a2)^2 / (SE1^2+SE2^2) tested on chi-sq distribution with 1 degree of freedom. +#' +#' @param fit Fitted \code{actmod} object with errors boostrapped (fit using \code{fitact} with \code{sample} argument != "none"). +#' @param times Numeric vector of radian times of day at which to compare activity levels. All pairwise comparisons are made. +#' @return A matrix with 4 columns: 1. differences between PDF values; 2. SEs of the differences; 3. Wald statistics; 4. p-values (H0 is no difference between estimates). Matrix rows give all possible pairwise comparisons, numbered in the order in which they appear in vector \code{times}. +#' @examples +#' data(BCItime) +#' tPaca <- 2*pi*BCItime$time[BCItime$species=="paca"] +#' fPaca <- fitact(tPaca, sample="data", reps=10) +#' plot(fPaca) +#' compareTimes(fPaca, c(5.5,6,0.5,1)) +#' @export +compareTimes <- function(fit, times) +{ if(class(fit)!="actmod") stop("fit input must be an actmod object") + if(!all(times>=0 & times<=2*pi)) stop("some times are <0 or >2*pi, expecting radian data") + + len <- length(times) + i <- rep(1:(len-1), (len-1):1) + j <- unlist(sapply(2:len, function(i) i:len)) + k <- findInterval(times, fit@pdf[,1]) + p <- (times-fit@pdf[,1][k]) / (fit@pdf[,1][k+1]-fit@pdf[,1][k]) + pdfs1 <- fit@pdf[,2][k] + pdfs2 <- fit@pdf[,2][k+1] + pdfs <- pdfs1 + p*(pdfs2-pdfs1) + sepdfs1 <- fit@pdf[,3][k] + sepdfs2 <- fit@pdf[,3][k+1] + sepdfs <- sepdfs1 + p*(sepdfs2-sepdfs1) + dif <- pdfs[i]-pdfs[j] + vardif <- sepdfs[i]^2 + sepdfs[j]^2 + W <- dif^2/vardif + prob <- 1-stats::pchisq(W,1) + res <- cbind(Difference=dif, SE=sqrt(vardif), W=W, p=prob) + dimnames(res)[[1]] <- paste(i,j,sep="v") + res +} + +#' Linear-circular kernel fit +#' +#' Fits a Von Mises kernel distribution describing a linear variable as a function of a circular predictor. +#' +#' @param x Numeric vector of radian values at which to evaluate the distribution. +#' @param circdat Numeric vector of radian data matched with \code{lindat}. +#' @param lindat Numeric vector of linear data matched with \code{circdat}. +#' @return A numeric vector of fitted \code{lindat} values matched with \code{x}. +#' @references Xu, H., Nichols, K. & Schoenberg, F.P. (2011) Directional kernel regression for wind and fire data. Forest Science, 57, 343-352. +#' @examples +#' data(BCIspeed) +#' i <- BCIspeed$species=="ocelot" +#' log_speed <- log(BCIspeed$speed[i]) +#' time <- BCIspeed$time[i]*2*pi +#' circseq <- seq(0,2*pi,pi/256) +#' trend <- lincircKern(circseq, time, log_speed) +#' plot(time, log_speed, xlim=c(0, 2*pi)) +#' lines(circseq, trend) +#' @export +lincircKern <- function(x,circdat,lindat) +{ if(length(lindat)!=length(circdat)) + stop("lindat and circdat lengths are unequal") + if(min(circdat)<0 | max(circdat)>2*pi) + stop("circdat values not between 0 and 2*pi, expecting radian data") + if(min(x)<0 | max(x)>2*pi) + stop("x values not between 0 and 2*pi, expecting radian values") + + hs <- 1.06 * min(stats::sd(lindat), (stats::quantile(lindat,0.75)-stats::quantile(lindat,0.25))/1.34) * + length(lindat)^-0.2 + bw <- 1/hs^2 + dx <- expand.grid(circdat,x) + dif <- abs(dx[,2]-dx[,1]) + i <- dif>pi + dif[i] <- 2*pi-dif[i] + prob <- matrix(circular::dvonmises(circular::circular(dif),circular::circular(0),bw), nrow=length(circdat)) + apply(prob,2,function(z) mean(z*lindat)/mean(z)) +} + +#' Linear-circular regression +#' +#' Fits a Von Mises kernel distribution describing a linear variable as a function +#' of a circular predictor, and boostraps the null distribution in order to evaluate +#' significance of radial variation in the linear variable. +#' +#' Deviation of \code{lindat} from the null expecation is assessed either visually +#' by the degree to which the fitted distribution departs from the null confidence +#' interval (use generic plot function), or quantitatively by column \code{p} of +#' slot \code{fit} in the resulting \code{lincircmod-class} object. +#' +#' @param circdat Numeric vector of radian data matched with \code{lindat}. +#' @param lindat Numeric vector of linear data matched with \code{circdat}. +#' @param pCI Single numeric value between 0 and 1 defining proportional confidence interval to return. +#' @param reps Integer number of bootstrap repetitions to perform. +#' @param res Resolution of fitted distribution and null confidence interval - specifically a single integer number of points on the circular scale at which to record distributions. +#' @return An object of type \code{\link{lincircmod-class}} +#' @references Xu, H., Nichols, K. & Schoenberg, F.P. (2011) Directional kernel regression for wind and fire data. Forest Science, 57, 343-352. +#' @examples +#' #Example with reps limited to increase speed +#' data(BCIspeed) +#' i <- BCIspeed$species=="ocelot" +#' sp <- log(BCIspeed$speed[i]) +#' tm <- BCIspeed$time[i]*2*pi +#' mod <- fitlincirc(tm, sp, reps=50) +#' plot(mod, CircScale=24, xaxp=c(0,24,4), xlab="Time", ylab="log(speed)") +#' legend(8,-3, c("Fitted speed", "Null CI"), col=1:2, lty=1:2) +#' @export +fitlincirc <- function(circdat, lindat, pCI=0.95, reps=10, res=512) +{ if(length(lindat)!=length(circdat)) + stop("lindat and circdat lengths are unequal") + if(min(circdat)<0 | max(circdat)>2*pi) + stop("circdat values not between 0 and 2*pi, expecting radian data") + if(max(circdat)<1) + warning("max(circdat) < 1, expecting radian data") + + n <- length(circdat) + x <- seq(0,2*pi,2*pi/res) + + bs <- pbapply::pbsapply(1:reps, function(i) + { j <- sample(1:n,n,TRUE) + lincircKern(x,circdat[j],lindat) + }) + nulllcl <- apply(bs,1,stats::quantile,(1-pCI)/2) + nullucl <- apply(bs,1,stats::quantile,(1+pCI)/2) + fit <- lincircKern(x,circdat,lindat) + p <- sapply(1:(res+1), function(i) + { f <- stats::ecdf(bs[i,]) + f(fit[i]) + }) + p[p>0.5] <- 1-p[p>0.5] + p <- 2*p + + methods::new("lincircmod", data=data.frame(circdat=circdat, lindat=lindat), + fit=data.frame(x=x, fit=fit, p=p, nullLCL=nulllcl, nullUCL=nullucl)) +} + +#' Plot activity distribution +#' +#' Plot an activity probability distribution from a fitted \code{actmod} object. +#' +#' @details When xunit=="clock", The underlying numeric range of the x-axis is [0,24] if centre=="day", +#' or [-12,12] if centre=="night". +#' @param x Object of class \code{actmod}. +#' @param xunit Character string defining x-axis unit. +#' @param yunit Character string defining y-axis unit. +#' @param data Character string defining whether to plot the data distribution and if so which style to use. +#' @param centre Character string defining whether to centre the plot on midday or midnight. +#' @param dline List of plotting parameters for data lines. +#' @param tline List of plotting parameters for trend line. +#' @param cline List of plotting parameters for trend confidence interval lines. +#' @param add Logical defining whether to create a new plot (default) or add to an existing plot. +#' @param xaxis List of plotting parameters to pass to axis command for x-axis plot (see axis for arguments). +#' @param ... Additional arguments passed to internal plot call affecting only the plot frame and y axis. Modify x axis through xaxis. +#' @examples +#' data(BCItime) +#' otm <- 2*pi*subset(BCItime, species=="ocelot")$time +#' btm <- 2*pi*subset(BCItime, species=="brocket")$time +#' omod <- fitact(otm) +#' bmod <- fitact(btm) +#' plot(omod, yunit="density", data="none") +#' plot(bmod, yunit="density", data="none", add=TRUE, tline=list(col="red")) +#' legend("topleft", c("Ocelot", "Brocket deer"), col=1:2, lty=1) +#' +#' mod <- fitact(otm, sample="data", reps=10) +#' plot(mod, dline=list(col="grey"), +#' tline=list(col="red", lwd=2), +#' cline=list(col="red", lty=3)) +#' +#' mod2 <- fitact(otm, bounds=c(pi*3/2, pi/2)) +#' plot(mod2, centre="night") +#' plot(mod2, centre="night", xlim=c(-6,6), xaxis=list(at=seq(-6,6,2))) +#' @export +plot.actmod <- function(x, xunit=c("clock","hours","radians"), yunit=c("frequency","density"), + data=c("histogram","rug","both","none"), centre=c("day","night"), + dline=list(lwd=ifelse(data=="rug", 0.1, 1)), tline=NULL, cline=list(lty=2), + add=FALSE, xaxis=NULL, ...){ + + #function sets up the plot parameters with appropriate scales and labels with flexibility for modifcation using ... + setup <- function(...){ + pprm <- list(...) + pprm <- append(list(x=xlim, y=ylim, type="n", xaxt="n"), pprm) + if(!("ylim" %in% names(pprm))) pprm <- append(pprm, list(ylim=ylim)) + if(!("ylab" %in% names(pprm))) + pprm <- switch(yunit, + "frequency"=append(pprm, list(ylab="Frequency")), + "density"=append(pprm, list(ylab="Density")) + ) + if(!("xlab" %in% names(pprm))) + pprm <- switch(xunit, + "clock"=append(pprm, list(xlab="Time")), + "hours"=append(pprm, list(xlab="Hours")), + "radians"=append(pprm, list(xlab="Radians")) + ) + do.call(graphics::plot, pprm) + xaxis <- append(list(side=1), xaxis[!(names(xaxis)=="side")]) + if(xunit=="clock"){ + if(!("at" %in% names(xaxis))){ + at <- switch(centre, + "day"=seq(0, 24, 6), + "night"=seq(-12, 12, 6) + ) + xaxis <- append(xaxis, list(at=at)) + } + if(!("labels" %in% names(xaxis))){ + at <- xaxis$at + hh <- trunc(at) + mm <- round(60*(at-trunc(at)), 0) + hh <- hh + ifelse(hh<0,24,0) - ifelse(mm<0, 1, 0) + mm <- mm + ifelse(mm<0, 60, 0) + lab <- paste0(ifelse(hh<10,"0",""), hh, ":", ifelse(mm<10,"0",""), mm) + xaxis <- append(xaxis, list(labels=lab)) + } + } + do.call(graphics::axis, xaxis) + } + ########## end setup function ############# + + data <- match.arg(data) + xunit <- match.arg(xunit) + yunit <- match.arg(yunit) + centre <- match.arg(centre) + if(!"lty" %in% names(cline)) cline <- append(cline, list(lty=2)) + + fit <- x + pdf <- fit@pdf + fdata <- fit@data + + if(centre=="night"){ + if(diff(range(pdf[,1]))==2*pi | min(pdf[,1]>=0)){ + if(diff(range(pdf[,1]))==2*pi) pdf <- pdf[-1,] + pdf[,1] <- pdf[,1] - ifelse(pdf[,1]>pi, 2*pi, 0) + pdf <- pdf[order(pdf[,1]), ] + } + fdata <- fdata-ifelse(fdata>pi,2*pi,0) + xlim <- c(-pi,pi) + }else xlim <- c(0,2*pi) + + x <- pdf[,"x"] + y <- pdf[,"y"] + if(ncol(pdf)==5) + { lcl <- pdf[,"lcl"] + ucl <- pdf[,"ucl"] + } else lcl <- ucl <- numeric(0) + + if(xunit=="radians") maxbrk <- 2*pi else { + xlim <- xlim*12/pi + x <- x*12/pi + fdata <- fdata*12/pi + maxbrk <- 24 + } + + if(yunit=="frequency"){ + y <- y*length(fdata)*pi/12 + lcl <- lcl*length(fdata)*pi/12 + ucl <- ucl*length(fdata)*pi/12 + }else{ + if(xunit %in% c("clock","hours")){ + y <- y*pi/12 + lcl <- lcl*pi/12 + ucl <- ucl*pi/12 + } + } + + if(data %in% c("histogram","both")){ + h <- switch(centre, + "day"=graphics::hist(fdata, breaks=seq(0,maxbrk,maxbrk/24), plot=F), + "night"=graphics::hist(fdata, breaks=seq(-maxbrk/2,maxbrk/2,maxbrk/24), plot=F)) + d <- switch(yunit, + "frequency"=h$counts, + "density"=h$density) + ylim <- c(0,max(y,d,ucl,na.rm=T)) + } else ylim <- c(0,max(y,ucl,na.rm=T)) + + #Plot data + if(!add) setup(...) + if(data %in% c("histogram","both")){ + if(!"lwd" %in% names(dline)) dline <- append(dline, list(lwd=1)) + do.call(graphics::lines, append(list(x=h$breaks, y=c(d,d[1]), type="s"), dline)) + } + if(data %in% c("rug","both")){ + if(!"lwd" %in% names(dline)) dline <- append(dline, list(lwd=0.1)) + for(i in 1:length(fdata)) + do.call(graphics::lines, append(list(x=rep(fdata[i],2), y=max(y,na.rm=T)*-c(0,0.03)), dline)) + } + + #Plot trend + i <- x < switch(centre, "day"=ifelse(xunit=="radians", pi, 12), 0) + do.call(graphics::lines, append(list(x=x[i], y=y[i]), tline)) + do.call(graphics::lines, append(list(x=x[!i], y=y[!i]), tline)) + + #Plot conf intervals + if(length(lcl)>0) + { do.call(graphics::lines, append(list(x=x[i], y=lcl[i]), cline)) + do.call(graphics::lines, append(list(x=x[i], y=ucl[i]), cline)) + do.call(graphics::lines, append(list(x=x[!i], y=lcl[!i]), cline)) + do.call(graphics::lines, append(list(x=x[!i], y=ucl[!i]), cline)) + } +} + +#' Plot linear-circular relationship +#' +#' Plot linear against circular data along with the fitted and null confidence limit distributions from a fitted \code{lincircmod} object. +#' +#' @param x Object of class \code{lincircmod}. +#' @param CircScale Single numeric value defining the plotting maximum of the circular scale. +#' @param tlim Numeric vector with two elements >=0 and <=1 defining the lower and upper limits at which to plot distributions; default plots the full range. +#' @param fcol,flty,ncol,nlty Define line colour (\code{col}) and type (\code{lty}) for fitted (\code{f}) and null (\code{n}) distributions; input types as for \code{col} and \code{lty}, see \code{\link{par}}. +#' @param ... Additional arguments passed to the inital plot construction, affecting axes and data plot symbols. +#' @export +plot.lincircmod <- function(x, CircScale=2*pi, tlim=c(0,1), fcol="black", flty=1, ncol="red", nlty=2, ...){ + + if(min(tlim)<0 | max(tlim)>1 | length(tlim)!=2) stop("tlim should contain two values >=0 and <=1") + + fit <- x + fit <- x@fit + dat <- x@data + xx <- fit$x*CircScale/(2*pi) + LinearData <- dat$lindat + CircularData <- dat$circdat*CircScale/(2*pi) + range <- tlim*CircScale + graphics::plot(CircularData, LinearData, ...) + if(range[1]=range[1] & xx<=range[2] + graphics::lines(xx[i], fit$fit[i], col=fcol, lty=flty) + graphics::lines(xx[i], fit$nullLCL[i], col=ncol, lty=nlty) + graphics::lines(xx[i], fit$nullUCL[i], col=ncol, lty=nlty) + } else + { i <- xx>=range[1] + graphics::lines(xx[i], fit$fit[i], col=fcol, lty=flty) + graphics::lines(xx[i], fit$nullLCL[i], col=ncol, lty=nlty) + graphics::lines(xx[i], fit$nullUCL[i], col=ncol, lty=nlty) + i <- xx<=range[2] + graphics::lines(xx[i], fit$fit[i], col=fcol, lty=flty) + graphics::lines(xx[i], fit$nullLCL[i], col=ncol, lty=nlty) + graphics::lines(xx[i], fit$nullUCL[i], col=ncol, lty=nlty) + } +} diff --git a/R/activity_code_1_1.r b/R/activity_code_1_1.r deleted file mode 100644 index a10b27c..0000000 --- a/R/activity_code_1_1.r +++ /dev/null @@ -1,669 +0,0 @@ -#' Animal activity statistics -#' -#' Provides functions to estimate and compare activity parameters from sensor data. -#' -#' @details Sensors that record active animals (eg camera traps) build up a record of -#' the distribution of activity over the course of the day. Records are more frequent -#' when animals are more active, and less frequent or absent when animals are inactive. -#' The area under the distribution of records thus contains information on the overall -#' level of activity in a sampled population. This package provides tools for plotting -#' activity distributions, quantifying the overall level of activity with error, and -#' statistically comparing distributions through bootstrapping. -#' -#' The core function is \code{fitact}, which creates an \code{actmod} object containing -#' the circular kernel PDF, and the activity level estimate derived from this. The -#' generic plot function for \code{actmod} objects plots the distribution. Functions -#' starting with \code{compare} make statistical comparisons between distributions or -#' activity estimates. Note that all time or other circular data should be in radians -#' (in the range 0 to 2*pi). -#' -#' @references Rowcliffe, M., Kays, R., Kranstauber, B., Carbone, C., Jansen, P.A. (2014) Quantifying animal activity level using camera trap data. Methods in Ecology and Evolution. -#' @seealso \code{\link{overlap}} -#' @docType package -#' @name activity -NULL - -#' Animal record time of day data -#' -#' Barro Colorado Island 2008 data: times of day at which animal records occured -#' (\code{time}), together with species (\code{species}). -#' -#' @format A dataframe with 17820 observations and 2 variables. -#' @source http://dx.doi.org/10.6084/m9.figshare.1160536 -#' @name BCItime -#' @docType data -NULL - -#' Animal speed data -#' -#' Barro Colorado Island 2008 data: speeds of animal passages past camera traps -#' (\code{speed}), together with species (\code{species}) and time of day (\code{time}) -#' for each record. -#' -#' @format A dataframe with 2204 observations and 3 variables. -#' @source http://dx.doi.org/10.6084/m9.figshare.1160536 -#' @name BCIspeed -#' @docType data -NULL - -#' Activity model class. -#' -#' An S4 class describing activity models fitted to time of observation data. -#' -#' @slot data Object of class \code{"numeric"}, the input data. -#' @slot wt Object of class \code{"numeric"}, weights applied to the data. -#' @slot bw Object of class \code{"numeric"}, kernel bandwidth. -#' @slot adj Object of class \code{"numeric"}, kernel bandwidth adjustment multiplier. -#' @slot pdf Object of class \code{"matrix"} describing fitted probability density function: -#' Column 1: Sequence of radian times at which PDF evaluated (specifically seq(0,2*pi,pi/256)). -#' Column 2: Corresponding circular kernel PDF values. -#' Additionally if errors boostrapped: -#' Column 3: PDF standard error. -#' Column 4: PDF lower confidence interval. Column 5: PDF upper confidence interval. -#' @slot act Object of class \code{"numeric"} giving activity level estimate and, if errors boostrapped, standard error and 95 percent confidence limits. -#' @method plot \code{signature(x = "actmod")}: Plots PDF, confidence interval if calculated, and optionally data distribution. -#' @export -setClass("actmod", - representation(data="numeric", wt="numeric", bw="numeric", adj="numeric", - pdf="matrix", act="numeric")) - -#' An S4 class describing linear-circular relationships. -#' -#' @slot data Object of class \code{"data.frame"}, the input data, with columns -#' \code{lindat} (linear data) and \code{circdat} (circular data). -#' @slot fit Object of class \code{"data.frame"}, summary of the model fit, with columns: -#' \code{x}: A regular ascending sequence from 0 to 2*pi at which other columns evaluated; -#' \code{fit}: The linear fitted values; -#' \code{p}: The two tailed probability of observing the fitted values under a random (null) circular distribution; -#' \code{nullLCL}: The lower confidence limit of the null distribution; -#' \code{nullUCL}: The upper confidence limit of the null distribution. -#' @method plot \code{signature(x = "actmod")}: Plots PDF, confidence interval if calculated, and optionally data distribution. -#' @export -setClass("lincircmod", representation(data="data.frame", fit="data.frame")) - -#' Impute empirical circular distribution. -#' -#' Imputes values at given points on an empirical circular distribution. -#' -#' Note that x is assumed circular, so first and last \code{y} values should be equal. Evaluation points \code{xvals} should also be within the range of \code{x}. -#' -#' @param xvals Numeric circular values at which to evaluate \code{y} -#' @param x Evenly spaced ascending numeric sequence of circular values -#' @param y Empirical numeric output distribution matched with \code{x} -#' @return A numeric vector of \code{y} values evaluated at \code{xvals} -#' @examples -#' #Abstract example -#' x <- seq(0,2*pi,length.out=11) -#' y <- c(0,1,2,3,4,5,4,3,2,1,0) -#' yfromx(0:6,x,y) -#' -#' #BCI data example -#' #Weighting ocelot activity pattern to correct for variation in speed -#' data(BCIspeed) -#' data(BCItime) -#' #Fit linear-circular model to log(speed) -#' i <- BCIspeed$species=="ocelot" -#' lcfit <- fitlincirc(BCIspeed$time[i]*2*pi, log(BCIspeed$speed[i]), reps=50) -#' #Fit weighted activity model using yfromx to create weights -#' j <- BCItime$species=="ocelot" -#' tdat <- BCItime$time[j]*2*pi -#' w <- 1/yfromx(tdat, lcfit@@fit$x, exp(lcfit@@fit$fit)) -#' mod <- fitact(tdat, wt=w, sample="none") -#' plot(mod) -#' #Oveplot unweighted model for comparison -#' mod2 <- fitact(tdat, sample="none") -#' plot(mod2, lcol=3, add=TRUE) -#' @export -yfromx <- function(xvals, x, y) -{ if(length(x)!=length(y)) stop("x and y lengths are unequal") - if(max(diff(x))-min(diff(x))>1e-10) stop("x is not an evenly spaced ascending sequence") - if(abs(y[1]-y[length(y)])>1e-10) stop("first and last y values should be equal") - if(min(xvals)max(x)) stop("some xvals are not within the range of x") - - unit <- x[2] - xvals[xvals==0] <- max(x) - i <- ceiling(xvals/unit) - y[i] + (y[i+1]-y[i])*(xvals-x[i])/unit -} - -#' Calculate circular kernel bandwidth. -#' -#' Uses an optimisation procedure to calculate the circular kernel bandwidth giving the best fit to the data. -#' -#' Mainly for internal use. -#' -#' @param dat Numeric data vector of radian times. -#' @param K Integer number of values of kappa over which to maximise (see references for details). -#' @return Single numeric bandwidth value. -#' @references Ridout, M.S. & Linkie, M. (2009) Estimating overlap of daily activity patterns from camera trap data. Journal of Agricultural Biological and Environmental Statistics, 14, 322-337. -#' @export -bwcalc <- function(dat,K=3) -{ if(!all(dat>=0 & dat<=2*pi)) warning("some dat values are <0 or >2*pi, expecting radian data") - if(max(dat)<1) warning("max(dat) < 1, expecting radian data") - - minfunc <- function(kap,k,dat) - { trigmom <- circular::trigonometric.moment(circular::circular(dat),k,center=T)$rho - (besselI(kap,k)/besselI(kap,0) - trigmom)^2 - } - kapk.calc <- function(k,dat) - optimise(minfunc,c(0,100),k,dat)$minimum - kap <- max(sapply(1:K, kapk.calc, dat)) - ((3*length(dat)*kap^2*besselI(2*kap,2)) / (4*pi^0.5*besselI(kap,0)^2))^(2/5) -} - -#' Circular kernel probability density function. -#' -#' Optionally weighted Von Mises kernel probability densities. -#' -#' If \code{bw} not provided it is calculated internally using \code{bw.calc}. The \code{adj} argument is used to adjust \code{bw} to facilitate exploration of fit flexibility. -#' -#' @param x Numeric vector of radian times at which to evaluate the PDF. -#' @param dat Numeric vector of radian time data to which the PDF is fitted. -#' @param wt A numeric vector of weights for each \code{dat} value. -#' @param bw Numeric value for kernel bandwidth. -#' @param adj Numeric kernel bandwidth multiplier. -#' @return Numeric vector of probability densities evaluated at \code{x}. -#' @seealso \code{\link{density.circular}, \link{bwcalc}} -#' @examples -#' #Example with made up input -#' tt <- runif(100,0,2*pi) -#' xx <- seq(0,2*pi, pi/256) -#' pdf <- dvmkern(xx, tt) -#' plot(xx, pdf, type="l") -#' @export -dvmkern <- function(x,dat,wt=NULL,bw=NULL,adj=1) -{ if(!all(dat>=0 & dat<=2*pi)) warning("some dat values are <0 or >2*pi, expecting radian data") - if(max(dat)<1) warning("max(dat) < 1, expecting radian data") - if(!all(x>=0 & x<=2*pi)) warning("some x values are <0 or >2*pi, expecting radian values") - if(!is.null(wt) & length(wt)!=length(dat)) stop("dat and wt have different lengths") - - if(is.null(bw)) bw <- bwcalc(dat) - if(is.null(wt)) wt <- rep(1,length(dat)) - dx <- expand.grid(dat,x) - dif <- abs(dx[,2]-dx[,1]) - i <- dif>pi - dif[i] <- 2*pi-dif[i] - prob <- circular::dvonmises(circular::circular(dif),circular::circular(0),bw*adj) - apply(matrix(prob*wt, nrow=length(dat)),2,sum)/sum(wt) -} - -#' Random circular kernel numbers. -#' -#' Random numbers drawn from a fitted Von Mises kernel distribution. -#' -#' @details If \code{fit} is a matrix, it should have two columns: -#' [,1] precisely seq(0,2*pi,pi/256) sequence of radian values at which pdf evaluated; -#' [,2] corresponding pdf values. -#' @param n Integer number of numbers to return. -#' @param fit An emprical distribution contained in either a matrix or a fitted \code{actmod} object (see Details). -#' @return A numeric vector of radian numbers. -#' @examples -#' #Matrix input -#' data(BCItime) -#' tm <- 2*pi*BCItime$time[BCItime$species=="paca"] -#' xx <- seq(0,2*pi, pi/256) -#' rn <- rvmkern(100, cbind(xx, dvmkern(xx, tm))) -#' -#' #actmod input -#' fit <- fitact(tm, sample="n") -#' rvmkern(100,fit) -#' @export -rvmkern <- function(n,fit) -{ fclass <- class(fit) - if(fclass=="actmod") fit <- fit@pdf else - if(fclass=="matrix") - { if(nrow(fit)!=513 | ncol(fit)!=2) - stop("fit matrix should have exactly 513 rows and 2 columns") - if(max(abs(seq(0,2*pi,pi/256)-fit[,1]))>0) - stop("first column of fit matrix should be seq(0,2*pi,pi/256)") - if(abs(1-sum(fit[,2])*2*pi/512)>0.01) - warning("second column of fit matrix doesn't look like a pdf and was rescaled to sum to 1") - } else - stop("fit should be of either matrix or actmod class") - cumpdf <- c(0,cumsum(fit[-513,2]*2*pi/512)) - cumpdf <- cumpdf/cumpdf[513] - rn <- runif(n) - i <- findInterval(rn,cumpdf) - fit[i,1] + ((rn-cumpdf[i])/(cumpdf[i+1]-cumpdf[i]))*2*pi/512 -} - -#' Fit activity model to time-of-day data -#' -#' Fits a circular kernel density to radian time-of-day data and estimates activity -#' level from this distribution. Optionally bootstraps the distribution, in which -#' case SEs and confidence limits are also stored for activity level and PDF. -#' -#' @details The bandwidth adjustment multiplier \code{adj} is provided to allow -#' exploration of the effect of adjusting the internally calculated bandwidth on -#' accuracy of activity level estimates. The alternative bootstrapping methods -#' defined by \code{sample} are: -#' \code{data}: sample from the data; -#' \code{model}: sample from the fitted probability density distribution; -#' \code{none}: no bootstrapping. -#' Confidence interval coverage seems to be better at large sample size -#' (greater than 100-200) using \code{"model"}, but better at small sample size -#' when using \code{"data"}. The reason for this needs further investigation. -#' @param dat A numeric vector of radian time-of-day data. -#' @param wt A numeric vector of weights for each \code{dat} value. -#' @param reps Number of boostrap iterations to perform. Ignored if sample=="none". -#' @param bw Numeric value for kernel bandwidth. If NULL, calculated internally. -#' @param adj Numeric bandwidth adjustment multiplier. -#' @param sample Character string defining sampling method for bootstrapping errors (see details). -#' @param show Logical whether or not to show a progress bar while bootstrapping. -#' @return An object of type \code{actmod} -#' @examples -#' #Fit without confidence limits -#' data(BCItime) -#' tdat <- 2*pi*BCItime$time[BCItime$species=="ocelot"] -#' mod1 <- fitact(tdat, sample="none") -#' plot(mod1) -#' -#' #Fit with confidence limits (limited reps to speed up) -#' mod2 <- fitact(tdat, reps=10) -#' plot(mod2) -#' -#' #Fit weighted function to correct for detection radius 1.21 times higher -#' #by day than by night, assuming day between pi/2 (6am) and pi*2/3 (6pm) -#' weight <- 1/ifelse(tdat>pi/2 & tdat=0 & dat<=2*pi)) warning("some dat values are <0 or >2*pi, expecting radian data") - if(max(dat)<1) warning("max(dat) < 1, expecting radian data") - if(!is.null(wt) & length(wt)!=length(dat)) stop("dat and wt have different lengths") - - sample <- match.arg(sample) - if(is.null(bw)) bw <- bwcalc(dat) - x <- seq(0,2*pi,pi/256) - pdf <- dvmkern(x, dat, wt, adj, bw) - act <- 1/(2*pi*max(pdf)) - if(sample=="none") - sepdf <- lclpdf <- uclpdf <- seact <- lclact <- uclact <- numeric(0) else - { if(sample=="model") samp <- matrix(rvmkern(reps*length(dat), cbind(x,pdf)), ncol=reps) else - if(sample=="data") samp <- matrix(sample(dat, reps*length(dat), replace=TRUE), ncol=reps) - if(show) - pdfs <- pbapply::pbapply(samp, 2, function(dat) dvmkern(x,dat,wt,adj=adj)) else - pdfs <- apply(samp, 2, function(dat) dvmkern(x,dat,wt,adj=adj)) - sepdf <- apply(pdfs,1,sd) - lclpdf <- apply(pdfs,1,quantile,probs=0.025) - uclpdf <- apply(pdfs,1,quantile,probs=0.975) - acts <- 1/(2*pi*apply(pdfs,2,max)) - seact <- sd(acts) - lclact <- quantile(acts,0.025) - uclact <- quantile(acts,0.975) - } - if(is.null(wt)) wt <- 1 - new("actmod", data=dat, wt=wt, bw=bw, adj=adj, - pdf=cbind(time=x, pdf=pdf, se=sepdf, lcl=lclpdf, ucl=uclpdf), - act=c(act=act, se=seact, lcl=lclact, ucl=uclact)) -} - -#' Compare circular distributions. -#' -#' Randomisation test for the probability that two sets of circular observations come from the same distribution. -#' -#' Calculates overlap index (see references) for the observed data samples, then generates a null distribution of overlap indices using data sampled randomly with replacement from the combined data. This randomised distribution is then used to estimate the probability that the observed overlap arose by chance. -#' -#' @param y1,y2 Numeric vectors of radian data. -#' @param reps Number of bootstrap iterations. -#' @param index Which of the three indices returned by overlap::overlapEst to use. -#' @return A named 2-element vector: Overlap = observed overlap index; p = probability observed index arose by chance. -#' @references Ridout, M.S. & Linkie, M. (2009) Estimating overlap of daily activity patterns from camera trap data. Journal of Agricultural Biological and Environmental Statistics, 14, 322-337. -#' @seealso \code{\link{overlapEst}} -#' @examples -#' #Example with bootstrap reps limited to speed up -#' data(BCItime) -#' tPaca <- 2*pi*BCItime$time[BCItime$species=="paca"] -#' tRat <- 2*pi*BCItime$time[BCItime$species=="rat"] -#' compareCkern(tPaca,tRat,reps=10) -#' @export -compareCkern <- function(y1,y2,reps=1000,index=c("Dhat4","Dhat5","Dhat1")) -{ index <- match.arg(index) - i <- switch(index, - "Dhat1"=1, - "Dhat4"=2, - "Dhat5"=3) - ovl <- overlap::overlapEst(y1,y2)[i] - names(ovl) <- NULL - samp <- matrix(sample(c(y1,y2), reps*length(c(y1,y2)), replace=TRUE), ncol=reps) - res <- pbapply::pbsapply(1:reps, function(j) - overlap::overlapEst(samp[1:length(y1),j], samp[length(y1)+(1:length(y2)),j])[i]) - fun <- ecdf(res) - c(Overlap=ovl, p=fun(ovl)) -} - -#' Compare activity level estimates -#' -#' Wald test for the statistical difference between two or more activitiy level estimates. -#' -#' Uses a Wald test to ask whether the difference between estimates a1 and a2 is -#' significantly different from 0: statistic W = (a1-a2)^2 / (SE1^2+SE2^2) tested -#' on chi-sq distribution with 1 degree of freedom. -#' -#' @param fits A list of fitted \code{actmod} objects -#' @return A matrix with 4 columns: 1. differences between estimates; 2. SEs of the differences; 3. Wald statistics; 4. p-values (H0 is no difference between estimates). Matrix rows give all possible pairwise comparisons, numbered in the order in which they entered in the list \code{fits}. -#' @examples -#' #Test whether paca have a sigificantly different activity level from rat. -#' #Bootstrap reps limited to speed up example. -#' data(BCItime) -#' tPaca <- 2*pi*BCItime$time[BCItime$species=="paca"] -#' tRat <- 2*pi*BCItime$time[BCItime$species=="rat"] -#' (fPaca <- fitact(tPaca, reps=10)) -#' (fRat <- fitact(tRat, reps=10)) -#' compareAct(list(fPaca,fRat)) -#' @export -compareAct <- function(fits) -{ if(class(fits)!="list" | !all(unlist(lapply(fits,class))=="actmod")) - stop("fits must be a list of actmod objects") - if(min(unlist(lapply(fits, function(x) length(x@act))))==1) - stop("all input model fits must be boostrapped") - - len <- length(fits) - i <- rep(1:(len-1), (len-1):1) - j <- unlist(sapply(2:len, function(i) i:len)) - acts <- unlist(lapply(fits, function(fit) fit@act[1])) - seacts <- unlist(lapply(fits, function(fit) fit@act[2])) - dif <- acts[i]-acts[j] - vardif <- seacts[i]^2 + seacts[j]^2 - W <- dif^2/vardif - prob <- 1-pchisq(W,1) - res <- cbind(Difference=dif, SE=sqrt(vardif), W=W, p=prob) - dimnames(res)[[1]] <- paste(i,j,sep="v") - res -} - -#' Compare activity across between times of day -#' -#' Uses a Wald test to statistically compare activity levels at given radian times of day for a fitted activity distribution. -#' -#' Bootrapping the activity model yields standard error estimates for the PDF. This function uses these SEs to compute a Wald statistic for the difference between PDF values (by inference activity levels) at given times of day: statistic W = (a1-a2)^2 / (SE1^2+SE2^2) tested on chi-sq distribution with 1 degree of freedom. -#' -#' @param fit Fitted \code{actmod} object with errors boostrapped (fit using \code{fitact} with \code{sample} argument != "none"). -#' @param times Numeric vector of radian times of day at which to compare activity levels. All pairwise comparisons are made. -#' @return A matrix with 4 columns: 1. differences between PDF values; 2. SEs of the differences; 3. Wald statistics; 4. p-values (H0 is no difference between estimates). Matrix rows give all possible pairwise comparisons, numbered in the order in which they appear in vector \code{times}. -#' @examples -#' data(BCItime) -#' tPaca <- 2*pi*BCItime$time[BCItime$species=="paca"] -#' fPaca <- fitact(tPaca, reps=10) -#' plot(fPaca, hrs=FALSE, frq=FALSE) -#' compareTimes(fPaca, c(5.5,6,0.5,1)) -#' @export -compareTimes <- function(fit, times) -{ if(class(fit)!="actmod") stop("fit input must be an actmod object") - if(!all(times>=0 & times<=2*pi)) stop("some times are <0 or >2*pi, expecting radian data") - - len <- length(times) - i <- rep(1:(len-1), (len-1):1) - j <- unlist(sapply(2:len, function(i) i:len)) - k <- findInterval(times, fit@pdf[,1]) - p <- (times-fit@pdf[,1][k]) / (fit@pdf[,1][k+1]-fit@pdf[,1][k]) - pdfs1 <- fit@pdf[,2][k] - pdfs2 <- fit@pdf[,2][k+1] - pdfs <- pdfs1 + p*(pdfs2-pdfs1) - sepdfs1 <- fit@pdf[,3][k] - sepdfs2 <- fit@pdf[,3][k+1] - sepdfs <- sepdfs1 + p*(sepdfs2-sepdfs1) - dif <- pdfs[i]-pdfs[j] - vardif <- sepdfs[i]^2 + sepdfs[j]^2 - W <- dif^2/vardif - prob <- 1-pchisq(W,1) - res <- cbind(Difference=dif, SE=sqrt(vardif), W=W, p=prob) - dimnames(res)[[1]] <- paste(i,j,sep="v") - res -} -#' Show activity level estimate -#' -#' Prints the \code{act} slot (activity level estimate) from an \code{actmod} object. -#' -#' @param object Object of class \code{actmod}. -#' @export -setMethod("show", "actmod", function(object) print(object@act)) - -#' Plot activity distribution -#' -#' Plot an activity probability distribution from a fitted \code{actmod} object. -#' -#' @param x Object of class \code{actmod}. -#' @param xunit Character string defining x-axis unit. -#' @param yunit Character string defining y-axis unit. -#' @param dat Character string defining whether to plot the data distribution and if so which style to use. -#' @param add Logical defining whether to create a new plot (default) or add the probability density to an existing plot (in which case no data are plotted). -#' @param dcol Numeric or character defining colour of data lines. -#' @param lcol Numeric or character defining colour of PDF lines. -#' @param ... Additional arguments passed to internal plot call affecting only axes. -#' @docType methods -#' @rdname plot-actmod-methods -#' @export -setMethod("plot", "actmod", - function(x, xunit=c("hours","radians"), yunit=c("frequency","density"), dat=c("histogram","rug","none"), add=FALSE, dcol=1, lcol=2, ...){ - - #setup function sets up the plot with appropriate scales and labels with flexibility for midifcation using ... - setup <- function(x,y,xunit,yunit){ - if(xunit=="hours"){ - Hours <- x - if(yunit=="frequency"){ - Frequency<-y - if("ylim" %in% names(list(...))){ - if("xaxp" %in% names(list(...))) - plot(Hours,Frequency, type="n", ...) else - plot(Hours,Frequency, type="n", xaxp=c(0,24,4), ...) - }else{ - if("xaxp" %in% names(list(...))) - plot(Hours,Frequency, type="n", ylim=c(0,max(y,d,ucl)), ...) else - plot(Hours,Frequency, type="n", xaxp=c(0,24,4), ylim=c(0,max(y,d,ucl)), ...) - } - }else{ - Density<-y - if("ylim" %in% names(list(...))){ - if("xaxp" %in% names(list(...))) - plot(Hours,Density, type="n", ...) else - plot(Hours,Density, type="n", xaxp=c(0,24,4), ...) - }else{ - if("xaxp" %in% names(list(...))) - plot(Hours,Density, type="n", ylim=c(0,max(y,d,ucl)), ...) else - plot(Hours,Density, type="n", xaxp=c(0,24,4), ylim=c(0,max(y,d,ucl)), ...) - } - } - }else{ - Radians <- x - if(yunit=="frequency"){ - Frequency<-y - if("ylim" %in% names(list(...))) - plot(Radians,Frequency, type="n", ...) else - plot(Radians,Frequency, type="n", ylim=c(0,max(y,d,ucl)), ...) - }else{ - Density<-y - if("ylim" %in% names(list(...))) - plot(Radians,Density, type="n", ...) else - plot(Radians,Density, type="n", ylim=c(0,max(y,d,ucl)), ...) - } - } - } - - xunit <- match.arg(xunit) - yunit <- match.arg(yunit) - dat <- match.arg(dat) - fit <- x - x <- fit@pdf[,1] - y <- fit@pdf[,2] - if(ncol(fit@pdf)==5) - { lcl <- fit@pdf[,4] - ucl <- fit@pdf[,5] - } else lcl <- ucl <- numeric(0) - data <- fit@data - if(xunit=="hours") - { x <- x*12/pi - data <- data*12/pi - maxbrk <- 24 - }else - { maxbrk <- 2*pi - } - h <- hist(data, breaks=seq(0,maxbrk,maxbrk/24), plot=F) - if(yunit=="frequency") - { y <- y*length(data)*pi/12 - lcl <- lcl*length(data)*pi/12 - ucl <- ucl*length(data)*pi/12 - d <- h$counts - }else - { if(xunit=="hours") - { y <- y*pi/12 - lcl <- lcl*pi/12 - ucl <- ucl*pi/12 - } - d <- h$density - } - if(!add) - { setup(range(x), c(0,max(y,d,ucl)), xunit, yunit) - if(dat=="histogram") - lines(h$breaks, c(d,d[1]), type="s", col=dcol) else - if(dat=="rug") - for(i in 1:length(data)) - lines(rep(data[i],2),max(y)*-c(0,0.03), lwd=0.1, col=dcol) - } - lines(x, y, col=lcol) - if(length(lcl)>0) - { lines(x, lcl, col=lcol, lty=2) - lines(x, ucl, col=lcol, lty=2) - } - } -) - -#' Linear-circular kernel fit -#' -#' Fits a Von Mises kernel distribution describing a linear variable as a function of a circular predictor. -#' -#' @param x Numeric vector of radian values at which to evaluate the distribution. -#' @param circdat Numeric vector of radian data matched with \code{lindat}. -#' @param lindat Numeric vector of linear data matched with \code{circdat}. -#' @return A numeric vector of fitted \code{lindat} values matched with \code{x}. -#' @references Xu, H., Nichols, K. & Schoenberg, F.P. (2011) Directional kernel regression for wind and fire data. Forest Science, 57, 343-352. -#' @examples -#' data(BCIspeed) -#' i <- BCIspeed$species=="ocelot" -#' sp <- log(BCIspeed$speed[i]) -#' tm <- BCIspeed$time[i]*2*pi -#' circseq <- seq(0,2*pi,pi/256) -#' trend <- lincircKern(circseq, tm, sp) -#' plot(circseq, trend, type="l") -#' @export -lincircKern <- function(x,circdat,lindat) -{ if(length(lindat)!=length(circdat)) - stop("lindat and circdat lengths are unequal") - if(min(circdat)<0 | max(circdat)>2*pi) - stop("circdat values not between 0 and 2*pi, expecting radian data") - if(min(x)<0 | max(x)>2*pi) - stop("x values not between 0 and 2*pi, expecting radian values") - - hs <- 1.06 * min(sd(lindat), (quantile(lindat,0.75)-quantile(lindat,0.25))/1.34) * - length(lindat)^-0.2 - bw <- 1/hs^2 - dx <- expand.grid(circdat,x) - dif <- abs(dx[,2]-dx[,1]) - i <- dif>pi - dif[i] <- 2*pi-dif[i] - prob <- matrix(circular::dvonmises(circular::circular(dif),circular::circular(0),bw), nrow=length(circdat)) - apply(prob,2,function(z) mean(z*lindat)/mean(z)) -} - -#' Linear-circular regression -#' -#' Fits a Von Mises kernel distribution describing a linear variable as a function -#' of a circular predictor, and boostraps the null distribution in order to evaluate -#' significance of radial variation in the linear variable. -#' -#' Deviation of \code{lindat} from the null expecation is assessed either visually -#' by the degree to which the fitted distribution departs from the null confidence -#' interval (use generic plot function), or quantitatively by column \code{p} of -#' slot \code{fit} in the resulting \code{lincircmod-class} object. -#' -#' @param circdat Numeric vector of radian data matched with \code{lindat}. -#' @param lindat Numeric vector of linear data matched with \code{circdat}. -#' @param pCI Single numeric value between 0 and 1 defining proportional confidence interval to return. -#' @param reps Integer number of bootstrap repetitions to perform. -#' @param res Resolution of fitted distribution and null confidence interval - specifically a single integer number of points on the circular scale at which to record distributions. -#' @return An object of type \code{\link{lincircmod-class}} -#' @references Xu, H., Nichols, K. & Schoenberg, F.P. (2011) Directional kernel regression for wind and fire data. Forest Science, 57, 343-352. -#' @examples -#' #Example with reps limited to increase speed -#' data(BCIspeed) -#' i <- BCIspeed$species=="ocelot" -#' sp <- log(BCIspeed$speed[i]) -#' tm <- BCIspeed$time[i]*2*pi -#' mod <- fitlincirc(tm, sp, reps=50) -#' plot(mod, CircScale=24, xaxp=c(0,24,4), -#' xlab="Time", ylab="log(speed m/s)") -#' legend(8,-3, c("Fitted speed", "Null CI"), col=1:2, lty=1:2) -#' @export -fitlincirc <- function(circdat,lindat,pCI=0.95,reps=1000,res=512) -{ if(length(lindat)!=length(circdat)) - stop("lindat and circdat lengths are unequal") - if(min(circdat)<0 | max(circdat)>2*pi) - stop("circdat values not between 0 and 2*pi, expecting radian data") - if(max(circdat)<1) - warning("max(circdat) < 1, expecting radian data") - - n <- length(circdat) - x <- seq(0,2*pi,2*pi/res) - - bs <- pbapply::pbsapply(1:reps, function(i) - { j <- sample(1:n,n,TRUE) - lincircKern(x,circdat[j],lindat) - }) - nulllcl <- apply(bs,1,quantile,(1-pCI)/2) - nullucl <- apply(bs,1,quantile,(1+pCI)/2) - fit <- lincircKern(x,circdat,lindat) - p <- sapply(1:(res+1), function(i) - { f <- ecdf(bs[i,]) - f(fit[i]) - }) - p[p>0.5] <- 1-p[p>0.5] - p <- 2*p - - new("lincircmod", data=data.frame(circdat=circdat, lindat=lindat), - fit=data.frame(x=x, fit=fit, p=p, nullLCL=nulllcl, nullUCL=nullucl)) -} - -#' Plot linear-circular relationship -#' -#' Plot linear against circular data along with the fitted and null confidence limit distributions from a fitted \code{lincircmod} object. -#' -#' @param x Object of class \code{lincircmod}. -#' @param CircScale Single numeric value defining the plotting maximum of the circular scale. -#' @param tlim Numeric vector with two elements >=0 and <=1 defining the lower and upper limits at which to plot distributions; default plots the full range. -#' @param fcol,flty,ncol,nlty Define line colour (\code{col}) and type (\code{lty}) for fitted (\code{f}) and null (\code{n}) distributions; input types as for \code{col} and \code{lty}, see \code{\link{par}}. -#' @param ... Additional arguments passed to the inital plot construction, affecting axes and data plot symbols. -#' @docType methods -#' @rdname plot-lincircmod-methods -#' @export -setMethod("plot", "lincircmod", - function(x, CircScale=2*pi, tlim=c(0,1), fcol="black", flty=1, ncol="red", nlty=2, ...) - { if(min(tlim)<0 | max(tlim)>1 | length(tlim)!=2) stop("tlim should contain two values >=0 and <=1") - - fit <- x - fit <- x@fit - dat <- x@data - xx <- fit$x*CircScale/(2*pi) - LinearData <- dat$lindat - CircularData <- dat$circdat*CircScale/(2*pi) - range <- tlim*CircScale - plot(CircularData, LinearData, ...) - if(range[1]=range[1] & xx<=range[2] - lines(xx[i], fit$fit[i], col=fcol, lty=flty) - lines(xx[i], fit$nullLCL[i], col=ncol, lty=nlty) - lines(xx[i], fit$nullUCL[i], col=ncol, lty=nlty) - } else - { i <- xx>=range[1] - lines(xx[i], fit$fit[i], col=fcol, lty=flty) - lines(xx[i], fit$nullLCL[i], col=ncol, lty=nlty) - lines(xx[i], fit$nullUCL[i], col=ncol, lty=nlty) - i <- xx<=range[2] - lines(xx[i], fit$fit[i], col=fcol, lty=flty) - lines(xx[i], fit$nullLCL[i], col=ncol, lty=nlty) - lines(xx[i], fit$nullUCL[i], col=ncol, lty=nlty) - } - } -) diff --git a/man/BCIspeed.Rd b/man/BCIspeed.Rd index dee3ed2..1339eb6 100644 --- a/man/BCIspeed.Rd +++ b/man/BCIspeed.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/activity_code_1_1.r +% Please edit documentation in R/activity_code.r \docType{data} \name{BCIspeed} \alias{BCIspeed} @@ -13,4 +13,3 @@ Barro Colorado Island 2008 data: speeds of animal passages past camera traps (\code{speed}), together with species (\code{species}) and time of day (\code{time}) for each record. } - diff --git a/man/BCItime.Rd b/man/BCItime.Rd index 5faa18f..cc8d717 100644 --- a/man/BCItime.Rd +++ b/man/BCItime.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/activity_code_1_1.r +% Please edit documentation in R/activity_code.r \docType{data} \name{BCItime} \alias{BCItime} @@ -12,4 +12,3 @@ http://dx.doi.org/10.6084/m9.figshare.1160536 Barro Colorado Island 2008 data: times of day at which animal records occured (\code{time}), together with species (\code{species}). } - diff --git a/man/activity.Rd b/man/activity.Rd index d836f0c..55d22a2 100644 --- a/man/activity.Rd +++ b/man/activity.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/activity_code_1_1.r +% Please edit documentation in R/activity_code.r \docType{package} \name{activity} \alias{activity} @@ -27,7 +27,3 @@ activity estimates. Note that all time or other circular data should be in radia \references{ Rowcliffe, M., Kays, R., Kranstauber, B., Carbone, C., Jansen, P.A. (2014) Quantifying animal activity level using camera trap data. Methods in Ecology and Evolution. } -\seealso{ -\code{\link{overlap}} -} - diff --git a/man/actmod-class.Rd b/man/actmod-class.Rd index 8ca75d3..343aa52 100644 --- a/man/actmod-class.Rd +++ b/man/actmod-class.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/activity_code_1_1.r +% Please edit documentation in R/activity_code.r \docType{class} \name{actmod-class} \alias{actmod-class} @@ -19,11 +19,11 @@ An S4 class describing activity models fitted to time of observation data. \item{\code{adj}}{Object of class \code{"numeric"}, kernel bandwidth adjustment multiplier.} \item{\code{pdf}}{Object of class \code{"matrix"} describing fitted probability density function: - Column 1: Sequence of radian times at which PDF evaluated (specifically seq(0,2*pi,pi/256)). + Column 1: A regular sequence of radian times at which PDF evaluated; range is [0, 2*pi] if unbounded, and sequence steps are range difference divided by 512. Column 2: Corresponding circular kernel PDF values. -Additionally if errors boostrapped: +Additionally if errors bootstrapped: Column 3: PDF standard error. - Column 4: PDF lower confidence interval. Column 5: PDF upper confidence interval.} + Column 4: PDF lower 95\% confidence limit. Column 5: PDF upper 95\% confidence limit.} \item{\code{act}}{Object of class \code{"numeric"} giving activity level estimate and, if errors boostrapped, standard error and 95 percent confidence limits.} }} diff --git a/man/bwcalc.Rd b/man/bwcalc.Rd index 82feac2..6d31d4c 100644 --- a/man/bwcalc.Rd +++ b/man/bwcalc.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/activity_code_1_1.r +% Please edit documentation in R/activity_code.r \name{bwcalc} \alias{bwcalc} \title{Calculate circular kernel bandwidth.} @@ -23,4 +23,3 @@ Mainly for internal use. \references{ Ridout, M.S. & Linkie, M. (2009) Estimating overlap of daily activity patterns from camera trap data. Journal of Agricultural Biological and Environmental Statistics, 14, 322-337. } - diff --git a/man/compareAct.Rd b/man/compareAct.Rd index 5ed61ff..676aca5 100644 --- a/man/compareAct.Rd +++ b/man/compareAct.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/activity_code_1_1.r +% Please edit documentation in R/activity_code.r \name{compareAct} \alias{compareAct} \title{Compare activity level estimates} @@ -24,10 +24,11 @@ on chi-sq distribution with 1 degree of freedom. #Test whether paca have a sigificantly different activity level from rat. #Bootstrap reps limited to speed up example. data(BCItime) -tPaca <- 2*pi*BCItime$time[BCItime$species=="paca"] +tPaca <- 2*pi*BCItime$time[BCItime$species=="ocelot"] tRat <- 2*pi*BCItime$time[BCItime$species=="rat"] -(fPaca <- fitact(tPaca, reps=10)) -(fRat <- fitact(tRat, reps=10)) +fPaca <- fitact(tPaca, sample="data", reps=10) +fRat <- fitact(tRat, sample="data", reps=10) +fPaca@act +fRat@act compareAct(list(fPaca,fRat)) } - diff --git a/man/compareCkern.Rd b/man/compareCkern.Rd index f2751c6..6a99262 100644 --- a/man/compareCkern.Rd +++ b/man/compareCkern.Rd @@ -1,38 +1,36 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/activity_code_1_1.r +% Please edit documentation in R/activity_code.r \name{compareCkern} \alias{compareCkern} \title{Compare circular distributions.} \usage{ -compareCkern(y1, y2, reps = 1000, index = c("Dhat4", "Dhat5", "Dhat1")) +compareCkern(fit1, fit2, reps = 999) } \arguments{ -\item{y1, y2}{Numeric vectors of radian data.} +\item{fit1, fit2}{Fitted activity models of class actmod created using function fitact.} \item{reps}{Number of bootstrap iterations.} - -\item{index}{Which of the three indices returned by overlap::overlapEst to use.} } \value{ -A named 2-element vector: Overlap = observed overlap index; p = probability observed index arose by chance. +A named 4-element vector: obs = observed overlap index; null = mean null overlap index; seNull = standard error of the null distribution; pNull = probability observed index arose by chance. } \description{ Randomisation test for the probability that two sets of circular observations come from the same distribution. } \details{ -Calculates overlap index (see references) for the observed data samples, then generates a null distribution of overlap indices using data sampled randomly with replacement from the combined data. This randomised distribution is then used to estimate the probability that the observed overlap arose by chance. +Calculates overlap index Dhat5 (see references) for the two fitted distributions, then generates a null distribution of overlap indices using data sampled randomly with replacement from the combined data. +This randomised distribution is then used to define an empirical probability distribution against which the probability that the observed overlap arose by chance is judged. +When one or both fitted models use weighted distributions, sampling probabilities are taken from the weights. If both models are weighted, the weights must therefore be on the same scale. } \examples{ -#Example with bootstrap reps limited to speed up +#Example with bootstrap reps limited to reduce run time data(BCItime) tPaca <- 2*pi*BCItime$time[BCItime$species=="paca"] tRat <- 2*pi*BCItime$time[BCItime$species=="rat"] -compareCkern(tPaca,tRat,reps=10) +fPaca <- fitact(tPaca) +fRat <- fitact(tRat) +compareCkern(fPaca,fRat,reps=10) } \references{ Ridout, M.S. & Linkie, M. (2009) Estimating overlap of daily activity patterns from camera trap data. Journal of Agricultural Biological and Environmental Statistics, 14, 322-337. } -\seealso{ -\code{\link{overlapEst}} -} - diff --git a/man/compareTimes.Rd b/man/compareTimes.Rd index 9939712..0bb0bb3 100644 --- a/man/compareTimes.Rd +++ b/man/compareTimes.Rd @@ -1,8 +1,8 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/activity_code_1_1.r +% Please edit documentation in R/activity_code.r \name{compareTimes} \alias{compareTimes} -\title{Compare activity across between times of day} +\title{Compare activity between times of day} \usage{ compareTimes(fit, times) } @@ -23,8 +23,7 @@ Bootrapping the activity model yields standard error estimates for the PDF. This \examples{ data(BCItime) tPaca <- 2*pi*BCItime$time[BCItime$species=="paca"] -fPaca <- fitact(tPaca, reps=10) -plot(fPaca, hrs=FALSE, frq=FALSE) +fPaca <- fitact(tPaca, sample="data", reps=10) +plot(fPaca) compareTimes(fPaca, c(5.5,6,0.5,1)) } - diff --git a/man/density2.Rd b/man/density2.Rd new file mode 100644 index 0000000..1fd5a0d --- /dev/null +++ b/man/density2.Rd @@ -0,0 +1,37 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/activity_code.r +\name{density2} +\alias{density2} +\title{Modified kernel density function} +\usage{ +density2(x, reps = 999, ...) +} +\arguments{ +\item{x}{numeric data vector} + +\item{reps}{bootstrap iterations for SE/interval calculation; set to NULL to suppress} + +\item{...}{Additional arguments passed to \code{stas::density}} +} +\value{ +A list with the same components as \code{stats::density} output plus: + \code{se}: standard error of the density + \code{lcl}, \code{ucl}: lower and upper 95\% confidence intervals of the density +} +\description{ +Modifies \code{stats::density} by: + Adding SE and 95\% confidence intervals for the density to the output; and + Truncating calculation (not just reporting) of density values on from and/or to. +} +\details{ +Truncation copes with cases where where no data are available outside truncation points. +Truncation is achieved by fitting the density to the data augmented by reflecting it +across each bound using the optimal bandwidth for the unaugmented data, and returning +the resulting densities for the region between the bounds. +} +\examples{ +data(BCItime) +tm <- subset(BCItime, species=="ocelot")$time +dens <- density2(tm, from=0.25, to=0.75) +plot(dens$x, dens$y, type="l") +} diff --git a/man/dvmkern.Rd b/man/dvmkern.Rd index 0dd32d6..4e3825f 100644 --- a/man/dvmkern.Rd +++ b/man/dvmkern.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/activity_code_1_1.r +% Please edit documentation in R/activity_code.r \name{dvmkern} \alias{dvmkern} \title{Circular kernel probability density function.} @@ -36,4 +36,3 @@ plot(xx, pdf, type="l") \seealso{ \code{\link{density.circular}, \link{bwcalc}} } - diff --git a/man/fitact.Rd b/man/fitact.Rd index 36aac54..f7db98e 100644 --- a/man/fitact.Rd +++ b/man/fitact.Rd @@ -1,11 +1,11 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/activity_code_1_1.r +% Please edit documentation in R/activity_code.r \name{fitact} \alias{fitact} \title{Fit activity model to time-of-day data} \usage{ -fitact(dat, wt = NULL, reps = 1000, bw = NULL, adj = 1, - sample = c("data", "model", "none"), show = TRUE) +fitact(dat, wt = NULL, reps = 999, bw = NULL, adj = 1, + sample = c("none", "data", "model"), bounds = NULL, show = TRUE) } \arguments{ \item{dat}{A numeric vector of radian time-of-day data.} @@ -20,45 +20,55 @@ fitact(dat, wt = NULL, reps = 1000, bw = NULL, adj = 1, \item{sample}{Character string defining sampling method for bootstrapping errors (see details).} +\item{bounds}{A two-element vector defining radian bounds at which to truncate.} + \item{show}{Logical whether or not to show a progress bar while bootstrapping.} } \value{ An object of type \code{actmod} } \description{ -Fits a circular kernel density to radian time-of-day data and estimates activity -level from this distribution. Optionally bootstraps the distribution, in which -case SEs and confidence limits are also stored for activity level and PDF. +Fits kernel density to radian time-of-day data and estimates activity level from this distribution. +Optionally: 1. bootstraps the distribution, in which case SEs and confidence limits are also +stored for activity level and PDF; 2. weights the distribution; 3. truncates the distribution at given times. } \details{ +When no \code{bounds} are given (default), a circular kernel distribution is fitted using \code{dvmkern}. +Otherwise, a normal kernel distribution is used, truncated at the values of \code{bounds}, using \code{density2}. + The bandwidth adjustment multiplier \code{adj} is provided to allow exploration of the effect of adjusting the internally calculated bandwidth on -accuracy of activity level estimates. The alternative bootstrapping methods -defined by \code{sample} are: +accuracy of activity level estimates. + +The alternative bootstrapping methods defined by \code{sample} are: + \code{none}: no bootstrapping. \code{data}: sample from the data; \code{model}: sample from the fitted probability density distribution; - \code{none}: no bootstrapping. Confidence interval coverage seems to be better at large sample size -(greater than 100-200) using \code{"model"}, but better at small sample size -when using \code{"data"}. The reason for this needs further investigation. +(greater than around 100) using \code{"model"}, but better at small sample size +when using \code{"data"}. } \examples{ #Fit without confidence limits data(BCItime) -tdat <- 2*pi*BCItime$time[BCItime$species=="ocelot"] -mod1 <- fitact(tdat, sample="none") +tm <- 2*pi*subset(BCItime, species=="brocket")$time +mod1 <- fitact(tm) plot(mod1) #Fit with confidence limits (limited reps to speed up) -mod2 <- fitact(tdat, reps=10) +mod2 <- fitact(tm, sample="data", reps=10) plot(mod2) -#Fit weighted function to correct for detection radius 1.21 times higher -#by day than by night, assuming day between pi/2 (6am) and pi*2/3 (6pm) -weight <- 1/ifelse(tdat>pi/2 & tdatpi/2 & tm