R/fit_rw_old.r

Defines functions fit_rw_old

Documented in fit_rw_old

##' Random Walk with time-varying autocorrelation
##'
##' Fit a random walk with time-varying autocorelation, via TMB,
##' to an ARGOS track and predict locations on a regular time step.
##'
##' The input track is given as a dataframe where each row is an
##' observed location and columns
##' \describe{
##' \item{'date'}{observation time (POSIXct, GMT),}
##' \item{'lon'}{observed longitude,}
##' \item{'lat'}{observed latitude,}
##' \item{'lc'}{ARGOS location class.}
##' }
##'
##' @title Random Walk with time-varying autocorrelation
##' @param d a dataframe of observations (see details)
##' @param tstep the time step to predict to (in days)
##' @param nu degrees of freedom parameter
##' @param gamma initial autocorrelation parameter
##' @param optim numerical optimizer
##' @param verbose report progress during minimization
##' @param amf ARGOS error multiplication factors
##' @param span loess parameter for location initial values
##' @return a list with components
##' \item{\code{predicted}}{a dataframe of estimated states; locations and gamma(t)}
##' \item{\code{fitted}}{a dataframe of fitted locations}
##' \item{\code{par}}{model parameter summmary}
##' \item{\code{data}}{input dataframe}
##' \item{\code{tmb}}{the TMB object}
##' \item{\code{opt}}{the object returned by the optimizer}
##' @useDynLib rwamr
##' @importFrom TMB MakeADFun sdreport
##' @importFrom tibble as_data_frame
##' @export
fit_rw_old <-
  function(d,
           optim = c("nlminb", "optim"),
           verbose = FALSE) {
    optim <- match.arg(optim)

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

    ## TMB - data and parameter lists
    if(d$error) {
      data <- list(y = cbind(d$data$lon.obs, d$data$lat.obs),
                   K = cbind(d$data$AMFlon, d$data$AMFlat),
                   nu = d$nu,
                   error = 1)
      parameters <- list(
        l_g = rep(0, e.rows),
        l_sigma = c(0,0),
#        l_rho = 0,
        l_sigma_g = 0,
        x = data$y,
        l_tau = c(0,0)#,
#        l_nu = 0
      )
      rdm <- c("x", "l_g")

    }
    else {
      data <- list(x = cbind(d$data$lon, d$data$lat), error = 0)
      parameters <- list(
        l_g = rep(0, e.rows),
        l_sigma = c(0, 0),
        l_rho = 0,
        l_sigma_g = 0
      )
      rdm <- "l_g"
    }

    ## 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")

    if(d$error == 0) {
      lg <- summary(rep, "random")
      fitted <- data_frame(
        date = d$data$date,
        l_g = c(NA, lg[, 1], NA),
        g = c(NA, plogis(lg[, 1]), NA),
        g.se = c(NA, lg[, 2], NA)
      )
    }
    else {
      rdm <- summary(rep, "random")
      x <- rdm[rownames(rdm) %in% "x", ]
      lg <- rdm[rownames(rdm) %in% "l_g", ]
#      x <- matrix(rdm[1:(nrow(d$data) * 4)], nrow(d$data), 4,
#                  dimnames = list(NULL, c("lon", "lat", "lon.se", "lat.se")))
#      lg <- matrix(rdm[(length(x) + 1):length(rdm)],
#                   ncol = 2, dimnames = list(NULL, c("lg", "lg.se")))
browser()
      fitted <- data_frame(
        date = d$data$date,
        lon = x[, 1],
        lat = x[, 2],
        lon.se = x[, 3],
        lat.se = x[, 4],
        l_g = c(NA, lg[, 1], NA),
        l_g.se = c(NA, lg[, 2], NA),
        g = c(NA, plogis(lg[, 1]), NA),
      )
    }



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

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

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