Model setup

The probability of being uninfected after a single, weekly exposure can written as following:

\[ s(i) = exp(-\beta_0 - \beta_{E}E_i) \] or \[ s(i) = exp(-\lambda(i)) \]

Instead of continuous time, we consider discretized time denoted by \(i \in \{1, ..., n+1\}\) for n+1 weeks of exposures (with n survived exposures). The likelihood follows for a single, infected participant (note that in discrete time we don’t apply an instantaneous hazard but use 1 - S(n) for the infectious week):

\[ L_j(n_j+1) = \prod_{i = 1}^{n_j} s_j(i) * (1-s_j(n_j+1))\]

for the j\(^{th}\) participant with a unique set of total exposures. We next setup the log-likelihood for the population (m participants) and use \(\Delta\) to denote an observed infection.

\[ \sum_{j = 1}^{m} log L_j(n_j) = \sum_{j = 1}^{m}\sum_{i = 1}^{n_j} log(s_j(i)) + \sum_{j = 1}^{m} \Delta_j * log(1-s_j(n_j+1))\]

The following assumptions are used: at-risk individuals are independent and the risk associated with a weekly exposures is unique from other exposures (i.e., non-infectious exposure weeks are exchangeable). From this formulation, we can find the maximum likelihood estimators for the parameters by minimizing the negative log likelihood. Both parameters are solved numerically.

The null model has the following form for weekly risk:

\[ s_0(i) = exp(-\beta_0) \]

I.e., risk is constant and not affected by household exposure. The null model has a simplified log-likelihood

\[ \sum_{j = 1}^{m} log L_j(n_j) = -\beta_0 \sum_{j = 1}^{m}n_j + log(1-exp(-\beta_0)) \sum_{j = 1}^{m} \Delta_j \]

with closed-form solution for the null risk

\[\hat{\beta_0} = -log(1 - \frac{\sum_{j = 1}^{m} \Delta_j}{\sum_{j = 1}^{m}n_j})\]

Run models

• Null model is just constant risk model.
• First set of fitted models are the marginal models
• Second set is the combined M + S mode

# note negative loglik is calculated
null_risk_mod = model_data %>% 
  group_by(virus, FamilyID) %>%
  summarize(
    infected = max(infectious_1wk),
    surv_weeks = max(c(0, infant_wks[which(infectious_1wk == 0)]))
  ) %>%
  group_by(virus) %>%
  summarize(
    null_beta = -log(1-sum(infected)/sum(surv_weeks)),
    null_loglik = null_beta * sum(surv_weeks) - log(1-exp(-null_beta)) * sum(infected)
  )

risk_mod = model_data_long %>%
  group_by(virus, idpar, cohort) %>%
  nest() %>%
  mutate(
    likdat = map(data, create_likdat),
    total = map_dbl(data, ~n_distinct(.x$FamilyID)), 
    total_infected = map_dbl(data, ~sum(.x$infectious_1wk)),
    fit_modNM = map(likdat,  ~optim(c(-12, -12), surv_logLik, likdat = .x)),
    fit_modBFGS = map(likdat,  ~optim(c(-12, -12), surv_logLik, likdat = .x,  method = "BFGS")),
    #fit_modE = map(likdat, ~optimize(surv_logLikE, c(-20, 0), likdat = .x)),
    fit_resNM = map(fit_modNM, tidy_fits),
    fit_resBFGS = map(fit_modBFGS, tidy_fits)
  ) 

# using both exposures
risk_mod_both = model_data %>%
  group_by(virus, cohort) %>%
  nest() %>%
  mutate(
    likdat = map(data, create_likdat2),
    total = map_dbl(data, ~n_distinct(.x$FamilyID)), 
    total_infected = map_dbl(data, ~sum(.x$infectious_1wk)),
    fit_modNM = map(likdat,  ~optim(c(-1, -3, -3), surv_logLik_2E, likdat = .x)),
    fit_modBFGS = map(likdat,  ~optim(c(-1, -3, -3), surv_logLik_2E, likdat = .x,  
                                      method = "BFGS")),
    fit_resNM = map(fit_modNM, tidy_fits2),
    fit_resBFGS = map(fit_modBFGS, tidy_fits2)
  ) 

For marginal models, BFGS finds equivalent or better loglik compared to NM. For combined model, the NM algorithm finds the equivalent or better.

risk_mod %>% 
  rename(NM = fit_resNM, BFGS = fit_resBFGS) %>%
  unnest(BFGS, .sep = "_") %>%
  unnest(NM, .sep = "_") %>% 
  select(virus, idpar, contains("beta"), contains("loglik")) %>%
  group_by(virus, idpar) %>%
  select(virus, idpar, BFGS_beta0, NM_beta0, BFGS_betaE, NM_betaE, BFGS_loglik, NM_loglik) %>%
  mutate_if(is.numeric, format, digits = 3) %>%
  kable(caption = "Comparing optimization results from marginal models (BFGS vs. Nelder-Mead (NM))", 
        digits = 3) %>%
  kable_styling(full_width = F)

`mutate_if()` ignored the following grouping variables:
Columns `virus`, `idpar`

risk_mod_both %>% 
  rename(NM = fit_resNM, BFGS = fit_resBFGS) %>%
  unnest(BFGS, .sep = "_") %>%
  unnest(NM, .sep = "_") %>% 
  select(virus, contains("beta"), contains("loglik"), -contains("betaE")) %>%
  group_by(virus) %>%
  select(virus, BFGS_beta0, NM_beta0, 
         BFGS_betaM, NM_betaM, 
         BFGS_betaS, NM_betaS, 
         BFGS_loglik, NM_loglik) %>%
  mutate_if(is.numeric, format, digits = 3) %>%
  kable(caption = "Comparing optimization results from full model (BFGS vs. Nelder-Mead (NM))",
        digits = 3) %>%
  kable_styling(full_width = F)

`mutate_if()` ignored the following grouping variables:
Column `virus`

Model results

marginal_results = risk_mod %>%
  left_join(null_risk_mod, by = c("virus")) %>%
  unnest(fit_resBFGS) %>%
  mutate(
    #null_betaE = map_dbl(fit_modE, ~exp(`$`(.x, "minimum"))),
    #betaE_loglik = map_dbl(fit_modE, ~`$`(.x, "objective")),
    LLR_stat = 2 * (null_loglik - loglik),
    pvalue = pchisq(LLR_stat, 1, lower.tail = F)
  ) %>%
  select(-contains("dat"), -contains("fit")) %>%
  arrange(virus)

risk_mod_loglik = marginal_results %>%
  subset(idpar != "HH") %>%
  select(virus, idpar, loglik, betaE) %>%
  gather(outcome, value, loglik, betaE) %>%
  unite(mod_res, outcome, idpar) %>%
  spread(mod_res, value)

combined_results = risk_mod_both %>%
  group_by(virus) %>%
  nest() %>%
  mutate(
    fit_res = if_else(virus == "CMV", map(data, unnest, fit_resBFGS),
                      map(data, unnest, fit_resNM))
  ) %>%
  unnest(fit_res) %>%
  select(-contains("fit"), -contains("data"), -contains("dat")) %>%
  left_join(risk_mod_loglik, by = c("virus")) %>%
  left_join(null_risk_mod, by = c("virus")) %>%
  mutate(
    LLR_stat_overall = 2 * (null_loglik - loglik),
    pvalue_overall = pchisq(LLR_stat_overall, 2, lower.tail = F),
    LLR_statM = 2 * (loglik_S - loglik),
    pvalueM = pchisq(LLR_statM, 1, lower.tail = F),
    LLR_statS = 2 * (loglik_M - loglik),
    pvalueS = pchisq(LLR_statS, 1, lower.tail = F)
    ) %>%
  arrange(virus)

options(knitr.kable.NA = '')

full_results = marginal_results %>%
  gather(parameter, estimate, null_beta, beta0, betaE) %>%
  bind_rows(
    combined_results %>% gather(parameter, estimate, beta0, betaM, betaS) %>%
      mutate(
        idpar = "CM",
        pvalue = if_else(parameter == "betaM", pvalue_overall, NA_real_),
      )
  ) %>%
  mutate(
    pvalue = if_else(parameter %in% c("null_beta", "beta0"), NA_real_, pvalue),
    idpar = if_else(parameter == "null_beta", "Constant", idpar),
    model = factor(idpar, 
                   levels = c("Constant", "S", "M", "HH", "CM"),
                   labels = c("Constant risk", "Secondary child", "Mother",
                              "Household sum", "Combined model")
    ),
    parameter = case_when(
      parameter == "null_beta" ~ "beta0",
      parameter == "betaE" ~ str_c("beta", idpar),
      TRUE ~ parameter
    ),
    estimate = if_else(estimate < 1e-25, 0, estimate)
  ) %>%
  select(virus, model, parameter, estimate, pvalue) %>%
  distinct() 


full_results %>%
  arrange(virus, model) %>%
  ungroup() %>%
  mutate(
    pvalue = clean_pvalues(pvalue, sig_alpha = 0)
  ) %>%
  mutate_if(is.numeric, format, digits = 3) %>%
  kable(digits = 3) %>%
  kable_styling(full_width = F) %>%
  collapse_rows(1:2, valign = "top")

id_calc = function(b0, bE, prob){
  (-log(1-prob) - b0)/(bE)
}

beta0_ests = full_results %>%
  subset(model != "Constant risk" & parameter == "beta0") %>%
  select(virus, model, estimate) %>%
  rename(beta0 = estimate)

full_results %>%
  subset(parameter != "beta0" & model != "Household sum" & virus == "HHV-6") %>%
  select(-pvalue) %>%
  group_by(virus, model) %>%
  left_join(beta0_ests, by = c("virus", "model")) %>%
  mutate(
    exposure_source = factor(case_when(
      parameter == "betaS" ~ "Secondary Child",
      parameter == "betaM" ~ "Mother"
    ), levels = c("Secondary Child", "Mother")),
    ID25 = id_calc(beta0, estimate, 0.25),
    ID50 = id_calc(beta0, estimate, 0.5),
    ID80 = id_calc(beta0, estimate, 0.8)
  ) %>%
  ungroup() %>%
  mutate_if(is.numeric, log10) %>%
  mutate(constant_risk = 100*(1-exp(-10^beta0))) %>%
  select(-estimate, -beta0, -parameter) %>%
  arrange(virus, model, exposure_source) %>%
  kable(digits = 2) %>%
  kable_styling(full_width = F) %>%
  collapse_rows(1:2, valign = "top") %>%
  add_header_above(c(" " = 3, "Infectious dose (ID)" = 3, ""))

exposure_max = exposure_data_long %>%
  subset(idpar != "HH") %>%
  group_by(virus, idpar) %>%
  summarize(
    max_exposure = max(exposure)
  )

risk_data = exposure_data_long %>%
  subset(idpar != "HH") %>%
  mutate(
    exposure_cat = floor(count) + 0.5
  ) %>%
  group_by(idpar, virus, exposure_cat) %>%
  summarize(
    total_exposures = n(),
    total_infected = sum(infectious_1wk),
    risk = mean(infectious_1wk)
  )

risk_grid = exposure_data %>%
  mutate(
    exposure_S = floor(S) + 0.5,
    exposure_M = floor(M) + 0.5
  ) %>%
  group_by(virus, exposure_S, exposure_M) %>%
  summarize(
    total_exposures = n(),
    total_infected = sum(infectious_1wk),
    risk = mean(infectious_1wk)
  )

# risk_prediction = marginal_results %>%
#   group_by(virus, idpar) %>%
#   left_join(exposure_max, by = c("virus", "idpar")) %>%
#   nest() %>%
#   mutate(
#     pred_res = map(data, ~with(.x, tibble(
#       max_exposure = log10(max_exposure),
#       exposure = seq(0, 5, length = 100),
#       risk_full = 1 - exp(-beta0 - betaE * 10^exposure)
#       #riskE = 1 - exp(-null_betaE * 10^exposure)
#     )))
#   ) %>%
#   unnest(pred_res) %>%
#   gather(risk_est, risk, risk_full)

risk_prediction_both = combined_results %>%
  gather(parameter, est, betaM, betaS) %>%
  mutate(idpar = str_remove(parameter, "beta")) %>%
  left_join(exposure_max, by = c("virus", "idpar")) %>%
  group_by(virus, idpar) %>%
  nest() %>%
  mutate(
    pred_res = map(data, ~with(.x, tibble(
      exposure = seq(0, log10(max_exposure), length = 100),
      risk = 1 - exp(-beta0 - est * 10^exposure)
    )))
  ) %>%
  unnest(pred_res) 

dr_plot = risk_data  %>%
  subset(virus == "HHV-6") %>%
  ggplot(aes(x = exposure_cat, y = risk)) +
  geom_bar(stat = "identity", fill = "grey50") +
  geom_line(data = subset(risk_prediction_both, virus == "HHV-6"),
            aes(x = exposure), size = 1.25) +
  scale_y_continuous("Estimated weekly risk of HHV-6 infection", breaks = 0:4/4) +
  scale_x_continuous(expression(paste("Log"[10], " HHV-6 VL (DNA copies/swab)")),
                     breaks = 0:6+0.5,
                     labels = paste(0:6, 1:7, sep = "-")) +
  facet_grid(.~idpar, labeller = as_labeller(c('S' = "Secondary Child", 'M' = "Mother"))) +
  theme(legend.position = "top")

exposure_heat = risk_grid %>%
  subset(virus == "HHV-6") %>%
  ggplot(aes(x = exposure_S, y = exposure_M, fill = total_exposures, 
             label = paste(total_infected, total_exposures, sep = "/"))) +
  geom_tile() +
  #geom_label(label.padding = unit(0.1, "lines"), label.size = 0, fill = "white") +
  geom_text(fontface = "bold", colour = "white", size = 3.5) +
  viridis::scale_fill_viridis("Total exposures", option = "viridis") +
  scale_y_continuous(expression(paste("Mother log"[10], " VL (DNA copies/swab)")), 
                     breaks = 0:6+0.5,
                     labels = paste(0:6, 1:7, sep = "-")) +
  scale_x_continuous(expression(paste("Secondary child log"[10], " VL (DNA copies/swab)")),
                     breaks = 0:6+0.5,
                     labels = paste(0:6, 1:7, sep = "-")) +
  theme_classic() +
  theme(legend.position = "bottom",
        legend.title = element_text(size = 8),
        legend.text = element_text(size = 8),
        legend.key.width = unit(0.75, "cm"))

  
risk_heat = crossing(
  virus = "HHV-6",
  exposure_S = seq(0, 
                   log10(subset(exposure_max, virus == "HHV-6" & idpar == "S")$max_exposure)+0.1,
                   length = 100),
  exposure_M = seq(0, 
                   log10(subset(exposure_max, virus == "HHV-6" & idpar == "M")$max_exposure)+0.1,
                   length = 100)
) %>%
  left_join(combined_results, by = "virus") %>%
  mutate(
    risk =  1 - exp(-beta0 - betaS * 10^exposure_S - betaM * 10^exposure_M)
  ) %>%
  ggplot(aes(x = exposure_S, y = exposure_M, fill = pmin(risk, 0.25))) +
  geom_tile() +
  scale_y_continuous(expression(paste("Mother log"[10], " VL (DNA copies/swab)")), 
                     breaks = 1:7) +
  scale_x_continuous(expression(paste("Secondary child log"[10], " VL (DNA copies/swab)")),
                     breaks = 1:7) +
  scale_fill_distiller("Weekly infection prob.", palette = "RdYlBu",  breaks= c(0:5/20),
                        labels = c(0:4/20, " >0.25")) +
  geom_text(data = subset(risk_grid, virus == "HHV-6"), aes(label = round(risk, 2)), 
            fontface = "bold", colour = "black") +
  theme_classic() +
  theme(legend.position = "bottom",
        legend.title = element_text(size = 8),
        legend.text = element_text(size = 8),
        legend.key.width = unit(0.75, "cm"))

plot_grid(dr_plot,
  plot_grid(exposure_heat, risk_heat, nrow = 1, labels = c("B", "C"), vjust = 1),
  nrow = 2, labels = "A")