Added drift models and namespace corrections

DESCRIPTION

@@ -1,18 +1,18 @@
 Package: crawl
 Type: Package
 Title: Fit Continuous-Time Correlated Random Walk Models to Animal Movement Data
-Version: 1.913
-Date: 2015-08-20
+Version: 1.920
+Date: 2015-08-21
 Author: Devin S. Johnson
 Maintainer: Devin S. Johnson <devin.johnson@noaa.gov>
 Imports:
     mvtnorm,
     sp,
     raster,
     Rcpp (>= 0.11.1),
-    RcppArmadillo (>= 0.4.200.0),
     methods
 LinkingTo: Rcpp, RcppArmadillo
+Suggests: ggplot2
 Description: The Correlated RAndom Walk Library of
     R functions was designed for fitting continuous-time
     correlated random walk (CTCRW) models with time indexed

NAMESPACE

@@ -19,7 +19,6 @@ export(flatten)
 export(intToPOSIX)
 export(mergeTrackStop)
 import(Rcpp)
-import(RcppArmadillo)
 import(mvtnorm)
 import(raster)
 import(sp)

R/RcppExports.R

@@ -1,23 +1,43 @@
 # This file was generated by Rcpp::compileAttributes
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
-makeT <- function(b, delta) {
-    .Call('crawl_makeT', PACKAGE = 'crawl', b, delta)
+CTCRWNLL_DRIFT <- function(y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P) {
+    .Call('crawl_CTCRWNLL_DRIFT', PACKAGE = 'crawl', y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P)
 }
 
-makeQ <- function(b, sig2, delta) {
-    .Call('crawl_makeQ', PACKAGE = 'crawl', b, sig2, delta)
+CTCRWNLL <- function(y, Hmat, beta, sig2, delta, noObs, active, a, P) {
+    .Call('crawl_CTCRWNLL', PACKAGE = 'crawl', y, Hmat, beta, sig2, delta, noObs, active, a, P)
 }
 
-CTCRWNLL <- function(y, Hmat, Qmat, Tmat, noObs, active, a, P) {
-    .Call('crawl_CTCRWNLL', PACKAGE = 'crawl', y, Hmat, Qmat, Tmat, noObs, active, a, P)
+CTCRWPREDICT_DRIFT <- function(y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P) {
+    .Call('crawl_CTCRWPREDICT_DRIFT', PACKAGE = 'crawl', y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P)
 }
 
-CTCRWPREDICT <- function(y, Hmat, Qmat, Tmat, noObs, active, a, P) {
-    .Call('crawl_CTCRWPREDICT', PACKAGE = 'crawl', y, Hmat, Qmat, Tmat, noObs, active, a, P)
+CTCRWPREDICT <- function(y, Hmat, beta, sig2, delta, noObs, active, a, P) {
+    .Call('crawl_CTCRWPREDICT', PACKAGE = 'crawl', y, Hmat, beta, sig2, delta, noObs, active, a, P)
 }
 
-CTCRWSAMPLE <- function(y, Hmat, Qmat, Tmat, noObs, active, a, P) {
-    .Call('crawl_CTCRWSAMPLE', PACKAGE = 'crawl', y, Hmat, Qmat, Tmat, noObs, active, a, P)
+CTCRWSAMPLE_DRIFT <- function(y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P) {
+    .Call('crawl_CTCRWSAMPLE_DRIFT', PACKAGE = 'crawl', y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P)
+}
+
+CTCRWSAMPLE <- function(y, Hmat, beta, sig2, delta, noObs, active, a, P) {
+    .Call('crawl_CTCRWSAMPLE', PACKAGE = 'crawl', y, Hmat, beta, sig2, delta, noObs, active, a, P)
+}
+
+makeT <- function(b, delta, active) {
+    .Call('crawl_makeT', PACKAGE = 'crawl', b, delta, active)
+}
+
+makeQ <- function(b, sig2, delta, active) {
+    .Call('crawl_makeQ', PACKAGE = 'crawl', b, sig2, delta, active)
+}
+
+makeT_drift <- function(b, b_drift, delta, active) {
+    .Call('crawl_makeT_drift', PACKAGE = 'crawl', b, b_drift, delta, active)
+}
+
+makeQ_drift <- function(b, b_drift, sig2, sig2_drift, delta, active) {
+    .Call('crawl_makeQ_drift', PACKAGE = 'crawl', b, b_drift, sig2, sig2_drift, delta, active)
 }
 

R/crawl-package.R

@@ -10,8 +10,8 @@
 #' \tabular{ll}{ 
 #' Package: \tab crawl\cr 
 #' Type: \tab Package\cr 
-#' Version: \tab 1.913\cr 
-#' Date: \tab August 20, 2015\cr 
+#' Version: \tab 1.920\cr 
+#' Date: \tab August 21, 2015\cr 
 #' License: \tab file LICENSE \cr 
 #' LazyLoad: \tab yes\cr 
 #' }
@@ -25,7 +25,7 @@
 #' @references Johnson, D., J. London, M. -A. Lea, and J. Durban (2008)
 #' Continuous-time correlated random walk model for animal telemetry data.
 #' Ecology 89(5) 1208-1215.
-#' @import Rcpp RcppArmadillo
+#' @import Rcpp
 #' @importFrom graphics layout
 #' @importFrom methods as slot 
 #' @importFrom stats approx filter model.frame model.matrix 

R/crwMLE.R

@@ -169,7 +169,7 @@ crwMLE = function(mov.model=~1, err.model=NULL, activity=NULL, drift=FALSE,
                      initial.state, theta, fixPar, method="L-BFGS-B", control=NULL, constr=list(lower=-Inf, upper=Inf), 
                      prior=NULL, need.hess=TRUE, initialSANN=list(maxit=200), attempts=1)
 {
-  if(drift) stop("At this time drift models are not supported with this function. Use 'crwMLE' for now.\n")
+  #if(drift) stop("At this time drift models are not supported with this function. Use 'crwMLE' for now.\n")
   st <- Sys.time()
 #  if (missing(Time.name) & !inherits(data,"STIDF")) stop("Argument 'Time.name' missing and NOT a spacetime STIDF object. Please specify")
 

R/crwN2ll.R

@@ -71,21 +71,17 @@ crwN2ll = function(theta, fixPar, y, noObs, delta, a,
     active <- ifelse(b==Inf, 0, 1)
     b <- ifelse(b==Inf, 0, b) 
   } else {active=rep(1,N)}
-  #if (driftMod) {
-  ### Change back for drift model
-  if(as.logical(FALSE)){
+  if (driftMod) {
     theta.drift <- par[(n.errX + n.errY + 2 * n.mov + 1):
                          (n.errX + n.errY + 2 * n.mov + 2)]
     b.drift <- exp(log(b) - log(1+exp(theta.drift[2])))
-    sig2.drift <- exp(log(sig2) + 2 * theta.drift[1]) #rep(exp(2*theta.drift[1]), length(sig2)) #
-    call.lik <- "crwdriftn2ll"
+    sig2.drift <- exp(log(sig2) + 2 * theta.drift[1]) 
+    ll <- CTCRWNLL_DRIFT( y=as.matrix(y), Hmat, b, b.drift, sig2, sig2.drift, delta, noObs, active, a,  P)$ll
   } else {
-    b.drift <- sig2.drift <- 0.0
-    call.lik <- CTCRWNLL
+    ll <- CTCRWNLL( y=as.matrix(y), Hmat, b, sig2, delta, noObs, active, a,  P)$ll
   }
-  movMats <- getQT(sig2, b, sig2.drift, b.drift, delta, driftMod=FALSE)
-
-  ll <- CTCRWNLL( y=as.matrix(y), Hmat, movMats[["Qmat"]], movMats[["Tmat"]], noObs, active, a,  P)$ll
+  #movMats <- getQT(sig2, b, sig2.drift, b.drift, delta, driftMod=FALSE)
+
   if(is.null(prior)) return(-2 * ll)
   else return(-2 * (ll + prior(theta)))
 }

R/crwPostIS.R

@@ -121,22 +121,16 @@
     active <- ifelse(b==Inf, 0, 1)
     b <- ifelse(b==Inf, 0, b) 
   } else active = rep(1,N)
-  #if (driftMod) {
-  ### Change back for drift model
-  if(as.logical(FALSE)){
+  if (driftMod) {
     theta.drift <- par[(n.errX + n.errY + 2 * n.mov + 1):
                          (n.errX + n.errY + 2 * n.mov + 2)]
     b.drift <- exp(log(b) - log(1+exp(theta.drift[2])))
-    sig2.drift <- exp(log(sig2) + 2 * theta.drift[1]) #rep(exp(2*theta.drift[1]), length(sig2)) #
-    call.lik <- "crwdriftn2ll"
+    sig2.drift <- exp(log(sig2) + 2 * theta.drift[1]) 
+    out=CTCRWSAMPLE_DRIFT(y, Hmat, b, b.drift, sig2, sig2.drift, delta, noObs, active, a, P)
   } else {
-    b.drift <- sig2.drift <- 0.0
-    #call.lik <- CTCRWNLL
+    out=CTCRWSAMPLE(y, Hmat, b, sig2, delta, noObs, active, a, P)
   }
-  movMats <- getQT(sig2, b, sig2.drift, b.drift, delta, driftMod)
 
-  out=CTCRWSAMPLE(y, Hmat, movMats$Qmat, movMats$Tmat, noObs, active, a, P)
-
   if(driftMod){
     colnames(out$sim) <- apply(expand.grid(c("mu","theta","gamma"), c("x","y")), 1, paste, collapse=".")
   }  else {

R/crwPredict.R

@@ -115,29 +115,21 @@ crwPredict=function(object.crwFit, predTime=NULL, flat=TRUE, ...)
   theta.mov <- par[(n.errX + n.errY + 1):(n.errX + n.errY + 2 * n.mov)]
   sig2 <- exp(2 * (mov.mf %*% theta.mov[1:n.mov]))
   b <- exp(mov.mf %*% theta.mov[(n.mov + 1):(2 * n.mov)])
-  #activity <- rep(1, N)
   if (!is.null(activity)) {
     theta.activ <- par[(n.errX + n.errY + 2 * n.mov + 1)]
     b <- b / ((activity) ^ exp(theta.activ))
     active <- ifelse(b==Inf, 0, 1)
     b <- ifelse(b==Inf, 0, b) 
   } else active = rep(1,N)
-  #if (driftMod) {
-  ### Change back for drift model
-  if(as.logical(FALSE)){
+  if (driftMod) {
     theta.drift <- par[(n.errX + n.errY + 2 * n.mov + 1):
                          (n.errX + n.errY + 2 * n.mov + 2)]
     b.drift <- exp(log(b) - log(1+exp(theta.drift[2])))
-    sig2.drift <- exp(log(sig2) + 2 * theta.drift[1]) #rep(exp(2*theta.drift[1]), length(sig2)) #
-    call.lik <- "crwdriftn2ll"
+    sig2.drift <- exp(log(sig2) + 2 * theta.drift[1])
+    out = CTCRWPREDICT_DRIFT(y, Hmat, b, b.drift, sig2, sig2.drift, delta, noObs, active, a, P)
   } else {
-    b.drift <- sig2.drift <- 0.0
-    #call.lik <- CTCRWNLL
+    out=CTCRWPREDICT(y, Hmat, b, sig2, delta, noObs, active, a, P)
   }
-  movMats <- getQT(sig2, b, sig2.drift, b.drift, delta, driftMod)
-
-  out=CTCRWPREDICT(y, Hmat, movMats$Qmat, movMats$Tmat, noObs, active, a, P) 
-
 
   pred <- data.frame(t(out$pred))
   if (driftMod) {
@@ -149,22 +141,23 @@ crwPredict=function(object.crwFit, predTime=NULL, flat=TRUE, ...)
   obsFit$outlier.chisq <- out$chisq
   obsFit$naive.p.val <- 1 - pchisq(obsFit$outlier.chisq, 2)
   if(getUseAvail){
-    warning("'getUseAvail' not implemented yet in this version of 'crawl' contact Devin to fix this! ")
-#     idx <- data$locType=="p"
-#     movMatsPred <- getQT(sig2[idx], b[idx], sig2.drift[idx], b.drift[idx], delta=c(diff(data[idx,tn]),1), driftMod)
-#     TmatP <- movMatsPred$Tmat
-#     QmatP <- movMatsPred$Qmat
-#     avail <- t(sapply(1:(nrow(TmatP)-1), makeAvail, Tmat=TmatP, Qmat=QmatP, predx=predx[idx,], predy=predy[idx,], 
-#                       vary=vary[,,idx], varx=varx[,,idx], driftMod=driftMod, lonadj=lonAdjVals[idx]))
-#     avail <- cbind(data[idx,tn][-1], avail)
-#     colnames(avail) <- c(tn, "meanAvail.x", "meanAvail.y", "varAvail.x", "varAvail.y")
-#     use <- cbind(data[idx,tn], predx[idx,1], predy[idx,1], varx[1,1,idx], vary[1,1,idx])[-1,]
-#     colnames(use) <- c(tn, "meanUse.x", "meanUse.y", "varUse.x", "varUse.y")
-#     UseAvail.lst <- list(use=use, avail=avail)
+    warning("'getUseAvail' not implemented yet in this version of 'crawl' contact maintainer to fix this! ")
+    #     idx <- data$locType=="p"
+    #     movMatsPred <- getQT(sig2[idx], b[idx], sig2.drift[idx], b.drift[idx], delta=c(diff(data[idx,tn]),1), driftMod)
+    #     TmatP <- movMatsPred$Tmat
+    #     QmatP <- movMatsPred$Qmat
+    #     avail <- t(sapply(1:(nrow(TmatP)-1), makeAvail, Tmat=TmatP, Qmat=QmatP, predx=predx[idx,], predy=predy[idx,], 
+    #                       vary=vary[,,idx], varx=varx[,,idx], driftMod=driftMod, lonadj=lonAdjVals[idx]))
+    #     avail <- cbind(data[idx,tn][-1], avail)
+    #     colnames(avail) <- c(tn, "meanAvail.x", "meanAvail.y", "varAvail.x", "varAvail.y")
+    #     use <- cbind(data[idx,tn], predx[idx,1], predy[idx,1], varx[1,1,idx], vary[1,1,idx])[-1,]
+    #     colnames(use) <- c(tn, "meanUse.x", "meanUse.y", "varUse.x", "varUse.y")
+    #     UseAvail.lst <- list(use=use, avail=avail)
     UseAvail.lst=NULL
   }
   else UseAvail.lst=NULL
-  speed = sqrt(apply(matrix(pred[,2:(2+driftMod)]), 1, sum)^2 + apply(matrix(pred[,(4+driftMod):(4+2*driftMod)]), 1, sum)^2)
+  speed = sqrt(apply(as.matrix(pred[,2:(2+driftMod)]), 1, sum)^2 + 
+                 apply(as.matrix(pred[,(4+driftMod):(4+2*driftMod)]), 1, sum)^2)
   out <- list(originalData=fillCols(data), alpha.hat=pred, 
               V.hat=var, speed=speed, loglik=out$ll, useAvail=UseAvail.lst)
   if (flat) {

demo/harborSeal.R

@@ -7,17 +7,18 @@ harborSeal$Argos_loc_class = factor(harborSeal$Argos_loc_class, levels=c("3","2"
 ## Project data ##
 library(rgdal)
 
-toProj = harborSeal[!is.na(harborSeal$latitude),]
+toProj = harborSeal[!is.na(harborSeal$latitude),c("Time","latitude","longitude")]
 coordinates(toProj) = ~longitude+latitude
 proj4string(toProj) <- CRS("+proj=longlat")
 toProj <- spTransform(toProj, CRS("+init=epsg:3338"))
 toProj = as.data.frame(toProj)
-harborSeal = merge(toProj, harborSeal, all=TRUE)
+colnames(toProj)[2:3] = c("x","y")
+harborSeal = merge(toProj, harborSeal, by="Time", all=TRUE)
 harborSeal = harborSeal[order(harborSeal$Time),]
 
-initial.cpp = list(
+initial = list(
   a=c(harborSeal$x[1],0,harborSeal$y[1],0),
-  P=diag(c(10000^2,54000^2,10000^2,5400^2))
+  P=diag(c(10000^2,5400^2,10000^2,5400^2))
 )
 
 ##Fit model as given in Johnson et al. (2008) Ecology 89:1208-1215
@@ -34,10 +35,10 @@ set.seed(123)
 fit1 <- crwMLE(
   mov.model=~1, err.model=list(x=~Argos_loc_class-1), activity=~I(1-DryTime),
   data=harborSeal, coord=c("x","y"), Time.name="Time", 
-  initial.state=initial.cpp, fixPar=fixPar, theta=c(rep(log(5000),3),log(3*3600), 3),
+  initial.state=initial, fixPar=fixPar, theta=c(rep(log(5000),3),log(3*3600), 0),
   constr=constr,
   control=list(maxit=2000, trace=1, REPORT=1),
-  initialSANN=list(maxit=200, trace=1, temp=20, tmax=20, REPORT=1)
+  initialSANN=list(maxit=200, trace=1, REPORT=1)
 )
 
 print(fit1)
@@ -47,9 +48,9 @@ p2=ggplot(aes(x=Time, y=mu.x), data=pred1) + geom_ribbon(aes(ymin=mu.x-2*se.mu.x
   geom_path(, col="red") + geom_point(aes(x=Time, y=x), col="blue", size=1)
 p3=ggplot(aes(x=Time, y=mu.y), data=pred1) + geom_ribbon(aes(ymin=mu.y-2*se.mu.y, ymax=mu.y+2*se.mu.y), fill="green", alpha=0.5)  + 
   geom_path(, col="red") + geom_point(aes(x=Time, y=y), col="blue", size=1)
-print(p1)
-print(p2)
-print(p3)
+suppressWarnings(print(p1))
+suppressWarnings(print(p2))
+suppressWarnings(print(p3))
 # ggsave("map.pdf", p1)
 # ggsave("xaxis.pdf", p2, width=10, height=2)
 # ggsave("yaxis.pdf", p3, width=10, height=2)

demo/harborSealSmooth.R

@@ -1,23 +1,23 @@
-require(crawl)
-require(splines)
+# library(crawl)
+library(ggplot2)
+library(splines)
+library(rgdal)
 data(harborSeal)
 head(harborSeal)
 harborSeal$Argos_loc_class = factor(harborSeal$Argos_loc_class, levels=c("3","2","1","0","A","B"))
 
 ## Project data ##
-library(sp)
-library(rgdal)
 
-toProj = harborSeal[!is.na(harborSeal$latitude),]
+toProj = harborSeal[!is.na(harborSeal$latitude),c("Time","latitude","longitude")]
 coordinates(toProj) = ~longitude+latitude
 proj4string(toProj) <- CRS("+proj=longlat")
 toProj <- spTransform(toProj, CRS("+init=epsg:3338"))
 toProj = as.data.frame(toProj)
-harborSeal = merge(toProj, harborSeal, all=TRUE)
+colnames(toProj)[2:3] = c("x","y")
+harborSeal = merge(toProj, harborSeal, by="Time", all=TRUE)
 harborSeal = harborSeal[order(harborSeal$Time),]
-#harborSeal$Time = harborSeal$Time*3600
 
-initial.cpp = list(
+initial = list(
   a=c(harborSeal$x[1],0,harborSeal$y[1],0),
   P=diag(c(10000^2,54000^2,10000^2,5400^2))
 )
@@ -43,8 +43,8 @@ fit1 <- crwMLE(
   constr=constr,
   theta = theta.start,
   prior=prior,
-  control=list(maxit=2000, trace=1, REPORT=1),
-  initialSANN=list(maxit=10000, temp=100, tmax=100, trace=1, REPORT=1)
+  control=list(maxit=2000, trace=1, REPORT=1)#,
+  #initialSANN=list(maxit=10000, temp=100, tmax=100, trace=1, REPORT=1)
 )
 
 print(fit1)

demo/northernFurSeal.R

@@ -10,34 +10,34 @@ proj4string(northernFurSeal) <- CRS("+proj=longlat")
 northernFurSeal <- spTransform(northernFurSeal, CRS("+proj=aea +lat_1=30 +lat_2=70 +lat_0=52 +lon_0=-170 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"))
 
 initial = list(
-  a=c(coordinates(northernFurSeal)[1,1],0,coordinates(northernFurSeal)[1,2],0),
-  P=diag(c(10000^2,54000^2,10000^2,5400^2))
+  a=c(coordinates(northernFurSeal)[1,1],0,0,coordinates(northernFurSeal)[1,2],0,0),
+  P=diag(c(10000^2,5400^2,5400^2,10000^2,5400^2,5400^2))
 )
 
-fixPar = c(log(250), log(500), log(1500), rep(NA,3), NA)
-displayPar( mov.model=~1, err.model=list(x=~Argos_loc_class-1),
+fixPar = c(log(250), log(500), log(1500), rep(NA,6))
+displayPar( mov.model=~1, err.model=list(x=~Argos_loc_class-1), drift = TRUE,
             data=northernFurSeal,fixPar=fixPar)
 constr=list(
-  lower=c(rep(log(1500),2), rep(-Inf,2)),
-  upper=rep(Inf,4)
+  lower=c(rep(log(1500),2), rep(-Inf,4)),
+  upper=rep(Inf,6)
 )
 
-ln.prior = function(theta){-abs(theta[4]-3)/0.5}
+ln.prior = function(theta){-abs(theta[4]+3)/0.5}
 
 #set.seed(123)
 fit1 <- crwMLE(
-  mov.model=~1, err.model=list(x=~Argos_loc_class-1), 
+  mov.model=~1, err.model=list(x=~Argos_loc_class-1), drift=TRUE,
   data=northernFurSeal, Time.name="Time", 
-  initial.state=initial, fixPar=fixPar, constr=constr, prior=ln.prior,
-  control=list(maxit=30, trace=1, REPORT=1),
-  initialSANN=list(maxit=200, trace=1, REPORT=1)
+  initial.state=initial, fixPar=fixPar, constr=constr,# prior=ln.prior,
+  control=list(maxit=300, trace=1, REPORT=1),
+  initialSANN=list(maxit=200, temp=100, tmax=20, trace=1, REPORT=1)
 )
 fit1
 ##Make hourly location predictions
 predTime <- seq(ceiling(min(northernFurSeal$Time)), floor(max(northernFurSeal$Time)), 1)
 predObj <- crwPredict(object.crwFit=fit1, predTime, speedEst=TRUE, flat=TRUE)
 head(predObj)
-crwPredictPlot(predObj, "map")
+crwPredictPlot(predObj, "map", asp=TRUE)
 
 ##Create simulation object with 100 parameter draws
 set.seed(123)
@@ -54,13 +54,13 @@ round(100/(1+(sd(w)/mean(w))^2))
 jet.colors <-
   colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                      "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
-crwPredictPlot(predObj, 'map')
+crwPredictPlot(predObj, 'map', asp=TRUE)
 
 ## Sample 20 tracks from posterior predictive distribution
 iter <- 20
 cols <- jet.colors(iter)
 for(i in 1:iter){
-  samp <- crwPostIS(simObj)
+  samp <- crwPostIS(simObj, fullPost = FALSE)
   lines(samp$alpha.sim[,'mu.x'], samp$alpha.sim[,'mu.y'],col=cols[i])
 }