---
title: "Estimating the prevalence of wasting"
author: Tomás Zaba
knitr:
  opts_chunk: 
    collapse: true
    comment: "#>"
vignette: >
  %\VignetteIndexEntry{Estimating the prevalence of wasting}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---

```{r}
#| label: load_library

library(mwana)
```

## Introduction 
This vignette demonstrates how to use the `mwana` package's functions to estimate the prevalence of wasting. The package allow users to estimate prevalence based on: 

  + Weight-for-height z-score (WFHZ) and/or edema;
  + Raw MUAC values and/or edema;
  + MUAC-for-age z-score (MFAZ) and/or edema, and
  + Combined prevalence. 

The prevalence functions in `mwana` were carefully conceived and designed to simplify the workflow of a nutrition data analyst, especially when dealing with data sets containing imperfections that require additional layers of analysis. Let's try to clarify this with two scenarios that I believe will remind you of the complexity involved: 

  + When analysing a multi-area data set, users will likely need to estimate the prevalence for each area individually. Afterward, they must extract the results and collate in a summary table to share.
  
  + When working with MUAC data, if age ratio test is rated as problematic, an additional tool is required to weight the prevalence and correct for age bias, thus the associated likely overestimation of the prevalence. In unfortunate cases where multiple areas face this issue, the workflow must be repeated several times, making the process cumbersome and highly error-prone 😬.

With `mwana` you no longer have to worry about this 🥳 as the functions are designed to deal with that. To demonstrate their use, we will use different data sets containing some imperfections alluded above:

  + `anthro.02` : a survey data with survey weights. Read more about this data with `?anthro.02`.
  + `anthro.03` : district-level SMART surveys with two districts whose WFHZ standard deviations are rated as problematic while the rest are within range.  Do `?anthro.03` for more details.
  + `anthro.04` : a community-based sentinel site data. The data has different characteristics that require different analysis approaches.

Now we can begin delving into each function.

### Estimation of the prevalence of wasting based on WFHZ {#sec-prevalence-wfhz}

To estimate the prevalence of wasting based on WFHZ we use the `mw_estimate_prevalence_wfhz()` function. The data set to supply must have been wrangled by `mw_wrangle_wfhz()`. 

As usual, we start off by inspecting our data set:

```{r}
#| label: inspect_anthro.02
#| echo: true
#| eval: false

tail(anthro.02)
```

```{r}
#| label: view_anthro.02
#| echo: false

tail(anthro.02)
```

We can see that the data set contains the required variables for a WFHZ prevalence analysis, including for a weighted analysis. This data set has already been wrangled, so we do not need to call the WFHZ wrangler in this case. We will begin the demonstration with an unweigthed analysis - typical of SMART surveys - and then we proceed to a weighted analysis.

#### Estimation of unweighted prevalence 

To achieve this we do: 
```{r}
#| label: unwt_wast_wfhz
#| echo: true
#| eval: false

anthro.02 |>
  mw_estimate_prevalence_wfhz(
    wt = NULL,
    edema = edema,
    .by = NULL
  )
```

This will return: 
```{r}
#| label: view_unwt_wast_wfhz
#| echo: false

mw_estimate_prevalence_wfhz(
  df = anthro.02,
  wt = NULL,
  edema = edema,
  .by = NULL
)
```

If for some reason the variable edema is not available in the data set, or it's there but not plausible, we can exclude it from the analysis by setting the argument `edema` to `NULL`:

```{r}
#| label: unwt_wast_wfhz_noedema
#| echo: true
#| eval: false

anthro.02 |>
  mw_estimate_prevalence_wfhz(
    wt = NULL,
    edema = NULL, # Setting edema to NULL
    .by = NULL
  )
```

And we get: 

```{r}
#| label: view_unwt_wast_wfhz_noedema
#| echo: false

mw_estimate_prevalence_wfhz(
  df = anthro.02,
  wt = NULL,
  edema = NULL,
  .by = NULL
)
```

If we inspect the `gam_n` and `gam_p` columns of this output table and the previous, we notice differences in the numbers. This occurs because edema cases were excluded in the second implementation. Note that you will observed a change if there are positive cases of edema in the data set; otherwise, setting `edema = NULL` will have no effect whatsoever. 

The above output summary does not show results by province. We can control that using the `.by` argument. In the above examples, it was set to `NULL`; now let's pass the name of the column containing the locations where the data was collected. In our case, the column is `province`: 

```{r}
#| label: unwt_wast_wfhz_province
#| echo: true
#| eval: false

anthro.02 |>
  mw_estimate_prevalence_wfhz(
    wt = NULL,
    edema = edema,
    .by = province # province is the variable's name holding data on where the survey was conducted.
  )
```

And _voila_ : 

```{r}
#| label: view_unwt_wast_wfhz_province
#| echo: false

anthro.02 |>
  mw_estimate_prevalence_wfhz(
    wt = NULL,
    edema = edema,
    .by = province
  )
```

A table with two rows is returned with each province's statistics.

#### Estimation of weighted prevalence 

To get the weighted prevalence, we make the use of the `wt` argument. We pass to it the column name containing the final survey weights. In our case, the column name is `wtfactor`:

```{r}
#| label: wt_wasting_wfhz
#| echo: true
#| eval: false

anthro.02 |>
  mw_estimate_prevalence_wfhz(
    wt = wtfactor, # Passing the wtfactor to wt
    edema = edema,
    .by = province
  )
```

And you get: 

```{r}
#| label: view_wt_wasting_wfhz
#| echo: false

anthro.02 |>
  mw_estimate_prevalence_wfhz(
    wt = wtfactor,
    edema = edema,
    .by = province
  )
```

:::{.callout-note}
## The work under the hood of `mw_estimate_prevalence_wfhz`

Under the hood, before starting the prevalence estimation, the function first checks the quality of the WFHZ standard deviation. If it is not rated as problematic, it proceeds with a complex sample-based analysis; otherwise, prevalence is estimated applying the PROBIT method. This is as you see in the body of the plausibility report generated by ENA. The `anthro.02` data set has no such issues, so you don't see `mw_estimate_prevalence_wfhz` in action on this regard. To see that, let's use the `anthro.03` data set.
:::

`anthro.03` contains problematic standard deviation in Metuge and Maravia districts, while the remaining districts are within range. 

Let's inspect our data set:

```{r}
#| label: anthro.3
#| echo: false

head(anthro.03)
```

Now let's apply the prevalence function. This is data is not wrangled, so we will have to wrangle it before passing to the prevalence function: 

```{r}
#| label: anthro.3_prev
#| echo: true
#| eval: false

anthro.03 |>
  mw_wrangle_wfhz(
    sex = sex,
    .recode_sex = TRUE,
    height = height,
    weight = weight
  ) |>
  mw_estimate_prevalence_wfhz(
    wt = NULL,
    edema = edema,
    .by = district
  )
```

The returned output will be: 

```{r}
#| label: view_anthro.3_prev
#| echo: false

anthro.03 |>
  mw_wrangle_wfhz(
    sex = sex,
    .recode_sex = TRUE,
    height = height,
    weight = weight
  ) |>
  mw_estimate_prevalence_wfhz(
    wt = NULL,
    edema = edema,
    .by = district
  )
```

Can you spot the differences? 😎 Yes, you're absolutely correct! While in Cahora-Bassa and Chiúta districts all columns are populated with numbers, in Metuge and Maravia, only the `gam_p`, `sam_p` and `mam_p` columns are filled with numbers, and everything else with `NA`. These are district where the PROBIT method was applied, while in Cahora-Bassa and Chiúta ditricts the standard complex sample analysis was done.

### Estimation of the prevalence of wasting based on MFAZ 

The prevalence of wasting based on MFAZ can be estimated using the `mw_estimate_prevalence_mfaz()` function. This function works and is implemented the same way as demonstrated in @sec-prevalence-wfhz, with the exception of the data wrangling that is based on MUAC. This was demonstrated in the [plausibility checks](https://nutriverse.io/mwana/articles/plausibility.html). In this way, to avoid redundancy, we will not demonstrate the workflow. 

### Estimation of the prevalence of wasting based on raw MUAC values {#sec-prevalence-muac}

This job is assigned to `mw_estimate_prevalence_muac()`. Once you call the function, before starting the prevalence estimation, it first evaluates the acceptability of the MFAZ standard deviation and the age ratio test. Yes, you read well, MFAZ's standard deviation, not on the raw values MUAC.

:::{.callout-important}
Although the acceptability is evaluated on the basis of MFAZ, the actual prevalence is estimated on the basis of the raw MUAC values. MFAZ is also used to detect outliers and flag them to be excluded from the prevalence analysis.
:::
 
The MFAZ standard deviation and the age ratio test results are used to control the prevalence analysis flow in this way: 

  + If the MFAZ standard deviation and the age ratio test are both not problematic, a standard complex sample-based prevalence is estimated. 
  + If the MFAZ standard deviation is not problematic but the age ratio test is problematic, the SMART MUAC tool age-weighting approach is applied. 
  + If the MFAZ standard deviation is problematic, even if age ratio is not problematic, no prevalence analysis is estimated, instead `NA` are thrown. 

When working with a multiple-area data set, these conditionals will still be applied according to each area's situation.

:::{.callout-note}
## How does it work on a multi-area data set

Fundamentally, the function performs the standard deviation and age ratio tests, evaluates their acceptability, and returns a summarized table by area. It then iterates over that summary table row by row checking the above conditionals. Based on the conditionals of each row (area), the function accesses the original data set, computes the prevalence accordingly, and returns the results. 
:::

To demonstrate this we will use the `anthro.04` data set. 

As usual, let's first inspect it: 
```{r}
#| label: anthro.04
#| echo: false

tail(anthro.04)
```

You see that this data has already been wrangled, so we will go straight to the prevalence estimation. 

:::{.callout-important}
As in ENA Software, make sure you run the plausibility check before you call the prevalence function. This is good to know about the acceptability of your data. If we do that with `anthro.04` we will see which province has issues, hence what we should be expecting to see in below demonstrations based on the conditionals stated above.
:::

```{r}
#| label: prev_muac
#| echo: true
#| eval: false

anthro.04 |>
  mw_estimate_prevalence_muac(
    wt = NULL,
    edema = edema,
    .by = province
  )
```

This will return: 

```{r}
#| label: view_prev_muac
#| echo: false

anthro.04 |>
  mw_estimate_prevalence_muac(
    wt = NULL,
    edema = edema,
    .by = province
  )
```

We see that in Province 1, all columns are filled with numbers; in Province 2, some columns are filled with numbers, while other columns are filled with `NA`s: this is where the age-weighting approach was applied. Lastly, in Province 3 a bunch of `NA` are filled everywhere - you know why 😉 . 

Alternatively, we can choose to apply the function that calculates the age-weighted prevalence estimates inside `mw_estimate_prevalence_muac()` directly on to our data set. This can be done by calling the `mw_estimate_smart_age_wt()` function. It worth noting that although possible, we recommend to use the main function. This is simply due the fact that if we decide to use the function independently, then we must, before calling it, check the acceptability of the standard deviation of MFAZ and of the age ratio test, and then evaluate if the conditions that fits the use `mw_estimate_smart_age_wt()` are there. We would have to do that ourselves. This introduces some kind of cumbersomeness in the workflow, and along with that, a risk of picking a wrong analysis workflow.

Nonetheless, if for any reason we decide to go for it anyway, then we would apply the function as demonstrated below. We will continue using the `anthro.04` data set. For this demonstration, we will just pull out the data set from `Province 2` where we already know that the conditions to apply `mw_estimate_smart_age_wt()` are met, and then we will pipe it in to the function:
```{r}
#| label: smart_wt
#| echo: true
#| eval: false

anthro.04 |>
  subset(province == "Province 2") |>
  mw_estimate_smart_age_wt(
    edema = edema,
    .by = NULL
  )
```

This returns the following: 
```{r}
#| label: smart_wt_view
#| echo: false

anthro.04 |>
  subset(province == "Province 2") |>
  mw_estimate_smart_age_wt()
```

#### Estimation of weighted prevalence

For this we go back `anthro.02` data set. 

We approach this task as follows: 

```{r}
#| label: wt_muac_prev
#| echo: true
#| eval: false

## Load library ----
library(dplyr)

## Compute prevalence ----
anthro.02 |>
  mw_wrangle_age(
    age = age,
    .decimals = 2
  ) |>
  mw_wrangle_muac(
    sex = sex,
    .recode_sex = FALSE,
    muac = muac,
    .recode_muac = TRUE,
    .to = "cm",
    age = "age"
  ) |>
  mutate(
    muac = recode_muac(muac, .to = "mm")
  ) |>
  mw_estimate_prevalence_muac(
    wt = wtfactor,
    edema = edema,
    .by = province
  )
```

This will return: 

```{r}
#| label: view_wt_muac_prev
#| echo: false
#| message: false

library(dplyr)

anthro.02 |>
  mw_wrangle_age(
    age = age,
    .decimals = 2
  ) |>
  mw_wrangle_muac(
    sex = sex,
    .recode_sex = FALSE,
    muac = muac,
    .recode_muac = TRUE,
    .to = "cm",
    age = "age"
  ) |>
  mutate(
    muac = recode_muac(muac, .to = "mm")
  ) |>
  mw_estimate_prevalence_muac(
    wt = wtfactor,
    edema = edema,
    .by = province
  )
```

:::{.callout-warning}
You may have noticed that in the above code block, we called the `recode_muac()` function inside `mutate()`. This is because after you use `mw_wrangle_muac()`, it puts the MUAC variable in centimeters. The `mw_estimate_prevalence_muac()` function was defined to accept MUAC in millimeters. Therefore, it must be converted to millimeters. 
:::

#### Estimation for non survey data

Thus far, the demonstration has been around survey data. However, it is also common in day-to-day practice to come across non survey data to analyse. Non survey data can be screenings or any other 
kind of community-based surveillance data. With this kind of data, the analysis workflow usually consists in a simple estimation of the point prevalence and the counts of the positive cases, without necessarily estimating the uncertainty. `mwana` provides a handy function for this task: `mw_estimate_prevalence_screening()`. Under the hood, this function works exactly the same way as `mw_estimate_prevalence_muac()`. The only difference is that it was designed to deal with non survey data. 

To demonstrate its usage, we will use the `anthro.04` data set.
```{r}
#| label: non-survey
#| echo: true
#| eval: false

anthro.04 |> 
  mw_estimate_prevalence_screening(
    muac = muac, 
    edema = edema, 
    .by = province
  )
```

The returned output is: 
```{r}
#| label: non-survey-view
#| echo: false

anthro.04 |> 
  mw_estimate_prevalence_screening(
    muac = muac, 
    edema = edema, 
    .by = province
  )
```

### Estimation of the combined prevalence of wasting

The estimation of the combined prevalence of wasting is a task attributed to the `mw_estimate_prevalence_combined()` function. The case-definition is based on the WFHZ, the raw MUAC values and edema. From the workflow standpoint, it combines the workflow demonstrated in @sec-prevalence-wfhz and in @sec-prevalence-muac.

To demonstrate it's implementation we will use the `anthro.01` data set.

Let's inspect the data:
```{r}
#| label: view_anthro.01
#| echo: false

head(anthro.01)
```

#### Data wrangling

Fundamentally, it combines the data wrangling workflow of WFHZ and MUAC: 
```{r}
#| label: combined_wrangling
#| echo: true
#| eval: false

## Load library ----
library(dplyr)

## Apply the wrangling workflow ----
anthro.01 |>
  mw_wrangle_age(
    dos = dos,
    dob = dob,
    age = age,
    .decimals = 2
  ) |>
  mw_wrangle_muac(
    sex = sex,
    .recode_sex = TRUE,
    muac = muac,
    .recode_muac = TRUE,
    .to = "cm",
    age = "age"
  ) |>
  mutate(
    muac = recode_muac(muac, .to = "mm")
  ) |>
  mw_wrangle_wfhz(
    sex = sex,
    weight = weight,
    height = height,
    .recode_sex = FALSE
  )
```

This is to get the `wfhz` and `flag_wfhz` the `mfaz` and `flag_mfaz` added to the data set. In the output below, we have just selected these columns: 

```{r}
#| label: view_combined_wrangling
#| echo: false
#| message: false

library(dplyr)

anthro.01 |>
  mw_wrangle_age(
    dos = dos,
    dob = dob,
    age = age,
    .decimals = 2
  ) |>
  mw_wrangle_muac(
    sex = sex,
    .recode_sex = TRUE,
    muac = muac,
    .recode_muac = TRUE,
    .to = "cm",
    age = "age"
  ) |>
  mutate(
    muac = recode_muac(muac, .to = "mm")
  ) |>
  mw_wrangle_wfhz(
    sex = sex,
    weight = weight,
    height = height,
    .recode_sex = FALSE
  ) |>
  select(area, wfhz, flag_wfhz, mfaz, flag_mfaz)
```

Under the hood, `mw_estimate_prevalence_combined()` applies the same analysis approach as in `mw_estimate_prevalence_wfhz` and in `mw_estimate_prevalence_muac()`. It checks the acceptability of the standard deviation of WFHZ and MFAZ and of the age ratio test. The following conditionals are checked and applied: 

  + If the standard deviation of WFHZ and of MFAZ, and the age ratio test are all concurrently not problematic, the standard complex sample-based estimation is applied.
  + If any of the above is rated problematic, the prevalence is not computed and `NA`s are thrown. 
  
In this function, a concept of "combined flags" is used.
  
:::{.callout-note}
## What is combined flag?

Combined flags consists of defining as flag any observation that is flagged in either `flag_wfhz` or `flag_mfaz` vectors. A new column `cflags` for combined flags is created and added to the data set. This ensures that all flagged observations from both WFHZ and MFAZ data are excluded from the prevalence analysis.
:::

| **flag_wfhz** | **flag_mfaz** | **cflags** |
| :---: | :---: | :---: |
| 1 | 0  | 1 |
| 0 | 1  | 1 |
| 0 | 0  | 0 |
: A glimpse of case-definition of combined flag {#tbl-Table 1}


Now that we understand what happens under the hood, we can now proceed to implement it: 

```{r}
#| label: cwasting
#| echo: true
#| eval: false

## Load library ----
library(dplyr)

## Apply the workflow ----
anthro.01 |>
  mw_wrangle_age(
    dos = dos,
    dob = dob,
    age = age,
    .decimals = 2
  ) |>
  mw_wrangle_muac(
    sex = sex,
    .recode_sex = TRUE,
    muac = muac,
    .recode_muac = TRUE,
    unit = "cm",
    .to = "age"
  ) |>
  mutate(
    muac = recode_muac(muac, .to = "mm")
  ) |>
  mw_wrangle_wfhz(
    sex = sex,
    weight = weight,
    height = height,
    .recode_sex = FALSE
  ) |>
  mw_estimate_prevalence_combined(
    wt = NULL,
    edema = edema,
    .by = area
  )
```

We get this: 

```{r}
#| label: view_cwasting
#| echo: false

library(dplyr)

anthro.01 |>
  mw_wrangle_age(
    dos = dos,
    dob = dob,
    age = age,
    .decimals = 2
  ) |>
  mw_wrangle_muac(
    sex = sex,
    .recode_sex = TRUE,
    muac = muac,
    .recode_muac = TRUE,
    .to = "cm",
    age = "age"
  ) |>
  mutate(
    muac = recode_muac(muac, .to = "mm")
  ) |>
  mw_wrangle_wfhz(
    sex = sex,
    weight = weight,
    height = height,
    .recode_sex = FALSE
  ) |>
  mw_estimate_prevalence_combined(
    wt = NULL,
    edema = edema,
    .by = area
  )
```

In district E `NA`s were returned because there were issues with the data. I leave it to you to figure out what was/were the issue/issues. 

:::{.callout-tip}
Consider running the plausibility checkers.
:::
