R/fit_rw.r

Defines functions fit_rw

Documented in fit_rw

##' Random Walk with time-varying autocorrelation
##'
##' Fit a random walk with time-varying autocorelation, via TMB,
##' to ssm-filtered track data.
##'
##' The input track is given as a data_frame where each row is a
##' location 'observed' without error and columns
##' \describe{
##' \item{'date'}{observation time (POSIXct, GMT),}
##' \item{'lon'}{observed longitude,}
##' \item{'lat'}{observed latitude,}
##' }
##'
##' @title Random Walk with time-varying autocorrelation
##' @param data a dataframe of observations (see details)
##' @param formula regression formula for the covariates; default: \code{~ 1} (no covariates)
##' @param tstep time step as a proportion of 24 h
##' @param nu fixed df for t-distributed measurement errors
##' @param g0 initial value for stationary gamma (move autocorrelation) on 0, 1
##' @param span degree of loess smoothing for location initial values
##' @param optim numerical optimizer
##' @param verbose report progress during minimization
##' @param amfCRAWL Argos error class multiplication factors
##' @return a list with components
##' \item{\code{fitted}}{a dataframe of estimated states: gamma (tva model); lon, lat, gamma (tvae model)}
##' \item{\code{par}}{model parameter summmary}
##' \item{\code{data}}{input data}
##' \item{\code{tmb}}{the TMB object}
##' \item{\code{opt}}{the object returned by the optimizer}
##' \item{\code{aic}}{AIC statistic: 2 * NLL + 2 * k}
##' @useDynLib rwamr
##' @importFrom TMB MakeADFun sdreport
##' @importFrom tibble data_frame
##' @importFrom stringr str_detect
##' @export
fit_rw <-
  function(data,
           formula = ~ 1,
           tstep = 24 / 24,
           nu = 10,
           g0 = 0.5,
           span = 0.4,
           optim = c("nlminb", "optim"),
           verbose = FALSE,
           amf = amfCRAWL()) {
    optim <- match.arg(optim)
    n <- names(data)

    ## specify model based on whether data contain Argos location class column
    ##  will need to modify to allow other track data types with errors
    if ("lc" %in% n || "lq" %in% n)
      model = "tvae"
    else {
      model = "tva"
    }
    ## check that lon, lat & date are all columns in the data
    if (!"lon" %in% n ||
        !"lat" %in% n)
      stop("\n'lon' or 'lat' columns are missing from data")
    if (!"date" %in% n)
      stop("\n'date' column is missing from data")

    # check that the formula is a formula
    is.formula <- function(x)
      tryCatch(
        inherits(x, "formula"),
        error = function(e) {
          FALSE
        }
      )

    if (!is.formula(formula))
      stop("\n'formula' must be specified as ~ x + ... or ~ 1 for no covariates")

    # check that there is no response variable in the formula
    if (attr(terms(formula), "response") != 0)
      stop("\n'formula' can not have a response variable")

    d <- data

    ## Ensure POSIXct dates
    d$date <-
      as.POSIXct(d$date, "%Y-%m-%d %H:%M:%S", tz = "GMT")
    d <- d[order(d$date), ]

    mm <- model.matrix(formula, d)

    if(model == "tva" && dim(mm)[2] > 1) model <- "tva_cov"
    else if(model == "tvae" && dim(mm)[2] > 1) {
      stop("\nCovariate models are not yet supported for location data with error\n")
      }

    switch(model,
           tva = {
             ## TMB - data and parameter lists
             data <- with(d,
                          list(
                            model = 0,
                            x = cbind(lon, lat),
                            dt = as.numeric(difftime(date[-1], date[-nrow(d)],
                                                     unit = "hours")) / 24
                          ))
             parameters <- list(lg = rep(0, dim(d)[1]),
                                l_sigma = c(0, 0),
                                l_s_g = 0
                                )

             rdm <- "lg"
           },
           tva_cov = {
             ## TMB - data and parameter lists
             data <- with(d,
                          list(
                            model = 1,
                            x = cbind(lon, lat),
                            dt = as.numeric(difftime(date[-1], date[-nrow(d)],
                                                     unit = "hours")) / 24,
                            m = mm
                          ))
             parameters <- list(
                                lg = rep(0, dim(d)[1]),
                                l_sigma = c(0, 0),
                                l_s_g = 0,
                                b = rep(0, dim(mm)[2])
                                )
             rdm <- "lg"
           },
           tvae = {
             tmb <- argos2tmb(d, tstep = tstep, amf = amf)

             ## Initialize parameters from loess smooths of the track
             fit.lon <- loess(
               lon ~ as.numeric(date),
               data = d,
               span = span,
               na.action = "na.exclude",
               control = loess.control(surface = "direct")
             )
             fit.lat <- loess(
               lat ~ as.numeric(date),
               data = d,
               span = span,
               na.action = "na.exclude",
               control = loess.control(surface = "direct")
             )

             ## Predict track for x starting values
             xs <-
               cbind(predict(fit.lon, newdata = data.frame(date = as.numeric(tmb$ts))),
                     predict(fit.lat, newdata = data.frame(date = as.numeric(tmb$ts))))

             ds <- xs[-1,] - xs[-nrow(xs),]
             es <- ds[-1,] - g0 * ds[-nrow(ds),]

             ## Estimate components of variance
             V <- cov(es)
             sigma <- sqrt(diag(V))
             rho <- V[1, 2] / prod(sqrt(diag(V)))
             tau <- c(sd(fit.lon$residuals / tmb$K[, 1]),
                      sd(fit.lat$residuals / tmb$K[, 2]))

             ## TMB - data and parameter lists
             data <- list(
               model = 2,
               y = tmb$y,
               idx = tmb$idx,
               w = tmb$w,
               K = tmb$K,
               nu = nu
             )
             parameters <- list(
               x = xs,
               lg = qlogis(rep(g0, nrow(es))),
               l_sigma = log(sigma),
               l_s_g = 0,
               l_tau = log(tau)
             )

             rdm <- c("x", "lg")
           })

    ## TMB - create objective function
    obj <-
      MakeADFun(
        data,
        parameters,
        map = NULL,
        random = rdm,
        DLL = "rwamr",
        silent = !verbose,
        hessian = TRUE
      )

    obj$env$inner.control$trace <- verbose
    obj$env$tracemgc <- verbose

    ## Minimize objective function
    opt <- suppressWarnings(switch(
      optim,
      nlminb = nlminb(obj$par, obj$fn, obj$gr),
      optim = do.call("optim", obj)
    ))

    ## Parameters, states and the fitted values
    rep <- sdreport(obj)
    fxd <- summary(rep, "report")
    rdm <- summary(rep, "random")

    switch(model,
           tva = {
             fitted <- data_frame(date = d$date,
                                  g = plogis(rdm[, 1]),
                                  g.se = rdm[, 2])
           },
           tva_cov = {
             fitted <- data_frame(date = d$date,
                                  g = plogis(rdm[, 1]),
                                  g.se = rdm[, 2])
           },
           tvae = {
             g_est <- rdm[rownames(rdm) %in% "lg", ]
             x <- rdm[rownames(rdm) %in% "x", ]
             fitted <- data_frame(
               date = tmb$ts,
               lon = x[1:(nrow(x) / 2), 1],
               lat = x[(nrow(x) / 2 + 1):nrow(x), 1],
               lon.se = x[1:(nrow(x) / 2), 2],
               lat.se = x[(nrow(x) / 2 + 1):nrow(x), 2],
               g = c(NA, plogis(g_est[, 1]), NA),
               lg.se = c(NA, g_est[, 2], NA)
             )
           })

    if (optim == "nlminb") {
      aic <- 2 * length(opt[["par"]]) + 2 * opt[["objective"]]
    }
    else {
      aic <- NA
    }

    list(
      fitted = fitted,
      par = fxd,
      data = d,
      tmb = obj,
      opt = opt,
      aic = aic
    )

  }
ianjonsen/bssm documentation built on July 3, 2017, 10:33 p.m.