ianjonsen / foieGras

Compare 95716f5 ... +7 ... ce6218f

No flags found

Use flags to group coverage reports by test type, project and/or folders.
Then setup custom commit statuses and notifications for each flag.

e.g., #unittest #integration

#production #enterprise

#frontend #backend

Learn more about Codecov Flags here.


@@ -16,8 +16,13 @@
Loading
16 16
  DATA_VECTOR(dt);         	        //  time diff in some appropriate unit. this should contain dt for both interp and obs positions.
17 17
  DATA_VECTOR(state0);              //  initial state
18 18
  DATA_IVECTOR(isd);                //  indexes observations (1) vs. interpolation points (0)
19 +
  DATA_IVECTOR(fidx);
20 +
  DATA_VECTOR(fdt);
21 +
  DATA_IVECTOR(pidx);
22 +
  DATA_VECTOR(pdt);
19 23
  DATA_IVECTOR(obs_mod);            //  indicates which obs error model to be used
20 24
  DATA_ARRAY_INDICATOR(keep, Y);    // for one step predictions
25 +
  DATA_SCALAR(se);                  // turn delta method on/off for sv (speeds up model fitting if SE's not req'd & this allows param SE's to be calculated regardless)
21 26
  
22 27
  // for KF observation model
23 28
  DATA_VECTOR(m);                 //  m is the semi-minor axis length
@@ -51,7 +56,8 @@
Loading
51 56
  Type jnll = 0.0;
52 57
  Type tiny = 1e-5;
53 58
  
54 -
  vector<Type> sv = dt.size();
59 +
  vector<Type> sf = fidx.size();
60 +
  vector<Type> sp = pidx.size();
55 61
  
56 62
  // Setup object for evaluating multivariate normal likelihood
57 63
  matrix<Type> cov(4,4);
@@ -85,8 +91,6 @@
Loading
85 91
    x_t(2) = (v(0,i) - v(0,i-1)); // /dt(i);
86 92
    x_t(3) = (v(1,i) - v(1,i-1)); // /dt(i);
87 93
    
88 -
    // 2-D velocity
89 -
    sv(i) = sqrt(pow(v(0,i) - v(0,i-1), 2) + pow(v(1,i) - v(1,i-1), 2)); 
90 94
    jnll += MVNORM<Type>(cov)(x_t); // Process likelihood
91 95
  }
92 96
  
@@ -145,6 +149,15 @@
Loading
145 149
    }
146 150
  }
147 151
  
152 +
// calculate 2-D vel along track separately for fitted (isd==1) vs predicted (isd==0) states
153 +
154 +
  for(int i = 1; i < fidx.size(); ++i) {
155 +
    sf(i) = sqrt(pow(mu(0, fidx(i)) - mu(0, fidx(i-1)), 2) + pow(mu(1, fidx(i)) - mu(1, fidx(i-1)), 2)) / fdt(i);
156 +
  }
157 +
  for(int i = 1; i < pidx.size(); ++i) {
158 +
    sp(i) = sqrt(pow(mu(0, pidx(i)) - mu(0, pidx(i-1)), 2) + pow(mu(1, pidx(i)) - mu(1, pidx(i-1)), 2)) / pdt(i);
159 +
  }
160 +
  
148 161
  SIMULATE {
149 162
    REPORT(Y);
150 163
  }
@@ -153,7 +166,13 @@
Loading
153 166
  ADREPORT(rho_o);
154 167
  ADREPORT(tau);
155 168
  ADREPORT(psi);
156 -
  ADREPORT(sv);
169 +
  if(se == 1) {
170 +
    ADREPORT(sf);
171 +
    ADREPORT(sp);
172 +
  } else if(se == 0) {
173 +
    REPORT(sf);
174 +
    REPORT(sp);
175 +
  }
157 176
  
158 177
  return jnll;
159 178
}

@@ -22,9 +22,8 @@
Loading
22 22
##' defaults are used when \code{method = "L-BFGS-B"}. Possible parameter names are same as \code{lower}
23 23
##' @param verbose integer; report progress during minimization: 0 = silent;
24 24
##' 1 = optimizer trace; 2 = parameter trace (default))
25 -
##' @param se logical; should standard errors for fixed effects be calculated (default = TRUE). 
26 -
##' Turning this off will speed up computation time at the expense of reporting uncertainty for 
27 -
##' fixed effects
25 +
##' @param se logical; should standard errors for speed estimates be calculated (default = FALSE). 
26 +
##' Turning this on will slow down computation time but provide SE's for speed-along-track calculations
28 27
##' @param ... control parameters for the chosen optimizer
29 28
##' @return Returns a list with components
30 29
##'   \item{\code{optim}}{the name of the numerical optimizer as a
@@ -53,7 +52,7 @@
Loading
53 52
           lower = NULL,
54 53
           upper = NULL,
55 54
           verbose = 1,
56 -
           se = TRUE,
55 +
           se = FALSE,
57 56
           ...) {
58 57
    optim <- match.arg(optim)
59 58
    method <- match.arg(method)

@@ -4,17 +4,21 @@
Loading
4 4
##' either a set of x,y (or lon,lat) coordinates from a \code{fG_ssm} fit with length
5 5
##' equal to the number of observations used in the SSM fit. 
6 6
##' @param x a compound \code{fG_ssm} model fit object (ignored if NULL)
7 -
##' @param reps number of replicate tracks to simulate from an \code{fG_ssm} model 
8 -
##' fit object (ignored if x is NULL)
9 7
##' @param what simulate fitted (typically irregular in time) or predicted 
10 8
##' (typically regular in time) locations 
11 -
##' @param cpf logical; should simulated tracks return to their start point (ie. a central-place forager)
9 +
##' @param reps number of replicate tracks to simulate from an \code{fG_ssm} model 
10 +
##' fit object
11 +
##' @param grad a rasterStack of x- and y-gradients as separate layers
12 +
##' @param beta a 2-element vector of parameters defining the potential function 
13 +
##' magnitude in x- and y-directions (ignored if \code{is.null(grad)}, 
14 +
##' ie. no potential function)
15 +
##' @param cpf logical; should simulated tracks return to their start piont (ie. a central-place forager)
12 16
##' @param sim_only logical, do not include \code{fG_ssm} estimated location in output 
13 17
##' (default is FALSE)
14 18
##' 
15 19
##' @examples 
16 -
##' fit <- fit_ssm(ellie, vmax = 4, model = "crw", time.step = 48, control = ssm_control(se = FALSE))
17 -
##' trs <- simfit(fit, reps = 2, what = "predicted")
20 +
##' fit <- fit_ssm(ellies, vmax = 4, model = "crw", time.step = 72)
21 +
##' trs <- simfit(fit, what = "predicted", reps = 2)
18 22
##' plot(trs)
19 23
##' 
20 24
##' @importFrom dplyr "%>%" 
@@ -24,19 +28,31 @@
Loading
24 28
##' @importFrom assertthat assert_that
25 29
##' @importFrom sf st_coordinates st_as_sf st_transform st_geometry<-
26 30
##' @importFrom stats rgamma runif
31 +
##' @importFrom raster extent extract nlayers
27 32
##' @importFrom CircStats rvm
28 33
##' @export
29 34
30 -
simfit <- function(x, 
31 -
                   what = c("fitted", "predicted"),
32 -
                   reps = 1,
33 -
                   cpf = FALSE,
34 -
                   sim_only = FALSE) {
35 -
35 +
simfit <-
36 +
  function(x,
37 +
           what = c("fitted", "predicted"),
38 +
           reps = 1,
39 +
           grad = NULL,
40 +
           beta = c(-250, -250),
41 +
           cpf = FALSE,
42 +
           sim_only = FALSE) {
43 +
    
44 +
  
36 45
  ################
37 46
  ## Check args ##
38 47
  ################
39 48
  what <- match.arg(what)
49 +
  if(!is.null(grad)) {
50 +
    if(!inherits(grad, "RasterStack")) 
51 +
      stop("grad must be NULL or a RasterStack with 2 layers")
52 +
    if(inherits(grad, "RasterStack") & nlayers(grad) != 2)
53 +
      stop("grad RasterStack must have 2 layers")
54 +
  }
55 +
  if(length(beta) != 2) stop("beta must be specified as a 2-element vector")
40 56
  
41 57
  assert_that(what %in% c("fitted", "predicted"),
42 58
              msg = "only `fitted` or `predicted` locations can be simulated
@@ -50,10 +66,10 @@
Loading
50 66
    model <- x$ssm[[k]]$pm
51 67
    switch(what,
52 68
           fitted = {
53 -
             loc <- grab(x[k, ], "fitted")
69 +
             loc <- grab(x[k,], "fitted", as_sf = FALSE)
54 70
           },
55 71
           predicted = {
56 -
             loc <- grab(x[k, ], "predicted")
72 +
             loc <- grab(x[k,], "predicted", as_sf = FALSE)
57 73
           })
58 74
    N <- nrow(loc)
59 75
    dts <- loc$date
@@ -73,137 +89,152 @@
Loading
73 89
           },
74 90
           rw = {
75 91
             Sigma <-
76 -
               diag(2) * c(x$ssm[[k]]$par["sigma_x", 1],
92 +
               diag(2) * c(x$ssm[[k]]$par["sigma_x", 1], 
77 93
                           x$ssm[[k]]$par["sigma_y", 1]) ^ 2
78 94
             Sigma[!Sigma] <-
79 -
               prod(Sigma[1, 1] ^ 0.5, Sigma[2, 2] ^ 0.5) * x$ssm[[k]]$par["rho_p", 1]
80 -
             vmin <- with(loc, c(min(diff(x), na.rm = TRUE),
95 +
               prod(Sigma[1, 1]^0.5, Sigma[2, 2]^0.5) * x$ssm[[k]]$par["rho_p", 1]
96 +
             vmin <- with(loc, c(min(diff(x), na.rm = TRUE), 
81 97
                                 min(diff(y), na.rm = TRUE)))
82 -
             vmax <- with(loc, c(max(diff(x), na.rm = TRUE),
98 +
             vmax <- with(loc, c(max(diff(x), na.rm = TRUE), 
83 99
                                 max(diff(y), na.rm = TRUE)))
84 100
           })
85 101
    
86 102
    ###############################
87 103
    ## Simulate movement process ##
88 104
    ###############################
105 +
    if(!is.null(grad)) {
106 +
      ex <- extent(grad[[1]])
107 +
    } else {
108 +
      ex <- c(-20077.51,20082.49,-19622.54,18437.46) ## approx extents of world mercator in km
109 +
    }
110 +
    wrap_x <- function(mu, x_rng) c((mu[1] - x_rng[1]) %% sum(abs(x_rng)) + x_rng[1], mu[2])
111 +
    reflect_y <- function(mu, y_rng) {
112 +
      if(mu[2] < y_rng[1]) {
113 +
        c(mu[1], y_rng[1] * 2 - mu[2])
114 +
      } else if(mu[2] > y_rng[2]) {
115 +
        c(mu[1], y_rng[2] * 2 - mu[2])
116 +
      } else {
117 +
        mu
118 +
      }
119 +
    }
120 +
    
89 121
    tmp <- lapply(1:reps, function(j) {
90 122
      switch(model,
91 123
             crw = {
92 124
               mu <- v <- matrix(NA, N, 2)
93 -
               mu[1, ] <- st_coordinates(loc$geometry)[1, ]
94 -
               v[1,] <- c(loc$u[1], loc$v[1])
125 +
               time <- seq(0, 1, length = N)
126 +
               mu[1,] <- c(loc$x[1], loc$y[1])
127 +
               v[1, ] <- c(loc$u[1], loc$v[1])
95 128
               for (i in 2:N) {
96 -
                 v[i, ] <- rtmvnorm(1,
97 -
                                    v[i - 1, ],
98 -
                                    sigma = Sigma * dt[i],
99 -
                                    lower = vmin,
100 -
                                    upper = vmax)
101 -
                 mu[i, ] <- mu[i - 1, ] + v[i, ] * dt[i]
102 -
                 ## keep within world Mercator y bounds (km)
103 -
                 #                     if(mu[i, 2] < -15496300) mu[i, 2] <- -15496300
104 -
                 #                     if(mu[i, 2] > 18764386) mu[i, 2] <- 18764386
129 +
                 v[i,] <- rtmvnorm(1,
130 +
                                   v[i - 1,],
131 +
                                   sigma = Sigma * dt[i],
132 +
                                   lower = vmin,
133 +
                                   upper = vmax)
134 +
                 ## wrap x values, reflect y values
135 +
                 mu1 <- wrap_x(mu[i-1,] + v[i,] * dt[i], ex[1:2])
136 +
                 mu1 <- reflect_y(mu1, ex[3:4])
137 +
                 if(!is.null(grad)) {
138 +
                   pv <- c(extract(grad[[1]], rbind(mu1))[1],
139 +
                           extract(grad[[2]], rbind(mu1))[1])
140 +
                   mu[i,] <- mu1 + pv * beta
141 +
                 } else {
142 +
                   mu[i,] <- mu1
143 +
                 }
105 144
               }
106 145
               df <- data.frame(
107 146
                 rep = j,
108 147
                 date = dts,
109 148
                 x = mu[, 1],
110 -
                 y = mu[, 2],
111 -
                 u = v[, 1],
112 -
                 v = v[, 2]
149 +
                 y = mu[, 2]
113 150
               )
114 151
               if(cpf) {
115 -
                 time <- seq(0, 1, length = N)
116 -
                 bb.x <- with(df, ((x - x[1]) - time * (x - x[1])[N]) + x[1])
117 -
                 bb.y <- with(df, ((y - y[1]) - time * (y - y[1])[N]) + y[1])
118 152
                 df <- df %>%
119 -
                   mutate(x = bb.x,
120 -
                          y = bb.y)
153 +
                   mutate(x = ((x - x[1]) - time *
154 +
                                 (x - x[1])[N]) + x[1]) %>%
155 +
                   mutate(y = ((y - y[1]) - time *
156 +
                                 (y - y[1])[N]) + y[1])
157 +
                 
121 158
               }
122 159
               df
123 160
             },
124 161
             rw = {
125 -
               mu <- matrix(NA, N, 2)
126 -
               mu[1,] <- st_coordinates(loc$geometry)[1, ]
127 -
               mu[2,] <-
128 -
                 rmvnorm(1, mu[1, ], sigma = Sigma * dt[2] ^ 2)
162 +
               mu <- matrix(NA, N, 2) 
163 +
               time <- seq(0, 1, length = N)
164 +
               mu[1, ] <- c(loc$x[1], loc$y[1])
165 +
               mu[2, ] <- rmvnorm(1, mu[1,], sigma = Sigma * dt[2]^2)
129 166
               for (i in 3:N) {
130 -
                 dxy <- rtmvnorm(
131 -
                   1,
132 -
                   mu[i - 1,] - mu[i - 2, ],
133 -
                   sigma = Sigma * dt[i] ^ 2,
134 -
                   lower = vmin,
135 -
                   upper = vmax
136 -
                 )
137 -
                 mu[i,] <- mu[i - 1,] + dxy
138 -
                 ## keep within world Mercator y bounds (km)
139 -
                 #                     if(mu[i, 2] < -15496300) mu[i, 2] <- -15496300
140 -
                 #                     if(mu[i, 2] > 18764386) mu[i, 2] <- 18764386
167 +
                 dxy <- rtmvnorm(100, mu[i - 1, ] - mu[i - 2,], 
168 +
                                 sigma = Sigma * dt[i]^2,
169 +
                                 lower = vmin,
170 +
                                 upper = vmax,
171 +
                                 algorithm = "gibbs", burn.in.samples = 100)
172 +
                 dxy <- dxy[which(!is.na(dxy))[1],]
173 +
                 ## wrap x values, reflect y values
174 +
                 mu1 <- wrap_x(mu[i-1,] + dxy, ex[1:2])
175 +
                 mu1 <- reflect_y(mu1, ex[3:4])
176 +
                 if(!is.null(grad)) {
177 +
                   pv <- c(extract(grad[[1]], rbind(mu1))[1],
178 +
                           extract(grad[[2]], rbind(mu1))[1])
179 +
                   mu[i,] <- mu1 + pv * beta
180 +
                 } else {
181 +
                   mu[i,] <- mu1
182 +
                 }
141 183
               }
142 184
               df <- data.frame(
143 185
                 rep = j,
144 186
                 date = dts,
145 187
                 x = mu[, 1],
146 188
                 y = mu[, 2]
147 189
               )
148 -
               ## convert to Brownian-Bridge to simulate a central-place forager
149 190
               if(cpf) {
150 -
                 n <- nrow(df)
151 -
                 time <- seq(0, 1, length = N)
152 -
                 df <- df %>% 
153 -
                   mutate(x = ((x - x[1]) - time * (x - x[1])[N]) + x[1],
154 -
                          y = ((y - y[1]) - time * (y - y[1])[N]) + y[1]
155 -
                          )
191 +
                 df <- df %>%
192 +
                   mutate(x = ((x - x[1]) - time * 
193 +
                                 (x - x[1])[N]) + x[1]) %>%
194 +
                   mutate(y = ((y - y[1]) - time * 
195 +
                                 (y - y[1])[N]) + y[1])
156 196
               }
157 197
               df
158 198
             })
159 -
    }) %>%
199 +
    }) %>% 
160 200
      do.call(rbind, .) %>%
161 201
      as_tibble()
162 -
    
202 +
163 203
    if (!sim_only) {
164 -
      loc <- grab(x[k,], what = what, as_sf = FALSE) %>%
204 +
      loc <- grab(x[k, ], what = what, as_sf = FALSE) %>%
165 205
        mutate(rep = 0)
166 -
      switch(model,
167 -
             crw = {
168 -
               loc <- loc %>% select(rep, date, x, y, u, v)
169 -
             },
170 -
             rw = {
171 -
               loc <- loc %>% select(rep, date, x, y)
172 -
             })
206 +
      loc <- loc %>% select(rep, date, x, y)
173 207
      tmp <- rbind(loc, tmp)
174 208
    }
175 209
    
176 210
    ## lon,lat
177 -
    tmp <-
178 -
      st_as_sf(tmp, coords = c("x", "y"), crs = "+proj=merc +units=km +datum=WGS84")
179 -
    xy <- st_coordinates(tmp) %>% as.data.frame()
180 -
    names(xy) <- c("x", "y")
181 -
    ll <- tmp %>%
211 +
    tmp1 <- try(st_as_sf(tmp, coords = c("x","y"), crs = "+proj=merc +units=km +datum=WGS84"), silent = TRUE)
212 +
    if(inherits(tmp1, "try-error")) {
213 +
      stop("oops something went wrong, try again", call. = FALSE)
214 +
    }
215 +
    
216 +
    xy <- st_coordinates(tmp1) %>% as.data.frame()
217 +
    names(xy) <- c("x","y")
218 +
    ll <- tmp1 %>% 
182 219
      st_transform(., crs = 4326) %>%
183 220
      st_coordinates(.) %>%
184 221
      as.data.frame(.)
185 -
    names(ll) <- c("lon", "lat")
186 -
    st_geometry(tmp) <- NULL
187 -
    cbind(tmp, xy, ll) %>%
222 +
    names(ll) <- c("lon","lat")
223 +
    st_geometry(tmp1) <- NULL
224 +
    cbind(tmp1, xy, ll) %>% 
188 225
      as_tibble() %>%
189 226
      select(rep, date, lon, lat, x, y, everything())
190 -
  })
191 -
227 +
  }) 
192 228
  
193 -
    
194 -
  d <- tibble(id = x$id,
195 -
              model = x$pmodel,
196 -
              sims = d)
229 +
  d <- tibble(id = x$id, model = x$pmodel, sims = d)
197 230
  switch(unique(x$pmodel),
198 -
         rw = {
231 +
         rw = { 
199 232
           class(d) <- append("fG_rws", class(d))
200 233
         },
201 234
         crw = {
202 235
           class(d) <- append("fG_crws", class(d))
203 236
         })
204 -
  
237 +
205 238
  class(d) <- append("fG_simfit", class(d))
206 239
  return(d)
207 -
208 240
}
209 -
                

@@ -261,21 +261,34 @@
Loading
261 261
               y = y - mean(y, na.rm = TRUE) / sd(y, na.rm = TRUE))
262 262
    }
263 263
    
264 +
    ## calculate fitted & predicted indices + delta t for proper speed estimates
265 +
    fidx <- which(d.all$isd)
266 +
    
267 +
    fdt <- as.numeric(difftime(d.all$date[fidx], lag(d.all$date[fidx]), 
268 +
                               units = "hours"))
269 +
    pidx <- which(!d.all$isd)
270 +
    pdt <- as.numeric(difftime(d.all$date[pidx], lag(d.all$date[pidx]), 
271 +
                               units = "hours"))
264 272
    
265 273
    data <- list(
266 274
      model_name = model,
267 275
      Y = rbind(d.all$x, d.all$y), 
268 276
      dt = dt,
269 277
      state0 = state0,
270 278
      isd = as.integer(d.all$isd),
279 +
      fidx = fidx-1, ## for C++ indexing
280 +
      fdt = fdt,
281 +
      pidx = pidx-1,
282 +
      pdt = pdt,
271 283
      obs_mod = as.integer(obs_mod),
284 +
      se = as.integer(control$se),
272 285
      m = d.all$smin,
273 286
      M = d.all$smaj,
274 287
      c = d.all$eor,
275 288
      K = cbind(d.all$emf.x, d.all$emf.y),
276 289
      GLerr = cbind(d.all$lonerr, d.all$laterr)
277 290
    )
278 -
    
291 +
279 292
    if(scale) d.all <- d.all.tmp
280 293
    
281 294
    ## TMB - create objective function
@@ -371,8 +384,8 @@
Loading
371 384
    
372 385
    ## if error then exit with limited output to aid debugging
373 386
    ## check if pdHess is FALSE at end and return warning
374 -
    rep <- try(sdreport(obj, skip.delta.method = !control$se)) #
375 -
   
387 +
    rep <- try(sdreport(obj)) #, skip.delta.method = !control$se
388 +
376 389
    options(warn = oldw) ## turn warnings back on
377 390
    
378 391
    if (!inherits(opt, "try-error") & !inherits(rep, "try-error")) {
@@ -386,6 +399,16 @@
Loading
386 399
      if("tau" %in% row.names(fxd)) {
387 400
        row.names(fxd)[which(row.names(fxd) == "tau")] <- c("tau_x","tau_y")
388 401
      }
402 +
      ## separate speed+se if present
403 +
      if("sf" %in% row.names(fxd)) {
404 +
        sf <- fxd[which(row.names(fxd) == "sf"), ]
405 +
        sf[1,] <- c(NA,NA)
406 +
        if("sp" %in% row.names(fxd)) {
407 +
          sp <- fxd[which(row.names(fxd) == "sp"), ]
408 +
          sp[1,] <- c(NA,NA)
409 +
        }
410 +
        fxd <- fxd[!row.names(fxd) %in% c("sf","sp"),]
411 +
      }
389 412
390 413
      switch(model,
391 414
             rw = {
@@ -419,7 +442,6 @@
Loading
419 442
               tmp <- summary(rep, "random")
420 443
               loc <- tmp[rownames(tmp) == "mu",]
421 444
               vel <- tmp[rownames(tmp) == "v",]
422 -
               ss <- fxd[which(row.names(fxd) == "sv"), ] ## speeds
423 445
               
424 446
               loc <-
425 447
                 cbind(loc[seq(1, dim(loc)[1], by = 2),],
@@ -432,8 +454,6 @@
Loading
432 454
                       vel[seq(2, dim(vel)[1], by = 2),]) %>%
433 455
                 as.data.frame(., row.names = 1:nrow(.))
434 456
               names(vel) <- c("u", "u.se", "v", "v.se")
435 -
               vel <- vel %>%
436 -
                 mutate(s = ss[, 1], s.se = ss[, 2])
437 457
               
438 458
               rdm <- bind_cols(loc, vel) %>%
439 459
                 mutate(
@@ -447,7 +467,7 @@
Loading
447 467
                          y = y  * sd(d.all$y, na.rm = TRUE) + mean(d.all$y, na.rm = TRUE))
448 468
               }
449 469
               rdm <- rdm %>%
450 -
                 select(id, date, x, y, x.se, y.se, u, v, u.se, v.se, s, s.se, isd)
470 +
                 select(id, date, x, y, x.se, y.se, u, v, u.se, v.se, isd)
451 471
             })
452 472
      
453 473
      ## coerce x,y back to sf object
@@ -458,24 +478,46 @@
Loading
458 478
      switch(model,
459 479
             rw = {
460 480
               rdm <- rdm %>% select(id, date, x.se, y.se, isd)
461 -
481 +
               ## fitted
482 +
               fv <- filter(rdm, isd) %>%
483 +
                 select(-isd)
484 +
               ##predicted
485 +
               if(all(!is.na(time.step))) {
486 +
                 pv <- filter(rdm, !isd) %>%
487 +
                   select(-isd)
488 +
               } else {
489 +
                 pv <- NULL
490 +
               }
462 491
             },
463 492
             crw = {
493 +
               if(control$se) {
494 +
                 sf.se <- sf[,2]
495 +
                 sf <- sf[,1]
496 +
                 if(all(!is.na(time.step))) {
497 +
                   sp.se <- sp[,2]
498 +
                   sp <- sp[,1]
499 +
                 }
500 +
               } else {
501 +
                 sf <- as.vector(obj$report()$sf)
502 +
                 sp <- as.vector(obj$report()$sp)
503 +
                 sf.se <- NA
504 +
                 sp.se <- NA
505 +
               }
464 506
               rdm <- rdm %>%
465 -
                 select(id, date, x.se, y.se, u, v, u.se, v.se, s, s.se, isd)
507 +
                 select(id, date, x.se, y.se, u, v, u.se, v.se, isd) 
508 +
               ## fitted
509 +
               fv <- filter(rdm, isd) %>%
510 +
                 mutate(s = sf, s.se = sf.se) %>%
511 +
                 select(-isd)
512 +
               ##predicted
513 +
               if(all(!is.na(time.step))) {
514 +
                 pv <- filter(rdm, !isd) %>%
515 +
                   mutate(s = sp, s.se = sp.se) %>%
516 +
                   select(-isd)
517 +
               } else {
518 +
                 pv <- NULL
519 +
               }
466 520
             })
467 -
468 -
      ## Fitted values (estimated locations at observation times)
469 -
      fv <- subset(rdm, isd) %>%
470 -
        select(-isd)
471 -
472 -
      ## Predicted values (estimated locations at regular time intervals, defined by `ts`)
473 -
      if(all(!is.na(time.step))) {
474 -
      pv <- subset(rdm, !isd) %>%
475 -
        select(-isd)
476 -
      } else {
477 -
        pv <- NULL
478 -
      }
479 521
      
480 522
      if (control$optim == "nlminb") {
481 523
        aic <- 2 * length(opt[["par"]]) + 2 * opt[["objective"]]

Learn more Showing 1 files with coverage changes found.

Changes in R/grab.R
-2
+2
Loading file...
Files Coverage
R -0.93% 78.40%
src 0.31% 80.57%
Project Totals (26 files) 78.65%
Loading