R/estc2.R

In jedalong/wildlifeTG: Time Geograhic Analysis of Wildlife Telemetry Data

# ---- roxygen documentation ----
#' @title Numerically estimate c2 parameter
#'
#' @description
#'   Numerically estimate the \code{c2} parameter for \code{fbtgUD} using a data-driven approach.
#' @details
#' THe estimation of the \code{c2} parameter for \code{fbtgUD} model takes an identical approach to that used with Brownian bridges, as originally proposed by Horne et al. (2007). A leave-one-out estimation technique is used, which essentially removes fixes and then estimates the fbtgUD surface between the two adjacent fixes  -- termed a segment -- and computes the probability associated with the missing fix. The \code{c2} value is returned that numerically maximizes the log-likelihood (for a given \code{timefun} -- see details in \code{fbtgUD}) for \code{length.out} evenly spaced values of \code{c2} between (user-defined) \code{min} and \code{max}. With \code{plot = TRUE} the user can verify that the chosen range of potential \code{c2} values is appropriate, and if not, retry using a differnt range. The process is computationally demanding and is highly dependent on the number of 'segments' used. The parameter \code{rand} can be used to adjust the number of segments used to minimize computational time. If \code{rand = NA} it removes every second fix, and estimates \code{c2} based on these n/2 segments (the default). Otherwise \code{rand} can be passed in as an integer, and \code{rand} randomly selected segments will be chosen to estimate \code{c2}. This is beneficial for trajectories with many fixes, where it might be useful to choose \code{rand <<< n/2}. Parrallelization is possible to further decrease computational time. This can be implemented by choosing an appropriate integer value for the \code{parallel} parameter (implemented using the \code{foreach} package). 
#'   
#' @param traj animal movement trajectory in the form of an \code{ltraj} object, see package \code{adehabitatLT}
#' @param tl a \code{TransitionLayer} object
#' @param timefun method for converting time into probability (see \code{fbtgUD}); one of:\cr 
#' -- \code{'inverse'} \eqn{= \frac{1}{c_2 * t}},\cr
#' -- \code{'inverse2'} \eqn{= \frac{1}{(c_2 * t)^2}},\cr
#' -- \code{'exp'} \eqn{= \exp{-c_2 * t}},\cr
#' -- \code{'norm'} \eqn{= \exp{-c_2 * t^2}},\cr
#' -- \code{'rootexp'} \eqn{= \exp{-c_2 * \sqrt{t}}},\cr
#' -- \code{'pareto'} \eqn{= \exp{-c_2 * \log{t}}},\cr
#' -- \code{'lognorm'}\eqn{= \exp{-c_2 * \log{t^2}}}.\cr
#' @param sigma location uncertainty parameter (see fbtgUD)
#' @param lower lower bound for \code{c2} testing range (default = 0) 
#' @param upper upper bound for \code{c2} testing range (default = 1)
#' @param rand if \code{NA} (the default) every second segment is evauluated (n/2), otherwise an integer indicating how many random segments to test.
#' @param niter used to define maximum number of iterations of golden-search routine
#' @param tolerance used to define precision of golden search routine (i.e., routine stops when the absolute difference between two consecutive test points is below this value). 
#' @param plot logical, whether or not to plot the log-likelihood curve.
#' 
#' @return
#'   This function returns a numerical value estimate for \code{c2} associated with the maximum of the log-likelihood.
#'
# @references
# @keywords 
#' @seealso fbtgUD
# @examples
#' @importFrom graphics abline points
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @export
# ---- End of roxygen documentation ----

estc2 <- function(traj,tl,timefun='exp',sigma=0,lower=0,upper=1,rand=NA,niter=10,tolerance=0.01,plot=TRUE){
  
  #use leave-one-out bootstrap to estimate c2 in a similar fashion to proposed by Horne et al. 2007 as is commonly used with Brownian bridge, can speed up by chosing smaller number of random segments to test.
  x <- ld(traj)
  n <- dim(x)[1]
  if (is.na(rand)){
    ii <- seq(1,(n-2),by=2)
  } else {
    ii <- sample(1:(n-2),rand)
  }

  tm <- transitionMatrix(tl)
  gr <- graph.adjacency(tm, mode="directed", weighted=TRUE)
  E(gr)$weight <- 1/E(gr)$weight
  
  A <- SpatialPoints(x[ii,c('x','y')])
  B <- SpatialPoints(x[ii+1,c('x','y')])
  C <- SpatialPoints(x[ii+2,c('x','y')])
  Ai <- cellFromXY(raster(tl),A)
  Bi <- cellFromXY(raster(tl),B)
  Ci <- cellFromXY(raster(tl),C)
  
  Tshort <- diag(costDistance(tl,A,C))
  Tai <- distances(gr,v=Ai,to=V(gr),mode='out')
  Tib <- distances(gr,v=Ci,to=V(gr),mode='in')
  t1 <- x$dt[ii]
  t2 <- x$dt[ii+1]
  tt <- t1+t2
  pz <- 0*ii
  
  ### Golden Search Routine
  
  #internal likelihood funciton
  c2func <- function(c2,ii,tl,sigma,timefun,x,Ai,Bi,Ci,Tshort,Tai,Tib,tt){
    for (i in 1:length(ii)){
      #Compute the likelihood for each segment
      pz[i] <- internalTS(t1[i],tt[i],Tai[i,],Tib[i,],tl,Tshort[i],sigma=sigma,timefun,c2,clipPPS=FALSE)[Bi[i]]
    }
    #Calculate the negative of the log-likelihood - we are using a minimizing golden search function
    LLpz <- -log(pz) 
    #REMOVE -INFs
    LLpz[is.infinite(LLpz)] <- NA
    LL <- sum(LLpz,na.rm=T)
    return(LL)
  }
  
  #### GOLDEN SEARCH ROUTINE ###
  #Progress Bar
  cat('Iterations: \n')
  
  golden.ratio = 2/(sqrt(5) + 1)
  
  ### Evaluate the function at the extremes
  fmin = c2func(lower,ii,tl,sigma,timefun,x,Ai,Bi,Ci,Tshort,Tai,Tib,tt)
  cat('1 \n')
  fmax = c2func(upper,ii,tl,sigma,timefun,x,Ai,Bi,Ci,Tshort,Tai,Tib,tt)
  cat('2 \n')
  
  ### Use the golden ratio to set the initial test points
  x1 = upper - golden.ratio*(upper - lower)
  x2 = lower + golden.ratio*(upper - lower)
  
  ### Evaluate the function at the first test points
  f1 <- c2func(x1,ii,tl,sigma,timefun,x,Ai,Bi,Ci,Tshort,Tai,Tib,tt)
  cat('3 \n')
  f2 <- c2func(x2,ii,tl,sigma,timefun,x,Ai,Bi,Ci,Tshort,Tai,Tib,tt)
  cat('4 \n')
  ### Output values storage
  c2.val <- c(lower,upper,x1,x2)
  LL.val <- c(fmin,fmax,f1,f2)
  
  #Search
  iteration = 0
  #Progress Bar
  
  while (iteration < (niter-4) & abs(upper - lower) > tolerance){
    iteration = iteration + 1
    if (f2 > f1){
      upper = x2
      x2 = x1
      f2 = f1
      x1 = upper - golden.ratio*(upper - lower)
      f1 = c2func(x1,ii,tl,sigma,timefun,x,Ai,Bi,Ci,Tshort,Tai,Tib,tt)
      c2.val <- c(c2.val,x1)
      LL.val <- c(LL.val,f1)
    } else {
      lower = x1
      x1 = x2
      f1 = f2
      x2 = lower + golden.ratio*(upper - lower)
      f2 = c2func(x2,ii,tl,sigma,timefun,x,Ai,Bi,Ci,Tshort,Tai,Tib,tt)
      c2.val <- c(c2.val,x2)
      LL.val <- c(LL.val,f2)
    }
    #update progress bar
    est.min = (lower + upper)/2
    cat(paste(iteration+4,' c2.est = ',est.min,'\n'))
  }
  #est.min = (lower + upper)/2

  if (plot){
    ord <- order(c2.val)
    plot(c2.val[ord],-LL.val[ord],xlab='theta',ylab='log-likelihood',pch=20)
    points(c2.val[ord],-LL.val[ord],type='l')
    abline(v=est.min,col='red')
  }

  return(est.min)
}