03_serological_analysis.rmd

---
title: "Serological Analysis"
date: "`r Sys.Date()`"
output:
  rmdformats::downcute:
    self_contained: true
    highlight: kate
    toc_depth: 3
    default_style: dark
    code_folding: hide
    code_download: true
    highlight_downlit: true
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
## Global options
knitr::opts_chunk$set(
  cache = TRUE,
  warning = FALSE,
  message = FALSE
)

library(tidyverse)
source("https://raw.githubusercontent.com/liwu1/glbh0034/main/glbh0034_functions_v3.R")
```

# Import data

```{r}
ffi_total <- readRDS("01_data/processed/ffi_total.rds")
```

## Format data

```{r}
ffi_tidy <- ffi_total %>%
  drop_na(age, gender) %>%
  pivot_longer(
    cols = c(pf_recent, pv_recent),
    names_to = "typeofplasmodium",
    values_to = "seropositive"
	) %>%
	mutate(
		seropositive = case_when(
			seropositive == "Positive" ~ 1,
			TRUE ~ 0
		)
	)
```

```{r}
ffi_tidy %>%
  rmarkdown::paged_table()
```

# Serological Analysis

## Logistic regression

```{r}
ffi_tidy	%>%
  group_by(ffi_is_district, typeofplasmodium) %>%
	summarise(
		broom::tidy(
     glm(
       seropositive ~ age + gender + 
			 				        ffi_is_fever_month + 
											ffi_is_mosq_net,
       family = "binomial"
     )
		)
	)  %>%
	mutate(
		p.value = scales::pvalue(p.value)
	) %>%
	print(n = Inf)
```

## Seroprevalence

```{r}
ffi_data_sero <- ffi_tidy %>%
  group_nest(ffi_is_district, typeofplasmodium) %>%
  mutate(
    data_sero = map(
      data,
      ~ create.data.object(
        .$age,
        .$seropositive
      )
    ),
    data_age_profile = map(
      data_sero,
      ~ create.age.profile(
        .,
        lag = 0.15,
        analysis = "overall"
      )
    ),
    plot_seroprevalence = pmap(
      list(
				data_age_profile,
				ffi_is_district,
				typeofplasmodium
			),
      ~ ..1$age.profiles %>%
        as_tibble() %>%
				ggplot(
					aes(
						x = age,
						y = sero.prev
					)
				) +
				geom_point(
					size = 4,
					color = "#0144D2",
					fill = "#0144D2",
					alpha = 0.5
				) +
				geom_linerange(
					aes(
						ymin = lower,
						ymax = upper
					),
					size = 0.8,
					color = "#0144D2"
				) +
				scale_y_continuous(
					limits = c(0, 1.1),
					breaks = c(seq(0, 1, 0.2))
				) +
				expand_limits(
					x = 0
				) +
				scale_x_continuous(
					breaks = seq(0, 50, 10)
				) +
				labs(
					x = "Age, years",
					y = "Sero-prevalence",
					title = paste0(..2, " - ", ..3)
				) +
				theme_classic(base_size = 14) +
				theme(
					axis.title = element_text(
						face = "bold"
					),
					plot.title = element_text(
						face = "bold",
						hjust = 0.5
					)
				)
    )
  )

ffi_data_sero$plot_seroprevalence
```


## Force of infection using sero-conversion rates

```{r error=TRUE}
force_infection <- ffi_tidy %>%
  group_nest(typeofplasmodium) %>%
  mutate(
		data_sero = map(
			data,
			~ create.data.object(
				.$age,
				.$seropositive,
				.$ffi_is_district
			)
		),
		# data_scr_fit = map(
		# 	data_sero,
		# 	~ simple.rcm.analysis(
		# 		.,
		# 		analysis = "split-unshared-rho", 
		# 		int.rho = c(0.001, 0.250)
		# 	)
		# )
	)
```

Try only with the first type of plasmodium:
```{r}
simple.rcm.analysis(force_infection$data_sero[[1]], analysis = "split-unshared-rho", int.rho = c(0.001, 0.05))
```