---
title: "QMRA (HEV)"
author: "Diego Garrido"
date: "`r format(Sys.time(), '%d %B, %Y')`"
format:
  html:
    code-fold: true
    code-tools: true
execute:
  echo: false
  tidy: true
  warning: false
  message: false
knitr: 
  opts_chunk:
    fig.align: center
bibliography:
  refs.bib
---

```{r}
library(dplyr)
library(ggplot2)
library(cowplot)
library(fitdistrplus)
```

## Data

### Original dataset

From `RiskHEV.xlsx$3-Dades_gomets` and `temperaturesmostreig gomets.xlsx`.

<details>

<summary>Details</summary>

```{r, echo = T}
tb = data.frame(
  spot = c(
    rep(c("llum", "porta", "bossa identificativa (pinça)", "contenidor manec ferro", "contenidor interior"), 4),
    "llum", "porta", "cinta contenidor","int", "ext-contenidor interior"),
  date = as.Date(c(
    rep("04/24/2024", 5), 
    rep("05/08/2024", 5), 
    rep("05/23/2024", 5), 
    rep("06/05/2024", 5), 
    rep("06/24/2024", 3), 
    rep("07/15/2024", 2)
  ), format = "%m/%d/%Y"),
  temp = c(
    rep(8.89, 5), 
    rep(12.6, 5), 
    rep(13.8, 5), 
    rep(17.06, 5), 
    rep(19.21, 3), 
    rep(22.2, 2)
  ),
  hev = c(
    4.77E+02, 3.14E+02, NA, 2.00E+03, 3.51E+04,
    NA, 4.90E-02, NA, 1.70E+00, NA,
    2.97E-01, NA, NA, 1.74E+02, 2.31E+02,
    1.29E+01, 4.17E-01, 1.20E-01, 1.34E+00, 6.40E-01,
    NA, 3.51E-01, 4.30E+02,
    NA, NA
  ))
```

```{r}
knitr::kable(format(tb, scientific = F))
```

</details>

### Data preprocessing

-   Remove spot: "*bossa identificativa (pinça)*"
-   Rename spots in english
-   Make groups: `S1`: "*door*" + "*light*"; `S2`: "*container_lid*" + "*grip_handle*"
-   Add helper columns: `group`, `datenum` (days from 2024-04-01), `dategroup`

```{r, echo = T}
tb %<>% 
  filter(spot != "bossa identificativa (pinça)") %>%
  mutate(
    spot = case_when(
      spot == "porta" ~ "door",
      spot == "llum" ~ "light",
      spot %in% c("cinta contenidor", "ext-contenidor interior", "contenidor interior", "int") ~ "container_lid",
      spot == "contenidor manec ferro" ~ "grip_handle",
      TRUE ~ spot
    ),
    group= case_when(
      spot %in% c("door", "light") ~ "S1",
      spot %in% c("container_lid", "grip_handle") ~ "S2"
    ),
    datenum = as.numeric(date - as.Date("2024-04-24")),
    dategroup = ifelse(datenum <= 29, "april-may", "june-july")
  )
```

```{r}
knitr::kable(format(tb, scientific = F))
```

### Descriptive stats

```{r, echo = T}
set.seed(123)
# ggplot(tb[!is.na(tb$hev), ], aes(y = log10(hev), x = spot)) +
#   geom_boxplot(col = "gray", outlier.size = 0) +
#   geom_jitter(width = 0.2) +
#   theme_classic(base_size = 18) +
#   theme(axis.text = element_text(angle = 45, hjust = 1))

tbs <- tb %>%
  filter(!is.na(hev)) %>%
  mutate(spot = forcats::fct_reorder(spot, hev, .fun = mean))

ggplot(tbs[!is.na(tbs$hev), ], aes(x = spot, y = log10(hev))) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(
    fun.data = mean_cl_normal, 
    geom = "errorbar", 
    width = 0.2
  ) +
  geom_jitter(width = 0.2, alpha = 0.6) +
  theme_classic(base_size = 18) +
  theme(axis.text = element_text(angle = 45, hjust = 1))

ggplot(tb[!is.na(tb$hev), ], aes(y = log10(hev), x = group)) +
  geom_boxplot(col = "gray", outlier.size = 0) +
  geom_jitter(width = 0.2) +
  theme_classic(base_size = 18) 

ggplot(tb[!is.na(tb$hev), ], aes(y = log10(hev), x = dategroup)) +
  geom_boxplot(col = "gray", outlier.size = 0) +
  geom_jitter(width = 0.2) +
  theme_classic(base_size = 18)
```

```{r, echo = T}
# fisher.test(with(tb, table(dategroup, is.na(hev))))

# ggplot(tb[!is.na(tb$hev), ], aes(x = datenum, y = log10(hev))) +
#   geom_point(alpha = 0.6) +                         
#   geom_smooth(method = "lm", se = TRUE,             
#               linewidth = 1, colour = "steelblue",
#               fill = "lightblue") +
#   theme_classic(base_size = 18) 
# 
# lm.fit <- lm(log10(hev)~datenum, data = tb)
# # plot(lm.fit)
# summary(lm.fit)

ggplot(tb[!is.na(tb$hev), ], aes(x = temp, y = log10(hev))) +
  geom_point(alpha = 0.6) +                         
  geom_smooth(method = "lm", se = TRUE,             
              linewidth = 1, colour = "steelblue",
              fill = "lightblue") +
  theme_classic(base_size = 18) 

lm.fit <- lm(log10(hev)~datenum, data = tb)
# plot(lm.fit)
summary(lm.fit)
```

-   `NA`s substituted by 0:

```{r, echo = T}
tb2 = tb
tb2$hev[is.na(tb2$hev)] = 1 # then log(1) = 0

# ggplot(tb2, aes(x = datenum, y = log10(hev))) +
#   geom_point(alpha = 0.6) +                         
#   geom_smooth(method = "lm", se = TRUE,
#               linewidth = 1, colour = "steelblue",
#               fill = "lightblue") +
#   theme_classic(base_size = 18)
# 
# lm.fit2 <- lm(log10(hev)~datenum, data = tb2)
# # plot(lm.fit)
# summary(lm.fit2)

ggplot(tb2, aes(x = temp, y = log10(hev))) +
  geom_point(alpha = 0.6) +                         
  geom_smooth(method = "lm", se = TRUE,
              linewidth = 1, colour = "steelblue",
              fill = "lightblue") +
  theme_classic(base_size = 18)

lm.fit2 <- lm(log10(hev)~datenum, data = tb2)
# plot(lm.fit)
summary(lm.fit2)
```

## Exposure assesment

### Surface model

$$d = C_{HEV} \times \frac{f_{conv}}{r_{eff}} \times f_{conv|surf} \times TE_{surf|hand} \times TE_{hand|mouth} \times f_{cont|mouth} \times n_{hours}$$

-   Concentration of HEV (surface, CG/cm$^2$), $C_{HEV} \sim$ Lognormal$(\mu, \sigma)$. We consider left-censored data.

<div style="margin-left: 2em;">
- Group `S1`:

<div style="margin-left: 4em;">
```{r, echo = T}
tbs1 = subset(tb, group == "S1")
tbs1cens = data.frame(left = tbs1$hev, right = tbs1$hev)
tbs1cens$right[is.na(tbs1cens$right)] = min(tbs1cens$right, na.rm = T)
fit1 = fitdistcens(tbs1cens, distr = "lnorm")
fit1
```

<mark>**Q1**. Detection limit is 4.56 GC/cm\^2. How is it possible that there are values below?</mark>
</div>

- Group `S2`:

<div style="margin-left: 4em;">
```{r, echo = T}
tbs2 = subset(tb, group == "S2")
tbs2cens = data.frame(left = tbs2$hev, right = tbs2$hev)
tbs2cens$right[is.na(tbs2cens$right)] = min(tbs2cens$right, na.rm = T)
fit2 = fitdistcens(tbs2cens, distr = "lnorm")
fit2
```
</div>
</div>


<div style="margin-left: 2em;">
<mark> **Q2**. Maybe we should do $U(0, max)$ for each spot. </mark>
</div>

-   Conversion factor (TCID50/GC), $f_{conv} = \frac{1}{700}$ [@McBride2013].
-   Recovery efficiency (%), $r_{eff} = 7.3$ [@Itarte2024].
-   Frequency of hand-to-surface contacts per work day ("*jornada laboral*") (events/day), $f_{conv|surf} \sim$ Uniform $(a, b)$ (from `RiskHEV.xlsx$3-Dades_gomets`):

<div style="margin-left: 2em;">
```{r}
# removed "bossa identificativa (pinça)"
contact = data.frame(spot = c("light", 
                              "door", 
                              "container_lid",
                              "grip_handle"),
                     a = c(6, 2, 0, 0),
                     b = c(8, 20, 4, 4))
knitr::kable(contact)
```
To integrate different spots within `S1` we may use a mixture between both distributions (also for HEV concentrations).

```{r}
which_dist <- sample(c(1, 2), size = 1, prob = c(0.5, 0.5)) # shared for HEV concentration
smpl = floor(ifelse(which_dist == 1, sample(6:8, size = 1), sample(2:20, size = 1)))
```
</diV>

-   Viral transfer efficiency from surface to bare hands (%), $TE_{surf|hand} \sim$ Normal $(\mu = 37, \sigma = 16)$ [@Lopez2013].

-   Viral transfer efficiency from hands to mouth (%), $TE_{hand|mouth} = 34$ [@Rusin2002].

-   Frequency of hand-to-mouth contacts (events/hour), $f_{cont|mouth}$. The distribution can be learnt from the data of @Wilson2021, available at https://doi.org/10.5518/798 ([data_face_contact.csv](https://archive.researchdata.leeds.ac.uk/680/1/data_face_contact.csv)).

<div style="margin-left: 2em;">

```{r, echo = T}
data_face_contact <- read.csv(url("https://archive.researchdata.leeds.ac.uk/680/1/data_face_contact.csv"))
f_contmouth = subset(data_face_contact, 
                     eatonly == "Non-eating" & 
                     age == "adult" & 
                     contactarea == "mouth")$frequency
pmf_contmouth <- table(f_contmouth) / length(f_contmouth)
# sum(pmf_contmouth)

# Non-parametric sampling
set.seed(123)
sampl <- sample(as.numeric(names(pmf_contmouth)), size = 1e6, replace = TRUE, prob = pmf_contmouth)
```

Distribution with $\mu\simeq$ `r round(mean(sampl),2)` and $\sigma \simeq$ `r round(sd(sampl),2)`.

</div>

- We thought of not including $n_{hours}$ as we consider daily doses. However, $f_{cont|mouth}$ is given in events/hour. <mark> **Q3.** So we need to define $n_{hours}$. We can use $U(4,6)$ as in @Itarte2024 </mark>.

## Dose-response model (Beta-Poisson)

### Exact model

$$ P_{inf}(d|\alpha,\beta) = 1 - {}_1F_1\left( \alpha, \alpha + \beta, -d \right) $$

### Approximation

$$ P_{inf}(d|\alpha,\beta) = 1 - \left(1 + \frac{d}{\beta} \right)^{-\alpha} $$ This approximation is valid when $\alpha << \beta$ and $\beta >> 1$ [@Haas2014].

### Reparametrization based on $N_{50}$

$$ P_{inf}(d |\alpha, N_{50}) = 1 - \left ( 1 + (2^{1/\alpha}-1)\cdot \frac{d}{N_{50}}\right ) ^{-\alpha}$$ where $\beta = \frac{N_{50}}{2^{1/\alpha}-1}$. A good description of the Beta-Poisson model is provided in @Teunis2000.

## Parameter estimation

Parameters derived for HEV from chimpanzee data via MLE in @Ruchusatsawat2022:

-   $\alpha = 216.89$
-   $N_{50} = 3.03 \cdot 10^7$ (intravenous inoculation)

95% CI for both parameters are shown in Figs. 2-3, but data (to the best of my knowledge) is not provided.

@Wilkins2025 refers to @Muller2017 and @AFFSA2009 and suggests to use $N_{50} = 10^{5.5}$ instead (based on consumption of pork-derived products). It seems more realistic to use this value in our scenario.

<mark>**Q4**. In @Itarte2024 we have an equation (eq. 4) to compute $P_{ill}$ (probability of illness) from $P_{inf}$ (probability of infection), but for HEV we do not have such equation (we lack the parameter estimates). In @Ruchusatsawat2022, they use a similar strategy (eqs. 6 and 11). We may be able to use this, but I am not sure whether it would be equivalent. Discuss in person.</mark>

## QMRA test (1)

Considering uniform distributions for $C_{HEV}$.

```{r, echo = T}
set.seed(123)
nsim = 1e4
res = list()
spots = c("door", "light", "grip_handle", "container_lid")
for (where in spots){
  C_HEV = runif(n = nsim, ifelse(is.na(min(subset(tb, spot==where)$hev)),
                                 0, min(subset(tb, spot==where)$hev)), 
                          max(subset(tb, spot==where)$hev, na.rm = T))
  f_conv = 1/700
  r_eff = 7.3
  range_vals <- switch(where,
                       door = 6:8,
                       light = 2:20,
                       grip_handle = 0:4,
                       container_lid = 0:4
  )
  f_conv_surf = sample(range_vals, replace = T, size = nsim)
  TE_surf_hand = rnorm(n = nsim, mean = 37, sd = 6)
  TE_hand_mouth = 34
  f_cont_mouth = sample(as.numeric(names(pmf_contmouth)), 
                        replace = T, size = nsim, prob = pmf_contmouth)
  n_hours = sample(4:6, size = nsim, replace = T)
  
  d = C_HEV*f_conv/r_eff*f_conv_surf*TE_surf_hand*TE_hand_mouth*f_cont_mouth*n_hours
  
  alpha = 216.89
  N50 = 10^5.5
  P_inf = 1 - (1 + (2^(1/alpha) - 1)*d/N50)^(-alpha)
  res[[where]] = P_inf
}

# For plotting
df = data.frame(
  spot = rep(spots, times = lengths(res)),
  probs = unlist(res, use.names = FALSE)
)
medians <- tapply(df$probs, df$spot, median)
df$spot <- factor(df$spot, levels = names(sort(medians)))

pl = ggplot(df, aes(x = spot, y = probs)) +
  geom_boxplot() +
  theme_minimal() +
  labs(y = expression(P[inf]))

plz = ggplot(df, aes(x = spot, y = probs)) +
  geom_boxplot() +
  theme_minimal() +
  coord_cartesian(ylim = c(0,0.05)) + 
  labs(y = expression(P[inf]))

plot_grid(pl, plz)
```

## References

::: {#refs}
:::