---
title: "Analysis"
author: "Hannah Rickman"
date: "2025-01-09"
output: html_document
editor_options: 
  markdown: 
    wrap: 72
---

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

## Libraries

```{r}
# Data manipulation and plotting
library(tidyverse)

# Descriptive tables
library(gtsummary)

# Bayesian models
library(brms)
library(broom)
library(broom.mixed)
library(emmeans)
library(car)  

# Spatial analysis and mapping
library(sf)
library(spdep)
library(tmap)

# Exploratory plots
library(GGally)
library(ggridges)
library(patchwork)

# Combining figures
library(gridExtra)
library(png)
library(grid)

# Blantyre data
library(mlwdata)
```

## Loading Files

```{r}
u5m <- readRDS("u5m_anon.rds") #Anonymised
cluster_sum_sf<- readRDS("cluster_sum_sf.rds") %>% # Contains shapefiles and cluster-level data
  #CNR currently per 100,000 person per year - change to per 1000 per year
  mutate(cnr_2020u_scaled = cnr_2020u/100,
  #Density is currently people per km2 - scale to per 10m x 10m
  dens_2020u_scaled = dens_2020u/10000)

#Adding cluster-level variables to u5m
cluster_data <- sf::st_drop_geometry(cluster_sum_sf)
u5m <- u5m %>% left_join(cluster_data, by = "cluster")
```

## Table 1.

```{r}
# Covariables for Table 1
table1_covs <- c(
                 "sex", 
                 "hivstatus_binary", "hiv_exp_bin" , "currentcough_ind" , "any_tbsx" , "bcg_ind" , 
         "tbrx_cat"  , "genhealth_cat" , "muacz_2", "hosp_amt_yn", "locschool_ind", "travel_yn",
         "clinic_catchment", "hh_members_n", "means_score_tertile", "poorest_third", "educ_3", "distance_to_clinic", 
         "tb_hhc_lifetime" , "qftp"
                 )

table1 <- u5m %>%
  tbl_summary(
    by = rec_location, # Stratify by "Clinic" vs "Community"
    include = c(age, all_of(table1_covs)),  # Use predefined vector of covariates
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})" ,# Median (IQR)
      all_categorical() ~ "{n}/{N} ({p}%)"              # Frequency (Percentage)
        ),
    missing = "no"
  ) %>%
  add_overall() %>%
  bold_labels() 

# Print Table 1
table1 %>% as.tibble()
```

## Supplementary data 2 - age-adjustedassociations between covariables and recruitment locations

```{r}
priors <- c(
  set_prior("normal(0, 2.5)", class = "b"),            # Fixed effects (age, sex)
  set_prior("normal(0, 5)", class = "Intercept")
)
```

```{r}
# Unadjusted analysis
tab1_supp_results_unadj <- vector("list", length(table1_covs))
names(tab1_supp_results_unadj) <- table1_covs

for (cov in table1_covs) {
  
  model_unadj <- brm(
    formula = as.formula(paste("as.factor(rec_location) ~", cov)),
    data = u5m,
    prior = priors,
    family = bernoulli(),
    silent = TRUE,
    refresh = 0
  )
  
  post_summary <- summary(model_unadj)$fixed
  df <- as_tibble(post_summary, rownames = "term") %>%
    filter(term != "Intercept") %>%
    mutate(
      OR = exp(Estimate),
      CI_low = exp(`l-95% CI`),
      CI_high = exp(`u-95% CI`),
      `OR (95% CrI)` = paste0(round(OR, 2), " (", round(CI_low, 2), "-", round(CI_high, 2), ")")
    ) %>%
    select(term, `OR (95% CrI)`) %>%
    rename(Covariate = term)
  
  tab1_supp_results_unadj[[cov]] <- df
}

# Age-adjusted analysis
tab1_supp_results_ageadj <- vector("list", length(table1_covs))
names(tab1_supp_results_ageadj) <- table1_covs

for (cov in table1_covs) {
  
  model_ageadj <- brm(
    formula = as.formula(paste("as.factor(rec_location) ~ age +", cov)),
    data = u5m,
    prior = priors,
    family = bernoulli(),
    silent = TRUE,
    refresh = 0
  )
  
  post_summary <- summary(model_ageadj)$fixed
  df <- as_tibble(post_summary, rownames = "term") %>%
    filter(term != "Intercept", term != "age") %>%
    mutate(
      OR = exp(Estimate),
      CI_low = exp(`l-95% CI`),
      CI_high = exp(`u-95% CI`),
      `Age-adjusted AOR (95% CrI)` = paste0(round(OR, 2), " (", round(CI_low, 2), "-", round(CI_high, 2), ")")
    ) %>%
    select(term, `Age-adjusted AOR (95% CrI)`) %>%
    rename(Covariate = term)
  
  tab1_supp_results_ageadj[[cov]] <- df
}


# Combine unadjusted results
unadj_df <- bind_rows(tab1_supp_results_unadj, .id = "Variable")  # Covariate = variable name (e.g., "sex")
# Rename coefficient name to "Term" for clarity
unadj_df <- unadj_df %>% rename(Value = Covariate)
# Combine age-adjusted results
ageadj_df <- bind_rows(tab1_supp_results_ageadj, .id = "Variable")  # Covariate = variable name (e.g., "sex")
# Rename coefficient name to "Term" for clarity
ageadj_df <- ageadj_df %>% rename(Value = Covariate)
# Now join by both variable and term
supp_tab_2 <- left_join(unadj_df, ageadj_df, by = c("Variable", "Value"))


# View results
supp_tab_2
```

## Calculating ARTI with 95% CIs for Table 1

```{r}
compute_risk_bayesian <- function(data, n_draws = 10000) {
  set.seed(123)
  n_positive <- sum(data$qftp == "Positive", na.rm = TRUE)
  n_total <- sum(!is.na(data$qftp))
  prevalence_draws <- rbeta(n_draws, n_positive + 1, n_total - n_positive + 1)

  age_values <- data$age[!is.na(data$age)]
  mean_age <- mean(age_values)
  se_age <- sd(age_values) / sqrt(length(age_values))
  age_draws <- rnorm(n_draws, mean = mean_age, sd = se_age)

  # Ensure age > 0 to avoid issues
  age_draws <- pmax(age_draws, 0.1)

  # Transform to risk
  risk_draws <- 1 - (1 - prevalence_draws)^(1 / age_draws)
  
  risk_ci <- quantile(risk_draws, probs = c(0.025, 0.975))
  risk_point <- mean(risk_draws)

  return(list(
    Risk_Estimate = risk_point,
    Risk_Draws = risk_draws,
    UI_Lower = risk_ci[1],
    UI_Upper = risk_ci[2]
  ))
}

group_results <- compute_risk_bayesian(u5m)
clinic_results <- compute_risk_bayesian(u5m %>% filter(rec_location == "Clinic"))
community_results <- compute_risk_bayesian(u5m %>% filter(rec_location == "Community"))

risk_diff_draws <- clinic_results$Risk_Draws - community_results$Risk_Draws
obs_risk_diff <- mean(risk_diff_draws)

prob_diff_gt_zero <- mean(risk_diff_draws > 0)  # Posterior probability that clinic > community
risk_diff_ci <- quantile(risk_diff_draws, probs = c(0.025, 0.975))

# Package results
risk_df <- tibble(
  Group = c("Overall", "Clinic", "Community"),
  Risk_Estimate = c(group_results$Risk_Estimate*100, clinic_results$Risk_Estimate*100, community_results$Risk_Estimate*100),
  UI_Lower = c(group_results$UI_Lower*100, clinic_results$UI_Lower*100, community_results$UI_Lower*100),
  UI_Upper = c(group_results$UI_Upper*100, clinic_results$UI_Upper*100, community_results$UI_Upper*100)
)

risk_diff_df <- tibble(
  Comparison = "Clinic vs Community",
  Observed_Risk_Difference = obs_risk_diff*100,
  UI_Lower = 100*risk_diff_ci[1],
  UI_Upper = 100*risk_diff_ci[2],
  Posterior_Prob_Clinic_Greater = prob_diff_gt_zero
)

risk_df
risk_diff_df
```

## Plotting graphically posterior distributions for ARTI from each area

```{r}
# Combine risk draws into one data frame for plotting
posterior_df <- tibble(
  Clinic = clinic_results$Risk_Draws*100,
  Household = community_results$Risk_Draws*100
) %>%
  pivot_longer(cols = everything(), names_to = "Group", values_to = "ARTI")

ggplot(posterior_df, aes(x = ARTI, fill = Group, color = Group)) +
  geom_density(alpha = 0.4) +
  geom_vline(data = risk_df %>% filter(Group!="Overall"), aes(xintercept = Risk_Estimate, color = Group),
             linetype = "dashed") +
  labs(
    title = "Posterior distributions of Annual Risk of Mtb Infection (ARTI)",
    x = "Estimated ARTI (%)",
    y = "Density"
  ) +
  theme_minimal() +
  scale_fill_brewer(palette = "Set2") +
  scale_color_brewer(palette = "Set2")
```

## Area-level COVS - covariance

```{r}
cluster_cov_labels2 <- c(
              dens_2020u = "Pop Density",
              scale_mean_pov = "Mean poverty score", 
              hiv_prev_mean = "HIV prevalence, %", 
              cnr_2020u = "CNR")
cluster_vars2<- c("dens_2020u", 
              "scale_mean_pov", 
              "hiv_prev_mean", 
              "cnr_2020u")
```

```{r, fig.width = 10, fig.height = 10}
# Function to compute and annotate R² directly on scatterplots
annotate_r2 <- function(data, mapping, ...) {
  x_var <- rlang::as_name(mapping$x)  # Extract variable names
  y_var <- rlang::as_name(mapping$y)

  model <- lm(reformulate(x_var, y_var), data = data)  # Fit model
  r2 <- summary(model)$r.squared  # Extract R²
  
  # Create scatterplot
  p <- ggplot(data, mapping) +
    geom_point(alpha = 0.6) +  # Scatterplot
    annotate("text", x = Inf, y = -Inf, label = paste0("R² = ", round(r2, 2)), 
             hjust = 1.1, vjust = -0.5, size = 4) +  # Annotate R² in bottom right
    theme_bw()
  
  return(p)
}

# Create ggpairs plot with R² annotations on scatterplots
cluster_data %>%
  ggpairs(
    columns = cluster_vars2,
    lower = list(continuous = wrap("points", alpha = 0.6)),  # Scatterplots in lower panel
    upper = list(continuous = annotate_r2),  # Annotate R² directly on scatterplots
    labeller = as_labeller(cluster_cov_labels2)
  ) +
  theme_bw()

```

## Table 2 summary statistics

```{r}
tab2_covs <- c("age", "sex",  "hiv_exp_bin", "muacz_2", "means_score_tertile", "educ_bin", "distance_to_clinic", "tb_hhc_lifetime", 
               "scale_mean_pov", "hiv_prev_mean", "cnr_2020u_scaled", "dens_2020u_scaled")
```

```{r}
table2 <- u5m %>%
  tbl_summary(
    by = qftp, # Stratify by QFTP
    include = all_of(tab2_covs), 
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})" ,# Median (IQR)
      all_categorical() ~ "{n}/{N} ({p}%)"              # Frequency (Percentage)
        ),
    missing = "no"
  ) %>%
  bold_labels() 

# table2 %>% as.tibble() %>% write.csv("table2.csv", na = "")
table2 %>% as.tibble()
```

## Univariable ORs for Table 2

```{r}
priors <- c(
  set_prior("normal(0, 2.5)", class = "b"),            # Fixed effects (age, sex)
  set_prior("normal(0, 5)", class = "Intercept"),      # Intercept
  set_prior("exponential(1)", class = "sd")           # Random effect (cluster variance)
)
```

```{r, eval = FALSE}
univar_models <- list()

univar_re_results <- lapply(tab2_covs, function(cov) {
  
  formula <- as.formula(
    paste("qftp ~ (1 | cluster) + ", cov)
  )
  
  model <- brm(
    formula = formula,
    data = u5m,
    family = bernoulli(link = "logit"),
    prior = priors,
    cores = 4,
    chains = 4,
    seed = 12345,
    iter = 4000,
    refresh = 0
  )
  
  univar_models[[cov]] <- model
  
  broom.mixed::tidy(
    model,
    conf.int = TRUE,
    exponentiate = TRUE
  ) %>%
    filter(term != "(Intercept)",
           term != "sd__(Intercept)") %>%
    mutate(
      AOR_CI = paste0(
        round(estimate, 2), " (",
        round(conf.low, 2), "-",
        round(conf.high, 2), ")"
      )
    ) %>%
    dplyr::select(term, AOR_CI)
})

univar_re_results_table <- bind_rows(univar_re_results)
```


## AORs for Table 2

```{r, eval = FALSE}
adj_models <- list()

adj_re_results <- lapply(tab2_covs, function(cov) {
  
  formula <- as.formula(
    paste("qftp ~ (1 | cluster) + age + sex +", cov)
  )
  
  model <- brm(
    formula = formula,
    data = u5m,
    family = bernoulli(link = "logit"),
    prior = priors,
    cores = 4,
    chains = 4,
    seed = 12345,
    iter = 4000,
    refresh = 0
  )
  
  adj_models[[cov]] <- model
  
  broom.mixed::tidy(
    model,
    conf.int = TRUE,
    exponentiate = TRUE
  ) %>%
    filter(term != "(Intercept)",
           term != "sd__(Intercept)") %>%
    mutate(
      AOR_CI = paste0(
        round(estimate, 2), " (",
        round(conf.low, 2), "-",
        round(conf.high, 2), ")"
      )
    ) %>%
    dplyr::select(term, AOR_CI)
})

adj_re_results_table <- bind_rows(adj_re_results)
```

## Sensitivity analysis additionally adjusting for recruitment location

```{r, eval = FALSE}
sens_models <- list()

sens_re_results <- lapply(tab2_covs, function(cov) {
  
  formula <- as.formula(
    paste("qftp ~ (1 | cluster) + age + sex + rec_location +", cov)
  )
  
  model <- brm(
    formula = formula,
    data = u5m,
    family = bernoulli(link = "logit"),
    prior = priors,
    cores = 4,
    chains = 4,
    seed = 12345,
    iter = 4000,
    refresh = 0
  )
  
  sens_models[[cov]] <- model
  
  broom.mixed::tidy(
    model,
    conf.int = TRUE,
    exponentiate = TRUE
  ) %>%
    filter(
      term != "(Intercept)",
      term != "sd__(Intercept)"
    ) %>%
    mutate(
      AOR_recloc_CI = paste0(
        round(estimate, 2), " (",
        round(conf.low, 2), "-",
        round(conf.high, 2), ")"
      )
    ) %>%
    dplyr::select(term, AOR_recloc_CI)
})

sens_re_results_table <- bind_rows(sens_re_results) %>%
  distinct(term, .keep_all = TRUE)
```

## ORs for poverty tertiles

```{r}
emm_univar <- emmeans(m_univar_means_score_tertile, ~ means_score_tertile)

contrast(
  emm_univar,
  method = "trt.vs.ctrl",
  ref = 3,
  type = "response"
)

emm_adj <- emmeans(m_adj_means_score_tertile, ~ means_score_tertile)

contrast(
  emm_adj,
  method = "trt.vs.ctrl",
  ref = 3,
  type = "response"
)

emm_sens <- emmeans(m_sens_means_score_tertile, ~ means_score_tertile)

contrast(
  emm_sens,
  method = "trt.vs.ctrl",
  ref = 3,
  type = "response"
)
```

## Combine for Table 2 and supp table 3

```{r}
tab2_ors <- univar_re_results_table %>% left_join(adj_re_results_table, by = c("term")) %>%
  left_join(sens_re_results_table, by = c("term"))

tab2_ors
# tab2_ors %>% write.csv("tab2_ors.csv")

```

## Correlation between cluster-level covariables, and Mtb immunoreativity prevalence

```{r, fig.width = 10, fig.height = 10}
prev_sum <- u5m %>%
  group_by(cluster) %>%
  summarise(qftp_prev = mean(qftp == "Positive") * 100, 
            cluster_n = n(), 
            cluster_pos = cluster_n * qftp_prev/100)

cluster_data <- cluster_data %>% left_join(prev_sum, by = c("cluster"))

cluster_cov_labels3 <- c(
              dens_2020u = "Pop Dens (2020)",
              scale_mean_pov = "Mean poverty score", 
              hiv_prev_mean = "HIV prevalence", 
              cnr_2020u = "CNR", 
              qftp_prev = "Prevalence QFTP+")
cluster_vars3<- c("dens_2020u", 
              "scale_mean_pov", 
              "hiv_prev_mean", 
              "cnr_2020u", 
              "qftp_prev")


# Create ggpairs plot with R² annotations on scatterplots
cluster_data %>%
  ggpairs(
    columns = cluster_vars3,
    lower = list(continuous = wrap("points", alpha = 0.6)),  # Scatterplots in lower panel
    upper = list(continuous = annotate_r2),  # Annotate R² directly on scatterplots
    labeller = as_labeller(cluster_cov_labels3)
  ) +
  theme_bw()

```

```{r}
cor_matrix <- cluster_data %>%
  dplyr::select(all_of(cluster_vars3)) %>%
  cor(use = "pairwise.complete.obs")
cor_matrix
```

```{r}
vif_values <- vif(lm(qftp_prev ~ ., data = cluster_data %>% dplyr::select(all_of(cluster_vars3))))
print(vif_values)
```

## Spatial variation in neighbourhood-level covariables (supp fig 5)

```{r}
hiv_plot<- cluster_sum_sf %>%
  ggplot(aes (fill = hiv_prev_mean)) +
  geom_sf() +
  scale_fill_gradient(
    low = "whitesmoke",    # Color for the lowest value
    high = "darkred", # Color for the highest value
    name = "HIV prevalence, %" # Optional: Label for the legend
  ) +
  guides(fill =guide_colorbar(direction = "horizontal")) +
  theme_void() +
  theme(legend.position = "bottom") 
  
dens_plot<- cluster_sum_sf %>%
  ggplot(aes (fill = dens_2020u_scaled)) +
  geom_sf() +
  scale_fill_gradient(
    low = "whitesmoke",    # Color for the lowest value
    high = "darkgreen", # Color for the highest value
    name = "Population density" # Optional: Label for the legend
  ) +
  guides(fill =guide_colorbar(direction = "horizontal")) +
  theme_void() +
  theme(legend.position = "bottom") 

pov_plot<- cluster_sum_sf %>%
  ggplot(aes (fill = scale_mean_pov)) +
  geom_sf() +
  scale_fill_gradient(
    low = "whitesmoke",    # Color for the lowest value
    high = "darkmagenta", # Color for the highest value
    name = "Median poverty score" # Optional: Label for the legend
  ) +
  guides(fill =guide_colorbar(direction = "horizontal")) +
  theme_void() +
  theme(legend.position = "bottom") 

hiv_plot + dens_plot + pov_plot +
  plot_annotation(tag_levels = "A")

```

## Model 1 - adjusted for baseline characteristics only (age + sex), no neighbourhood-level COVs

```{r, eval = FALSE}
m1 <- brm(
  qftp ~ age + sex + (1|cluster),
  data = u5m,
  prior = priors,
  family = bernoulli(link = "logit"),
  chains = 4, iter = 4000, cores = 4, 
  seed = 12345
)
```

```{r}
m1 %>% summary()
```


### Model checking
```{r}
plot(m1, ask = FALSE)
```

```{r}
pp_check(m1)  # Posterior predictive check
```

```{r}
plot(conditional_effects(m1), ask = FALSE)
```

```{r}
mcmc_plot(m1, type = "dens_overlay", prob = 0.95)
```

### ICC - looking how much variation is cluster-level

```{r}
# Extract standard deviation of cluster-level intercept
sd_cluster <- posterior_summary(m1, pars = "sd_cluster__Intercept")[, "Estimate"]

# Compute variance of the random intercept
var_cluster <- sd_cluster^2

# Residual variance for logistic regression (logit link)
var_residual <- (pi^2) / 3  # approx 3.29

# Compute ICC
icc <- var_cluster / (var_cluster + var_residual)

# Print ICC
print(paste("ICC:", round(icc, 4)))
```


```{r}
m1_ranef <- ranef(m1)$cluster

# Create a data frame with the posterior mean and SD of the random intercept for each cluster.
m1_ranef_df <- data.frame(
  cluster = rownames(m1_ranef),
  m1_re_mean = m1_ranef[, "Estimate", "Intercept"],   # posterior mean for the random intercept
  m1_re_sd   = m1_ranef[, "Est.Error", "Intercept"]     # posterior standard deviation (SD) for the random intercept
) %>%
    mutate(
    m1_re_or_mean  = exp(m1_re_mean),  
    m1_re_or_lower = exp(m1_re_mean - 1.96 * m1_re_sd),  # 95% CI lower bound
    m1_re_or_upper = exp(m1_re_mean + 1.96 * m1_re_sd),   # 95% CI upper bound
    cluster_num = as.numeric(gsub("c", "", cluster)),
  cluster_rank = case_when(cluster_num >14 & cluster_num <48 ~ cluster_num - 12, 
                             cluster_num >= 48 ~ cluster_num - 16,
                             cluster_num <=14 ~ cluster_num)
    )


m1_ranef_df %>% dplyr::select(cluster, cluster_rank, m1_re_or_mean, m1_re_or_lower, m1_re_or_upper) %>% arrange(m1_re_or_mean)

```

```{r}
m1_post <- posterior_samples(m1, pars = "r_cluster")

# Reshape to long format for plotting
m1_post_long <- m1_post %>%
  pivot_longer(cols = everything(), names_to = "cluster", values_to = "log_odds") %>%
  mutate(
    cluster_num = as.numeric(gsub("r_cluster\\[c|,Intercept\\]", "", cluster)),  # Clean cluster names
    cluster_rank = case_when(cluster_num >14 & cluster_num <48 ~ cluster_num - 12, 
                             cluster_num >= 48 ~ cluster_num - 16,
                             cluster_num <=14 ~ cluster_num),
    cluster_rank = as.factor(cluster_rank),
    odds_ratio = exp(log_odds)  
  )
m1_post_long 
```

```{r}
cluster_sum_sf <- cluster_sum_sf %>% left_join(m1_ranef_df, by = c("cluster", "cluster_num", "cluster_rank"))
```


### Density distributions for ARTI posterior distributions (supplementary data 4)

```{r}
ggplot() + 
  # Density ridgelines for full posterior distributions
  geom_density_ridges(
    data = m1_post_long, 
    aes(x = odds_ratio, y = fct_reorder(cluster_rank, odds_ratio, .fun = median
                                        ),
        ), 
    fill = "lightseagreen",
    alpha = 0.5,
    scale = 1, rel_min_height = 0.01
  ) +
  # Add Mean OR Points (from summary)
  geom_point(
    data = m1_ranef_df, 
    aes(x = m1_re_or_mean, y = as.factor(cluster_rank)), 
    color = "black", size = 2
  ) +
  
  # Add 95% CI Error Bars (from summary)
  geom_errorbarh(
    data = m1_ranef_df, 
    aes(xmin = m1_re_or_lower, xmax = m1_re_or_upper, y = as.factor(cluster_rank)), 
    height = 0.2, color = "black"
  ) +
  geom_vline(
    xintercept = 1, 
    linetype = "dashed", 
    colour = "darkred"
  ) +
  scale_x_log10() +  # Log-scale for ORs
  labs(
    title = "Neighbourhood-level odds ratios for Mtb immunoreactivity",
    x = "Estimated odds ratio (log scale)",
    y = "Neighbourhood",
    caption = "Black points: mean OR; black bars: 95% CrI"
  ) +
  
  theme_bw() +
  coord_cartesian(xlim = c(0.1,10))
```


## Sensitivity analysis: adjusting for recruitment type

```{r, eval = FALSE}
m1_sens <- brm(
  qftp ~ age + sex + rec_location + (1|cluster),
  data = u5m,
  prior = priors,
  family = bernoulli(link = "logit"),
  chains = 4, iter = 4000, cores = 4, 
  seed = 12345
)
```

```{r}
m1_sens %>% summary()
```

### Model checking

```{r}
plot(m1_sens, ask = FALSE)
```

```{r}
pp_check(m1_sens)  # Posterior predictive check
```

```{r}
plot(conditional_effects(m1_sens), ask = FALSE)
```

```{r}
mcmc_plot(m1_sens, type = "dens_overlay", prob = 0.95)
```


```{r}
m1_sens_ranef <- ranef(m1_sens)$cluster


m1_sens_ranef_df <- data.frame(
  cluster = rownames(m1_sens_ranef),
  m1_sens_re_mean = m1_sens_ranef[, "Estimate", "Intercept"],   
  m1_sens_re_sd   = m1_sens_ranef[, "Est.Error", "Intercept"]     
) %>%
    mutate(
    m1_sens_re_or_mean  = exp(m1_sens_re_mean),  
    m1_sens_re_or_lower = exp(m1_sens_re_mean - 1.96 * m1_sens_re_sd),  # 95% CI lower bound
    m1_sens_re_or_upper = exp(m1_sens_re_mean + 1.96 * m1_sens_re_sd),   # 95% CI upper bound
    cluster_num = as.numeric(gsub("c", "", cluster)),
  cluster_rank = case_when(cluster_num >14 & cluster_num <48 ~ cluster_num - 12, 
                             cluster_num >= 48 ~ cluster_num - 16,
                             cluster_num <=14 ~ cluster_num)
    )


m1_sens_ranef_df %>% dplyr::select(cluster, cluster_rank, m1_sens_re_or_mean, m1_sens_re_or_lower, m1_sens_re_or_upper) %>% arrange(m1_sens_re_or_mean)

```


## Model 2 - adding area-level covariables

```{r, eval = FALSE}
m2_mod <- brm(
  qftp ~ age + sex + 
        poorest_third + hiv_prev_mean + dens_2020u_scaled +
        (1 | cluster),
  data = u5m,
  prior = priors,
  family = bernoulli(link = "logit"),
  seed = 12345,
  chains = 4, iter = 4000, cores = 4, 
  control=list(adapt_delta=0.9)
)

```

```{r}
plot(m2_mod, ask = FALSE)
```

```{r}
m2_mod %>% summary()
```

```{r}
pp_check(m2_mod)  # Posterior predictive check
```

```{r}
plot(conditional_effects(m2_mod), ask = FALSE)
```


## Comparing models

```{r}
m1_tidy <- 
  broom::tidy(m1, conf.int = TRUE, exponentiate = TRUE) %>%
  mutate(simple = paste0(
    round(estimate, 2), " (", round(conf.low, 2), "-", round(conf.high, 2), ")"
  )) %>%
  dplyr::select(term, simple)

m2_mod_tidy<- 
  broom::tidy(m2_mod, conf.int = TRUE, exponentiate = TRUE) %>%
  mutate(extended = paste0(
    round(estimate, 2), " (", round(conf.low, 2), "-", round(conf.high, 2), ")"
  )) %>%
  dplyr::select(term, extended)

tidy_models <- 
m1_tidy %>% right_join(m2_mod_tidy, by = c("term"))

tidy_models
```
### LOO and WAIC 

```{r}
# LOO
loo_m1 <- loo(m1)
loo_m2_mod <- loo(m2_mod)
loo_m1_sens <- loo(m1_sens)

loo_comparison <- loo_compare(
  loo_m1,
  loo_m2_mod,
  loo_m1_sens
)

# Alternatively, WAIC
waic_m1 <- waic(m1)
waic_m2_mod <- waic(m2_mod)
waic_m1_sens <- waic(m1_sens)

waic_comparison <- loo_compare(
  waic_m1,
  waic_m2_mod,
  waic_m1_sens
)

# Extract model fit statistics
model_fit <- tibble(
  term = c("LOOIC", "WAIC"),
  
  simple = c(
    round(loo_m1$estimates["looic", "Estimate"], 2),
    round(waic_m1$estimates["waic", "Estimate"], 2)
  ),
  
  extended = c(
    round(loo_m2_mod$estimates["looic", "Estimate"], 2),
    round(waic_m2_mod$estimates["waic", "Estimate"], 2)
  ),
  
  sensitivity = c(
    round(loo_m1_sens$estimates["looic", "Estimate"], 2),
    round(waic_m1_sens$estimates["waic", "Estimate"], 2)
  )
)

model_fit
```

```{r}
post_samples_m1 <- as_draws_df(m1)
# Standard deviation of the residual (fixed)
sd_residual <- sqrt(pi^2 / 3)

# Posterior samples of the random effect standard deviation (sd(Intercept))
sd_cluster_samples_m1 <- post_samples_m1$`sd_cluster__Intercept`

# Compute ICC for each posterior sample
icc_samples_m1 <- sd_cluster_samples_m1^2 / (sd_cluster_samples_m1^2 + sd_residual^2)

# Calculate 95% credible interval
icc_ci_m1 <- quantile(icc_samples_m1, probs = c(0.025, 0.5, 0.975))

#Now for M2
post_samples_m2 <- as_draws_df(m2_mod)
sd_cluster_samples_m2 <- post_samples_m2$`sd_cluster__Intercept`
icc_samples_m2 <- sd_cluster_samples_m2^2 / (sd_cluster_samples_m2^2 + sd_residual^2)

# Calculate 95% credible interval
icc_ci_m2 <- quantile(icc_samples_m2, probs = c(0.025, 0.5, 0.975))

icc_tib<- 
icc_ci_m1 %>% bind_rows(icc_ci_m2) %>% 
  mutate(model= c("simple", "extended"), 
         icc = paste0(
    round(`50%`, 3), " (", round(`2.5%`, 3), "-", round(`97.5%`, 3), ")"
         )) %>% 
  dplyr::select(model, icc) %>%
  pivot_wider(names_from = model, values_from = icc) %>% mutate(term = "ICC")
icc_tib
```

### Model comparison table

```{r}
# Merge with model results
models_table <- tidy_models %>%
  bind_rows(model_fit) %>%
  filter(term !="(Intercept") %>%
  bind_rows(icc_tib)

models_table
```

## Figure 2 maps

```{r}
m2_mod_ranef <- ranef(m2_mod)$cluster

m2_mod_ranef_df <- data.frame(
  cluster = rownames(m2_mod_ranef),
  m2_mod_re_mean = m2_mod_ranef[, "Estimate", "Intercept"],   
  m2_mod_re_sd   = m2_mod_ranef[, "Est.Error", "Intercept"]     
) %>% 
  mutate(m2_mod_re_or_mean = exp(m2_mod_re_mean)) # exponentiate for OR
cluster_sum_sf <- cluster_sum_sf%>% left_join(m2_mod_ranef_df, by = c("cluster"))
cluster_sum_sf %>% dplyr::select(cluster, contains("re_mean"))
```

### Draw maps with OSM backgrounds
Background

```{r}
timasamala_area<-st_union(cluster_sum_sf) 
timasamala_buffer <- timasamala_area %>%
  st_buffer(dist = 0.008) %>% st_bbox()

ti_osm_buffer<- readRDS("ti_osm_buffer.rds")
```

```{r}
clinics<-blantyre_clinics%>% st_set_crs(4326) #in mlwdata
timasamala_clinics <- blantyre_clinics %>%
  filter(clinic=="Bangwe Health Centre"|
           clinic=="Limbe Health Centre"|
           clinic=="Ndirande Health Centre") %>%
  mutate(clinic = str_remove(clinic, " Health Centre"))
qech <- blantyre_clinics %>%
  filter(clinic=="Queen Elizabeth hospital") %>%
  mutate(clinic="QECH")
```

```{r}
m1_osm <-tm_shape(ti_osm_buffer) +
      tm_rgb() +
      tm_shape(cluster_sum_sf) + 
      tm_polygons(col = "m1_re_or_mean", 
              palette = c("blue", "white", "red"),  # Manually setting colors
              title = "OR",  # Change legend label
              style = "cont",
              breaks = c(0.5, 1, 2),       # Matching ggplot2 breakpoints
              contrast = 1,                         # Ensures full color contrast
              auto.palette.mapping = FALSE,         # Ensures manual colors are used
              legend.is.portrait = TRUE,
              legend.reverse = TRUE,      
              alpha = 0.9) +
      tm_layout(legend.position = c("right","top"),
                legend.outside = FALSE,
                         legend.bg.color = "white", 
                title = "A",
                asp = 0) +
  tm_shape(timasamala_clinics) +
  tm_dots(col="black", size=0.5) +
  tm_text("clinic", col="black", bg.color = "white", bg.alpha = 0.75,
          just = "left", xmod=0.3, ymod=-0.3) + 
  tm_scale_bar(breaks = c(0, 1, 2, 3, 4)) +
  tm_add_legend('symbol', col="black", size=0.5,
                      labels="Clinic") 
  
m2_osm <-tm_shape(ti_osm_buffer) +
      tm_rgb() +
      tm_shape(cluster_sum_sf) + 
      tm_polygons(col = "m2_mod_re_or_mean", 
              palette = c("blue", "white", "red"),  # Manually setting colors
              title = "OR",  # Change legend label
              style = "cont",
              breaks = c(0.5, 1, 2),       # Matching ggplot2 breakpoints
              contrast = 1,                         # Ensures full color contrast
              auto.palette.mapping = FALSE,         # Ensures manual colors are used
              legend.is.portrait = TRUE,
              legend.reverse = TRUE,      
              alpha = 0.9) +
      tm_layout(
  legend.position = c("right", "top"),
  legend.outside = FALSE,
  legend.bg.color = "white",
  title = "B",  
  asp = 0) +
  tm_shape(timasamala_clinics) +
  tm_dots(col="black", size=0.5) +
  tm_text("clinic", col="black", bg.color = "white", bg.alpha = 0.75,
          just="left", xmod=0.3, ymod=-0.3) + 
  tm_scale_bar(breaks = c(0, 1, 2, 3, 4)) +
  tm_add_legend('symbol', col="black", size=0.5,
                      labels="Clinic") 


tmap <- tmap_arrange(m1_osm, m2_osm, ncol = 2, outer.margins = FALSE)
tmap
```

## Inspect for residual spatial autocorrelation using Moran's I

### Construct weights matrix

```{r}
nb <- poly2nb(cluster_sum_sf, queen = TRUE, snap = 0.001)   # Define neighbors

#Adjacency matrix 
adj_matrix <- nb2mat(nb, style = "B", zero.policy = TRUE) 
listw <- nb2listw(nb, style = "W", zero.policy = TRUE) 

```

### Add residuals from models to cluster_sum_sf

```{r}
ind_residuals_m1 <- data.frame(
  cluster = m1$data$cluster,  # Ensure cluster ID is available
  resid = residuals(m1)  # Extract individual-level residuals
)

cluster_residuals_m1 <- ind_residuals_m1 %>%
  group_by(cluster) %>%
  summarise(mean_resid_m1 = mean(resid.Estimate, na.rm = TRUE))  # Aggregate residuals
cluster_sum_sf <- cluster_sum_sf %>%
  left_join(cluster_residuals_m1, by = "cluster")

```

```{r}
moran_test_m1 <- moran.test(cluster_sum_sf$mean_resid_m1, listw)
print(moran_test_m1)
```

### And now for M2
Add residuals from models to cluster_sum_sf

```{r}
ind_residuals_m2 <- data.frame(
  cluster = m2_mod$data$cluster,  # Ensure cluster ID is available
  resid = residuals(m2_mod)  # Extract individual-level residuals
)

cluster_residuals_m2 <- ind_residuals_m2 %>%
  group_by(cluster) %>%
  summarise(mean_resid_m2 = mean(resid.Estimate, na.rm = TRUE))  # Aggregate residuals
cluster_sum_sf <- cluster_sum_sf %>%
  left_join(cluster_residuals_m2, by = "cluster")

moran_test_m2 <- moran.test(cluster_sum_sf$mean_resid_m2, listw)
print(moran_test_m2)
```

## Calculating ARTIs

-   With 50% male and female\
-   Fix age at 3\

```{r}
# Create new dataset with both Male & Female for each cluster
new_data <- expand.grid(
  cluster = unique(u5m$cluster),
  sex = c("Male", "Female"),
  age = 3  # Fix age at 3
)

# Extract posterior draws for both sexes
m1_post_draws <- posterior_linpred(m1, newdata = new_data, re_formula = ~(1|cluster))

# Convert posterior draws to a tibble (long format)
m1_post_df <- as_tibble(m1_post_draws) %>%
  mutate(draw = row_number()) %>%  # Add draw ID
  pivot_longer(-draw, names_to = "obs_id", values_to = "logit_qftp") %>%  # Long format
  mutate(
    cluster = new_data$cluster[as.numeric(gsub("V", "", obs_id))],  # Match cluster
    sex = new_data$sex[as.numeric(gsub("V", "", obs_id))]  # Match sex
  ) %>%
  dplyr::select(-obs_id)  

# Convert logit to probability scale
m1_post_df <- m1_post_df %>%
  mutate(prob_qftp = plogis(logit_qftp))

# Average Male & Female predictions within each draw & cluster
m1_post_avg <- m1_post_df %>%
  group_by(cluster) %>%
  summarise(mean_prob_qftp = mean(prob_qftp), 
            median_prob_qftp = median(prob_qftp),  
            qftp_ci_lo = quantile(prob_qftp, 0.025),  # 2.5% CI
            qftp_ci_hi = quantile(prob_qftp, 0.975),
            .groups = "drop") %>%
  mutate(mean_arti = 1-(1-mean_prob_qftp)^(1/3), 
         m1_arti_ci_lo = 1-(1-qftp_ci_lo)^(1/3),
         m1_arti_ci_hi = 1-(1-qftp_ci_hi)^(1/3)
         )

# View first few rows
m1_post_avg <- m1_post_avg %>%
  mutate(m1_qftp = mean_prob_qftp, 
         m1_arti = mean_arti, 
         cluster_num = as.numeric(gsub("c", "", cluster)),
    cluster_rank = rank(cluster_num)) %>%
  dplyr::select(cluster_rank, cluster, m1_qftp, m1_arti, m1_arti_ci_lo, m1_arti_ci_hi)

m1_post_avg %>% arrange(m1_arti) %>% dplyr::select(cluster_rank, starts_with("m1_arti"))

```

## Sensitivity analysis

```{r}
m1_sens_ranef_df <- data.frame(
  cluster = rownames(m1_sens_ranef),
  m1_sens_re_mean = m1_sens_ranef[, "Estimate", "Intercept"],
  m1_sens_re_sd   = m1_sens_ranef[, "Est.Error", "Intercept"]
) 

m1_sens_ranef_df %>%
  dplyr::select(
    cluster,
    m1_sens_re_or_mean,
    m1_sens_re_or_lower,
    m1_sens_re_or_upper
  ) %>%
  arrange(m1_sens_re_or_mean)
```

### Map of neighbourhoods using output of sensitivity analysis

```{r}
# Join to spatial data
cluster_sum_sf_sens <- cluster_sum_sf %>%
  left_join(
    m1_sens_ranef_df,
    by = c("cluster", "cluster_num")
  )


# Plot the sensitivity model
m1_sens_osm <- tm_shape(ti_osm_buffer) +
  tm_rgb() +
  tm_shape(cluster_sum_sf_sens) + 
  tm_polygons(
    col = "m1_sens_re_or_mean", 
    palette = c("blue", "white", "red"),
    title = "OR",
    style = "cont",
    breaks = c(0.5, 1, 2),
    contrast = 1,
    auto.palette.mapping = FALSE,
    legend.is.portrait = TRUE,
    legend.reverse = TRUE,
    alpha = 0.9
  ) +
  tm_layout(
    legend.position = c("right", "top"),
    legend.outside = FALSE,
    legend.bg.color = "white", 
    title = "A",
    asp = 0
  ) +
  tm_shape(timasamala_clinics) +
  tm_dots(col = "black", size = 0.5) +
  tm_text(
    "clinic",
    col = "black",
    bg.color = "white",
    bg.alpha = 0.75,
    just = "left",
    xmod = 0.3,
    ymod = -0.3
  ) + 
  tm_scale_bar(breaks = c(0, 1, 2, 3, 4)) +
  tm_add_legend(
    "symbol",
    col = "black",
    size = 0.5,
    labels = "Clinic"
  )

m1_sens_osm
```

### Ranking

```{r}
# Random effects from original model
m1_re <- tibble(
  cluster = rownames(m1_ranef),
  m1_ranef = m1_ranef[, "Estimate", "Intercept"]
)

# Random effects from sensitivity model
sens_re <- tibble(
  cluster = rownames(m1_sens_ranef),
  m1_sens_ranef = m1_sens_ranef[, "Estimate", "Intercept"]
)

# Compare
re_compare <- m1_re %>%
  left_join(sens_re, by = "cluster") %>%
  mutate(
    OR_m1 = exp(m1_ranef),
    OR_sens = exp(m1_sens_ranef),
    difference = m1_sens_ranef - m1_ranef,
    OR_ratio = exp(m1_sens_ranef - m1_ranef)
  ) %>%
  mutate(
    rank_m1 = rank(m1_ranef),
    rank_sens = rank(m1_sens_ranef)
  )

cor(
  re_compare$rank_m1,
  re_compare$rank_sens,
  method = "spearman"
)
```

### Compare sensitivity and base model

```{r}
rho <- cor(
  re_compare$m1_ranef,
  re_compare$m1_sens_ranef,
  method = "spearman",
  use = "complete.obs"
)

ggplot(re_compare, aes(x = m1_ranef, y = m1_sens_ranef)) +
  geom_point() +
  geom_abline(
    slope = 1,
    intercept = 0,
    linetype = "dashed"
  ) +
  annotate(
    "text",
    x = -Inf,
    y = Inf,
    label = paste0("Spearman \u03c1 = ", round(rho, 4)),
    hjust = -0.1,
    vjust = 1.5,
    size = 5
  ) +
  labs(
    x = "Cluster random effect: main model",
    y = "Cluster random effect: sensitivity model"
  ) +
  theme_minimal()
```

## Figure 2

### Figure 2A

```{r}
cluster_sum_sf <- cluster_sum_sf %>% left_join(m1_post_avg, by=c("cluster"))

#Add ranks
cluster_sum_sf <-cluster_sum_sf %>% 
  mutate(m1_arti_rank = rank(m1_arti), 
         cnr_rank = rank(cnr_2020u)) 

#Select highest ARTI areas
highlight_arti <- cluster_sum_sf %>%
  filter(m1_arti_rank >= 26)

#Draw map
cnr_osm <-tm_shape(ti_osm_buffer) +
      tm_rgb() +
      tm_shape(cluster_sum_sf) + 
      tm_polygons(col = "cnr_2020u", 
              palette = c("lemonchiffon", "orange", "red"),  # Manually setting colors
              title = "CNR",  # Change legend label
              style = "cont",
              contrast = 1,                         # Ensures full color contrast
              auto.palette.mapping = FALSE,         # Ensures manual colors are used
              legend.is.portrait = TRUE,
              legend.reverse = TRUE) +
  tm_shape(highlight_arti) + 
  tm_polygons(
    border.col = "blue",  # Set border color to blue
    alpha = 0.5,  # Adjust transparency if needed
    lwd = 2.5  # Line width of the borders
  ) +
  tm_add_legend(
    type = "fill",  # Add a fill legend
    col="white",
    labels = "Highest ARTI",  # Legend label
    size = 3,       # Adjust the size of the legend symbol
    border.col = "blue"  # Set the border color of the legend box to blue
  ) +
  tm_shape(timasamala_clinics) +
  tm_dots(col="black", size=0.5) +
  tm_text("clinic", col="black", bg.color = "white", bg.alpha = 0.75,
          just = "left", xmod=0.3, ymod=-0.3) + 
  tm_scale_bar(breaks = c(0, 1, 2, 3, 4)) +
  tm_add_legend('symbol', col="black", size=0.5,
                      labels="Clinic") +
  tm_layout(legend.position = c("right","top"),
                legend.outside = FALSE,
                         legend.bg.color = "white",
              title = "A"
) 

```

### Figure 2B

```{r}
# Top 8 by cnr_2020u (red box)
top_cnr <- cluster_sum_sf %>%
  top_n(8, cnr_2020u)

red_box <- data.frame(
  ymin = min(top_cnr$m1_arti - 0.0015),
  ymax = max(top_cnr$m1_arti + 0.0015),
  xmin = min(top_cnr$cnr_2020u - 5),
  xmax = max(top_cnr$cnr_2020u + 5),
  box_type = "Highest CNR neighbourhoods"  # Add a column to categorize the red box
)

# Top 8 by m1_arti (blue box)
top_arti <- cluster_sum_sf %>%
  top_n(8, m1_arti)

blue_box <- data.frame(
  ymin = min(top_arti$m1_arti - 0.0005),
  ymax = max(top_arti$m1_arti + 0.0005),
  xmin = min(top_arti$cnr_2020u - 15),
  xmax = max(top_arti$cnr_2020u + 15),
  box_type = "Highest ARTI neighbourhoods"  # Add a column to categorize the blue box
)

# Combine both red and blue box data
boxes <- rbind(red_box, blue_box)

# Add both rectangles to plot
g_cnr_arti2 <- cluster_sum_sf %>%
  ggplot(aes(
    x = cnr_2020u, 
    y = m1_arti, 
    label = cluster
  )) + 
  geom_point(alpha = 0.5, size = 3) +
  geom_rect(data = boxes, 
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, color = box_type),  # Map color to the box type
            fill = NA, linetype = "dashed", size = 1, inherit.aes = FALSE) +
  scale_color_manual(
    values = c("Highest CNR neighbourhoods" = "red", "Highest ARTI neighbourhoods" = "blue"),  # Map colors to box types
    labels = c("Highest ARTI neighbourhoods", "Highest CNR neighbourhoods")  # Labels for the legend
  ) +
  theme_bw() + 
  labs(
    x = "CNR, per 100,000 per year", 
    y = "ARTI", 
    title = "B"
  ) +
  scale_y_continuous(labels = scales::percent, limits = c(0, NA)) +  # Set lower limit of y-axis to 0
  scale_x_continuous(limits = c(0, NA)) +  # Set lower limit of x-axis to 0
  theme(
    legend.title = element_blank(),   # Remove legend title
    legend.position="bottom", 
    aspect.ratio = 1, 
    axis.text = element_text(size = 14),  # Increase axis text size
    axis.title = element_text(size = 16),  # Increase axis title size
    legend.text = element_text(size = 14),  # Increase legend text size
  ) 

```

### Combined figure

```{r}
cnr_osm_grob <- grid.grabExpr(print(cnr_osm))

# Combine the tmap image with the ggplot object

fig2 <- grid.arrange(
  cnr_osm_grob,  # tmap plot as an image
  g_cnr_arti2,  # CNR ggplot 
  ncol = 2,  
  widths = c(0.73, 1)  # Control the relative width of the plots 
)
plot(fig2)

```