calculatingCombinedGAM.Rmd

---
title: Notes on calculating combined GAM estimates within the Rapid Assessment Method (RAM)
author: "Mark Myatt"
date: "15 February 2020" #'`r format(Sys.Date(), "%d %B %Y")`'
fontsize: 12pt
geometry: margin=2cm
documentclass: article
classoption: a4paper
bibliography: bibliography.bib
lot: FALSE
lof: FALSE
link-citations: TRUE
links-as-notes: FALSE
colorlinks: TRUE
linkcolor: blue
citecolor: blue
urlcolor: blue
---

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

if(!require(ggplot2)) install.packages("ggplot2")
if(!require(magrittr)) install.packages("magrittr")

themeSettings <- theme_bw() + 
  theme(panel.border = element_rect(colour = "#993300",
                                    size = 0.5),
        panel.grid.major = element_line(linetype = 1, 
                                        size = 0.2, 
                                        colour = "gray90"),
        panel.grid.minor = element_line(linetype = 0),
        strip.background = element_rect(colour = "#993300",
                                        fill = "#993300"),
        strip.text = element_text(colour = "white"),
        legend.key = element_rect(linetype = 0),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.ticks = element_line(colour = "#993300", size = 0.5))

```

\newpage

# Background

A June 2019 article in Field Exchange advocates for the reporting of combined GAM estimate (i.e., a single prevalence estimate for GAM cases identified by WHZ and by MUAC). By January 2020, the SMART ENA software has included the calculation and reporting of this single GAM prevalence estimate.

This document reports on our current thining on how this single prevalence estimate can be calculated within the Rapid Assessment Method (RAM). Additionally, we explore how oedema cases can also be added in calculating this single GAM prevalence estimate.

# Proposed approach

PROBIT gives a probability so we look to combining two probabilities:

$$ P(GAM_{\text{MUAC}} ~ \cup ~ GAM_{\text{WHZ}}) ~ = ~ P(GAM_{\text{MUAC}}) ~ + ~ P(GAM_{\text{WHZ}}) $$

However, the problem is that we do not have `independent` probabilities. We overestimate because the intersection gets counted twice. Therefore we need:

$$ P(GAM_{\text{MUAC}} ~ \cup ~ GAM_{\text{WHZ}}) ~ = ~ P(GAM_{\text{MUAC}}) ~ + ~ P(GAM_{\text{WHZ}}) ~ - ~ P(GAM_{\text{MUAC}} ~ \cap ~ GAM_{\text{WHZ}}) $$

We have the first two terms but not the third. We can estimate the third term from a 2 by 2 table:

```{r twoByTwo, echo = FALSE, eval = TRUE}
twoByTwo <- data.frame(rbind(c("a", "b"), c("c", "d")))
row.names(twoByTwo) <- c("$\\text{MUAC} < 125$", "$\\text{MUAC} \\geq 125$")

knitr::kable(x = twoByTwo,
             booktabs = TRUE,
             row.names = TRUE,
             col.names = c("$\\text{WHZ} < -2$", "$\\text{WHZ} \\geq -2$"),
             align = "c",
             escape = FALSE,
             format = "latex") %>%
  kableExtra::kable_styling(latex_options = c("striped", "HOLD_position"), 
                            font_size = 12) %>%  
  kableExtra::row_spec(row = 0, bold = TRUE) %>%   
  kableExtra::row_spec(row = 1:2, monospace = TRUE) %>%
  kableExtra::column_spec(column = 1:3, monospace = TRUE) %>%
  kableExtra::column_spec(column = 1, bold = TRUE)  
```

and

$$ P(GAM_{\text{MUAC}} ~ \cap ~ GAM_{\text{WHZ}}) ~ = ~ \frac{a}{a ~ + ~ b ~ + ~ c ~ + ~ d} $$

We have a small sample size so the estimate will lack precision but I think that being "clever" and using something like:

$$ P(GAM_{\text{MUAC}} ~ \cap ~ GAM_{\text{WHZ}}) ~ = ~ P(GAM_{\text{MUAC}}) ~ \times ~ P(GAM_{\text{WHZ}}) $$

will not work as it assumes independence.

We can try to move forward with this hybrid method.

\newpage

# Worked example for combined GAM by WHZ and GAM by MUAC

We try the proposed approach in R using a dataset from Uganda.

```{r example1, echo = TRUE, eval = TRUE}
## Read dataset from Uganda
x <- read.table(file = "data/ugan01.csv", header = TRUE, sep = ",")

## First 10 rows of the Uganda dataset
```

```{r example1a, echo = FALSE, eval = TRUE}
head(x, 10)
```

&nbsp;

We then create case definitions as follows:

```{r example2, echo = TRUE, eval = TRUE}
## Case definitions
x$gamWHZ <- ifelse(x$whz < -2, 1, 2)               ## GAM by WHZ
x$gamMUAC <- ifelse(x$muac < 125, 1, 2)            ## GAM by MUAC
x$cGAM <- ifelse(x$whz < -2 | x$muac < 125, 1, 2)  ## GAM by WHZ and MUAC

## First 10 rows of the updated Uganda dataset contining case definitions
```

```{r example2a, echo = FALSE, eval = TRUE}
head(x, 10)
```

## Combined GAM estimation using the classical approach

Calculating for prevalence using the classical approach we get:

```{r example3, echo = TRUE, eval = TRUE}
## Classic prevalence for GAM by MUAC
round(prop.table(table(x$gamMUAC))[1] * 100, 2)

## Classic prevalence for GAM by WHZ
round(prop.table(table(x$gamWHZ))[1] * 100, 2)

## Classic prevalence for GAM by WHZ and MUAC
round(prop.table(table(x$cGAM))[1] * 100, 2)
```

&nbsp;

We can test whether GAM cases by MUAC and GAM cases by WHZ are independent.

```{r example4, echo = TRUE, eval = FALSE}
## Test if the two case definitions are independent
chisq.test(table(x$gamMUAC, x$gamWHZ))
```

&nbsp;

The chi-square test has a *p-value* of `r chisq.test(table(x$gamMUAC, x$gamWHZ))$p.value` indicating that the two case definitions are not independent.

\newpage

## Combined GAM estimation using proposed hybrid approach

We then proceed with our proposed hybrid approach using simple PROBIT prevalence estimation.

```{r example6, echo = TRUE, eval = TRUE}
## Simple PROBIT prevalence for GAM by MUAC
pMUAC <- pnorm(125, mean(x$muac), sd(x$muac))
```

```{r example6a, echo = FALSE, eval = TRUE}
round(pMUAC * 100, 2)
```

&nbsp;

```{r example7, echo = TRUE, eval = TRUE}
## Simple PROBIT prevalence for GAM by WHZ
pWHZ <- pnorm(-2, mean(x$whz), sd(x$whz))
```

```{r example7a, echo = FALSE, eval = TRUE}
round(pWHZ * 100, 2)
```

&nbsp;

```{r example8, echo = TRUE, eval = TRUE}
## Estimate the UNION probability
pUNION <- table(x$gamMUAC, x$gamWHZ)[1,1] / sum(table(x$gamMUAC, x$gamWHZ))
```

```{r example8a, echo = FALSE, eval = TRUE}
round(pUNION * 100, 2)
```

&nbsp;

```{r example9, echo = TRUE, eval = TRUE}
## Estimate of GAM by MUAC and WHZ by PROBIT
round((pMUAC + pWHZ - pUNION) * 100, 2)
```

\newpage

# Extending the approach to include oedema cases

Using the same concepts, we can extend the approach to include oedema cases. For this we need to take into account not only the intersection between GAM by MUAC and GAM by WHZ but also the intersection of GAM by MUAC and oedema cases, the intersection of GAM by WHZ and oedema cases, and the intersection between GAM by MUAC, GAM by WHZ and oedema cases as follows:

$$\begin{aligned}
P(GAM_{\text{MUAC}} ~ \cup ~ & GAM_{\text{WHZ}} ~ \cup ~ SAM_{\text{oedema}}) ~ = \\
& P(GAM_{\text{MUAC}}) ~ + ~ P(GAM_{\text{WHZ}}) ~ + ~ P(SAM_\text{oedema}) \\ 
~ & - ~ [P(GAM_{\text{MUAC}} ~ \cap ~ GAM_{\text{WHZ}}) ~ - ~ P(GAM_{\text{MUAC}} ~ \cap ~ GAM_{\text{WHZ}} ~ \cap ~ SAM_{\text{oedema}})] \\
~ & - ~ [P(GAM_{\text{MUAC}} ~ \cap ~ SAM_{\text{oedema}}) ~ - ~ P(GAM_{\text{MUAC}} ~ \cap ~ GAM_{\text{WHZ}} ~ \cap ~ SAM_{\text{oedema}})] \\
~ & - ~ [P(GAM_{\text{WHZ}} ~ \cap ~ SAM_{\text{oedema}}) ~ - ~ P(GAM_{\text{MUAC}} ~ \cap ~ GAM_{\text{WHZ}} ~ \cap ~ SAM_{\text{oedema}})] \\
~ & - ~ P(GAM_{\text{MUAC}} ~ \cap ~ GAM_{\text{WHZ}} ~ \cap ~ SAM_{\text{oedema}})
\end{aligned}$$

Simplifying this we get:

$$\begin{aligned}
P(GAM_{\text{MUAC}} ~ \cup ~ & GAM_{\text{WHZ}} ~ \cup ~ SAM_{\text{oedema}}) ~ = \\
& P(GAM_{\text{MUAC}}) ~ + ~ P(GAM_{\text{WHZ}}) ~ + ~ P(SAM_\text{oedema}) \\ 
~ & - ~ P(GAM_{\text{MUAC}} ~ \cap ~ GAM_{\text{WHZ}}) \\
~ & - ~ P(GAM_{\text{MUAC}} ~ \cap ~ SAM_{\text{oedema}} \\
~ & - ~ P(GAM_{\text{WHZ}} ~ \cap ~ SAM_{\text{oedema}}) \\ 
~ & + ~ (2 \times P(GAM_{\text{MUAC}} ~ \cap ~ GAM_{\text{WHZ}} ~ \cap ~ SAM_{\text{oedema}}))
\end{aligned}$$

# Worked example for combined GAM by WHZ, GAM by MUAC and oedema

We try the proposed approach in R using the same dataset from Uganda used in the previous example. We already have an oedema variable that is coded into 1 for an oedema case and 2 for a non-oedema case. We just need to re-define the `cGAM` variable to include oedema cases.

\newpage

```{r example10, echo = TRUE, eval = TRUE}
## Case definition for combined GAM including oedema cases
x$cGAM <- ifelse(x$whz < -2 | x$muac < 125 | x$oedema == 1, 1, 2)

## First 10 rows of the updated Uganda dataset contining case definitions
```

```{r example10a, echo = FALSE, eval = TRUE}
head(x, 10)
```

## Combined GAM estimation using the classical approach

Calculating for oedema and combine GAM prevalence using the classical estimator we get:

```{r example11, echo = TRUE, eval = TRUE}
## Classic prevalence for oedema cases
round(prop.table(table(x$oedema))[1] * 100, 2)

## Classic prevalence for GAM by WHZ and MUAC and oedema cases
round(prop.table(table(x$cGAM))[1] * 100, 2)
```

## Combined GAM estimation using proposed hybrid approach

We then proceed with our proposed hybrid approach using simple PROBIT prevalence estimation to include oedema cases.

```{r example12, echo = TRUE, eval = TRUE}
## Estimate the probability of oedema cases
pOedema <- prop.table(table(x$oedema))[1]
```

```{r example13, echo = TRUE, eval = TRUE}
## Estimate union probability for WHZ and oedema
pOedemaWHZ <- table(x$gamWHZ, x$oedema)[1,1] / sum(table(x$gamWHZ, x$oedema))
```

```{r example13a, echo = FALSE, eval = TRUE}
round(pOedemaWHZ * 100, 2)
```

&nbsp;

```{r example14, echo = TRUE, eval = TRUE}
## Estimate union probability for MUAC and oedema
pOedemaMUAC <- table(x$gamMUAC, x$oedema)[1,1] / sum(table(x$gamMUAC, x$oedema))
```

```{r example14a, echo = FALSE, eval = TRUE}
round(pOedemaMUAC * 100, 2)
```

&nbsp;

```{r example15, echo = TRUE, eval = TRUE}
## Create variable for GAM meeting WHZ, MUAC and oedema criteria
x$gamUnion <- ifelse(x$oedema == 1 & x$gamWHZ == 1 & x$gamMUAC == 1, 1, 2)

## Estimate union probability for MUAC, WHZ and oedema
pGAMunion <- prop.table(table(x$gamUnion))[1]
```

```{r example15a, echo = FALSE, eval = TRUE}
round(pGAMunion * 100, 2)
```

&nbsp;

```{r example16, echo = TRUE, eval = TRUE}
## Estimate combined GAM including oedema cases
round((pMUAC + pWHZ + pOedema - pUNION - pOedemaWHZ - 
         pOedemaMUAC + (2 * pGAMunion)) * 100, 2)
```