Create weekly categories

Plot checks that recoding was done right (no overlap on y-axis across steps)

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). This could be done with group_by and nest.

First plot displays extent of left censoring Second plot shows where intepolation occured.

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

  out = temp_data %>%
  group_by(FamilyID, unique_id, momhiv, virus, idpar, 
           infant_wks, first_infant_week, 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, 
                                 first_infant_week, 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.fill(infectious_1wk,  0),
      final_exposure = na.fill(final_exposure,  0),
      count = na.approx(count, rule = 2)
      )
  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)))
  if(any(out$interpolate_idpar != "")){
    testthat::expect_equal(unique(out$infectious_1wk[out$interpolate_idpar != ""]), 
                             expected = 0, 
                           info=paste("Check infectious_1wk interpolation",
                                      unique(out$unique_id)))
    testthat::expect_equal(unique(out$final_exposure[out$interpolate_idpar != ""]), 
                             expected = 0, 
                           info=paste("Check infectious_1wk interpolation",
                                      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")
  

all_exposures %>%
  select(FamilyID, virus, idpar, first_infant_week) %>%
  distinct() %>%
  group_by(virus, first_infant_week, idpar) %>%
  summarize(total = n()) %>%
  ggplot(aes(x = factor(first_infant_week), y = total)) +
  geom_histogram(stat = "identity") +
  geom_text(aes(label = total), vjust = 1) +
  facet_grid(idpar~virus)

Warning: Ignoring unknown parameters: binwidth, bins, pad

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

all_exposures %>%
  group_by(FamilyID, infant_wks, virus) %>%
  summarize(
    interpolate = str_c(interpolate_idpar, collapse = "")
  ) %>%
  arrange(FamilyID, virus) %>%
  ggplot(aes(y = infant_wks, x = FamilyID, fill = factor(interpolate))) +
  geom_tile() +
  scale_fill_manual("", values = c("black", "red", "blue", "gray"),
                    breaks = c("", "S", "MS", "M"),
                    labels = c("", "Interpolate - S", 
                               "Interpolate - M,S", "interpolate - M")) +
  coord_flip() +
  labs(y = "infant weeks post-dob pre-infection") +
  theme_bw() +
  theme(legend.position = "top", axis.text.y = element_text(size = 7)) +
  facet_wrap(~virus, nrow = 2, strip.position = "right", scales = "free_y")

Past versions of interpolate-missing-weeks-2.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.6.0 (2019-04-26)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/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] forcats_0.4.0   stringr_1.4.0   dplyr_0.8.1     purrr_0.3.2    
 [5] readr_1.3.1     tidyr_0.8.3     tibble_2.1.3    ggplot2_3.1.1  
 [9] tidyverse_1.2.1 zoo_1.8-6      

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