add dcrwrice and dcrwvm distributions for correlated steps and turn a…

…ngles; fix bug in calcuating initial velocity terms in crw models when there are multiple individuals

DESCRIPTION

@@ -63,7 +63,8 @@ Suggests:
     gdistance,
     methods,
     mgcv,
-    optimx
+    optimx,
+    pracma
 RoxygenNote: 7.3.2
 Encoding: UTF-8
 NeedsCompilation: yes

NAMESPACE

@@ -171,8 +171,10 @@ importFrom(sp,spDistsN1)
 importFrom(sp,spTransform)
 importFrom(stats,acf)
 importFrom(stats,aggregate)
+importFrom(stats,approxfun)
 importFrom(stats,as.formula)
 importFrom(stats,dbinom)
+importFrom(stats,density)
 importFrom(stats,dnbinom)
 importFrom(stats,dnorm)
 importFrom(stats,dunif)
@@ -188,6 +190,7 @@ importFrom(stats,optim)
 importFrom(stats,optimHess)
 importFrom(stats,optimise)
 importFrom(stats,pbinom)
+importFrom(stats,pchisq)
 importFrom(stats,plogis)
 importFrom(stats,pnbinom)
 importFrom(stats,predict)

R/CIreal.R

@@ -142,7 +142,7 @@ CIreal.default <- function(m,alpha=0.95,covs=NULL,parms=NULL)
     } else {
       par <- as.vector(t(m$mle[[i]]))
     }
-    if(!(inputs$dist[[i]] %in% angledists) | (inputs$dist[[i]] %in% angledists & m$conditions$estAngleMean[[i]] & !m$conditions$Bndind[[i]])) {
+    if(!(inputs$dist[[i]] %in% angledists) | (inputs$dist[[i]] %in% angledists & m$conditions$estAngleMean[[i]] & !m$conditions$Bndind[[i]]) | (inputs$dist[[i]]=="crwvm")) {
       Par[[i]] <- get_CI(m$mod$estimate,par,m,parindex[[i]]+1:parCount[[i]],fullDM[[i]],DMind[[i]],p$bounds[[i]],Sigma,m$conditions$circularAngleMean[[i]],inputs$consensus[[i]],nbStates,alpha,tmpParNames,m$stateNames,nc[[i]],meanind[[i]],m$conditions$workBounds[[i]],inputs$dist[[i]])
     } else {
       if(!m$conditions$estAngleMean[[i]])

R/RcppExports.R

@@ -235,6 +235,32 @@ dt_rcpp <- function(x, df, ncp) {
     .Call(`_momentuHMM_dt_rcpp`, x, df, ncp)
 }
 
+#' Correlated random walk Rice distribution
+#'
+#' Probability density function of Rice distribution under the CTCRW model of Johnson et al. (2008) (written in C++)
+#'
+#' @param x numeric data vector of length \code{n + n} where the first \code{n} entries correspond to the step lengths and the next \code{n} entries to the corresponding previous step lengths (the first of which is NA and ignored).
+#' @param beta correlation parameter 
+#' @param sigma speed parameter
+#'
+#' @return Vector of densities
+dcrwrice_rcpp <- function(x, beta, sigma) {
+    .Call(`_momentuHMM_dcrwrice_rcpp`, x, beta, sigma)
+}
+
+#' Correlated random walk von Mises density function
+#'
+#' Probability density function of the Von Mises distribution under the CTCRW model of Johnson et al. (2008) (written in C++)
+#'
+#' @param x numeric data vector of length \code{n + n + n} where the first \code{n} entries correspond to angles (von Mises distribution), the next \code{n} entries to the corresponding step lengths, and the last \code{n} entries to the corresponding previous step lengths.
+#' @param beta correlation parameter 
+#' @param sigma speed parameter
+#'
+#' @return Vector of densities
+dcrwvm_rcpp <- function(x, beta, sigma) {
+    .Call(`_momentuHMM_dcrwvm_rcpp`, x, beta, sigma)
+}
+
 #' Matrix Exponential
 #'
 #' This function computes the exponential of a square matrix \code{A}, defined as the sum from \code{r=0} to infinity of \code{A^r/r!}, using the default method of \code{\link[expm]{expm}}.

R/addSmoothGradient.R

@@ -26,7 +26,7 @@ addSmoothGradient <- function(data,#CT=FALSE,Time.name="time",Time.unit='auto',
   if(!is.momentuHMMData(data)) stop("data must be a momentuHMMData object (as returned by prepData or simData) ")
   if(inherits(data,"ctds")) stop("data cannot be a 'ctds' object (as returned by prepCTDS")
 
-  if(isTRUE(inherits(data,"CT"))){
+  if(isTRUE(attr(data,"CT"))){
     CT <- TRUE
     Time.name <- attr(data,"Time.name")
     Time.unit <- attr(data,"Time.unit")
@@ -70,7 +70,6 @@ addSmoothGradient <- function(data,#CT=FALSE,Time.name="time",Time.unit='auto',
   rx <- diff(data[iInd,coordNames[1]])
   ry <- diff(data[iInd,coordNames[2]])
   if(isTRUE(CT)){
-    if(!isTRUE(inherits(data,"CT"))) warning("'data' appears to be in discrete time -- should 'CT' be FALSE?")
     if(is.null(Time.name)) stop("'Time.name' cannot be NULL when CT=TRUE")
     else {
       if(!Time.name %in% names(data)) stop("Time.name not found in 'data'") 
@@ -84,7 +83,6 @@ addSmoothGradient <- function(data,#CT=FALSE,Time.name="time",Time.unit='auto',
       }
     }
   } else {
-    if(isTRUE(inherits(data,"CT"))) warning("'data' appears to be in continuous time -- should 'CT' be TRUE?")
     dt <- rep(1,length(rx))
   }
   locs <- utils::head(data[iInd,coordNames],-1)

R/allProbs.R

@@ -108,11 +108,21 @@ allProbs <- function(m)
         genData <- c(data[[paste0(i,".x")]],data[[paste0(i,".y")]])
       else if(dist[[i]]=="mvnorm3" || dist[[i]]=="rw_mvnorm3")
         genData <- c(data[[paste0(i,".x")]],data[[paste0(i,".y")]],data[[paste0(i,".z")]])
+      genInd <- which(!is.na(genData[1:nbObs]))
+    } else if(dist[[i]]=="crwrice" | dist[[i]]=="crwvm"){
+      genData <- cbind(data$step,data$angle)
+      step_tm1 <- c(NA,data$step[-length(data$step)])
+      if(dist[[i]]=="crwrice"){
+        genData[c(1,1+cumsum(table(m$data$ID)[-(length(unique(m$data$ID)))])),1] <- NA
+        genInd <- which(!is.na(genData[,1]))
+      } else genInd <- which(!is.na(data[[i]]))
+      step_tm1[-genInd] <- NA
+
     } else {
       genData <- data[[i]]
+      genInd <- which(!is.na(genData[1:nbObs]))
     }
 
-    genInd <- which(!is.na(genData[1:nbObs]))
     sp <- par[[i]]
     genProb[[i]] <- matrix(1,nbObs,nbStates)
 
@@ -121,7 +131,11 @@ allProbs <- function(m)
       genFun <- Fun[[i]]
 
       # Constitute the lists of state-dependent parameters for the step and angle
-      genArgs <- list(genData[which(!is.na(genData))])
+      if(dist[[i]]=="crwrice"){
+          genArgs <- list(genData[genInd,1,drop=FALSE])
+      } else if(dist[[i]]=="crwvm"){
+          genArgs <- list(genData[genInd,2,drop=FALSE])
+      } else genArgs <- list(genData[which(!is.na(genData))])
 
       zeromass <- 0
       onemass <- 0
@@ -167,6 +181,12 @@ allProbs <- function(m)
         genArgs[[2]] <- shape
         genArgs[[3]] <- 1/scale # dgamma expects rate=1/scale
       }
+      if(dist[[i]]=="crwrice"){
+        genArgs[[j+2]] <- step_tm1[genInd]
+      } else if(dist[[i]]=="crwvm") {
+        genArgs[[j+2]] <- genData[genInd,1] # add steps argument
+        genArgs[[j+3]] <- step_tm1[genInd] # add steps_tm1 argument
+      }
       if(zeroInflation[[i]] | oneInflation[[i]]) {
         if(zeroInflation[[i]] & !oneInflation[[i]]){
           genProb[[i]][genInd,state] <- ifelse(genData[genInd]==0,
@@ -182,7 +202,9 @@ allProbs <- function(m)
           genProb[[i]][genInd,state][genData[genInd]>0 & genData[genInd]<1] <- (1.-zeromass[genData[genInd]>0 & genData[genInd]<1]-onemass[genData[genInd]>0 & genData[genInd]<1]) * do.call(genFun,genArgs)[genData[genInd]>0 & genData[genInd]<1] # if gen !=0 and gen!=1
         }
       }
-      else genProb[[i]][genInd,state] <- do.call(genFun,genArgs)
+      else {
+        genProb[[i]][genInd,state] <- do.call(genFun,genArgs)
+      }
     }
     zeroInd <- which(apply(genProb[[i]],1,function(x) all(x < .Machine$double.xmin)))
     if(length(zeroInd)){

R/checkInputs.R

@@ -55,6 +55,8 @@ checkInputs<-function(nbStates,dist,Par,estAngleMean,circularAngleMean,zeroInfla
 
   ndist <- dist
   for(i in distnames){
+    if(dist[[i]]=="crwrice" & i!="step") stop("'crwrice' distribution observations must be named 'step'")
+    if(dist[[i]]=="crwvm" & i!="angle") stop("'crwvm' distribution observations must be named 'angle'")
     if(dist[[i]]=="vmConsensus"){
       consensus[[i]] <- TRUE
       estAngleMean[[i]] <- TRUE
@@ -69,6 +71,27 @@ checkInputs<-function(nbStates,dist,Par,estAngleMean,circularAngleMean,zeroInfla
   circularAngleMean<-circularAngleMean[distnames]
   consensus<-consensus[distnames]
 
+  crwST <- FALSE
+  if("angle" %in% distnames){ 
+    if(dist[["angle"]] %in% angledists & ("step" %in% distnames)){
+      if((dist[["angle"]]=="crwvm" & dist[["step"]]!="crwrice") | (dist[["angle"]]!="crwvm" & dist[["step"]]=="crwrice")) stop("'crwvm' angle distributions must be paired with 'crwrice' step distribution")
+      if(dist[["angle"]]=="crwvm" & dist[["step"]]=="crwrice") {
+        crwST <- TRUE
+        if(!is.null(DM)){
+          if(!is.null(DM$angle)){
+            if(!isTRUE(all.equal(DM$step,DM$angle,check.attributes=FALSE))) stop("DM$step and DM$angle should be the same for 'crwrice' and 'crwvm' distributions")
+          }
+        }
+        if(isTRUE(zeroInflation$step)){
+          if(!isTRUE(all.equal(Par$step[1:length(Par$angle)],Par$angle,check.attributes=FALSE))) stop("Par0$step and Par0$angle 'beta' and 'sigma' parameters should be the same for 'crwrice' and 'crwvm' distributions")
+        } else if(!isTRUE(all.equal(Par$step,Par$angle,check.attributes=FALSE))) stop("Par0$step and Par0$angle should be the same for 'crwrice' and 'crwvm' distributions")
+        estAngleMean$angle <- TRUE
+        circularAngleMean$angle <- FALSE
+        consensus$angle <- FALSE
+      }
+    } else if(dist[["angle"]]=="crwvm") stop("'crwvm' angle distributions must be paired with 'crwrice' step distribution")
+  }
+
   if(!is.null(stateNames)){
     if(!is.character(stateNames)) stop("stateNames must be a character vector")
     if(length(stateNames)!=nbStates) stop("stateNames must have length ",nbStates)
@@ -118,7 +141,7 @@ checkInputs<-function(nbStates,dist,Par,estAngleMean,circularAngleMean,zeroInfla
 
   checkNames(distnames,estAngleMean=estAngleMean,circularAngleMean=circularAngleMean,zeroInflation=zeroInflation,oneInflation=oneInflation,DM=DM,userBounds=userBounds)
 
-  return(list(p=p,dist=ndist,estAngleMean=estAngleMean,circularAngleMean=circularAngleMean,consensus=consensus,DM=DM))
+  return(list(p=p,dist=ndist,estAngleMean=estAngleMean,circularAngleMean=circularAngleMean,consensus=consensus,DM=DM,crwST=crwST))
 }
 
 checkNames <- function(distnames,...){