Create weekly categories

wk_cuts = 0:ceiling(max(primary_exposures_idpar$infant_days)/7) * 7
wk_labels = head(wk_cuts, -1)/7
all_exposures_raw$infant_wks = cut(all_exposures_raw$infant_days, include.lowest = T, ordered_result = T,
                                 breaks = wk_cuts, labels = wk_labels)

all_exposures_raw$final_infant_wk = as.numeric(as.character(cut(all_exposures_raw$final_infant_day, 
                                                   include.lowest = T, ordered_result = T,
                                 breaks = 0:100 * 7, labels = F))) - 1

testthat::expect_equal(min(all_exposures_raw$final_infant_wk) ,
                       expected = 0, 
                       info = "check infant_wk rescale")

all_exposures_raw$infant_wks = as.numeric(as.character(all_exposures_raw$infant_wks))

ggplot(arrange(all_exposures_raw, infant_days), aes(y = infant_wks, x = infant_days)) +
  geom_tile()

Past versions of setup-weeks-1.png

Fill in missing weeks

• All interpolation is done using the interpolation of the log viral load (zoo:na.approx).
• map_df was used so that the data is summarized by a refactored infant_wk so complete can be used to find missing weeks for a giving exposure set (which varies by infant and exposure source).

all_exposures = map_df(unique(all_exposures_raw$unique_id), function(uid){
  temp_data = subset(ungroup(all_exposures_raw), unique_id == uid)
  
  # refactor levels for complete()
  temp_data$infant_wks = factor(temp_data$infant_wks, 
                                levels = min(temp_data$infant_wks):max(temp_data$infant_wks)) 

  out = temp_data %>%
  group_by(FamilyID, unique_id, momhiv, virus, idpar, infant_wks, final_infant_wk) %>%
  summarize(
    count = max(count),
    infected = unique(infected),
    infectious_1wk = max(infectious_1wk),
    final_exposure = max(final_exposure)
    ) %>% 
    ungroup() %>%
    complete(infant_wks, nesting(FamilyID, momhiv, virus, idpar, final_infant_wk, 
                               infected, unique_id)) %>%
    arrange(infant_wks) %>%
    mutate(
      interpolate_idpar = if_else(is.na(count), unique(temp_data$idpar), ""),
      infant_wks = as.numeric(as.character(infant_wks)),
      infectious_1wk = na.locf(infectious_1wk), # should only be zero, tested after
      final_exposure = na.locf(final_exposure), # same as above
      count = na.approx(count)
      )
  if(nrow(temp_data) == 1) return(out)

  testthat::expect_equal(n_distinct(diff(out$infant_wks)), expected = 1, 
                         info=paste("Check infant_wks interpolation worked (common interval)", unique(out$unique_id)))
  testthat::expect_equal(unique(diff(out$infant_wks)), expected = 1, 
                         info=paste("Check infant_wks interpolation worked (interval = one)", unique(out$unique_id)))
  
  out
  
}) %>%
  mutate(
    count = if_else(count >= lower_limit - 1, count, 0) # the 1 is a small tolerance factor
  )


testthat::expect_equal(all_exposures %>% group_by(unique_id) %>% 
  summarize(test = sum(infectious_1wk), test2 = sum(final_exposure)) %>%
  filter(test > 1 | test2 > 1) %>% nrow(), expected = 0, 
  info = "Verifying infectious_1wk and final_exposure after interpolation")
  
ggplot(arrange(all_exposures, unique_id), aes(y = infant_wks, x = unique_id)) +
  geom_tile() +
  coord_flip()

Past versions of interpolate-missing-weeks-1.png

all_exposures_wide = all_exposures %>%
  group_by(FamilyID, virus, infant_wks) %>%
  mutate(
    interpolate_idpar = str_trim(str_c(sort(unique(interpolate_idpar)), collapse = " "))
  ) %>%
  ungroup() %>%
  reshape2::dcast(FamilyID + virus + infant_wks + infectious_1wk + final_infant_wk +
                                   infected + momhiv + final_exposure + interpolate_idpar ~ idpar,
                                     data = ., value.var = "count") %>%
  mutate(
    HH = log10(10^M + 10^S),
    HH = if_else(HH <= lower_limit, 0, HH)
    ) %>% 
  ungroup()

exposure_data = all_exposures_wide %>%
  filter(!is.na(S)  & !is.na(M)) %>%
  group_by(FamilyID, virus) %>%
  mutate(
    obs_infected = infected * max(infectious_1wk), 
    final_wk = max(infant_wks),
    outcome_time = ifelse(obs_infected, final_infant_wk, final_wk + 1)
  ) %>%
  ungroup() %>%
  left_join(age_data, by = "FamilyID")


# extra step because empty string patterns are not supported

tmp_chr = "-"
exposure_data_long = exposure_data %>%
  gather(idpar, count, S, M, HH) %>%
  mutate(
    interpolate_idpar_tmp = if_else(interpolate_idpar == "", tmp_chr, interpolate_idpar),
    interpolated = if_else(interpolate_idpar != "" & idpar == "HH", T,
                           str_detect(interpolate_idpar_tmp, idpar))
  ) %>%
  select(-interpolate_idpar_tmp)

testthat::expect_equal(exposure_data %>% group_by(FamilyID, virus) %>% 
  summarize(test = sum(infectious_1wk), test2 = sum(final_exposure)) %>%
  filter(test > 1 | test2 > 1) %>% nrow(), expected = 0, 
  info = "Verifying final exposure has at most one infectious dose per infant")

testthat::expect_equal(exposure_data_long %>% group_by(FamilyID, idpar, virus) %>% 
  summarize(test = sum(infectious_1wk), test2 = sum(final_exposure)) %>%
  filter(test > 1 | test2 > 1) %>% nrow(), expected = 0, 
  info = "Verifying final exposure has at most one infectious dose per infant")

# save the data
if(save_data) {
  write_csv(exposure_data, "data/exposure_data.csv")
  write_csv(exposure_data_long, "data/exposure_data_long.csv")
}

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14.4

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] bindrcpp_0.2.2  forcats_0.3.0   stringr_1.4.0   dplyr_0.7.8    
 [5] purrr_0.3.0     readr_1.3.1     tidyr_0.8.2     tibble_2.0.1   
 [9] ggplot2_3.1.0   tidyverse_1.2.1 zoo_1.8-4      

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.5 xfun_0.4         reshape2_1.4.3   haven_2.0.0     
 [5] lattice_0.20-38  testthat_2.0.1   colorspace_1.4-0 generics_0.0.2  
 [9] htmltools_0.3.6  yaml_2.2.0       rlang_0.3.1      pillar_1.3.1    
[13] withr_2.1.2      glue_1.3.0       modelr_0.1.2     readxl_1.2.0    
[17] bindr_0.1.1      plyr_1.8.4       munsell_0.5.0    gtable_0.2.0    
[21] workflowr_1.2.0  cellranger_1.1.0 rvest_0.3.2      evaluate_0.12   
[25] labeling_0.3     knitr_1.21       broom_0.5.1      Rcpp_1.0.0      
[29] backports_1.1.3  scales_1.0.0     jsonlite_1.6     fs_1.2.6        
[33] hms_0.4.2        digest_0.6.18    stringi_1.2.4    grid_3.5.1      
[37] rprojroot_1.3-2  cli_1.0.1        tools_3.5.1      magrittr_1.5    
[41] lazyeval_0.2.1   crayon_1.3.4     whisker_0.3-2    pkgconfig_2.0.2 
[45] xml2_1.2.0       lubridate_1.7.4  assertthat_0.2.0 rmarkdown_1.8   
[49] httr_1.4.0       rstudioapi_0.9.0 R6_2.3.0         nlme_3.1-137    
[53] git2r_0.24.0     compiler_3.5.1