---
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
bibliography:
  refs.bib
---

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

## Data

### Original dataset

From `RiskHEV.xlsx$3-Dades_gomets`.

<mark>**Q1**. Confirm that "*cinta contenedor*" and "*ext-contenidor interior*" are actually "*contenidor manec ferro*"; and that "*int*" is "*contenidor interior*"</mark>

```{r, echo = T}
tb = data.frame(
  spot = c(
    rep(c("llum", "porta", "bossa identificativa (pinça)", "contenidor manec ferro", "contenidor interior"), 4),
    "llum", "porta", "contenidor manec ferro","contenidor interior", "contenidor manec ferro"),
  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"),
  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))
```

### Data preprocessing

- Make groups: `S1`: "*porta*" + "*llum*"; `S2`: "*contenidor[...]*"
- Remove spot: "*bossa identificativa (pinça)*"
- Add helper columns: `group`, `datenum` (days from 2024-04-01), `dategroup`

```{r, echo = T}
tb %<>% 
  filter(spot != "bossa identificativa (pinça)") %>%
  mutate(
    group= case_when(
      spot %in% c("porta", "llum") ~ "S1",
      grepl("contenidor", spot) ~ "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, 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))
ggplot(tb, 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, aes(y = log10(hev), x = group)) +
  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))

ggplot(tb, aes(y = log10(hev), x = dategroup)) +
  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))
```
```{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)
```

- `NA`s to min detection:
```{r, echo = T}
tb2 = tb
tb2$hev[is.na(tb2$hev)] = min(tb2$hev, na.rm = T)

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)
```

<mark>**Q2**. It would be interesting to add temperature data to this descriptive analysis. </mark>

## 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}$$

We do not include $n_{hours}$ as we consider daily doses.

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

-Group `S1`:
```{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>**Q3**. Detection limit (common or specific to door and light?)</mark>

-Group `S2`:
```{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
```

- Conversion factor (TCID50/GC), $f_{conv} = \frac{1}{700}$ [@McBride2013].
- Recovert 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`):

```{r}
# removed "bossa identificativa (pinça)"
contact = data.frame(spot = c("llum", 
                              "porta", 
                              "contenidor manec ferro",
                              "contenidor interior"),
                     a = c(6, 2, 0, 0),
                     b = c(8, 20, 4, 4))
knitr::kable(contact)
```
<mark>**Q4**. How to integrate different spots within `S1`?</mark>

- Viral transfer efficiency from surface to bare hands (%), $TE_{surf|hand} \sim$ Normal $(\mu = 37, \sigma = 16)$ [@Lopez2013]. <mark>**Q5**. Make sure we can use this distribution/parameters here too.</mark>

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

- Frequency of hand-to-mouth contacts (events/hour), $f_{cont|mouth}$ [@Wilson2021]. Here we have parameters $\mu \simeq 2.8, \sigma \simeq 2.7$. <mark>**Q6**. But which distribution? Normal or Lognormal? in addition, this is events/hour. Shouldn't we consider events/day? </mark>

### Splashing model

<mark>**Q7**. I assume this model is only for `S2`. Not sure about which equation should I employ here. In @Itarte2024 (see section 2.4.2.1, eq. 1) and in `RiskHEV.xlsx$2-Model_splashing` there are different models. Not sure either about $V_{leachate}$ (units do not match with $C_{HEV}$: ml $vs$ GC/cm$^2$.</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).

<mark>**Q8**. Discuss which value of $N_{50}$ should we use</mark>. 

<mark>**Q9**. 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. We should discuss it.</mark>

## References

<div id="refs"></div>