---
title: "HIV self-testing positivity rate and linkage to confirmatory testing and care: a telephone survey in Côte d’Ivoire, Mali and Senegal"
author: "Arsène Kouassi Kra et al"
date: today
toc: true
format:
  html:
    code-tools: true
    code-fold: true
    self-contained: true
    toc: false
    dpi: 300
    fig-responsive: true
execute:
  warning: false
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, message=FALSE}
# Loading packages
library(tidyverse)
library(gtsummary)
library(labelled)
library(scales)
library(ggrepel)
library(RColorBrewer)
library(utils)
library(patchwork)

# Data import
data <- readr::read_csv("data.csv", show_col_types = FALSE)
```

## Table 1. Self-reported test result and the reported number of lines among participants of phase 1, and positivity rates according to different hypotheses

```{r fig.height=6, fig.width=12, message=FALSE,warning=FALSE}
# Creating "test_result_number_lines_reported" variable
data <- data |>
  mutate(
    test_result_number_lines_reported =
      case_when(
        (HIVST_reported_lines == "2 lines" & HIVST_reported_result == "reactive") ~ "2 lines / reactive",
        (HIVST_reported_lines == "2 lines" & HIVST_reported_result == "non reactive") ~ "2 lines / non-reactive",
        (HIVST_reported_lines == "2 lines" & (HIVST_reported_result == "DK" | HIVST_reported_result == "R")) ~ "2 lines / DK-R",
        (HIVST_reported_lines == "1 line" & HIVST_reported_result == "reactive") ~ "1 line / reactive",
        (HIVST_reported_lines == "1 line" & HIVST_reported_result == "non reactive") ~ "1 line / non-reactive",
        (HIVST_reported_lines == "1 line" & (HIVST_reported_result == "DK" | HIVST_reported_result == "R")) ~ "1 line / DK-R", (HIVST_reported_lines == "0 line" & HIVST_reported_result == "non reactive") ~ "0 line / non-reactive", (HIVST_reported_lines == "0 line" & HIVST_reported_result == "reactive") ~ "0 line / reactive", (HIVST_reported_lines == "0 line" & (HIVST_reported_result == "DK" | HIVST_reported_result == "R")) ~ "0 line / DK-R", ((HIVST_reported_lines == "0 line" | HIVST_reported_lines == "1 line") & HIVST_reported_result == "invalid") ~ "0-1 line / invalid",
        (HIVST_reported_lines == "2 lines" & HIVST_reported_result == "invalid") ~ "2 lines / invalid",
        ((HIVST_reported_lines == "DK" | HIVST_reported_lines == "R") & HIVST_reported_result == "invalid") ~ "DK-R / invalid",
        ((HIVST_reported_lines == "DK" | HIVST_reported_lines == "R") & HIVST_reported_result == "reactive") ~ "DK-R / reactive",
        ((HIVST_reported_lines == "DK" | HIVST_reported_lines == "R") & HIVST_reported_result == "non reactive") ~ "DK-R / non-reactive",
        ((HIVST_reported_lines == "DK" | HIVST_reported_lines == "R") & (HIVST_reported_result == "DK" | HIVST_reported_result == "R")) ~ "DK-R / DK-R"
      )
  ) |>
  set_variable_labels(test_result_number_lines_reported = "Reported number of lines / self-interpreted HIVST result")

# Reorganizing the labels of the variable "test_result_number_lines_reported"
data$test_result_number_lines_reported <- data$test_result_number_lines_reported |>
  fct_relevel(
    "2 lines / reactive",
    "1 line / non-reactive",
    "0-1 line / invalid",
    "1 line / reactive",
    "2 lines / non-reactive", "0 line / non-reactive",
    "0 line / DK-R", "1 line / DK-R", "2 lines / DK-R", "DK-R / reactive",
    "DK-R / DK-R", "DK-R / non-reactive"
  )

# Generating Table 1
tbl1 <-
  data |>
  tbl_summary(
    include = test_result_number_lines_reported,
    label = test_result_number_lines_reported ~ ""
  ) |>
  modify_header(label = "**Phase 1 participants**") |>
  as_gt() |>
  gt::tab_row_group(
    label = gt::md("*Partial reponse (P)*"),
    rows = 8:13
  ) |>
  gt::tab_row_group(
    label = gt::md("*Inconsistant reponse (I)*"),
    rows = 5:7
  ) |>
  gt::tab_row_group(
    label = gt::md("*Consistant reponse (C)*"),
    rows = 2:4
  )

tbl1
```

## Figure 3. Positivity rates based on self-interpreted HIVST results or the reported number of visible lines, by distribution channel, gender and country, among participants of the first survey phase in Côte d'Ivoire, Mali, and Senegal (2021).

```{r}
#| fig-width: 12
#| fig-height: 8
#| column: page


## Fig 3.a BASED ON SELF-REPORTED RESULTS

# Creating "delivery_channel_by_sex" variable
data <- data |>
  mutate(
    delivery_channel_by_sex = interaction(sex, delivery_channel_grouped) |>
      fct_recode() |>
      fct_relevel(
        "man.MSM-based channels", "woman.MSM-based channels", "man.FSW-based channels",
        "woman.FSW-based channels", "man.Other delivery channels", "woman.Other delivery channels"
      )
  )

# Recoding variable "delivery_channel_by_sex"
d <- data |>
  mutate(
    delivery_channel_by_sex = delivery_channel_by_sex |>
      fct_recode(
        "man\nMSM-\nbased\nchannels" = "man.MSM-based channels",
        "woman\nMSM-\nbased\nchannels" = "woman.MSM-based channels",
        "man\nFSW-\n based\nchannels" = "man.FSW-based channels",
        "woman\nFSW-\nbased\nchannels" = "woman.FSW-based channels",
        "man\nother\nchannels" = "man.Other delivery channels",
        "woman\nother\nchannels" = "woman.Other delivery channels"
      )
  )

# Computing numerators and denominators for the 3 hypothesis
dh1 <- d |>
  group_by(delivery_channel_by_sex, country) |>
  summarise(
    hypothesis = "Low",
    n = sum(HIVST_reported_result %in% c("reactive")),
    N = n()
  )

dh2 <- d |>
  group_by(delivery_channel_by_sex, country) |>
  summarise(
    hypothesis = "Central",
    n = sum(HIVST_reported_result %in% c("reactive")),
    N = sum(!HIVST_reported_result %in% c("DK", "R"))
  )

dh3 <- d |>
  group_by(delivery_channel_by_sex, country) |>
  summarise(
    hypothesis = "High",
    n = sum(HIVST_reported_result %in% c("reactive", "DK", "R")),
    N = n()
  )

# Computing positivity rates based on self interpreted results
d_self_interpreted <- bind_rows(dh1, dh2, dh3) |>
  mutate(
    hypothesis = factor(hypothesis, levels = c("Low", "Central", "High")),
    positivity_rate = n / N,
    label = scales::percent(positivity_rate, suffix = "", accuracy = .1)
  ) |>
  rowwise() |>
  mutate(
    ci_low = prop.test(n, N)$conf.int[1],
    ci_high = prop.test(n, N)$conf.int[2]
  ) |>
  mutate(
    ylabel = ci_high + .01
  )

# Dropping cells with denominator below 25
to_drop <- d_self_interpreted$N < 25
d_self_interpreted$positivity_rate[to_drop] <- NA
d_self_interpreted$ci_low[to_drop] <- NA
d_self_interpreted$ci_high[to_drop] <- NA
d_self_interpreted$label[to_drop] <- ""
d_self_interpreted$label[to_drop & d_self_interpreted$hypothesis == "Central"] <- "*"
d_self_interpreted$ylabel[to_drop] <- 0

# Generating Figure 3.a.
g_self_interpreted <-
  d_self_interpreted |>
  ggplot() +
  aes(
    x = delivery_channel_by_sex,
    y = positivity_rate,
    ymin = ci_low,
    ymax = ci_high,
    color = hypothesis,
    shape = hypothesis,
  ) +
  geom_errorbar(
    position = position_dodge(width = 0.9),
    color = "#555555",
    width = .2
  ) +
  geom_point(
    stat = "identity",
    position = position_dodge(width = 0.9),
    size = 2
  ) +
  geom_text(
    mapping = aes(label = label, y = ylabel),
    position = position_dodge(width = 0.9),
    color = "black",
    size = 2.5,
    vjust = 0,
    face = "bold"
  ) +
  facet_grid(cols = vars(country)) +
  labs(title = "", x = "", y = "") +
  ggtitle("a. Positivity rates based on self-interpreted HIVST result") +
  scale_y_continuous(labels = scales::percent, limits = c(0, .5)) +
  scale_colour_manual(values = c("#ACD39E", "#5AAE61", "#1B7837")) +
  theme_classic() +
  theme(
    panel.grid.major.y = element_line(colour = "#DDDDDD", linetype = "dotted"),
    legend.title = element_blank(),
    legend.position = "bottom"
  )

## Fig 3.b BASED ON THE NUMBER OF LINES

# Computing numerators and denominators for the 3 hypothesis
dh1 <- d |>
  group_by(delivery_channel_by_sex, country) |>
  summarise(
    hypothesis = "Low",
    n = sum(HIVST_reported_lines %in% c("2 lines")),
    N = n()
  )

dh2 <- d |>
  group_by(delivery_channel_by_sex, country) |>
  summarise(
    hypothesis = "Central",
    n = sum(HIVST_reported_lines %in% c("2 lines")),
    N = sum(!HIVST_reported_result %in% c("DK", "R"))
  )

dh3 <- d |>
  group_by(delivery_channel_by_sex, country) |>
  summarise(
    hypothesis = "High",
    n = sum(HIVST_reported_lines %in% c("2 lines", "DK", "R")),
    N = n()
  )

# Computing positivity rates based on the number of reported lines
d_reported_lines <- bind_rows(dh1, dh2, dh3) |>
  mutate(
    hypothesis = factor(hypothesis, levels = c("Low", "Central", "High")),
    positivity_rate = n / N,
    label = scales::percent(positivity_rate, suffix = "", accuracy = .1)
  ) |>
  rowwise() |>
  mutate(
    ci_low = prop.test(n, N)$conf.int[1],
    ci_high = prop.test(n, N)$conf.int[2]
  ) |>
  mutate(
    ylabel = ci_high + .01
  )

# Dropping cells with denominator below 25
to_drop <- d_reported_lines$N < 25
d_reported_lines$positivity_rate[to_drop] <- NA
d_reported_lines$ci_low[to_drop] <- NA
d_reported_lines$ci_high[to_drop] <- NA
d_reported_lines$label[to_drop] <- ""
d_reported_lines$label[to_drop & d_reported_lines$hypothesis == "Central"] <- "*"
d_reported_lines$ylabel[to_drop] <- 0

# Generating Figure 3.b.
g_reported_lines <-
  d_reported_lines |>
  ggplot() +
  aes(
    x = delivery_channel_by_sex,
    y = positivity_rate,
    ymin = ci_low,
    ymax = ci_high,
    color = hypothesis,
    shape = hypothesis,
  ) +
  geom_errorbar(
    position = position_dodge(width = 0.9),
    color = "#555555",
    width = .2
  ) +
  geom_point(
    stat = "identity",
    position = position_dodge(width = 0.9),
    size = 2
  ) +
  geom_text(
    mapping = aes(label = label, y = ylabel),
    position = position_dodge(width = 0.9),
    color = "black",
    size = 2.5,
    vjust = 0,
    face = "bold"
  ) +
  facet_grid(cols = vars(country)) +
  labs(title = "", x = "", y = "") +
  ggtitle("b. Positivity rates based on the reported number of lines") +
  scale_y_continuous(labels = scales::percent, limits = c(0, .5)) +
  scale_colour_manual(values = c("#F67E4B", "#DD3D2D", "#A50026")) +
  theme_classic() +
  theme(
    panel.grid.major.y = element_line(colour = "#DDDDDD", linetype = "dotted"),
    legend.title = element_blank(),
    legend.position = "bottom"
  )

# Combining fig 3.a and fig 3.b
wrap_plots(
  g_self_interpreted,
  g_reported_lines,
  nrow = 2
) &
  theme(
    axis.text = element_text(size = 8.5),
    legend.text = element_text(size = 9),
    panel.spacing = unit(0.5, "lines"),
    strip.text = element_text(size = 9, face = "bold"),
    plot.title = element_text(size = 10, face = "bold")
  )
```

## Figure 4: Elements for the flow chart of the participant selection process for Phase 2 of the survey

```{r}
n_eligible_for_phase_2 <- data |>
  filter((HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines")) |>
  nrow()

n_eligible_agreed_recontacted <- data |>
  filter((HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines") & (recontact_phase2 == "yes")) |>
  nrow()

n_eligible_refusal_recontacted <- data |>
  filter((HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines") & (recontact_phase2 == "no")) |>
  nrow()

n_unreachable_time_phase2 <- data |>
  filter((HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines") & (recontact_phase2 == "yes") & (status_phase2 == "unavaible" | is.na(status_phase2))) |>
  nrow()

n_successfully_recontacted <- data |>
  filter((HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines") & (recontact_phase2 == "yes") & (status_phase2 != "unavaible" & !is.na(status_phase2))) |>
  nrow()

n_refused_partcipate <- data |>
  filter((HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines") & (recontact_phase2 == "yes") & (status_phase2 != "unavaible" & !is.na(status_phase2)) & status_phase2 == "refusal") |>
  nrow()

n_accepted_partcipate <- data |>
  filter((HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines") & (recontact_phase2 == "yes") & (status_phase2 != "unavaible" & !is.na(status_phase2)) & status_phase2 != "refusal") |>
  nrow()

n_disconneted <- data |>
  filter((HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines") & (recontact_phase2 == "yes") & (status_phase2 != "unavaible" & !is.na(status_phase2)) & (status_phase2 != "refusal" & status_phase2 == "disconnected")) |>
  nrow()

n_dropped <- data |>
  filter((HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines") & (recontact_phase2 == "yes") & (status_phase2 != "unavaible" & !is.na(status_phase2)) & status_phase2 != "refusal" & status_phase2 == "dropped out before the end") |>
  nrow()

n_completed_quest_phase2 <- data |>
  filter((HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines") & (recontact_phase2 == "yes") & status_phase2 == "questionnaires completed") |>
  nrow()
```

-   `r n_eligible_for_phase_2` participants who reported a reactive test and/or two lines in phase 1

-   `r n_eligible_agreed_recontacted` participants who agreed to be recalled several months later

-   `r n_eligible_refusal_recontacted` participants refusal to be recalled several months later

-   `r n_unreachable_time_phase2` unreachable at the time of phase 2

-   `r n_successfully_recontacted` successfully recontacted for phase 2

-   `r n_refused_partcipate` refused to participate

-   `r n_accepted_partcipate` accepted to participate in phase 2 survey

-   `r n_disconneted` disconnected and unreachable

-   `r n_dropped` dropped out before the end

-   `r n_completed_quest_phase2` completed questionnaire phase 2

## Table 2. Linkage to confirmatory testing, proportion being confirmed HIV positive and treatment initiation, by reported number of lines and self-interpreted HIVST result among phase 2 eligible participants who completed their questionnaire

```{r fig.height=6, fig.width=12, message=FALSE,warning=FALSE}
# Selection of individuals eligible for phase 2 who completed their phase 2 questionnaire
phase2 <- data |>
  filter(
    HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines",
    recontact_phase2 == "yes",
    status_phase2 == "questionnaires completed"
  ) |>
  mutate(all = "ALL")
phase2$test_result_number_lines_reported <- droplevels(phase2$test_result_number_lines_reported)

# Computing the numbers of individuals who completed phase 2 questionnaire
tbl_completed <-
  phase2 |>
  tbl_summary(
    include = c(all, test_result_number_lines_reported),
    statistic = ~"{n}"
  ) |>
  modify_header(stat_0 = "**n**") |>
  modify_footnote(update = everything() ~ NA)

# Computing the proportion who linked to confirmatory testing
tbl_linked <-
  phase2 |>
  tbl_summary(
    by = confirmatory_test,
    include = c(all, test_result_number_lines_reported),
    percent = "row",
    digits = ~0
  ) |>
  add_ci(style_fun = ~ label_percent(accuracy = 1, suffix = "")) |>
  modify_column_hide(c(stat_1, ci_stat_1)) |>
  modify_header(stat_2 ~ "**n (%)**") |>
  modify_footnote(update = everything() ~ NA)

# Computing the proportion who were confirmed HIV-positive
tbl_confirmed <-
  phase2 |>
  filter(confirmatory_test == "yes") |>
  mutate(confirmed = confirmatory_test_result == "positive") |>
  tbl_summary(
    by = confirmed,
    include = c(all, test_result_number_lines_reported),
    percent = "row",
    digits = ~0
  ) |>
  add_ci(style_fun = ~ label_percent(accuracy = 1, suffix = "")) |>
  modify_column_hide(c(stat_1, ci_stat_1)) |>
  modify_header(stat_2 ~ "**n (%)**") |>
  modify_footnote(update = everything() ~ NA)

# Computing the proportion who initiated ART
tbl_initiated <-
  phase2 |>
  filter(confirmatory_test_result == "positive") |>
  mutate(initiated = consulted_health_prof == "yes") |>
  mutate(test_result_number_lines_reported = fct_drop(test_result_number_lines_reported)) |>
  tbl_summary(
    by = initiated,
    include = c(all, test_result_number_lines_reported),
    percent = "row",
    digits = ~0
  ) |>
  add_ci(style_fun = ~ label_percent(accuracy = 1, suffix = "")) |>
  modify_column_hide(c(stat_1, ci_stat_1)) |>
  modify_header(stat_2 ~ "**n (%)**") |>
  modify_footnote(update = everything() ~ NA)

# Merging tables
tbl_merge(
  list(tbl_completed, tbl_linked, tbl_confirmed, tbl_initiated),
  tab_spanner = c(
    "**Completed phase 2**",
    "**Linked to confirmatory testing**",
    "**Confirmed HIV positive**",
    "**Initiated ART**"
  )
)
```

## Table S1. Eligibility and participation in phase 2 survey by sociodemographic characteristics, distribution channel, HIV testing history, the reported number of lines and the self-interpreted HIVST result

```{r fig.height=6, fig.width=12, message=FALSE,warning=FALSE}
# Reordoring value labels
data <- data |>
  mutate(
    marital_status = marital_status |>
      fct_relevel(
        "single", "divorced / separated / widowed", "living with partner / married"
      ),
    educational_level = educational_level |>
      fct_relevel(
        "none / primary", "secondary", "higher"
      )
  )


# Improving variable labels
data <- data |>
  set_variable_labels(
    country = "Country",
    delivery_channel_by_sex = "Sex and distribution channel",
    age_group = "Age group",
    marital_status = "Marital status",
    educational_level = "Educational level",
    first_time_tester = "Firt-time tester",
    test_result_number_lines_reported = "Reported number of lines & self-interpreted HIVST result"
  )

# Computing Table S1.
data |>
  tbl_summary(
    include =
      c(
        country, delivery_channel_by_sex, age_group,
        marital_status, educational_level, first_time_tester, test_result_number_lines_reported
      ),
    by = final_status_phas,
    type = list(c(first_time_tester) ~ "categorical")
  ) |>
  add_overall(TRUE) |>
  add_p(
    test = list(
      test_result_number_lines_reported ~ "chisq.test",
      delivery_channel_by_sex ~ "chisq.test"
    )
  )
```

## Table S2.a: Positivity rates based on self-interpreted HIVST result, by distribution channel, gender and country, among phase 1 participants

```{r fig.height=6, fig.width=12, message=FALSE,warning=FALSE}
# Creation of a positivity rate table according to the low hypothesis
tblS2a_low <- data |>
  to_factor() |>
  tbl_custom_summary(
    include = c(country),
    by = delivery_channel_by_sex,
    stat_fns = ~ proportion_summary(variable = "HIVST_reported_result", value = c("reactive")),
    statistic = ~"{prop}% ({n}/{N})",
    digits = ~ list(
      function(x) {
        style_percent(x, digits = 0)
      }, 0, 0
    ),
    overall_row = TRUE
  ) |>
  add_overall(last = FALSE)

# Creation of a positivity rate table according to the central hypothesis
tblS2a_central <- data |>
  filter(HIVST_reported_result != "DK", HIVST_reported_result != "R") |>
  to_factor() |>
  tbl_custom_summary(
    include = c(country),
    by = delivery_channel_by_sex,
    stat_fns = ~ proportion_summary(variable = "HIVST_reported_result", value = c("reactive")),
    statistic = ~"{prop}% ({n}/{N})",
    digits = ~ list(
      function(x) {
        style_percent(x, digits = 0)
      }, 0, 0
    ),
    overall_row = TRUE
  ) |>
  add_overall(last = FALSE)

# Creation of a positivity rate table according to the high hypothesis
tblS2a_high <- data |>
  to_factor() |>
  tbl_custom_summary(
    include = c(country),
    by = delivery_channel_by_sex,
    stat_fns = ~ proportion_summary(variable = "HIVST_reported_result", value = c("reactive", "DK", "R")),
    statistic = ~"{prop}% ({n}/{N})",
    digits = ~ list(
      function(x) {
        style_percent(x, digits = 0)
      }, 0, 0
    ),
    overall_row = TRUE
  ) |>
  add_overall(last = FALSE)

# Merging the different tables
tblS2a <- tbl_stack(list(tblS2a_low, tblS2a_central, tblS2a_high),
  group_header = c(
    "Low: positivity rate based on self reported test result and DK-R are considered not reactive",
    "Central: positivity rate based on self reported test result and DK-R are excluded",
    "High: positivity rate based on self reported test result and DK-R are considered reactive"
  )
)

tblS2a |>
  modify_header(label ~ "**Variable**") |>
  modify_spanning_header(
    c("stat_1", "stat_2") ~ "**MSM-based channels**",
    c("stat_3", "stat_4") ~ "**FSW-based channels**",
    c("stat_5", "stat_6") ~ "**Others delivery channels**",
    stat_0 ~ "**Total**"
  ) |>
  modify_header(
    stat_1 ~ "**Man**", stat_2 ~ "**Woman**",
    stat_3 ~ "**Man**", stat_4 ~ "**Woman**",
    stat_5 ~ "**Man**", stat_6 ~ "**Woman**",
    stat_0 ~ " "
  ) |>
  modify_footnote(update = everything() ~ NA)
```

## Table S2.b: Positivity rates based on the reported number of lines, by distribution channel, gender and country, among phase 1 participants

```{r fig.height=6, fig.width=12, message=FALSE,warning=FALSE}
# Creation of a positivity rate table according to the low hypothesis
tblS2b_low <- data |>
  to_factor() |>
  tbl_custom_summary(
    include = c(country),
    by = delivery_channel_by_sex,
    stat_fns = ~ proportion_summary(variable = "HIVST_reported_lines", value = c("2 lines")),
    statistic = ~"{prop}% ({n}/{N})",
    digits = ~ list(
      function(x) {
        style_percent(x, digits = 0)
      }, 0, 0
    ),
    overall_row = TRUE
  ) |>
  add_overall(last = FALSE)

# Creation of a positivity rate table according to the central hypothesis
tblS2b_central <- data |>
  filter(HIVST_reported_lines != "DK", HIVST_reported_lines != "R") |>
  to_factor() |>
  tbl_custom_summary(
    include = c(country),
    by = delivery_channel_by_sex,
    stat_fns = ~ proportion_summary(variable = "HIVST_reported_lines", value = c("2 lines")),
    statistic = ~"{prop}% ({n}/{N})",
    digits = ~ list(
      function(x) {
        style_percent(x, digits = 0)
      }, 0, 0
    ),
    overall_row = TRUE
  ) |>
  add_overall(last = FALSE)

# Creation of a positivity rate table according to the high hypothesis
tblS2b_high <- data |>
  to_factor() |>
  tbl_custom_summary(
    include = c(country),
    by = delivery_channel_by_sex,
    stat_fns = ~ proportion_summary(variable = "HIVST_reported_lines", value = c("2 lines", "DK", "R")),
    statistic = ~"{prop}% ({n}/{N})",
    digits = ~ list(
      function(x) {
        style_percent(x, digits = 0)
      }, 0, 0
    ),
    overall_row = TRUE
  ) |>
  add_overall(last = FALSE)


# Merging of the different tables
tblS2b <- tbl_stack(list(tblS2b_low, tblS2b_central, tblS2b_high),
  group_header = c(
    "Low: positivity rate based on self reported number of lines and DK-R are considered not 2 lines",
    "Central: positivity rate based on self reported number of lines and DK-R are excluded",
    "High: positivity rate based on self reported number of lines and DK-R are considered 2 lines"
  )
)

tblS2b |>
  modify_header(label ~ "**Variable**") |>
  modify_spanning_header(
    c("stat_1", "stat_2") ~ "**MSM-based channels**",
    c("stat_3", "stat_4") ~ "**FSW-based channels**",
    c("stat_5", "stat_6") ~ "**Others delivery channels**",
    stat_0 ~ "**Total**"
  ) |>
  modify_header(
    stat_1 ~ "**Man**", stat_2 ~ "**Woman**",
    stat_3 ~ "**Man**", stat_4 ~ "**Woman**",
    stat_5 ~ "**Man**", stat_6 ~ "**Woman**",
    stat_0 ~ " "
  ) |>
  modify_footnote(update = everything() ~ NA)
```

## Table S3.a: Positivity rates based on self-interpreted HIVST result, by age group and country, among phase 1 participants

```{r fig.height=6, fig.width=12, message=FALSE,warning=FALSE}
# Creation of a positivity rate table according to the low hypothesis
tblS3a_low <- data |>
  to_factor() |>
  tbl_custom_summary(
    include = c(country),
    by = age_group,
    stat_fns = ~ proportion_summary(variable = "HIVST_reported_result", value = c("reactive")),
    statistic = ~"{prop}% ({n}/{N})",
    digits = ~ list(
      function(x) {
        style_percent(x, digits = 0)
      }, 0, 0
    ),
    overall_row = TRUE
  ) |>
  add_overall(last = FALSE)

# Creation of a positivity rate table according to the central hypothesis
tblS3a_central <- data |>
  filter(HIVST_reported_result != "DK", HIVST_reported_result != "R") |>
  to_factor() |>
  tbl_custom_summary(
    include = c(country),
    by = age_group,
    stat_fns = ~ proportion_summary(variable = "HIVST_reported_result", value = c("reactive")),
    statistic = ~"{prop}% ({n}/{N})",
    digits = ~ list(
      function(x) {
        style_percent(x, digits = 0)
      }, 0, 0
    ),
    overall_row = TRUE
  ) |>
  add_overall(last = FALSE)

# Creation of a positivity rate table according to the high hypothesis
tblS3a_high <- data |>
  to_factor() |>
  tbl_custom_summary(
    include = c(country),
    by = age_group,
    stat_fns = ~ proportion_summary(variable = "HIVST_reported_result", value = c("reactive", "DK", "R")),
    statistic = ~"{prop}% ({n}/{N})",
    digits = ~ list(
      function(x) {
        style_percent(x, digits = 0)
      }, 0, 0
    ),
    overall_row = TRUE
  ) |>
  add_overall(last = FALSE)

# Merging of the different tables
tblS3a <- tbl_stack(list(tblS3a_low, tblS3a_central, tblS3a_high),
  group_header = c(
    "Low: positivity rate based on self reported test result and DK-R are considered not reactive",
    "Central: positivity rate based on self reported test result and DK-R are excluded",
    "High: positivity rate based on self reported test result and DK-R are considered reactive"
  )
)

tblS3a
```

## Table S3.b: Positivity rates based on the reported number of lines, by age group and country, among phase 1 participants

```{r fig.height=6, fig.width=12, message=FALSE,warning=FALSE}
# Creation of a positivity rate table according to the low hypothesis
tblS3b_low <- data |>
  to_factor() |>
  tbl_custom_summary(
    include = c(country),
    by = age_group,
    stat_fns = ~ proportion_summary(variable = "HIVST_reported_lines", value = c("2 lines")),
    statistic = ~"{prop}% ({n}/{N})",
    digits = ~ list(
      function(x) {
        style_percent(x, digits = 0)
      }, 0, 0
    ),
    overall_row = TRUE
  ) |>
  add_overall(last = FALSE)

# Creation of a positivity rate table according to the central hypothesis
tblS3b_central <- data |>
  filter(HIVST_reported_lines != "DK", HIVST_reported_lines != "R") |>
  to_factor() |>
  tbl_custom_summary(
    include = c(country),
    by = age_group,
    stat_fns = ~ proportion_summary(variable = "HIVST_reported_lines", value = c("2 lines")),
    statistic = ~"{prop}% ({n}/{N})",
    digits = ~ list(
      function(x) {
        style_percent(x, digits = 0)
      }, 0, 0
    ),
    overall_row = TRUE
  ) |>
  add_overall(last = FALSE)

# Creation of a positivity rate table according to the high hypothesis
tblS3b_high <- data |>
  to_factor() |>
  tbl_custom_summary(
    include = c(country),
    by = age_group,
    stat_fns = ~ proportion_summary(variable = "HIVST_reported_lines", value = c("2 lines", "DK", "R")),
    statistic = ~"{prop}% ({n}/{N})",
    digits = ~ list(
      function(x) {
        style_percent(x, digits = 0)
      }, 0, 0
    ),
    overall_row = TRUE
  ) |>
  add_overall(last = FALSE)


# Merging of the different tables
tblS3b <- tbl_stack(list(tblS3b_low, tblS3b_central, tblS3b_high),
  group_header = c(
    "Low: positivity rate based on self reported number of lines and DK-R are considered not 2 lines",
    "Central: positivity rate based on self reported number of lines and DK-R are excluded",
    "High: positivity rate based on self reported number of lines and DK-R are considered 2 lines"
  )
)

tblS3b
```

## Table S4: Main reason for not linking to confirmatory testing among phase 2 participants who did not link to confirmatory testing, by reported number of lines and self-interpreted HIVST result

```{r fig.height=6, fig.width=12, message=FALSE,warning=FALSE}
data <- data |>
  set_variable_labels(
    reason_non_confirmatory_test = "Main reason for not  `\n` linking to confirmatory testing among phase 2  `\n`  participants who did not link to confirmatory testing"
  )

data |>
  filter(
    status_phase2 == "questionnaires completed",
    HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines",
    recontact_phase2 == "yes", confirmatory_test == "no"
  ) |>
  mutate(
    test_result_number_lines_reported = droplevels(test_result_number_lines_reported)
  ) |>
  tbl_summary(
    include = c(reason_non_confirmatory_test),
    by = test_result_number_lines_reported, sort = reason_non_confirmatory_test ~ "frequency"
  ) |>
  add_overall()
```

## Table S5: Time between HIVST and confirmatory testing among phase 2 participants who did link to confirmatory testing, by reported number of lines and self-interpreted HIVST result

```{r fig.height=6, fig.width=12, message=FALSE,warning=FALSE}
data <- data |>
  mutate(
    time_to_confirmation = time_to_confirmation |>
      fct_relevel(
        "last than a week after the HIVST", "between 1 and 2 weeks",
        "between 3 and 4 weeks", "between 1 and 2 months after", "3 months later or more"
      )
  )
data <- data |>
  set_variable_labels(
    time_to_confirmation = "Time between HIVST and confirmatory testing `\n` among phase 2 participants who did link to confirmatory testing"
  )

data |>
  filter(
    recontact_phase2 == "yes",
    HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines",
    confirmatory_test_result == "positive",
    confirmatory_test == "yes",
    status_phase2 == "questionnaires completed",
    consulted_health_prof == "yes",
  ) |>
  mutate(
    test_result_number_lines_reported = droplevels(test_result_number_lines_reported)
  ) |>
  tbl_summary(
    include = c(time_to_confirmation),
    by = test_result_number_lines_reported
  ) |>
  add_overall()
```

## Table S6: Place of confirmatory testing among phase 2 participants who did link to confirmatory testing, by reported number of lines and self-interpreted HIVST result

```{r fig.height=6, fig.width=12, message=FALSE,warning=FALSE}
data <- data |>
  set_variable_labels(
    place_confirmatory_test = "Place of confirmatory testing among phase 2 participants who did link to confirmatory testing"
  )

data |>
  filter(
    recontact_phase2 == "yes",
    HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines",
    status_phase2 == "questionnaires completed",
    confirmatory_test == "yes"
  ) |>
  mutate(
    test_result_number_lines_reported = droplevels(test_result_number_lines_reported)
  ) |>
  tbl_summary(
    include = c(place_confirmatory_test),
    by = test_result_number_lines_reported
  ) |>
  add_overall()
```

## Table S7. Time between phase 1 and phase 2 interviews among phase 2 participants who did link to confirmatory testing, by reported number of lines and self-interpreted HIVST result

```{r fig.height=6, fig.width=12,message=FALSE,warning=FALSE}
# Computing time between phase 1 and phase 2
data$date_phase1 <- ymd_hms(data$date_phase1)
data$date_phase2 <- ymd_hms(data$date_phase2)
data <-
  data |>
  mutate(
    time_between_phase1_2 = time_length(date_phase2 - date_phase1, "month")
  )
data$time_between_phase1_2 <- cut(data$time_between_phase1_2,
  include.lowest = TRUE,
  right = TRUE,
  dig.lab = 4,
  breaks = c(3, 4, 6, 9),
  labels = c("less than 4 months", "Between 4 and 6 months", "More than 6 months")
)

# labeling the "place_confirmatory_test" variable
data <- data |>
  set_variable_labels(
    time_between_phase1_2 = "Time between phase 1 and phase 2 interviews among phase 2 participants who did link to confirmatory testing"
  )

# Generating the S7 table
data |>
  filter(
    recontact_phase2 == "yes",
    HIVST_reported_result == "reactive" | HIVST_reported_lines == "2 lines",
    status_phase2 == "questionnaires completed"
  ) |>
  mutate(
    test_result_number_lines_reported = droplevels(test_result_number_lines_reported)
  ) |>
  tbl_summary(
    include = c(time_between_phase1_2),
    by = test_result_number_lines_reported
  ) |>
  add_overall()
```
