Cell counts

spss <- getSpssData()
Age <- zap_label(spss$LFT)
spss <- as_factor(spss)
setDT(spss)
complete_days[, username2 := as.numeric(username)]
spss[, spssusername := as.numeric(username)]
spss[, LFT := Age]
complete_days[spss, on = .(username2 == spssusername), `:=`(Sex = Geslacht,
                                        Age = LFT,
                                        Urbanicity = STEDBUURT)]

splitdat <- complete_days[, .(`Trip Time` = unique(total_trip_time)), .(person_day, Age, Sex, Urbanicity, os, device_id)]
splitdat[, Age := as.numeric(Age)]
splitdat[, Sex := factor(Sex, levels = c(1, 2), labels = c("Male", "Female"))]
splitdat[, Urban := factor(Urbanicity, levels = levels(splitdat$Urbanicity), labels = c("Very high", "High", "Moderate", "Slight", "Rural", "Unknown"))]
splitdat[, .(mean = mean(`Trip Time`), min = min(`Trip Time`), max = max(`Trip Time`), sd = sd(`Trip Time`), med = median(`Trip Time`), N = uniqueN(device_id)), .(Urban, Age = cut(Age, breaks = c(15, 18, 35, 55, 77)), Sex)][order(Urban, Age, Sex)][N > 4] %>% 
  kable(digits = 0, caption = "Total travel time split by characteristics") %>% 
  kable_styling(full_width = FALSE)

We are only able to look at the main effects of these categories. Table .. might be indicative of differential effects, but with so few people in each group, finding effects will be difficult.

splitdat[, .(mean = mean(`Trip Time`), min = min(`Trip Time`), max = max(`Trip Time`), sd = sd(`Trip Time`), med = median(`Trip Time`), N = uniqueN(device_id)), .(Urban)][order(Urban)][N > 4] %>% 
  kable(digits = 0, caption = "By urbanicity") %>% 
  kable_styling(full_width = FALSE) -> urbtab
splitdat[, .(mean = mean(`Trip Time`), min = min(`Trip Time`), max = max(`Trip Time`), sd = sd(`Trip Time`), med = median(`Trip Time`), N = uniqueN(device_id)), .(Age = cut(Age, breaks = c(15, 18, 35, 55, 77)))][order(Age)][N > 4] %>% 
  kable(digits = 0, caption = "By age") %>% 
  kable_styling(full_width = FALSE)  -> agetab
splitdat[, .(mean = mean(`Trip Time`), min = min(`Trip Time`), max = max(`Trip Time`), sd = sd(`Trip Time`), med = median(`Trip Time`), N = uniqueN(device_id)), .(Sex)][order(Sex)][N > 4] %>% 
  kable(digits = 0, caption = "By sex") %>% 
  kable_styling(full_width = FALSE)  -> sextab

cat(c('<table><tr valign="top"><td>', urbtab, '</td><td>', agetab, '</td><td>', sextab, '</td><tr></table>'),
    sep = '')

If we consider the number of people who sent us at least some location data, our N per category grows significantly. We can consider this the total possible N to achieve without considering non-response.

data[, username := as.numeric(username)]
data[spss, on = .(username == spssusername), `:=`(Sex = Geslacht, Age = LFT, Urbanicity = STEDBUURT)]
data[, Age := as.numeric(Age)]
data[, AgeC := cut(Age, breaks = c(14, 18, 35, 55, 91))]
# data[, Sex := factor(Sex, levels = c(1, 2), labels = c("Male", "Female"))]
data[, Urban := factor(Urbanicity, levels = levels(splitdat$Urbanicity), labels = c("Very high", "High", "Moderate", "Slight", "Rural", "Unknown"))]
data[, .(`Cell count` = uniqueN(username)), .(Urban, AgeC, Sex)][order(Urban, AgeC, Sex)] %>% kable() %>% kable_styling(full_width = F)

Linear interpolation

data2 <- interpolatePoints(data, min_time = minutes(3), max_time = minutes(30), min_dist = -1)

data2[, newttn := shift(timestamp, n = -1) - timestamp, person_day]
data2[, gap := TRUE]
data2[newttn < minutes(5), gap := FALSE]

sg_complete <-
  copy(data2[gap == FALSE, per_day_coverage := sum(newttn), person_day][per_day_coverage > minutes(1368)])

sg_complete_pd <- sg_complete[, unique(person_day)]

for (i in seq_along(sg_complete_pd)) {
  this_person_day <- sg_complete_pd[i]
  res <-
    stopFinder(copy(sg_complete[person_day == this_person_day]),
               thetaD = 200,
               thetaT = minutes(5))
  
  mergingCycle(
    res,
    thetaD = 200,
    small_track_action = "exclude",
    max_time = minutes(5),
    max_locs = Inf,
    max_dist = 200
  )
  
  events <- returnStateEvents(res)
  
  if (any(events$state == "moving")) {
    segs <- getSegments(tdtr(
      res[state == "moving"],
      col_names = list(
        entity_id_col = "device_id",
        timestamp_col = "timestamp",
        latitude_col
        = "latitude",
        longitude_col = "longitude"
      ),
      group_col = "move_id",
      max_error = 30
    ),
    group = TRUE)
    
    set(
      events,
      j = c(
        "total_travel_dist",
        "total_trip_time",
        "avg_trip_time",
        "n_trips"
      ),
      value = list(
        sum(segs[["segdist"]]),
        events[state == "moving", sum(time_length(as.duration(end_time - begin_time), "minutes"))],
        events[state == "moving", mean(time_length(as.duration(end_time - begin_time), "minutes"))],
        max(events[["move_id"]], na.rm = TRUE)
      )
    )
    events[, avg_dist_per_trip := total_travel_dist / n_trips]
    
    # events[, total_travel_dist := sum(segs$segdist)]
    # events[, n_trips := max(move_id, na.rm = TRUE)]
    # events[, total_trip_time := sum(end_time[state == "moving"] - begin_time[state == "moving"])]
    # events[, avg_trip_time := mean(end_time[state == "moving"] - begin_time[state == "moving"])]
    
  } else {
    set(
      events,
      j = c(
        "total_travel_dist",
        "avg_dist_per_trip",
        "total_trip_time",
        "avg_trip_time",
        "n_trips"
      ),
      value = list(0)
    )
  }
  if (any(events$state == "stopped")) {
    events[, n_stops := max(stop_id, na.rm = TRUE)]
    
  } else {
    set(events,
        j = c("n_stops"),
        value = list(0))
  }
  
  sg_complete[person_day == this_person_day,
                c(
                  "total_travel_dist",
                  "avg_dist_per_trip",
                  "total_trip_time",
                  "avg_trip_time",
                  "n_trips",
                  "n_stops",
                  "rog"
                ) :=
                  c(unique(events[, .(
                    total_travel_dist,
                    avg_dist_per_trip,
                    total_trip_time,
                    avg_trip_time,
                    n_trips,
                    n_stops
                  )]),
                  radiusOfGyrationDT(latitude, longitude, timestamp))]
}
sg_complete[, day_of_week := na.exclude(day_of_week)[1], person_day]

tabdat_sg_day <-
  unique(sg_complete[, .(
    day_of_week,
    n_trips,
    n_stops,
    total_travel_dist,
    avg_dist_per_trip,
    total_trip_time,
    avg_trip_time,
    rog
  ), person_day])

Using only interpolation, we are able to use 604 days across 160 persons. We are interested in what can be gained from also filling in the long gaps. DTWBI looks at the pattern of movement behavior over time, aggregating it into time series with the goal of finding similar patterns. When there exists a long gap within a respondent’s day, in order to generate informative data against which we can compare, we need a buffer of sufficient length on either side of the missing data. The length of this buffer can vary – the shorter the buffer, the less specific we are when selecting imputation candidates. For our purposes, a buffer of at least one hour is necessary.

Problems

Because respondents may have multiple gaps within a 24-hour period and because some gaps may occur at the beginning or end of the period, not all days can be imputed with DTWBI. We require that the total sum of covered time before and after the gap is at least 110 minutes. To reduce our workload, if we can calculate that we will not have at least 95% of a full day covered (1368 minutes), we will not impute that day. It might be possible in the future to include these days so that they can become reference queries for imputation, even if they will lack sufficient information to be considered themselves.

data2[gap == FALSE, new_per_day_coverage := sum(newttn), person_day]
data2[, new_complete := FALSE]
data2[new_per_day_coverage >= minutes(1368), new_complete := TRUE]
data2[new_complete == FALSE, candtw := bufferTest(.SD, hours(1), minutes(1368)), person_day]

Requiring a buffer on both sides for DTWBI only provides us with 187 new days. This is suboptimal, and it occurs because we have prematurely divided our dataset into days. If instead, we maintained the data sets as uninterrupted streams until after the imputation, we could often get the data on both sides of, for example, night missingness. We have two options in this case: 1. DTWBI for beginning/end missingness, which incorporates match.buffers on only one side, or 2. perform DTWBI before the data are split on days, then split after imputation.

# data2[, lastttn := (hms::as_hms("04:00:00") - hms::as_hms(timestamp[.N])), person_day]
# data2[lastttn < 0, lastttn := lastttn + as.duration(hours(24))]
# data2[is.na(newttn), newttn := lastttn, person_day]

A problem we run into with solution 2 is that many of the long gaps (which could be handled with DTWBI) are interspersed with single measurements, which keeps us from having the full length of matching buffer. Essentially having those extra pieces of data keep the model from working. The figure below illustrates this issue. The black line shows the gap time, and the color on either side of it shows the one hour match buffer. The best solution here is likely to combine the intersecting gaps into a single gap, making note of the distance that is noted (if any) when we receive the single points flanked by long gaps. Then handle as necessary.

data3 <- copy(data2[, .(device_id, latitude, longitude, accuracy, timestamp, timestamp_date, adj_date, day_of_week, person_day)])
data3[, tfp := timestamp - shift(timestamp, n = 1), device_id]
# data3[, tfp := shift(ttn)]
data3[, newid := rleid(device_id, !is.na(tfp) & ifelse(tfp > hours(12), tfp, 0))]
data3[, person_block := paste(device_id, newid, sep = "_")]

data3[, ttn := shift(timestamp,n = -1) - timestamp, person_block]
data3[, gap := TRUE]
data3[ttn < minutes(5), gap := FALSE]

test <- data3[person_block == "100_1"]
gaps <- test[gap == TRUE, .(person_block, gapid = .I, gaptime = timestamp, ttn, start = timestamp - hours(1), end = timestamp + ttn + hours(1))]
ggplot(gaps, aes(x = gaptime, y = gapid)) +
  geom_linerange(aes(xmin = start, xmax = end, color = factor(gapid))) +
  geom_linerange(aes(xmin = gaptime, xmax = gaptime + ttn), fill = "black", size = 1) + 
  scale_x_datetime("Time")+
  scale_y_continuous("Gap number")+
  scale_color_discrete("Gap number")+
  theme_minimal() +
  theme(legend.position = "none")

gaps <- data3[gap == TRUE, .(gapid = 1:.N, gaptime = timestamp, start = timestamp - hours(1), end = timestamp + ttn + hours(1)), person_block]
gaps[, newgap := mergeGap(end, start), person_block]
gaps[, `:=`(newstart = min(start), newend = max(end, na.rm = TRUE)), .(person_block, newgap)]
gaps[start == newstart, forebuffer_end := gaptime]
gaps[end == newend, endbuffer_start := gaptime]
gaps[, `:=`(endbuffer_start = na.exclude(endbuffer_start)[1], forebuffer_end = na.exclude(forebuffer_end)[1]), .(person_block, newgap)]
data3[gaps,
     on = .(person_block, timestamp >= newstart, timestamp <= newend),
     `:=`(newgap = i.newgap, newstart = i.newstart, newend = i.newend, forebuffer_end = i.forebuffer_end, endbuffer_start = i.endbuffer_start)]

data3[, extra_point := FALSE]
data3[!is.na(newgap) & gap == FALSE, extra_point := timestamp %between% c(forebuffer_end[1], endbuffer_start[1]), .(person_block, newgap)]
data3[!is.na(newgap), candtw := testBufferSingleGap(.SD), .(person_block, newgap)]
data3[!is.na(newgap), atleast24h := testBufferAfterMerge(.SD, min.day.minutes = minutes(1368)), .(person_block)]
data3[gap == TRUE | extra_point == TRUE, interdist := sum(geodist_vec(longitude, latitude, sequential = TRUE, measure = "haversine")), .(person_block, newgap)]
data3[!is.na(newgap), interdist := na.exclude(interdist)[1], .(person_block, newgap)]
data3[, newgaplength := (endbuffer_start[1] - forebuffer_end[1]), .(person_block, newgap)]
data3[newgaplength > hours(12), candtw := FALSE]
data3[, aftermergeid := paste(person_block, rleid(newgap), sep = "_"), .(person_block)]

If we merge the long gaps, we end up with 2354 long gaps that can be imputed. DTWBI requires that we create uniform discrete time, aggregating the travel behavior within to that level. We have chosen 15 minutes to strike a balance between the periods being short enough to represent differential travel behavior and long enough to ensure that they encompass multiple location measurements when the minimum gap is set to 5 minutes in length.

setorder(data3, device_id, timestamp)
data3[, hasmissings := any(gap == TRUE), aftermergeid]
full <- copy(data3[hasmissings == FALSE | candtw == TRUE]) # we lose about 1m records and 4000 person_blocks containing not enough information
full[, time_elapsed := timestamp - timestamp[1], aftermergeid]
full[, ftchunk := floor(as.period(time_elapsed)/minutes(15))]
full[, fifteens := sprintf("%04d", ftchunk)]
full[, grouping := paste(aftermergeid, fifteens, sep = "-")]

setorder(full, aftermergeid, ftchunk)
started.at=proc.time()
res <- tdtr(full, col_names = list(entity_id_col = "aftermergeid", timestamp_col = "timestamp", latitude_col
                                   = "latitude", longitude_col = "longitude"), group_col = "grouping", max_error = 30)
cat("Finished in",timetaken(started.at),"\n") ## 13 minutes
saveRDS(res, "./data/tdtrdata_221810.RDS")
segs <- getSegments(res, group = TRUE)
saveRDS(segs, "./data/segs_221810.RDS")

We add in the missing time series entries, setting the distance to NA.

segs <- readRDS("./data/segs_221810.RDS")
setDT(segs)
databyfifteen <- segs[, .(tsstart = min(seg_start_time),
                          tsend = max(seg_end_time),
                          dist = sum(segdist),
                          tsstartlat = seg_start_lat[1],
                          tsendlat = seg_end_lat[.N],
                          tsstartlon = seg_start_lon[1],
                          tsendlon = seg_end_lon[.N]),
                      .(entity_id, .id)]
databyfifteen[, c("person_block", "fifteens") := tstrsplit(.id, split = "-")]
databyfifteen[, ftchunk := as.numeric(fifteens)]
databyfifteen[, tsstart := as_datetime(tsstart)]
databyfifteen[, tsend := as_datetime(tsend)]
databyfifteen[, first_timestamp := (min(tsstart)), entity_id]

missingfifteens <- databyfifteen[, seq(from = 0, to = max(ftchunk))[!seq(from = 0, to = max(ftchunk)) %in% ftchunk], entity_id]
missingfifteens[, `:=`(lastbefore = min(V1) - 1,
                       firstafter = max(V1) + 1), entity_id]
newrows <- databyfifteen[missingfifteens, on = .(entity_id, ftchunk == lastbefore),
                         .(entity_id,
                           .id = paste(entity_id, sprintf("%04d", i.V1), sep = "-"),
                           tsstart = first_timestamp + minutes(i.V1 * 15) ,
                           tsend = first_timestamp + minutes((i.V1 + 1) * 15),
                           dist = NA,
                           tsstartlat = x.tsendlat,
                           tsendlat = NA_real_,
                           tsstartlon = x.tsendlon,
                           tsendlon = NA_real_,
                           person_block,
                           fifteens = sprintf("%04d", i.V1),
                           ftchunk = i.V1,
                           first_timestamp,
                           firstafter)]
newrows[databyfifteen, on = .(entity_id, firstafter == ftchunk), `:=`(tsendlat = i.tsstartlat, tsendlon = i.tsstartlon)]
set(newrows, j = "firstafter", value = NULL)

databyfifteen <- rbindlist(list(
  databyfifteen, newrows
))

setorder(databyfifteen, entity_id, ftchunk)
databyfifteen[, first_timestamp := NULL]
databyfifteen[, time := hms::as_hms(tsstart)]
databyfifteen[, maxftchunk := max(ftchunk), entity_id]

This gives us 4928 complete or imputable slices of varying length. The next step is imputation using DTWBI with the Burst parameters. First we separate out the ‘queries’ which are the sets containing long gap missing data, and we order those by the amount of missingness.

n_imps <- 10
prob.power <- 4
window.size = hours(3)
match.window = 4

impdata <-
    copy(rbindlist(replicate(n_imps, databyfifteen, simplify = FALSE), idcol = "imp"))
impdata[, orig_dist := dist]

queries <- impdata[imp == 1, .SD[any(is.na(dist))], entity_id]
query_order <- queries[, sum(is.na(dist)), entity_id][order(V1, decreasing = FALSE), entity_id]

library(doParallel)
detectCores()

## [1] 16

registerDoParallel(10)

allimps <- foreach(
    this_imp = seq_len(n_imps),
    .combine = rbind,
    .inorder = FALSE,
    .multicombine = TRUE
  ) %dopar% {
    for (rechain_nr in seq_len(2)) {
      cat("Chain number", rechain_nr, "\n", file = "")
      # for (i in seq_along(query_order)) {
      for (i in 2) {
        refs <-
          impdata[imp == this_imp &
                    entity_id != query_order[i], # Select only those for this imputation chain
                  .SD[!any(is.na(dist))],                   # That are complete (includes previously imputed)
                  entity_id]                                     # Grouped by entity_id so we get only those with no NAs when we use any()
        this_query <-
          queries[entity_id == query_order[i]] # Select the query for innermost loop
        miss       <-
          this_query[, is.na(dist)]       # Extract the vector of missingness
        obs        <-
          !miss                           # The inverse observed vector

        # query <- this_query
        impdata <- imputeMidmissing(
          impdata  = impdata,
          this_imp = this_imp,
          query    = this_query,
          obs      = obs,
          refs     = refs,
          n_imps   = 1,
          prob.power   = prob.power,
          window       = window.size,
          match.window = match.window
        )

      } # End query iteration
    }   # End chaining iteration
    impdata[imp == this_imp]
  }     # End %dopar% over n_imps



  allimps[, `:=`(
    window.size = window.size,
    prob.power = prob.power,
    match.window = match.window
  )]
  
allimps[is.na(orig_dist) & !is.na(dist)]