From c19dcaf80c9c0d455c7ea9b777cd945e6989b784 Mon Sep 17 00:00:00 2001 From: Brett McClintock Date: Thu, 8 Aug 2024 17:02:05 -0700 Subject: [PATCH] add dcrwrice and dcrwvm distributions for correlated steps and turn angles; fix bug in calcuating initial velocity terms in crw models when there are multiple individuals --- DESCRIPTION | 3 +- NAMESPACE | 3 + R/CIreal.R | 2 +- R/RcppExports.R | 26 ++++++ R/addSmoothGradient.R | 4 +- R/allProbs.R | 28 +++++- R/checkInputs.R | 25 +++++- R/fitHMM.R | 43 ++++++--- R/fitTMB.R | 101 ++++++++++++++------- R/getDM.R | 33 +++++-- R/get_fixParIndex.R | 51 +++++++---- R/mixtureProbs.R | 5 ++ R/momentuHMM_utils.R | 168 +++++++++++++++++++++++++++++++++-- R/n2w.R | 4 +- R/nLogLike.R | 18 +++- R/parDef.R | 11 +++ R/plot_momentuHMM.R | 40 ++++++++- R/pseudoRes.R | 31 ++++++- R/simData.R | 23 ++++- R/simHierData.R | 13 ++- R/w2n.R | 4 +- man/dcrwrice_rcpp.Rd | 21 +++++ man/dcrwvm_rcpp.Rd | 21 +++++ man/nLogLike.Rd | 5 +- src/RcppExports.cpp | 26 ++++++ src/TMB/include/dist.hpp | 4 + src/TMB/include/dist_def.hpp | 78 ++++++++++++++++ src/combine.h | 10 ++- src/densities.h | 80 +++++++++++++++++ src/momentuHMM_init.c | 6 +- src/nLogLike.cpp | 74 +++++++++++++-- 31 files changed, 853 insertions(+), 108 deletions(-) create mode 100644 man/dcrwrice_rcpp.Rd create mode 100644 man/dcrwvm_rcpp.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 6f7b474..9182b19 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -63,7 +63,8 @@ Suggests: gdistance, methods, mgcv, - optimx + optimx, + pracma RoxygenNote: 7.3.2 Encoding: UTF-8 NeedsCompilation: yes diff --git a/NAMESPACE b/NAMESPACE index b817715..aeb9740 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/CIreal.R b/R/CIreal.R index 4bd6f4e..3028a0d 100644 --- a/R/CIreal.R +++ b/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]]) diff --git a/R/RcppExports.R b/R/RcppExports.R index f5a3c9b..b14f40e 100644 --- a/R/RcppExports.R +++ b/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}}. diff --git a/R/addSmoothGradient.R b/R/addSmoothGradient.R index 32b30cc..984c466 100644 --- a/R/addSmoothGradient.R +++ b/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) diff --git a/R/allProbs.R b/R/allProbs.R index 929f24a..866ac05 100644 --- a/R/allProbs.R +++ b/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)){ diff --git a/R/checkInputs.R b/R/checkInputs.R index 61f1b12..ba3a90b 100644 --- a/R/checkInputs.R +++ b/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,...){ diff --git a/R/fitHMM.R b/R/fitHMM.R index bfaa3c3..9bfb13a 100644 --- a/R/fitHMM.R +++ b/R/fitHMM.R @@ -740,6 +740,19 @@ fitHMM.momentuHMMData <- function(data,nbStates,dist, inputs <- checkInputs(nbStates,dist,Par0,estAngleMean,circularAngleMean,zeroInflation,oneInflation,DM,userBounds,stateNames,checkInflation = TRUE) p <- inputs$p + crwST <- inputs$crwST + if(isTRUE(crwST)){ + stepNA <- which(is.na(data$step)) + angleNA <- which(is.na(data$angle)) + if(length(stepNA)){ + if(!all(stepNA %in% (aInd+table(data$ID)-1))) stop("sorry, 'crwrice' distribution observations for 'step' cannot include missing values (except for the last observation for each individual)") + } + if(!all(is.na(data$angle[aInd]))) stop("The initial turn angle for each individual must be NA") + #if(length(angleNA)){ + # if(!all(angleNA %in% c(1,nrow(data)))) stop("sorry, 'crwvm' distribution observations for 'angle' cannot include missing values") + #} + } + DMinputs<-getDM(data,inputs$DM,inputs$dist,nbStates,p$parNames,p$bounds,Par0,zeroInflation,oneInflation,inputs$circularAngleMean,TMB=isTRUE(optMethod=="TMB")) fullDM <- DMinputs$fullDM DMind <- DMinputs$DMind @@ -754,8 +767,8 @@ fitHMM.momentuHMMData <- function(data,nbStates,dist, #################################### ## Prepare initial values for nlm ## #################################### - fixParIndex <- get_fixParIndex(Par0,beta0,delta0,fixPar,distnames,inputs,p,nbStates,DMinputs,recharge,nbG0covs,nbRecovs,workBounds,mixtures,newformula,formulaStates,nbCovs,betaCons,betaRef,deltaCons,stationary,nbCovsDelta,formulaDelta,formulaPi,nbCovsPi,data,newdata,kappa,optMethod) - + fixParIndex <- get_fixParIndex(Par0,beta0,delta0,fixPar,distnames,inputs,p,nbStates,DMinputs,recharge,nbG0covs,nbRecovs,workBounds,mixtures,newformula,formulaStates,nbCovs,betaCons,betaRef,deltaCons,stationary,nbCovsDelta,formulaDelta,formulaPi,nbCovsPi,data,newdata,kappa,optMethod,crwST) + parCount<- lapply(fullDM,ncol) for(i in distnames[!unlist(lapply(inputs$circularAngleMean,isFALSE))]){ parCount[[i]] <- length(unique(gsub("cos","",gsub("sin","",colnames(fullDM[[i]]))))) @@ -894,7 +907,7 @@ fitHMM.momentuHMMData <- function(data,nbStates,dist, curmod$minimum <- nLogLike(optPar,nbStates,newformula,p$bounds,p$parSize,data,inputs$dist,covs, inputs$estAngleMean,inputs$circularAngleMean,inputs$consensus,zeroInflation,oneInflation, stationary,fullDM,DMind,p$Bndind,knownStates,unlist(fixParIndex$fixPar),fixParIndex$wparIndex, - nc,meanind,covsDelta,workBounds,prior,betaCons,fixParIndex$betaRef,deltaCons,optInd,recovs,g0covs,mixtures,covsPi,hierRecharge,aInd,isTRUE(list(...)$CT),dtIndex,kappa) + nc,meanind,covsDelta,workBounds,prior,betaCons,fixParIndex$betaRef,deltaCons,optInd,recovs,g0covs,mixtures,covsPi,hierRecharge,aInd,isTRUE(list(...)$CT),dtIndex,kappa,crwST) curmod$estimate <- numeric() } else if(optMethod=="TMB"){ if(!is.null(userBounds)) stop("sorry, userBounds cannot be specified when optMethod='TMB'") @@ -902,7 +915,7 @@ fitHMM.momentuHMMData <- function(data,nbStates,dist, if(any(unlist(inputs$consensus))) stop("sorry, consensus models are not currently supported when optMethod='TMB'") if(any(unlist(inputs$circularAngleMean))) stop("sorry, circular-circular regression models are not currently supported when optMethod='TMB'") if(!is.null(hierRecharge)) stop("sorry, recharge models are not currently supported when optMethod='TMB'") - withCallingHandlers(curmod <- tryCatch(fitTMB(data,inputs$dist,nbStates,p,inputs$estAngleMean,oneInflation,zeroInflation,inputs$DM,DMinputs,newformula,fixParIndex,stationary,prior,knownStates,betaCons,control,formulaDelta,covsDelta,workBounds,fit,isTRUE(list(...)$CT),dtIndex,kappa,hessian),error=function(e) e),warning=h) + withCallingHandlers(curmod <- tryCatch(fitTMB(data,inputs$dist,nbStates,p,inputs$estAngleMean,oneInflation,zeroInflation,inputs$DM,DMinputs,newformula,fixParIndex,stationary,prior,knownStates,betaCons,control,formulaDelta,covsDelta,workBounds,fit,isTRUE(list(...)$CT),dtIndex,kappa,hessian,crwST),error=function(e) e),warning=h) if(!inherits(curmod,"error") & (fitCount+1)>retryFits){ if(!is.null(prior)) { tmbPrior <- curmod$prior @@ -915,7 +928,7 @@ fitHMM.momentuHMMData <- function(data,nbStates,dist, withCallingHandlers(curmod <- tryCatch(stats::nlm(nLogLike,optPar,nbStates=nbStates,formula=newformula,bounds=p$bounds,parSize=p$parSize,data=data,dist=inputs$dist,covs=covs, estAngleMean=inputs$estAngleMean,circularAngleMean=inputs$circularAngleMean,consensus=inputs$consensus,zeroInflation=zeroInflation,oneInflation=oneInflation, stationary=stationary,fullDM=fullDM,DMind=DMind,Bndind=p$Bndind,knownStates=knownStates,fixPar=unlist(fixParIndex$fixPar),wparIndex=fixParIndex$wparIndex, - nc=nc,meanind=meanind,covsDelta=covsDelta,workBounds=workBounds,prior=prior,betaCons=betaCons,betaRef=fixParIndex$betaRef,deltaCons=deltaCons,optInd=optInd,recovs=recovs,g0covs=g0covs,mixture=mixtures,covsPi=covsPi,recharge=hierRecharge,aInd=aInd,CT=isTRUE(list(...)$CT),dtIndex=dtIndex,kappa=kappa, + nc=nc,meanind=meanind,covsDelta=covsDelta,workBounds=workBounds,prior=prior,betaCons=betaCons,betaRef=fixParIndex$betaRef,deltaCons=deltaCons,optInd=optInd,recovs=recovs,g0covs=g0covs,mixture=mixtures,covsPi=covsPi,recharge=hierRecharge,aInd=aInd,CT=isTRUE(list(...)$CT),dtIndex=dtIndex,kappa=kappa,crwST=crwST, fscale=fscale,print.level=print.level,gradtol=gradtol, stepmax=stepmax,steptol=steptol, iterlim=iterlim,hessian=ifelse(is.null(nlmPar$hessian) & is.null(prior),TRUE,ifelse(!is.null(prior),FALSE,ifelse(is.null(nlmPar$hessian),TRUE,nlmPar$hessian)))),error=function(e) e),warning=h) @@ -925,7 +938,7 @@ fitHMM.momentuHMMData <- function(data,nbStates,dist, curmod$hessian <- tryCatch(numDeriv::hessian(nLogLike,curmod$estimate,nbStates=nbStates,formula=newformula,bounds=p$bounds,parSize=p$parSize,data=data,dist=inputs$dist,covs=covs, estAngleMean=inputs$estAngleMean,circularAngleMean=inputs$circularAngleMean,consensus=inputs$consensus,zeroInflation=zeroInflation,oneInflation=oneInflation, stationary=stationary,fullDM=fullDM,DMind=DMind,Bndind=p$Bndind,knownStates=knownStates,fixPar=unlist(fixParIndex$fixPar),wparIndex=fixParIndex$wparIndex, - nc=nc,meanind=meanind,covsDelta=covsDelta,workBounds=workBounds,prior=NULL,betaCons=betaCons,betaRef=fixParIndex$betaRef,deltaCons=deltaCons,optInd=optInd,recovs=recovs,g0covs=g0covs,mixture=mixtures,covsPi=covsPi,recharge=hierRecharge,aInd=aInd,CT=isTRUE(list(...)$CT),dtIndex=dtIndex,kappa=kappa),error=function(e) e) + nc=nc,meanind=meanind,covsDelta=covsDelta,workBounds=workBounds,prior=NULL,betaCons=betaCons,betaRef=fixParIndex$betaRef,deltaCons=deltaCons,optInd=optInd,recovs=recovs,g0covs=g0covs,mixture=mixtures,covsPi=covsPi,recharge=hierRecharge,aInd=aInd,CT=isTRUE(list(...)$CT),dtIndex=dtIndex,kappa=kappa,crwST=crwST),error=function(e) e) if(!inherits(curmod$hessian,"error")) dimnames(curmod$hessian) <- list(names(optPar),names(optPar)) } } @@ -933,14 +946,14 @@ fitHMM.momentuHMMData <- function(data,nbStates,dist, withCallingHandlers(curmod <- tryCatch(stats::optim(optPar,nLogLike,gr=NULL,nbStates,newformula,p$bounds,p$parSize,data,inputs$dist,covs, inputs$estAngleMean,inputs$circularAngleMean,inputs$consensus,zeroInflation,oneInflation, stationary,fullDM,DMind,p$Bndind,knownStates,unlist(fixParIndex$fixPar),fixParIndex$wparIndex, - nc,meanind,covsDelta,workBounds,prior,betaCons,fixParIndex$betaRef,deltaCons,optInd,recovs,g0covs,mixtures,covsPi,hierRecharge,aInd,isTRUE(list(...)$CT),dtIndex,kappa, + nc,meanind,covsDelta,workBounds,prior,betaCons,fixParIndex$betaRef,deltaCons,optInd,recovs,g0covs,mixtures,covsPi,hierRecharge,aInd,isTRUE(list(...)$CT),dtIndex,kappa,crwST, method=optMethod,control=control,hessian=ifelse(is.null(prior),hessian,FALSE)),error=function(e) e),warning=h) if(!inherits(curmod,"error")){ if(!is.null(prior) & hessian){ curmod$hessian <- tryCatch(stats::optimHess(optPar,curmod$estimate,gr=NULL,nbStates,newformula,p$bounds,p$parSize,data,inputs$dist,covs, inputs$estAngleMean,inputs$circularAngleMean,inputs$consensus,zeroInflation,oneInflation, stationary,fullDM,DMind,p$Bndind,knownStates,unlist(fixParIndex$fixPar),fixParIndex$wparIndex, - nc,meanind,covsDelta,workBounds,NULL,betaCons,fixParIndex$betaRef,deltaCons,optInd,recovs,g0covs,mixtures,covsPi,hierRecharge,aInd,isTRUE(list(...)$CT),dtIndex,kappa, + nc,meanind,covsDelta,workBounds,NULL,betaCons,fixParIndex$betaRef,deltaCons,optInd,recovs,g0covs,mixtures,covsPi,hierRecharge,aInd,isTRUE(list(...)$CT),dtIndex,kappa,crwST, control=control),error=function(e) e) } } @@ -1033,8 +1046,7 @@ fitHMM.momentuHMMData <- function(data,nbStates,dist, } mod$wpar <- mod$estimate wpar <- mod$estimate <- expandPar(mod$wpar,optInd,fixParIndex$fixPar,fixParIndex$wparIndex,betaCons,deltaCons,nbStates,nbCovsDelta,stationary,nbCovs,nbRecovs+nbG0covs,mixtures,nbCovsPi,isTRUE(optMethod=="TMB"),fixParIndex$Par0,fixParIndex$beta0$beta,fixParIndex$delta0) - mle <- w2n(wpar,p$bounds,p$parSize,nbStates,nbCovs,inputs$estAngleMean,inputs$circularAngleMean,inputs$consensus,stationary,fullDM,DMind,nrow(data),inputs$dist,p$Bndind,nc,meanind,covsDelta,workBounds,covsPi,isTRUE(optMethod=="TMB")) - + if(fit & !is.null(mod$hessian)){ Sigma <- tryCatch(MASS::ginv(mod$hessian),error=function(e) e) mod$Sigma <- Sigma @@ -1043,6 +1055,9 @@ fitHMM.momentuHMMData <- function(data,nbStates,dist, mod$Sigma <- matrix(0,length(mod$estimate),length(mod$estimate)) mod$Sigma[(1:length(mod$estimate))[-optInd],(1:length(mod$estimate))[-optInd]] <- Sigma } + if(crwST){ + mod$Sigma[parindex[["angle"]]+1:parCount$angle,parindex[["angle"]]+1:parCount$angle] <- mod$Sigma[parindex[["step"]]+1:parCount$angle,parindex[["step"]]+1:parCount$angle] + } } else { warning("ginv of the hessian failed -- ",Sigma) } @@ -1053,12 +1068,14 @@ fitHMM.momentuHMMData <- function(data,nbStates,dist, mod$minimum <- nLogLike(optPar,nbStates,newformula,p$bounds,p$parSize,data,inputs$dist,covs, inputs$estAngleMean,inputs$circularAngleMean,inputs$consensus,zeroInflation,oneInflation, stationary,fullDM,DMind,p$Bndind,knownStates,unlist(fixParIndex$fixPar),fixParIndex$wparIndex, - nc,meanind,covsDelta,workBounds,prior,betaCons,fixParIndex$betaRef,deltaCons,optInd,recovs,g0covs,mixtures,covsPi,hierRecharge,aInd,isTRUE(list(...)$CT),dtIndex,kappa) + nc,meanind,covsDelta,workBounds,prior,betaCons,fixParIndex$betaRef,deltaCons,optInd,recovs,g0covs,mixtures,covsPi,hierRecharge,aInd,isTRUE(list(...)$CT),dtIndex,kappa,crwST) mod$estimate <- wpar mod$wpar <- optPar - mle <- w2n(wpar,p$bounds,p$parSize,nbStates,nbCovs,inputs$estAngleMean,inputs$circularAngleMean,inputs$consensus,stationary,fullDM,DMind,nrow(data),inputs$dist,p$Bndind,nc,meanind,covsDelta,workBounds,covsPi,isTRUE(optMethod=="TMB")) } + if(crwST) wpar[parindex[["angle"]]+1:parCount$angle] <- mod$estimate[parindex[["angle"]]+1:parCount$angle] <- wpar[parindex[["step"]]+1:parCount$angle] + mle <- w2n(wpar,p$bounds,p$parSize,nbStates,nbCovs,inputs$estAngleMean,inputs$circularAngleMean,inputs$consensus,stationary,fullDM,DMind,nrow(data),inputs$dist,p$Bndind,nc,meanind,covsDelta,workBounds,covsPi,isTRUE(optMethod=="TMB")) + #################### ## Prepare output ## @@ -1180,7 +1197,7 @@ fitHMM.momentuHMMData <- function(data,nbStates,dist, # conditions of the fit conditions <- list(dist=dist,zeroInflation=zeroInflation,oneInflation=oneInflation, - estAngleMean=inputs$estAngleMean,circularAngleMean=inputs$circularAngleMean,stationary=stationary,formula=formula,userBounds=userBounds,workBounds=workBounds,bounds=p$bounds,Bndind=p$Bndind,DM=DM,fullDM=fullDM,DMind=DMind,fixPar=fixParIndex$ofixPar,wparIndex=fixParIndex$wparIndex,formulaDelta=formulaDelta,betaCons=betaCons,betaRef=fixParIndex$betaRef,deltaCons=deltaCons,optInd=optInd,recharge=recharge,mvnCoords=mvnCoords,mixtures=mixtures,formulaPi=formulaPi,fit=fit,initialValues=initialValues,kappa=kappa,optMethod=optMethod) + estAngleMean=inputs$estAngleMean,circularAngleMean=inputs$circularAngleMean,stationary=stationary,formula=formula,userBounds=userBounds,workBounds=workBounds,bounds=p$bounds,Bndind=p$Bndind,DM=DM,fullDM=fullDM,DMind=DMind,fixPar=fixParIndex$ofixPar,wparIndex=fixParIndex$wparIndex,formulaDelta=formulaDelta,betaCons=betaCons,betaRef=fixParIndex$betaRef,deltaCons=deltaCons,optInd=optInd,recharge=recharge,mvnCoords=mvnCoords,mixtures=mixtures,formulaPi=formulaPi,fit=fit,initialValues=initialValues,kappa=kappa,optMethod=optMethod,crwST=crwST) if(isTRUE(list(...)$CT)) conditions$dtIndex <- dtIndex diff --git a/R/fitTMB.R b/R/fitTMB.R index 45186bf..afa0b05 100644 --- a/R/fitTMB.R +++ b/R/fitTMB.R @@ -5,7 +5,7 @@ #' @importFrom stats update # #' @importFrom mgcv gam # #' @importFrom optimx optimx -fitTMB <- function(data,dist,nbStates,p,estAngleMean,oneInflation,zeroInflation,DM,DMinputs,formula,fixParIndex,stationary,prior,knownStates,betaCons,control,formulaDelta,covsDelta,workBounds,fit,CT,dtIndex,kappa,hessian) { +fitTMB <- function(data,dist,nbStates,p,estAngleMean,oneInflation,zeroInflation,DM,DMinputs,formula,fixParIndex,stationary,prior,knownStates,betaCons,control,formulaDelta,covsDelta,workBounds,fit,CT,dtIndex,kappa,hessian,crwST=FALSE) { for(pkg in c("Matrix","optimx","mgcv")){ if (!requireNamespace(pkg, quietly = TRUE)) { @@ -28,6 +28,10 @@ fitTMB <- function(data,dist,nbStates,p,estAngleMean,oneInflation,zeroInflation, distcode <- as.vector(sapply(distnames,function(x) dist_code(dist[[x]],zeroInflation[[x]],oneInflation[[x]])))#as.vector(sapply(self$obs()$dists(), function(d) d$code())) names(distcode) <- distnames + aInd <- NULL + for(i in 1:nbAnimals) + aInd <- c(aInd,which(data$ID==unique(data$ID)[i])[1]) + if(CT){ dt <- data$dt } else dt <- rep(1,nrow(data)) @@ -117,6 +121,7 @@ fitTMB <- function(data,dist,nbStates,p,estAngleMean,oneInflation,zeroInflation, } fixParIndex$fixPar[[i]] <- reorderFix(fixParIndex$fixPar[[i]]) } + lag <- 0 if(is.list(DM[[i]]) | is.null(DM[[i]])) { DM[i] <- make_formulas(DM[i],i,p$parNames[i],nbStates) if(dist[[i]] %in% rwdists){ @@ -124,6 +129,12 @@ fitTMB <- function(data,dist,nbStates,p,estAngleMean,oneInflation,zeroInflation, for(j in names(DM[[i]])[which(grepl("mean",names(DM[[i]])))]){ for(s in names(DM[[i]][[j]])){ formTerms <- stats::terms(DM[[i]][[j]][[s]],keep.order=TRUE) + lag <- get_crwlag(DM[[i]][[j]][[s]],lag) + if(lag){ + for(jj in all.vars(DM[[i]][[j]][[s]],data)){ + attr(data[[jj]],"aInd") <- aInd + } + } if(!length(attr(formTerms,"term.labels")) && !attr(formTerms,"intercept")){ DM[[i]][[j]][[s]] <- ~1 #DM[[i]][[paste0(j,"_tm1")]][[s]] <- stats::as.formula(paste0("~0+",i,gsub("mean","",j),"_tm1")) @@ -131,7 +142,8 @@ fitTMB <- function(data,dist,nbStates,p,estAngleMean,oneInflation,zeroInflation, } else if(grepl("_tm1",j)){ tInd <- which(!attr(formTerms,"term.labels") %in% paste0(i,gsub("mean","",j))) if(length(tInd)){ - DM[[i]][[j]][[s]] <- stats::as.formula(paste0("~0+",paste0(i,gsub("mean","",j)),paste0("+I(",attr(formTerms,"term.labels")[tInd]," * dt)",collapse=""))) + tmpTerms <- gsub(":","*",attr(formTerms,"term.labels")[tInd]) + DM[[i]][[j]][[s]] <- stats::as.formula(paste0("~0+",paste0(i,gsub("mean","",j)),paste0("+I(",tmpTerms," * dt)",collapse=""))) } } } @@ -449,29 +461,32 @@ fitTMB <- function(data,dist,nbStates,p,estAngleMean,oneInflation,zeroInflation, mod$estimate <- tmb_obj$par } # reorder rwdist parms - parInd <- which(names(mod$estimate)=="coeff_fe_obs") - convInd <- lapply(distnames,function(x) match(DMnames[[x]],names(Par0[[x]])[which(!is.na(fixParIndex$fixPar[[x]]) & !duplicated(fixParIndex$fixPar[[x]]))],nomatch=NA)) - names(convInd) <- distnames - k <- 0 - for(i in distnames){ - if(any(is.na(convInd[[i]]))) convInd[[i]] <- convInd[[i]][-which(is.na(convInd[[i]]))] - if(!all(is.na(fixParIndex$fixPar[[i]]))){ - maxf <- max(convInd[[i]],na.rm=TRUE) - convInd[[i]] <- k + convInd[[i]] - k <- k + maxf + if(!isTRUE(crwST)){ + parInd <- which(names(mod$estimate)=="coeff_fe_obs") + convInd <- lapply(distnames,function(x) match(DMnames[[x]],names(Par0[[x]])[which(!is.na(fixParIndex$fixPar[[x]]) & !duplicated(fixParIndex$fixPar[[x]]))],nomatch=NA)) + names(convInd) <- distnames + k <- 0 + for(i in distnames){ + if(any(is.na(convInd[[i]]))) convInd[[i]] <- convInd[[i]][-which(is.na(convInd[[i]]))] + if(!all(is.na(fixParIndex$fixPar[[i]]))){ + maxf <- max(convInd[[i]],na.rm=TRUE) + convInd[[i]] <- k + convInd[[i]] + k <- k + maxf + } } - } - convInd <- unlist(convInd) - mod$estimate[parInd] <- mod$estimate[convInd] - if(hessian){ - if(max(convInd) 1) { - tmp <- cbind(tmp, obs_var[,1:(wh[i] - 1)]) - tmpnms <- c(tmpnms, colnames(obs_var)[1:(wh[i] - 1)]) + if (wh[i] > 1) { + tmp <- obs_var[,1:(ind - 1)] + tmpnms <- c(tmpnms, colnames(obs_var)[1:(ind - 1)]) } tmp <- cbind(tmp, v) tmpnms <- c(tmpnms, rep(names(wh[i]), ncol(v))) - if (wh < ncol(obs_var)) { - tmp <- cbind(tmp, obs_var[,(wh[i] + 1):ncol(obs_var)]) - tmpnms <- c(tmpnms, colnames(obs_var)[(wh[i] + 1):ncol(obs_var)]) + if (ind < ncol(obs_var)) { + tmp <- cbind(tmp, obs_var[,(ind + 1):ncol(obs_var),drop=FALSE]) + tmpnms <- c(tmpnms, colnames(obs_var)[(ind + 1):ncol(obs_var)]) } obs_var <- tmp colnames(obs_var) <- tmpnms diff --git a/R/getDM.R b/R/getDM.R index 71fe639..4d31d7c 100644 --- a/R/getDM.R +++ b/R/getDM.R @@ -13,6 +13,13 @@ getDM<-function(data,DM,dist,nbStates,parNames,bounds,Par,zeroInflation,oneInfla lanInd <- vector('list',length(dist)) names(lanInd) <- distnames + aInd <- NULL + nbAnimals <- length(unique(data$ID)) + if(nbAnimals){ + for(i in 1:nbAnimals) + aInd <- c(aInd,which(data$ID==unique(data$ID)[i])[1]) + } else aInd <- 1 + for(i in distnames){ if(dist[[i]]=="ctds" && #i==attr(data,"ctdsData") && (is.null(DM[[i]]) || !is.list(DM[[i]]))) stop("DM$",i," must be specified as a list with a formula for parameter 'lambda'") @@ -56,7 +63,12 @@ getDM<-function(data,DM,dist,nbStates,parNames,bounds,Par,zeroInflation,oneInfla for(j in names(DM[[i]])){ tmpCov[[j]] <- vector('list',nbStates) for(state in 1:nbStates){ - if(wlag) lag <- get_crwlag(formulaStates[[j]][[state]],lag) + lag <- get_crwlag(formulaStates[[j]][[state]],lag) + if(lag){ + for(jj in all.vars(formulaStates[[j]][[state]],data)){ + attr(data[[jj]],"aInd") <- aInd + } + } tmpCov[[j]][[state]]<-stats::model.matrix(formulaStates[[j]][[state]],data) if(!isFALSE(circularAngleMean[[i]])){ if(j=="mean"){ @@ -169,12 +181,16 @@ getDM<-function(data,DM,dist,nbStates,parNames,bounds,Par,zeroInflation,oneInfla tmpcovj <- gsub(factorcovs[factorvar][j],tmpcov[j],tmpcovj) } tform <- stats::formula(paste("~ 0 + ",tmpcovj)) - if(wlag) lag <- get_crwlag(tform,lag) + lag <- get_crwlag(tform,lag) + if(lag){ + for(jj in all.vars(tform,data)){ + attr(data[[jj]],"aInd") <- aInd + } + } tmpcovs<-stats::model.matrix(tform,data) tmpcovs<-tmpcovs[,which(gsub(" ","",colnames(tmpcovs)) %in% gsub(" ","",cov))] covs<-cbind(covs,tmpcovs) } else { - if(wlag) lag <- get_crwlag(form,lag) #if(!grepl("langevin\\(",cov)){ # still need to figure out solution for objects in formula that are not in data for langevin models tmpCol <- tryCatch(stats::get_all_vars(form,data),error=function(e) e)#rownames(attr(stats::terms(formula),"factors"))#attr(stats::terms(formula),"term.labels")#seq(1,ncol(data))[-match(c("ID","x","y",distnames),names(data),nomatch=0)] if(inherits(tmpCol,"error")){ @@ -192,6 +208,12 @@ getDM<-function(data,DM,dist,nbStates,parNames,bounds,Par,zeroInflation,oneInfla } rm(tmpCol) #} + lag <- get_crwlag(form,lag) + if(lag){ + for(jj in all.vars(form,data)){ + attr(data[[jj]],"aInd") <- aInd + } + } tmpcovs<-stats::model.matrix(form,data)[,2] covs<-cbind(covs,tmpcovs) } @@ -269,7 +291,7 @@ getDM<-function(data,DM,dist,nbStates,parNames,bounds,Par,zeroInflation,oneInfla return(list(fullDM=simpDM,DMind=DMind,lag=lag,langInd = any(!is.null(unlist(lanInd))))) } -get_crwlag <- function(form,clag){ +get_crwlag <- function(form,clag=0){ lag <- 0 if(length(attr(stats::terms(form, specials="crw"),"specials")$crw)){ if (!requireNamespace("prodlim", quietly = TRUE)) { @@ -278,7 +300,8 @@ get_crwlag <- function(form,clag){ } Terms <- stats::terms(form,specials="crw") tm1 <- attr(prodlim::strip.terms(Terms[attr(Terms,"specials")$crw],specials="crw",arguments=list(crw=list("lag"=1,"dt"=1))),"term.labels") - lag <- as.numeric(attr(prodlim::strip.terms(Terms[attr(Terms,"specials")$crw],specials="crw",arguments=list(crw=list("lag"=1,"dt"=1))),"stripped.arguments")$crw[[tm1]]$lag) + strTerm <- attr(prodlim::strip.terms(Terms[attr(Terms,"specials")$crw],specials="crw",arguments=list(crw=list("lag"=1,"dt"=1))),"stripped.arguments")$crw + lag <- max(as.numeric(unlist(strTerm)[which( names(unlist(strTerm)) %in% paste0(tm1,".lag"))])) } max(c(clag,lag)) } diff --git a/R/get_fixParIndex.R b/R/get_fixParIndex.R index 2900f3c..f30c4f4 100644 --- a/R/get_fixParIndex.R +++ b/R/get_fixParIndex.R @@ -1,6 +1,17 @@ -get_fixParIndex <- function(Par0,beta0,delta0,fixPar,distnames,inputs,p,nbStates,DMinputs,recharge,nbG0covs,nbRecovs,workBounds,mixtures,newformula,formulaStates,nbCovs,betaCons,betaRef,deltaCons,stationary,nbCovsDelta,formulaDelta,formulaPi,nbCovsPi,data,newdata,kappa=Inf,optMethod){ +get_fixParIndex <- function(Par0,beta0,delta0,fixPar,distnames,inputs,p,nbStates,DMinputs,recharge,nbG0covs,nbRecovs,workBounds,mixtures,newformula,formulaStates,nbCovs,betaCons,betaRef,deltaCons,stationary,nbCovsDelta,formulaDelta,formulaPi,nbCovsPi,data,newdata,kappa=Inf,optMethod,crwST=FALSE){ + if(isTRUE(crwST)){ + if(!is.null(fixPar)){ + if(!is.null(fixPar$angle)){ + if(!is.null(fixPar$step)) { + if(!isTRUE(all.equal(fixPar$step,fixPar$angle))) stop("fixPar$step and fixPar$angle should be the same for 'crwrice' and 'crwvm' distributions") + } + } else fixPar$angle <- fixPar$step + } + } + if(optMethod!="TMB"){ + if(nbStates>1){ if(!is.null(betaCons)){ if(!is.matrix(betaCons)) stop("betaCons must be a matrix") @@ -149,8 +160,12 @@ get_fixParIndex <- function(Par0,beta0,delta0,fixPar,distnames,inputs,p,nbStates parindex <- c(0,cumsum(unlist(lapply(Par0,length)))) names(parindex) <- c(distnames,"beta") ofixPar <- fixPar + for(i in distnames){ - if(!is.null(fixPar[[i]])){ + if(isTRUE(crwST) & i=="angle"){ + tmp <- 1:length(Par0[[i]]) + wparIndex <- c(wparIndex,parindex[[i]]+tmp) + } else if(!is.null(fixPar[[i]])){ if(length(fixPar[[i]])!=length(Par0[[i]])) stop("fixPar$",i," must be of length ",length(Par0[[i]])) tmp <- which(!is.na(fixPar[[i]])) Par0[[i]][tmp]<-fixPar[[i]][tmp] @@ -357,7 +372,10 @@ get_fixParIndex <- function(Par0,beta0,delta0,fixPar,distnames,inputs,p,nbStates names(parindex) <- c(distnames,"beta") ofixPar <- fixPar for(i in distnames){ - if(!is.null(fixPar[[i]])){ + if(isTRUE(crwST) & i=="angle"){ + tmp <- 1:length(Par0[[i]]) + wparIndex <- c(wparIndex,parindex[[i]]+tmp) + } else if(!is.null(fixPar[[i]])){ if(length(fixPar[[i]])!=length(Par0[[i]])) stop("fixPar$",i," must be of length ",length(Par0[[i]])) tmp <- sort(which(duplicated(fixPar[[i]]) | is.na(fixPar[[i]]))) wparIndex <- c(wparIndex,parindex[[i]]+tmp) @@ -396,19 +414,21 @@ get_fixParIndex <- function(Par0,beta0,delta0,fixPar,distnames,inputs,p,nbStates fixPar <- fixPar[parVec] ofixPar <- ofixPar[parVec] for(i in names(fixPar)){ - tmp <- c(fixPar[[i]]) - if(!all(is.na(tmp))){ - if(!all(tmp[which(!duplicated(tmp) & !is.na(tmp))]==(1:length(tmp[which(!duplicated(tmp) & !is.na(tmp))])))){ - k <- 1 - newtmp <- tmp - for(j in unique(tmp)){ - if(!is.na(j)){ - newtmp[which(tmp==j)] <- k - k <- k + 1 + if(!(isTRUE(crwST) && i=="angle")){ + tmp <- c(fixPar[[i]]) + if(!all(is.na(tmp))){ + if(!all(tmp[which(!duplicated(tmp) & !is.na(tmp))]==(1:length(tmp[which(!duplicated(tmp) & !is.na(tmp))])))){ + k <- 1 + newtmp <- tmp + for(j in unique(tmp)){ + if(!is.na(j)){ + newtmp[which(tmp==j)] <- k + k <- k + 1 + } } - } - stop("fixPar$",i," entries that are not fixed (i.e. NA) must be ordered sequentially starting from 1.\n Should it be c(",paste0(newtmp,collapse=", "),")?") - } + stop("fixPar$",i," entries that are not fixed (i.e. NA) must be ordered sequentially starting from 1.\n Should it be c(",paste0(newtmp,collapse=", "),")?") + } + } } } k <- 0 @@ -419,6 +439,7 @@ get_fixParIndex <- function(Par0,beta0,delta0,fixPar,distnames,inputs,p,nbStates k <- k + maxf } } + if(isTRUE(crwST)) fixPar$angle <- fixPar$step delta0 <- logdelta0 } diff --git a/R/mixtureProbs.R b/R/mixtureProbs.R index 0461195..5260c5a 100644 --- a/R/mixtureProbs.R +++ b/R/mixtureProbs.R @@ -227,6 +227,11 @@ get_mixProbs <- function(optPar,mod,mixture){ } par[["pi"]] <- matrix(1,1,1) + #if(isTRUE(mod$conditions$crwST)) { + # dist$step <- "crwST" + # dist$angle <- NULL + #} + mixProbs <- lnum <- la <- numeric(mixtures) for(mix in 1:mixtures){ par$beta <- beta[(mix-1)*(nbCovs+1)+1:(nbCovs+1),,drop=FALSE] diff --git a/R/momentuHMM_utils.R b/R/momentuHMM_utils.R index 1239206..c65ebd7 100644 --- a/R/momentuHMM_utils.R +++ b/R/momentuHMM_utils.R @@ -1,10 +1,10 @@ -momentuHMMdists<-sort(c('gamma','weibull','exp','lnorm','beta','pois','wrpcauchy','vm','norm','bern','vmConsensus','mvnorm2','mvnorm3','rw_mvnorm2','rw_mvnorm3','rw_norm','cat','negbinom','logis','t','ctds')) +momentuHMMdists<-sort(c('gamma','weibull','exp','lnorm','beta','pois','wrpcauchy','vm','norm','bern','vmConsensus','mvnorm2','mvnorm3','rw_mvnorm2','rw_mvnorm3','rw_norm','cat','negbinom','logis','t','ctds','crwrice','crwvm')) moveHMMdists<-sort(c('gamma','weibull','exp','lnorm','wrpcauchy','vm')) -angledists<-sort(c('wrpcauchy','vm','vmConsensus')) -stepdists<-sort(c('gamma','weibull','exp','lnorm')) +angledists<-sort(c('wrpcauchy','vm','vmConsensus','crwvm')) +stepdists<-sort(c('gamma','weibull','exp','lnorm','crwrice')) singleParmdists<-sort(c('exp','pois','bern')) -nonnegativedists<-sort(c('gamma','weibull','exp','lnorm','pois','negbinom')) -zeroInflationdists<-sort(c('gamma','weibull','exp','lnorm','beta')) +nonnegativedists<-sort(c('gamma','weibull','exp','lnorm','pois','negbinom','crwrice')) +zeroInflationdists<-sort(c('gamma','weibull','exp','lnorm','beta','crwrice')) oneInflationdists<-sort(c('beta')) integerdists<-sort(c('bern','pois','cat','negbinom','ctds')) mvndists <- c('mvnorm2','mvnorm3','rw_mvnorm2','rw_mvnorm3') @@ -17,7 +17,7 @@ meansListNoTime<-c("numeric","integer") modeList<- c("factor","logical") plotArgs <- c("cex","cex.main","cex.lab","cex.axis","cex.legend","lwd","asp","legend.pos") fitMethods<-c("nlm","TMB","Nelder-Mead","SANN") -badNames <- c("beta", "delta", "pi", "g0", "theta") +badNames <- c("beta", "delta", "pi", "g0", "theta","dt") #' @importFrom stats dbinom dbern <- function (x, prob, log = FALSE) @@ -55,6 +55,151 @@ rnegbinom<- function (n, mu, size) return(stats::rnbinom(n, size = size, mu = mu)) } +dcrwrice <- function(x, beta, sigma, step_tm1 = 0, log = FALSE){ + + n <- length(x) + if(length(beta)==1) beta <- rep(beta,n) + if(length(sigma)==1) sigma <- rep(sigma,n) + if(n>1 & length(beta)!=n) stop("'beta' should be of length ",n) + if(n>1 & length(sigma)!=n) stop("'sigma' should be of length ",n) + + mu <- step_tm1*(1-exp(-beta))/beta + var <- sigma*sigma/(beta*beta) * (1. - 2. / beta * (1.-exp(-beta))+1. / (2.*beta)*(1-exp(-2*beta))); + sd <- sqrt(var) + xabs <- abs(x * mu/var) + logdens <- log(x) - 2 * log(sd) + (-(x^2 + mu^2)/(2 * var)) + log(besselI(xabs, nu = 0, expon.scaled = TRUE)) + xabs + if (log) + logdens + else exp(logdens) +} + +#' @importFrom stats pchisq +pcrwrice <- function (q, beta, sigma, step_tm1 = 0, lower.tail = TRUE, log.p = FALSE) +{ + n <- length(q) + if(length(beta)==1) beta <- rep(beta,n) + if(length(sigma)==1) sigma <- rep(sigma,n) + if(n>1 & length(beta)!=n) stop("'beta' should be of length ",n) + if(n>1 & length(sigma)!=n) stop("'sigma' should be of length ",n) + + mu <- step_tm1*(1-exp(-beta))/beta + var <- sigma*sigma/(beta*beta) * (1. - 2. / beta * (1.-exp(-beta))+1. / (2.*beta)*(1-exp(-2*beta))); + sd <- sqrt(var) + a <- mu/sd + b <- q/sd + m <- 1 + + out <- stats::pchisq(b^2, df = 2 * m, ncp = a^2, lower.tail = lower.tail, log.p = log.p) + return(out) +} + +#' @importFrom stats rnorm +rcrwrice <- function(n, beta, sigma, step_tm1=0) +{ + steps <- numeric(n+1) + steps[1] <- step_tm1 + + sd <- sqrt(sigma*sigma/(beta*beta) * (1. - 2. / beta * (1.-exp(-beta))+1. / (2.*beta)*(1-exp(-2*beta)))) + theta <- 1 + + for(i in 2:(n+1)){ + mu <- steps[i-1]*(1-exp(-beta))/beta + + X <- rnorm(1, mean = mu * cos(theta), sd = sd) + Y <- rnorm(1, mean = mu * sin(theta), sd = sd) + steps[i] <- sqrt(X^2 + Y^2) + } + return(steps[-1]) +} + +#' @importFrom stats integrate approxfun +intdcrwrice <- function(step, beta, sigma, xmin, xmax, stepDens, log = FALSE){ + n <- length(step) + out <- numeric(n) + for(i in 1:n){ + out[i] <- integrate(function(step_tm1, step, beta, sigma, log = FALSE){ + + mu <- step_tm1*(1-exp(-beta))/beta + var <- sigma*sigma/(beta*beta) * (1. - 2. / beta * (1.-exp(-beta))+1. / (2.*beta)*(1-exp(-2*beta))); + sd <- sqrt(var) + xabs <- abs(step[i] * mu/var) + logdens <- log(step[i]) - 2 * log(sd) + (-(step[i]^2 + mu^2)/(2 * var)) + log(besselI(xabs, nu = 0, expon.scaled = TRUE)) + xabs + logdens <- logdens + log(stepDens(step_tm1)) + if (log) + logdens + else exp(logdens) + + },lower=xmin,upper=xmax,step=step,beta=beta,sigma=sigma)$value + } + if(log) out <- log(out) + return(out) +} + +dcrwvm <- function(x, beta, sigma, steps, steps_tm1, log = FALSE) +{ + n <- length(x) + + if(missing(steps)) stop("'steps' must be provided") + if(missing(steps_tm1)) stop("'steps_tm1' must be provided") + if(length(steps)!=n) stop("'steps' must be of length ",n) + if(length(steps_tm1)!=n) stop("'steps_tm1' must be of length ",n) + if(length(beta)==1) beta <- rep(beta,n) + if(length(sigma)==1) sigma <- rep(sigma,n) + + var <- sigma*sigma/(beta*beta) * (1. - 2. / beta * (1.-exp(-beta))+1. / (2.*beta)*(1-exp(-2*beta))); + kappa <- steps_tm1*steps*exp(-beta)/var + b <- besselI(kappa, nu = 0, expon.scaled = TRUE) + logdens <- -log(2*pi*b) + kappa * (cos(x)-1) + if(log) logdens + else exp(logdens) +} + +intdcrwvm <- function(angle, beta, sigma, xmin, xmax, stepDens, log = FALSE){ + n <- length(angle) + out <- numeric(n) + for(i in 1:n){ + out[i] <- pracma::integral2(function(step_tm1, step, angle, beta, sigma, log = FALSE){ + + var <- sigma*sigma/(beta*beta) * (1. - 2. / beta * (1.-exp(-beta))+1. / (2.*beta)*(1-exp(-2*beta))); + kappa <- step_tm1*step*exp(-beta)/var + b <- besselI(kappa, nu = 0, expon.scaled = TRUE) + logdens <- -log(2*pi*b) + kappa * (cos(angle)-1) + logdens <- logdens + log(stepDens(step_tm1)) + log(stepDens(step)) + if(log) logdens + else exp(logdens) + },xmin=xmin,xmax=xmax,ymin=xmin,ymax=xmax,angle=angle[i],beta=beta,sigma=sigma,singular = TRUE)$Q + } + return(out) +} + +#' @importFrom CircStats rvm +rcrwvm <- function(n, beta, sigma, steps) +{ + if(missing(steps)) steps <- rcrwrice(2,beta,sigma) + else if(length(steps)!=2) stop("'steps' must be of length 2") + + var <- sigma*sigma/(beta*beta) * (1. - 2. / beta * (1.-exp(-beta))+1. / (2.*beta)*(1-exp(-2*beta))); + kappa <- steps[-1]*steps[-length(steps)]*exp(-beta)/var + return(CircStats::rvm(n,0,kappa)) +} + +pcrwvm <- function(x, beta, sigma, step_tm1 = 0, log = FALSE) +{ + n <- nrow(x) + steps <- x[,1] + angles <- x[,2] + + if(length(beta)==1) beta <- rep(beta,n) + if(length(sigma)==1) sigma <- rep(sigma,n) + + var <- sigma*sigma/(beta*beta) * (1. - 2. / beta * (1.-exp(-beta))+1. / (2.*beta)*(1-exp(-2*beta))); + kappa <- c(step_tm1*steps[1]*exp(-beta[1])/var[1],steps[-1]*steps[-n]*exp(-beta[-1])/var[-1]) + b <- besselI(kappa, nu = 0, expon.scaled = TRUE) + logdens <- -log(2*pi*b) + kappa * (cos(angles)-1) + if(log) logdens + else exp(logdens) +} + dmvnorm2 <- dmvnorm3 <- drw_mvnorm2 <- drw_mvnorm3 <- function(x,mean,sigma){ dmvnorm_rcpp(x,mean,sigma) } @@ -163,12 +308,17 @@ crw <- function(x_tm1,lag=1,dt){ call. = FALSE) } } - if(missing(dt)) dt <- rep(1,length(x_tm1)) + if(missing(dt) || is.function(dt)) dt <- rep(1,length(x_tm1)) else if(length(dt)==1) dt <- rep(dt,length(x_tm1)) else if(length(dt)!=length(x_tm1)) stop("argument 'dt' in 'crw' function must be of length 1 or ",length(x_tm1)) - diff <- dplyr::lag(x_tm1,n=lag-1,default=x_tm1[1])-dplyr::lag(x_tm1,n=lag,default=x_tm1[1]) - vel <- diff/dplyr::lag(dt,n=lag,default=0) + aInd <- c(attr(x_tm1,"aInd"),length(x_tm1)+1) + nbAnimals <- length(aInd) - 1 + diff <- vel <- numeric(length(x_tm1)) + for(i in 1:nbAnimals){ + diff[aInd[i]:(aInd[i+1]-1)] <- dplyr::lag(x_tm1[aInd[i]:(aInd[i+1]-1)],n=lag-1,default=x_tm1[aInd[i]])-dplyr::lag(x_tm1[aInd[i]:(aInd[i+1]-1)],n=lag,default=x_tm1[aInd[i]]) + vel[aInd[i]:(aInd[i+1]-1)] <- diff[aInd[i]:(aInd[i+1]-1)]/dplyr::lag(dt[aInd[i]:(aInd[i+1]-1)],n=lag,default=0) + } if(any(!is.finite(vel))){ vel[which(!is.finite(vel))] <- 0 } diff --git a/R/n2w.R b/R/n2w.R index 4dd2901..9238682 100644 --- a/R/n2w.R +++ b/R/n2w.R @@ -54,7 +54,7 @@ n2w <- function(par,bounds,beta,delta=NULL,nbStates,estAngleMean,DM,Bndind,dist, for(i in names(par)){ p <- par[[i]] if(is.null(DM[[i]])){ - if(estAngleMean[[i]] & Bndind[[i]]){ + if(estAngleMean[[i]] & Bndind[[i]] & dist[[i]]!="crwvm"){ if(length(which(par[[i]][1:nbStates]<=bounds[[i]][1:nbStates,1] | par[[i]][1:nbStates]>bounds[[i]][1:nbStates,2]))>0) stop(paste0("Check the parameter bounds for ",i," (the initial parameters should be ", "strictly between the bounds of their parameter space). The angle mean should be in (-pi,pi].")) @@ -66,7 +66,7 @@ n2w <- function(par,bounds,beta,delta=NULL,nbStates,estAngleMean,DM,Bndind,dist, stop(paste0("Check the parameter bounds for ",i," (the initial parameters should be ", "strictly between the bounds of their parameter space).")) } - if(estAngleMean[[i]] & Bndind[[i]] & !TMB){ + if(estAngleMean[[i]] & Bndind[[i]] & !TMB & dist[[i]]!="crwvm"){ bounds[[i]][,1] <- -Inf bounds[[i]][which(bounds[[i]][,2]!=1),2] <- Inf diff --git a/R/nLogLike.R b/R/nLogLike.R index fb4d0a7..c9adc28 100644 --- a/R/nLogLike.R +++ b/R/nLogLike.R @@ -45,6 +45,7 @@ #' @param CT logical indicating whether to fit discrete-time approximation of a continuous-time model #' @param dtIndex time difference index for calculating transition probabilities of hierarchical continuous-time models #' @param kappa maximum allowed value for the row sums of the off-diagonal elements in the state transition rate matrix, such that the minimum value for the diagonal elements is \code{-kappa}. Default: \code{Inf}. Setting less than \code{Inf} can help avoid numerical issues during optimization, in which case the transition rate parameters \code{beta} are on the logit scale (instead of the log scale). Ignored unless \code{CT=TRUE}. +#' @param crwST logical indicating whether to fit a correlated step and turn random walk model #' #' @return The negative log-likelihood of the parameters given the data. #' @@ -79,12 +80,16 @@ nLogLike <- function(optPar,nbStates,formula,bounds,parSize,data,dist,covs, estAngleMean,circularAngleMean,consensus,zeroInflation,oneInflation, - stationary=FALSE,fullDM,DMind,Bndind,knownStates,fixPar,wparIndex,nc,meanind,covsDelta,workBounds,prior=NULL,betaCons=NULL,betaRef,deltaCons=NULL,optInd=NULL,recovs=NULL,g0covs=NULL,mixtures=1,covsPi,recharge=NULL,aInd=aInd,CT=FALSE,dtIndex=NULL,kappa=Inf) + stationary=FALSE,fullDM,DMind,Bndind,knownStates,fixPar,wparIndex,nc,meanind,covsDelta,workBounds,prior=NULL,betaCons=NULL,betaRef,deltaCons=NULL,optInd=NULL,recovs=NULL,g0covs=NULL,mixtures=1,covsPi,recharge=NULL,aInd=aInd,CT=FALSE,dtIndex=NULL,kappa=Inf,crwST=FALSE) { # check arguments - distnames<-names(dist) - + if(crwST) { + odist <- dist + dist$step <- "crwST" + dist$angle <- NULL + } + #covs <- stats::model.matrix(formula,data) nbCovs <- ncol(covs)-1 if(!is.null(recovs)) { @@ -96,6 +101,11 @@ nLogLike <- function(optPar,nbStates,formula,bounds,parSize,data,dist,covs, wpar <- expandPar(optPar,optInd,fixPar,wparIndex,betaCons,deltaCons,nbStates,ncol(covsDelta)-1,stationary,nbCovs,nbRecovs+nbG0covs,mixtures,ncol(covsPi)-1) par <- tryCatch(w2n(wpar,bounds,parSize,nbStates,nbCovs,estAngleMean,circularAngleMean,consensus,stationary,fullDM,DMind,nrow(data),dist,Bndind,nc,meanind,covsDelta,workBounds,covsPi),error=function(e) e) + if(crwST) { + dist <- odist + par$angle <- par$step[1:(2*nbStates),] + } + if(inherits(par,"error")){ return(NaN) } else { @@ -136,6 +146,8 @@ nLogLike <- function(optPar,nbStates,formula,bounds,parSize,data,dist,covs, par[which(unlist(lapply(dist,is.null)))]<-matrix(NA) } + distnames<-names(dist) + if(stationary) par$delta <- c(NA) if(nbStates==1) { diff --git a/R/parDef.R b/R/parDef.R index f9bfb4a..26f9828 100644 --- a/R/parDef.R +++ b/R/parDef.R @@ -63,6 +63,17 @@ parDef <- function(dist,nbStates,estAngleMean,zeroInflation,oneInflation,DM,user tmpbounds <- matrix(rep(c(0,1),parSize[[i]] * nbStates),ncol=2,byrow=TRUE) parNames[[i]]<-paste0("prob",1:parSize[[i]]) }, + "crwrice"={ + parSize[[i]] <- 2 + zeroInflation[[i]] + tmpbounds <- matrix(rep(c(0,Inf),parSize[[i]] * nbStates),ncol=2,byrow=TRUE) + parNames[[i]]<-c("beta","sigma") + }, + "crwvm"={ + estAngleMean[[i]] <- FALSE + parSize[[i]] <- 2 + tmpbounds <- matrix(rep(c(0,Inf),parSize[[i]] * nbStates),ncol=2,byrow=TRUE) + parNames[[i]]<-c("beta","sigma") + }, "ctds"={ parSize[[i]] <- dimCat-1 tmpbounds <- matrix(rep(c(0,Inf),parSize[[i]] * nbStates),ncol=2,byrow=TRUE) diff --git a/R/plot_momentuHMM.R b/R/plot_momentuHMM.R index 5c965f9..1af4e06 100644 --- a/R/plot_momentuHMM.R +++ b/R/plot_momentuHMM.R @@ -50,7 +50,7 @@ #' @export #' @importFrom graphics legend lines segments arrows layout image contour barplot plot.new #' @importFrom grDevices adjustcolor gray hcl colorRampPalette recordPlot -#' @importFrom stats as.formula +#' @importFrom stats as.formula density approxfun #' @importFrom CircStats circ.mean # @importFrom scatterplot3d scatterplot3d #' @importFrom MASS kde2d @@ -354,6 +354,11 @@ plot.momentuHMM <- function(x,animals=NULL,covs=NULL,ask=TRUE,breaks="Sturges",h call. = FALSE) } dcat <- dctds <- extraDistr::dcat + } else if(Fun[[i]]=="dcrwvm"){ + if (!requireNamespace("pracma", quietly = TRUE)) { + stop("Package \"pracma\" needed for correlated random walk von Mises ('crwvm') distribution. Please install it.", + call. = FALSE) + } } } @@ -596,8 +601,24 @@ plot.momentuHMM <- function(x,animals=NULL,covs=NULL,ask=TRUE,breaks="Sturges",h } # (weighted by the proportion of each state in the Viterbi states sequence) if(m$conditions$zeroInflation[[i]] | m$conditions$oneInflation[[i]]){ + if(inputs$dist[[i]]=="crwrice"){ + genFun <- "intdcrwrice" + genArgs[[4]] <- min(m$data[[i]],na.rm=TRUE) + genArgs[[5]] <- max(m$data[[i]],na.rm=TRUE) + d <- density(m$data[[i]],na.rm=TRUE) + genArgs[[6]] <- stats::approxfun(d$x,d$y) + message("Integrating over step length sample distribution for 'step' state ",state,"...") + } genDensities[[state]] <- cbind(grid,(1-zeroMass[[i]][state]-oneMass[[i]][state])*w[[i]][state]*do.call(genFun,genArgs)) } else if(infInd) { + if(inputs$dist[[i]]=="crwvm"){ + genFun <- "intdcrwvm" + genArgs[[4]] <- min(m$data$step,na.rm=TRUE) + genArgs[[5]] <- max(m$data$step,na.rm=TRUE) + d <- density(m$data$step,na.rm=TRUE) + genArgs[[6]] <- stats::approxfun(d$x,d$y) + message("Integrating over step length sample distribution for 'angle' state ",state,"...") + } genDensities[[state]] <- cbind(grid,(1-zeroMass$step[state])*w[[i]][state]*do.call(genFun,genArgs)) } else if(inputs$dist[[i]] %in% mvndists){ if(inputs$dist[[i]]=="mvnorm2" || inputs$dist[[i]]=="rw_mvnorm2"){ @@ -605,6 +626,23 @@ plot.momentuHMM <- function(x,animals=NULL,covs=NULL,ask=TRUE,breaks="Sturges",h genDensities[[state]] <- list(x=genArgs[[1]][1:100], y=genArgs[[1]][101:200], z=w[[i]][state]*dens) } } else { + if(inputs$dist[[i]]=="crwrice" | inputs$dist[[i]]=="crwvm"){ + if(inputs$dist[[i]]=="crwrice"){ + genFun <- "intdcrwrice" + genArgs[[4]] <- min(m$data[[i]],na.rm=TRUE) + genArgs[[5]] <- max(m$data[[i]],na.rm=TRUE) + d <- density(m$data[[i]],na.rm=TRUE) + genArgs[[6]] <- stats::approxfun(d$x,d$y) + message("Integrating over step length sample distribution for 'step' state ",state,"...") + } else { + genFun <- "intdcrwvm" + genArgs[[4]] <- min(m$data$step,na.rm=TRUE) + genArgs[[5]] <- max(m$data$step,na.rm=TRUE) + d <- density(m$data$step,na.rm=TRUE) + genArgs[[6]] <- stats::approxfun(d$x,d$y) + message("Integrating over step length sample distribution for 'angle' state ",state,"...") + } + } genDensities[[state]] <- cbind(grid,w[[i]][state]*do.call(genFun,genArgs)) } diff --git a/R/pseudoRes.R b/R/pseudoRes.R index 7695617..e3b67c5 100644 --- a/R/pseudoRes.R +++ b/R/pseudoRes.R @@ -203,11 +203,20 @@ pseudoRes <- function(m, ncores = 1) d2<-function(q,mean,sigma){ stats::pchisq(stats::mahalanobis(q,mean,convertSigma(matrix(sigma,length(mean),length(mean)),length(mean))),df=length(mean)) } + genInd <- which(!is.na(genData[1:nbObs])) + } else if(dist[[j]]=="crwrice" | dist[[j]]=="crwvm"){ + genData <- data[[j]] + step_tm1 <- c(NA,data$step[-length(data$step)]) + if(dist[[j]]=="crwrice"){ + genData[c(1,1+cumsum(table(m$data$ID)[1:(length(unique(m$data$ID))-1)]))] <- NA + genInd <- which(!is.na(genData)) + } else genInd <- which(!is.na(data[[j]])) + step_tm1[-genInd] <- NA } else { genData <- data[[j]] + genInd <- which(!is.na(genData[1:nbObs])) } - genInd <- which(!is.na(genData[1:nbObs])) zeroInflation <- m$conditions$zeroInflation[[j]] oneInflation <- m$conditions$oneInflation[[j]] @@ -217,7 +226,8 @@ pseudoRes <- function(m, ncores = 1) if(!(dist[[j]] %in% angledists)){ - genArgs <- list(genData[which(!is.na(genData))]) + if(dist[[j]] %in% mvndists) genArgs <- list(genData[which(!is.na(genData))]) + else genArgs <- list(genData[genInd]) zeromass <- 0 onemass <- 0 @@ -271,6 +281,10 @@ pseudoRes <- function(m, ncores = 1) genArgs[[3]] <- 1/scale # dgamma expects rate=1/scale } + if(dist[[j]]=="crwrice"){ + genArgs[[k+2]] <- step_tm1[genInd] + } + if(zeroInflation | oneInflation) { if(zeroInflation & !oneInflation){ pgenMat[genInd,state] <- ifelse(genData[genInd]==0, @@ -299,7 +313,7 @@ pseudoRes <- function(m, ncores = 1) } else { - genpiInd <- which(genData!=pi & !is.na(genData)) + genpiInd <- union(which(genData!=pi),genInd) genArgs <- list(Fun[[j]],-pi,genData[1]) # to pass to function "integrate" below @@ -307,7 +321,16 @@ pseudoRes <- function(m, ncores = 1) genArgs[[3]]<-genData[i] for(k in 1:(nrow(genPar)/nbStates)) genArgs[[k+3]] <- genPar[(k-1)*nbStates+state,i] - + if(dist[[j]]=="crwvm"){ + mu <- 0 + beta <- genPar[state,i] + sigma <- genPar[nbStates+state,i] + var <- sigma*sigma/(beta*beta) * (1. - 2. / beta * (1.-exp(-beta))+1. / (2.*beta)*(1-exp(-2*beta))) + kappa <- step_tm1[i] * m$data$step[i] *exp(-beta)/var + genArgs[[4]] <- 0 + genArgs[[5]] <- kappa + genArgs[[1]] <- "dvm" + } pgenMat[i,state] <- do.call(integrate,genArgs)$value } } diff --git a/R/simData.R b/R/simData.R index 5429e1b..2febb61 100644 --- a/R/simData.R +++ b/R/simData.R @@ -403,6 +403,8 @@ simData <- function(nbAnimals=1,nbStates=2,dist, rmvnorm3 <- rmvnorm3 rvm <- CircStats::rvm rwrpcauchy <- CircStats::rwrpcauchy + rcrwrice <- rcrwrice + rcrwvm <- rcrwvm trMatrix_rcpp <- trMatrix_rcpp ctPar <- ctPar if (requireNamespace("extraDistr", quietly = TRUE)){ @@ -661,6 +663,7 @@ simData <- function(nbAnimals=1,nbStates=2,dist, zi <- FALSE if(!is.null(zeroInflation$step)) zi <- zeroInflation$step if(is.null(dist$angle)) dist$angle<-"none" + if(nbStates==1) states <- FALSE # work around bug in moveHMM data <- moveHMM::simData(nbAnimals, nbStates, dist$step, dist$angle, Par$step, Par$angle, beta, covs, nbCovs, zi, obsPerAnimal, model, states) attr(data,"class") <- "data.frame" data$ID <- as.factor(data$ID) @@ -730,8 +733,10 @@ simData <- function(nbAnimals=1,nbStates=2,dist, estAngleMean <- vector('list',length(distnames)) names(estAngleMean) <- distnames for(i in distnames){ - if(dist[[i]] %in% angledists) estAngleMean[[i]]<-TRUE - else estAngleMean[[i]]<-FALSE + if(dist[[i]] %in% angledists) { + estAngleMean[[i]]<-TRUE + if(dist[[i]]=="crwvm") estAngleMean[[i]]<-FALSE + } else estAngleMean[[i]]<-FALSE } maxRate <- Inf @@ -1072,7 +1077,7 @@ simData <- function(nbAnimals=1,nbStates=2,dist, allStates <- NULL allSpatialcovs<-NULL - #make sure 'step' preceeds 'angle' + #make sure 'step' precedes 'angle' if(all(c("step","angle") %in% distnames)){ distnames<-c("step","angle",distnames[!(distnames %in% c("step","angle"))]) } @@ -1808,10 +1813,22 @@ simData <- function(nbAnimals=1,nbStates=2,dist, } else { for(j in 1:(parSize[[i]]-zeroInflation[[i]]-oneInflation[[i]])) genArgs[[i]][[j+1]] <- subPar[[i]][(j-1)*nbStates+Z[k]] + if(inputs$dist[[i]]=="crwrice"){ + if(k>1) genArgs[[i]][[j+2]] <- genData[[i]][k-1] # previous step + else genArgs[[i]][[j+2]] <- 0 + } else if(inputs$dist[[i]]=="crwvm"){ + if(k>1) genArgs[[i]][[j+2]] <- genData[["step"]][c(k-1,k)] # previous step + } } if(inputs$dist[[i]] %in% angledists){ + if(dist[[i]]=="crwvm"){ + if(k>1 && !all(genArgs[[i]][[j+2]]>0)) { + genArgs[[i]][[2]] <- genArgs[[i]][[3]] <- 12 # uniform turn angle if previous steps are zero + genArgs[[i]][[4]] <- NULL + } + } genData[[i]][k] <- do.call(Fun[[i]],genArgs[[i]]) if(genData[[i]][k] > pi) genData[[i]][k] <- genData[[i]][k]-2*pi if(genData[[i]][k] < -pi) genData[[i]][k] <- genData[[i]][k]+2*pi diff --git a/R/simHierData.R b/R/simHierData.R index f0e8fc0..3e69291 100644 --- a/R/simHierData.R +++ b/R/simHierData.R @@ -100,6 +100,8 @@ simHierData <- function(nbAnimals=1,hierStates,hierDist, rmvnorm3 <- rmvnorm3 rvm <- CircStats::rvm rwrpcauchy <- CircStats::rwrpcauchy + rcrwrice <- rcrwrice + rcrwvm <- rcrwvm trMatrix_rcpp <- trMatrix_rcpp ctPar <- ctPar if (requireNamespace("extraDistr", quietly = TRUE)){ @@ -414,7 +416,10 @@ simHierData <- function(nbAnimals=1,hierStates,hierDist, estAngleMean <- vector('list',length(distnames)) names(estAngleMean) <- distnames for(i in distnames){ - if(dist[[i]] %in% angledists) estAngleMean[[i]]<-TRUE + if(dist[[i]] %in% angledists) { + estAngleMean[[i]]<-TRUE + if(dist[[i]]=="crwvm") estAngleMean[[i]]<-FALSE + } else estAngleMean[[i]]<-FALSE } @@ -1509,6 +1514,12 @@ simHierData <- function(nbAnimals=1,hierStates,hierDist, } else { for(j in 1:(parSize[[i]]-zeroInflation[[i]]-oneInflation[[i]])) genArgs[[i]][[j+1]] <- subPar[[i]][(j-1)*nbStates+Z[k]] + if(inputs$dist[[i]]=="crwrice"){ + if(k>1) genArgs[[i]][[j+2]] <- genData[[i]][k-1] # previous step + else genArgs[[i]][[j+2]] <- 0 + } else if(inputs$dist[[i]]=="crwvm"){ + if(k>1) genArgs[[i]][[j+2]] <- genData[["step"]][c(k-1,k)] # previous step + } } if(inputs$dist[[i]] %in% angledists){ diff --git a/R/w2n.R b/R/w2n.R index be2224c..fee92a1 100644 --- a/R/w2n.R +++ b/R/w2n.R @@ -133,7 +133,7 @@ w2n <- function(wpar,bounds,parSize,nbStates,nbCovs,estAngleMean,circularAngleMe tmpwpar<-wpar[parindex[[i]]+1:parCount[[i]]] - if(estAngleMean[[i]] & Bndind[[i]] & !TMB){ + if(estAngleMean[[i]] & Bndind[[i]] & !TMB & dist[[i]]!="crwvm"){ bounds[[i]][,1] <- -Inf bounds[[i]][which(bounds[[i]][,2]!=1),2] <- Inf @@ -148,7 +148,7 @@ w2n <- function(wpar,bounds,parSize,nbStates,nbCovs,estAngleMean,circularAngleMe parlist[[i]] <- tryCatch(w2nDM(tmpwpar,bounds[[i]],fullDM[[i]],DMind[[i]],nbObs,circularAngleMean[[i]],consensus[[i]],nbStates,0,nc[[i]],meanind[[i]],workBounds[[i]],dist[[i]],TMB),error=function(e) e) if(inherits(parlist[[i]],"error")) stop("\nFailed to calculate '",i,"' natural scale parameters: ",parlist[[i]],"\n check initial values, workBounds, userBounds, nlmPar or control, etc.") - if((dist[[i]] %in% angledists) & !estAngleMean[[i]]){ + if((dist[[i]] %in% angledists) & !estAngleMean[[i]] & dist[[i]]!="crwvm"){ tmp<-matrix(0,nrow=(parSize[[i]]+1)*nbStates,ncol=nbObs) tmp[nbStates+1:nbStates,]<-parlist[[i]] parlist[[i]] <- tmp diff --git a/man/dcrwrice_rcpp.Rd b/man/dcrwrice_rcpp.Rd new file mode 100644 index 0000000..618768d --- /dev/null +++ b/man/dcrwrice_rcpp.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{dcrwrice_rcpp} +\alias{dcrwrice_rcpp} +\title{Correlated random walk Rice distribution} +\usage{ +dcrwrice_rcpp(x, beta, sigma) +} +\arguments{ +\item{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).} + +\item{beta}{correlation parameter} + +\item{sigma}{speed parameter} +} +\value{ +Vector of densities +} +\description{ +Probability density function of Rice distribution under the CTCRW model of Johnson et al. (2008) (written in C++) +} diff --git a/man/dcrwvm_rcpp.Rd b/man/dcrwvm_rcpp.Rd new file mode 100644 index 0000000..0b993a2 --- /dev/null +++ b/man/dcrwvm_rcpp.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{dcrwvm_rcpp} +\alias{dcrwvm_rcpp} +\title{Correlated random walk von Mises density function} +\usage{ +dcrwvm_rcpp(x, beta, sigma) +} +\arguments{ +\item{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.} + +\item{beta}{correlation parameter} + +\item{sigma}{speed parameter} +} +\value{ +Vector of densities +} +\description{ +Probability density function of the Von Mises distribution under the CTCRW model of Johnson et al. (2008) (written in C++) +} diff --git a/man/nLogLike.Rd b/man/nLogLike.Rd index c65c7f4..e92cf36 100644 --- a/man/nLogLike.Rd +++ b/man/nLogLike.Rd @@ -42,7 +42,8 @@ nLogLike( aInd = aInd, CT = FALSE, dtIndex = NULL, - kappa = Inf + kappa = Inf, + crwST = FALSE ) } \arguments{ @@ -127,6 +128,8 @@ model (if any).} \item{dtIndex}{time difference index for calculating transition probabilities of hierarchical continuous-time models} \item{kappa}{maximum allowed value for the row sums of the off-diagonal elements in the state transition rate matrix, such that the minimum value for the diagonal elements is \code{-kappa}. Default: \code{Inf}. Setting less than \code{Inf} can help avoid numerical issues during optimization, in which case the transition rate parameters \code{beta} are on the logit scale (instead of the log scale). Ignored unless \code{CT=TRUE}.} + +\item{crwST}{logical indicating whether to fit a correlated step and turn random walk model} } \value{ The negative log-likelihood of the parameters given the data. diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 82c483f..aebc492 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -294,6 +294,32 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// dcrwrice_rcpp +arma::colvec dcrwrice_rcpp(NumericVector x, arma::mat beta, arma::mat sigma); +RcppExport SEXP _momentuHMM_dcrwrice_rcpp(SEXP xSEXP, SEXP betaSEXP, SEXP sigmaSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< arma::mat >::type beta(betaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type sigma(sigmaSEXP); + rcpp_result_gen = Rcpp::wrap(dcrwrice_rcpp(x, beta, sigma)); + return rcpp_result_gen; +END_RCPP +} +// dcrwvm_rcpp +arma::colvec dcrwvm_rcpp(NumericVector x, arma::mat beta, arma::mat sigma); +RcppExport SEXP _momentuHMM_dcrwvm_rcpp(SEXP xSEXP, SEXP betaSEXP, SEXP sigmaSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< arma::mat >::type beta(betaSEXP); + Rcpp::traits::input_parameter< arma::mat >::type sigma(sigmaSEXP); + rcpp_result_gen = Rcpp::wrap(dcrwvm_rcpp(x, beta, sigma)); + return rcpp_result_gen; +END_RCPP +} // expmatrix_rcpp arma::mat expmatrix_rcpp(arma::mat x, double kappa, bool check); RcppExport SEXP _momentuHMM_expmatrix_rcpp(SEXP xSEXP, SEXP kappaSEXP, SEXP checkSEXP) { diff --git a/src/TMB/include/dist.hpp b/src/TMB/include/dist.hpp index b95754d..8164b6a 100644 --- a/src/TMB/include/dist.hpp +++ b/src/TMB/include/dist.hpp @@ -71,6 +71,10 @@ std::unique_ptr> dist_generator(const int& code) { return(std::unique_ptr>(new RWcovMultivariateNormal)); case 28: return(std::unique_ptr>(new RWMultivariateNormal)); + case 29: + return(std::unique_ptr>(new CRWRice)); + case 30: + return(std::unique_ptr>(new CRWVonMises)); default: return(std::unique_ptr>(new Normal)); } diff --git a/src/TMB/include/dist_def.hpp b/src/TMB/include/dist_def.hpp index be47a72..88919af 100644 --- a/src/TMB/include/dist_def.hpp +++ b/src/TMB/include/dist_def.hpp @@ -1426,4 +1426,82 @@ class Dirichlet : public Dist { } }; + +template +class CRWRice : public Dist { +public: + // Constructor + CRWRice() {}; + // Link function + vector link(const vector& par, const int& n_states) { + vector wpar(par.size()); + // mean and sd + wpar = log(par); + return(wpar); + } + // Inverse link function + matrix invlink(const vector& wpar, const int& n_states) { + int n_par = wpar.size()/n_states; + matrix par(n_states, n_par); + // mean + for (int i = 0; i < n_states; ++i) par(i, 0) = exp(wpar(i)); + // sd + for (int i = 0; i < n_states; ++i) par(i, 1) = exp(wpar(i + n_states)); + return(par); + } + // Probability density/mass function + Type pdf(const vector& x, const vector& par, const Type& delta_t, const bool& logpdf) { + Type beta = par(0); + Type sigma = par(1); + Type step_tm1 = x(0); + Type step = x(1); + Type mu = step_tm1*(Type(1.)-exp(-beta))/beta; + Type var = sigma*sigma/(beta*beta) * (Type(1.) - Type(2.) / beta * (Type(1.)-exp(-beta))+Type(1.) / (Type(2.)*beta)*(Type(1.)-exp(-Type(2.)*beta))); + Type sd = sqrt(var); + Type xabs = fabs(step*mu/var); + Type val = log(step) - Type(2.0) * log(sd) + + (-(step*step+mu*mu)/(Type(2.0)*var)) + + log(exp(-xabs)*besselI(xabs,Type(0))) + xabs; + if(!logpdf) val = exp(val); + return(val); + } +}; + +template +class CRWVonMises : public Dist { +public: + // Constructor + CRWVonMises() {}; + // Link function + vector link(const vector& par, const int& n_states) { + vector wpar(par.size()); + // mean and sd + wpar = log(par); + return(wpar); + } + // Inverse link function + matrix invlink(const vector& wpar, const int& n_states) { + int n_par = wpar.size()/n_states; + matrix par(n_states, n_par); + // mean + for (int i = 0; i < n_states; ++i) par(i, 0) = exp(wpar(i)); + // sd + for (int i = 0; i < n_states; ++i) par(i, 1) = exp(wpar(i + n_states)); + return(par); + } + // Probability density/mass function + Type pdf(const vector& x, const vector& par, const Type& delta_t, const bool& logpdf) { + Type beta = par(0); + Type sigma = par(1); + Type step_tm1 = x(0); + Type step = x(1); + Type angle = x(2); + Type var = sigma*sigma/(beta*beta) * (Type(1.) - Type(2.) / beta * (Type(1.)-exp(-beta))+Type(1.) / (Type(2.)*beta)*(Type(1.)-exp(-Type(2.)*beta))); + Type kappa = step_tm1*step*exp(-beta)/var; + Type b = exp(-kappa) * besselI(kappa, Type(0)); + Type val = -log(Type(2.) * M_PI * b) + kappa * (cos(angle)-Type(1.)); + if(!logpdf) val = exp(val); + return(val); + } +}; #endif diff --git a/src/combine.h b/src/combine.h index 01c10a9..7c3fe3d 100644 --- a/src/combine.h +++ b/src/combine.h @@ -72,4 +72,12 @@ arma::mat cbindsigma3(arma::mat x, arma::mat y, arma::mat z, arma::mat xy, arma: out(4,_) = as(wrap(xz)); out(5,_) = as(wrap(yz)); return as(out); -} \ No newline at end of file +} + +// // [[Rcpp::export]] +//arma::mat concat(arma::mat x, arma::mat y){ +// NumericMatrix out(1,x.size() + y.size()); +// std::copy(x.begin(),x.end(), out.begin()); +// std::copy(y.begin(),y.end(), out.begin() + x.size()); +// return as(out); +//} \ No newline at end of file diff --git a/src/densities.h b/src/densities.h index eda32a7..1710f47 100644 --- a/src/densities.h +++ b/src/densities.h @@ -521,6 +521,86 @@ arma::colvec dt_rcpp(NumericVector x, arma::mat df, arma::mat ncp) return res; } +//' 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 +// [[Rcpp::export]] +arma::colvec dcrwrice_rcpp(NumericVector x, arma::mat beta, arma::mat sigma) +{ + double mu; + double var; + double sd; + double xabs; + double step_tm1; + unsigned int xs = (unsigned int) beta.n_cols; + arma::colvec res(xs); + res(0) = 0.; + + for(unsigned int i=0;i #include "expm.h" -// Mostly generated with tools::package_native_routine_registration_skeleton (R 4.0) +// Mostly generated with tools::package_native_routine_registration_skeleton(getwd(),,,FALSE) (R 4.0) // additionally requires: // #include "expm.h" in header // expm = (void (*) (double*, int, double*, precond_type)) R_GetCCallable("expm", "expm"); in R_init_momentuHMM @@ -18,6 +18,8 @@ extern SEXP _momentuHMM_combine(SEXP, SEXP); extern SEXP _momentuHMM_dbern_rcpp(SEXP, SEXP, SEXP); extern SEXP _momentuHMM_dbeta_rcpp(SEXP, SEXP, SEXP); extern SEXP _momentuHMM_dcat_rcpp(SEXP, SEXP, SEXP); +extern SEXP _momentuHMM_dcrwrice_rcpp(SEXP, SEXP, SEXP); +extern SEXP _momentuHMM_dcrwvm_rcpp(SEXP, SEXP, SEXP); extern SEXP _momentuHMM_dexp_rcpp(SEXP, SEXP, SEXP); extern SEXP _momentuHMM_dgamma_rcpp(SEXP, SEXP, SEXP); extern SEXP _momentuHMM_dlnorm_rcpp(SEXP, SEXP, SEXP); @@ -46,6 +48,8 @@ static const R_CallMethodDef CallEntries[] = { {"_momentuHMM_dbern_rcpp", (DL_FUNC) &_momentuHMM_dbern_rcpp, 3}, {"_momentuHMM_dbeta_rcpp", (DL_FUNC) &_momentuHMM_dbeta_rcpp, 3}, {"_momentuHMM_dcat_rcpp", (DL_FUNC) &_momentuHMM_dcat_rcpp, 3}, + {"_momentuHMM_dcrwrice_rcpp", (DL_FUNC) &_momentuHMM_dcrwrice_rcpp, 3}, + {"_momentuHMM_dcrwvm_rcpp", (DL_FUNC) &_momentuHMM_dcrwvm_rcpp, 3}, {"_momentuHMM_dexp_rcpp", (DL_FUNC) &_momentuHMM_dexp_rcpp, 3}, {"_momentuHMM_dgamma_rcpp", (DL_FUNC) &_momentuHMM_dgamma_rcpp, 3}, {"_momentuHMM_dlnorm_rcpp", (DL_FUNC) &_momentuHMM_dlnorm_rcpp, 3}, diff --git a/src/nLogLike.cpp b/src/nLogLike.cpp index f546d26..90dbc16 100644 --- a/src/nLogLike.cpp +++ b/src/nLogLike.cpp @@ -216,6 +216,8 @@ double nLogLike_rcpp(int nbStates, arma::mat covs, DataFrame data, CharacterVect funMap["bern"] = dbern_rcpp; funMap["beta"] = dbeta_rcpp; funMap["cat"] = dcat_rcpp; + funMap["crwrice"] = dcrwrice_rcpp; + funMap["crwvm"] = dcrwvm_rcpp; funMap["ctds"] = dcat_rcpp; funMap["exp"] = dexp_rcpp; funMap["gamma"] = dgamma_rcpp; @@ -259,8 +261,8 @@ double nLogLike_rcpp(int nbStates, arma::mat covs, DataFrame data, CharacterVect unsigned int nDists = (unsigned int) dist.size(); - for(unsigned int k=0;k(dataNames[k]); + for(unsigned int d=0;d(dataNames[d]); genDist = as(dist[genname]); if(genDist=="mvnorm2" || genDist=="rw_mvnorm2"){ L = List::create(as(data[genname+".x"]) ,as(data[genname+".y"])); @@ -270,6 +272,12 @@ double nLogLike_rcpp(int nbStates, arma::mat covs, DataFrame data, CharacterVect L = List::create(as(data[genname+".x"]) ,as(data[genname+".y"]),as(data[genname+".z"])); //NumericVector genData = combine(L); mvn = true; + } else if(genDist=="crwrice"){ + L = List::create(as(data["step"]),as(data["step"])); + mvn = false; + } else if(genDist=="crwvm"){ + L = List::create(as(data["angle"]),as(data["step"]),as(data["step"])); + mvn = false; } else { L = List::create(as(data[genname])); mvn = false; @@ -289,18 +297,63 @@ double nLogLike_rcpp(int nbStates, arma::mat covs, DataFrame data, CharacterVect else zeroInd = 0; + unsigned int k=0; // animal index + // remove the NAs from step (impossible to subset a vector with NAs) genData = combine(L,NAvalue); arma::uvec noNAs = arma::find_finite(as(L[0])); - + arma::colvec tmp(nbObs); + if(genDist=="crwrice") { + for(int i=0; i < nbObs; i++){ + if(k0.) genData[nbObs+i+1] = genData[i]; + else if(genData[i]==0){ + genData[nbObs+i+1] = DBL_MIN; + genData[nbObs+i] = NAvalue; + } + } + } + std::copy(genData.begin(),genData.begin()+nbObs,tmp.begin()); + noNAs = arma::find(tmp>=0.); + L[0] = tmp; + //k=0; + //for(int i=0; i < nbObs; i++){ + // Rprintf("k %d i %d step %f step_tm1 %f tmp[i] %f \n",k,i,genData[i],genData[nbObs+i],tmp[i]); + // if((k+1(L[0])); + for(int i : noNAs){ + if(i>0) { + genData[2*nbObs+i] = genData[nbObs+i-1]; + } + } + for(int i : NAs){ + genData[nbObs+i] = NAvalue; + genData[2*nbObs+i] = NAvalue; + } + } + // extract zero-mass and one-mass parameters if necessary if(genzeroInflation || genoneInflation) { if(genzeroInflation){ zeromass = genPar.rows(genPar.n_rows-oneInd-nbStates,genPar.n_rows-oneInd-1); //genPar(arma::span(genPar.n_rows-1),arma::span(),arma::span()); - noZeros = arma::find(as(genData)>0); - nbZeros = arma::find(as(genData)==0); + if(genDist=="crwrice"){ + noZeros = arma::find(as(L[0])>0); + nbZeros = arma::find(as(L[0])==0); + } else { + noZeros = arma::find(as(genData)>0); + nbZeros = arma::find(as(genData)==0); + } } if(genoneInflation){ @@ -355,7 +408,7 @@ double nLogLike_rcpp(int nbStates, arma::mat covs, DataFrame data, CharacterVect zerom = zeromass.row(state); // compute probability of non-zero observations - genProb.elem(noZeros) = (1. - zerom.elem(noZeros)) % funMap[genDist](genData[genData>0.],genArgs1.elem(noZeros),genArgs2.elem(noZeros)); + genProb.elem(noZeros) = (1. - zerom.elem(noZeros)) % funMap[genDist](genData[genData>0.],genArgs1.cols(noZeros),genArgs2.cols(noZeros)); // compute probability of zero observations genProb.elem(nbZeros) = zerom.elem(nbZeros); @@ -365,7 +418,7 @@ double nLogLike_rcpp(int nbStates, arma::mat covs, DataFrame data, CharacterVect onem = onemass.row(state); // compute probability of non-one observations - genProb.elem(noOnes) = (1. - onem.elem(noOnes)) % funMap[genDist](genData[(genData!=NAvalue) & (genData<1.)],genArgs1.elem(noOnes),genArgs2.elem(noOnes)); + genProb.elem(noOnes) = (1. - onem.elem(noOnes)) % funMap[genDist](genData[(genData!=NAvalue) & (genData<1.)],genArgs1.cols(noOnes),genArgs2.cols(noOnes)); // compute probability of one observations genProb.elem(nbOnes) = onem.elem(nbOnes); @@ -376,7 +429,7 @@ double nLogLike_rcpp(int nbStates, arma::mat covs, DataFrame data, CharacterVect onem = onemass.row(state); // compute probability of non-zero and non-one observations - genProb.elem(noZerosOnes) = (1. - zerom.elem(noZerosOnes) - onem.elem(noZerosOnes)) % funMap[genDist](genData[(genData>0) & (genData<1)],genArgs1.elem(noZerosOnes),genArgs2.elem(noZerosOnes)); + genProb.elem(noZerosOnes) = (1. - zerom.elem(noZerosOnes) - onem.elem(noZerosOnes)) % funMap[genDist](genData[(genData>0) & (genData<1)],genArgs1.cols(noZerosOnes),genArgs2.cols(noZerosOnes)); // compute probability of zero observations genProb.elem(nbZeros) = zerom.elem(nbZeros); @@ -385,6 +438,11 @@ double nLogLike_rcpp(int nbStates, arma::mat covs, DataFrame data, CharacterVect genProb.elem(nbOnes) = onem.elem(nbOnes); } else { + //k=0; + //for(unsigned int i : noNAs){ + // Rprintf("k %d i %d step %f step_tm1 %f beta %f sigma %f tmp[i] %f \n",k,i,genData[i],genData[nbObs+i],genArgs1(i),genArgs2(i),tmp[i]); + // if((k+1