---
title: "iron inadequacy"
author: "R Goto"
date: "2023-12-03"
output: html_document
---
## 7.1. Settings
```{r eval=FALSE}
library(tidyverse)
library(summarytools)
library(survey)
library(srvyr)
```

## 7.2. Dataset
```{r eval=FALSE}
nafet <- read_csv("nafet.csv")  # estimated micronutrient intake per day per AFE in no fortification
nafetF <- read_csv("nafetF.csv")  # estimated micronutrient intake per day per AFE in status quo and full fortification scenarios
dataAw <- read_csv("dataAw.csv")  # data including weighting factors
```
options(scipen = 10, digits=3)

## 7.3. Full probability approach for estimaing iron inadequecy (see Allen et al. 2006, pages 156-159)
### 7.3.1. No fortification
```{r eval=FALSE}
nafetAw <- left_join(nafet, dataAw, by = "y4_hhid")

# in total
fe_prob <- nafetAw %>%
  mutate(feprob = case_when(
    feafet <= 15 ~ 1,
    feafet <= 16.7 & feafet > 15 ~ 0.96,
    feafet <= 18.7 & feafet > 16.7 ~ 0.93,
    feafet <= 21.4 & feafet > 18.7 ~ 0.85,
    feafet <= 23.6 & feafet > 21.4 ~ 0.75,
    feafet <= 25.7 & feafet > 23.6 ~ 0.65,
    feafet <= 27.8 & feafet > 25.7 ~ 0.55,
    feafet <= 30.2 & feafet > 27.8 ~ 0.45,
    feafet <= 33.2 & feafet > 30.2 ~ 0.35,
    feafet <= 37.3 & feafet > 33.2 ~ 0.25,
    feafet <= 45.0 & feafet > 37.3 ~ 0.15,
    feafet <= 53.5 & feafet > 45.0 ~ 0.08,
    feafet <= 63.0 & feafet > 53.5 ~ 0.04,
    feafet > 63 ~ 0)) %>% 
  count(feprob) %>% 
  mutate(feprev=feprob*n) %>% 
  summarise(sumprov=sum(feprev)) # e.g. calculate the value 2458.37/3290=0.7472 (74.7%)

# in rural/urban
fe_probur <- nafetAw %>%
  mutate(feprob = case_when(
    feafet <= 15 ~ 1,
    feafet <= 16.7 & feafet > 15 ~ 0.96,
    feafet <= 18.7 & feafet > 16.7 ~ 0.93,
    feafet <= 21.4 & feafet > 18.7 ~ 0.85,
    feafet <= 23.6 & feafet > 21.4 ~ 0.75,
    feafet <= 25.7 & feafet > 23.6 ~ 0.65,
    feafet <= 27.8 & feafet > 25.7 ~ 0.55,
    feafet <= 30.2 & feafet > 27.8 ~ 0.45,
    feafet <= 33.2 & feafet > 30.2 ~ 0.35,
    feafet <= 37.3 & feafet > 33.2 ~ 0.25,
    feafet <= 45.0 & feafet > 37.3 ~ 0.15,
    feafet <= 53.5 & feafet > 45.0 ~ 0.08,
    feafet <= 63.0 & feafet > 53.5 ~ 0.04,
    feafet > 63 ~ 0)) %>% 
  group_by(clustertype) %>% 
  count(feprob) %>% 
  mutate(feprev=feprob*n) %>% 
  summarise(sumprov=sum(feprev)) 
# e.g. rural 1399.41/1958=71.5%, urban 1058.96/1332=79.5%

# in strata
fe_probdo <- nafetAw %>%
  mutate(feprob = case_when(
    feafet <= 15 ~ 1,
    feafet <= 16.7 & feafet > 15 ~ 0.96,
    feafet <= 18.7 & feafet > 16.7 ~ 0.93,
    feafet <= 21.4 & feafet > 18.7 ~ 0.85,
    feafet <= 23.6 & feafet > 21.4 ~ 0.75,
    feafet <= 25.7 & feafet > 23.6 ~ 0.65,
    feafet <= 27.8 & feafet > 25.7 ~ 0.55,
    feafet <= 30.2 & feafet > 27.8 ~ 0.45,
    feafet <= 33.2 & feafet > 30.2 ~ 0.35,
    feafet <= 37.3 & feafet > 33.2 ~ 0.25,
    feafet <= 45.0 & feafet > 37.3 ~ 0.15,
    feafet <= 53.5 & feafet > 45.0 ~ 0.08,
    feafet <= 63.0 & feafet > 53.5 ~ 0.04,
    feafet > 63 ~ 0)) %>% 
  group_by(domain) %>% 
  count(feprob) %>% 
  mutate(feprev=feprob*n) %>% 
  summarise(sumprov=sum(feprev))
# e.g. Dar es salaam 431.72/539=80.1%, Mainland urban 387.22/529=73.2%, Mainland rural 1221.15/1757=69.5%, Zanzibar 418.28/465=90.0%

```
 
### 7.3.2. Status quo
```{r eval=FALSE}
names(nafetF)
nafetFAw <- left_join(nafetF, dataAw, by = "y4_hhid")

# in total
fe_probpf <- nafetFAw %>%
  mutate(feprob = case_when(
    feafetpf <= 15 ~ 1,
    feafetpf <= 16.7 & feafetpf > 15 ~ 0.96,
    feafetpf <= 18.7 & feafetpf > 16.7 ~ 0.93,
    feafetpf <= 21.4 & feafetpf > 18.7 ~ 0.85,
    feafetpf <= 23.6 & feafetpf > 21.4 ~ 0.75,
    feafetpf <= 25.7 & feafetpf > 23.6 ~ 0.65,
    feafetpf <= 27.8 & feafetpf > 25.7 ~ 0.55,
    feafetpf <= 30.2 & feafetpf > 27.8 ~ 0.45,
    feafetpf <= 33.2 & feafetpf > 30.2 ~ 0.35,
    feafetpf <= 37.3 & feafetpf > 33.2 ~ 0.25,
    feafetpf <= 45.0 & feafetpf > 37.3 ~ 0.15,
    feafetpf <= 53.5 & feafetpf > 45.0 ~ 0.08,
    feafetpf <= 63.0 & feafetpf > 53.5 ~ 0.04,
    feafetpf > 63 ~ 0)) %>% 
  count(feprob) %>% 
  mutate(feprev=feprob*n) %>% 
  summarise(sumprov=sum(feprev))
  
# by urban/rural
fe_probpfur <- nafetFAw %>%
  mutate(feprob = case_when(
    feafetpf <= 15 ~ 1,
    feafetpf <= 16.7 & feafetpf > 15 ~ 0.96,
    feafetpf <= 18.7 & feafetpf > 16.7 ~ 0.93,
    feafetpf <= 21.4 & feafetpf > 18.7 ~ 0.85,
    feafetpf <= 23.6 & feafetpf > 21.4 ~ 0.75,
    feafetpf <= 25.7 & feafetpf > 23.6 ~ 0.65,
    feafetpf <= 27.8 & feafetpf > 25.7 ~ 0.55,
    feafetpf <= 30.2 & feafetpf > 27.8 ~ 0.45,
    feafetpf <= 33.2 & feafetpf > 30.2 ~ 0.35,
    feafetpf <= 37.3 & feafetpf > 33.2 ~ 0.25,
    feafetpf <= 45.0 & feafetpf > 37.3 ~ 0.15,
    feafetpf <= 53.5 & feafetpf > 45.0 ~ 0.08,
    feafetpf <= 63.0 & feafetpf > 53.5 ~ 0.04,
    feafetpf > 63 ~ 0)) %>% 
  group_by(clustertype) %>% 
  count(feprob) %>% 
  mutate(feprev=feprob*n) %>% 
  summarise(sumprov=sum(feprev)) 

# by strata
fe_probpfdo <- nafetFAw %>%
  mutate(feprob = case_when(
    feafetpf <= 15 ~ 1,
    feafetpf <= 16.7 & feafetpf > 15 ~ 0.96,
    feafetpf <= 18.7 & feafetpf > 16.7 ~ 0.93,
    feafetpf <= 21.4 & feafetpf > 18.7 ~ 0.85,
    feafetpf <= 23.6 & feafetpf > 21.4 ~ 0.75,
    feafetpf <= 25.7 & feafetpf > 23.6 ~ 0.65,
    feafetpf <= 27.8 & feafetpf > 25.7 ~ 0.55,
    feafetpf <= 30.2 & feafetpf > 27.8 ~ 0.45,
    feafetpf <= 33.2 & feafetpf > 30.2 ~ 0.35,
    feafetpf <= 37.3 & feafetpf > 33.2 ~ 0.25,
    feafetpf <= 45.0 & feafetpf > 37.3 ~ 0.15,
    feafetpf <= 53.5 & feafetpf > 45.0 ~ 0.08,
    feafetpf <= 63.0 & feafetpf > 53.5 ~ 0.04,
    feafetpf > 63 ~ 0)) %>% 
  group_by(domain) %>% 
  count(feprob) %>% 
  mutate(feprev=feprob*n) %>% 
  summarise(sumprov=sum(feprev)) 
```

### 7.3.3. full fortification
```{r eval=FALSE}
# in total
fe_probff <- nafetFAw %>%
  mutate(feprob = case_when(
    feafetff <= 15 ~ 1,
    feafetff <= 16.7 & feafetff > 15 ~ 0.96,
    feafetff <= 18.7 & feafetff > 16.7 ~ 0.93,
    feafetff <= 21.4 & feafetff > 18.7 ~ 0.85,
    feafetff <= 23.6 & feafetff > 21.4 ~ 0.75,
    feafetff <= 25.7 & feafetff > 23.6 ~ 0.65,
    feafetff <= 27.8 & feafetff > 25.7 ~ 0.55,
    feafetff <= 30.2 & feafetff > 27.8 ~ 0.45,
    feafetff <= 33.2 & feafetff > 30.2 ~ 0.35,
    feafetff <= 37.3 & feafetff > 33.2 ~ 0.25,
    feafetff <= 45.0 & feafetff > 37.3 ~ 0.15,
    feafetff <= 53.5 & feafetff > 45.0 ~ 0.08,
    feafetff <= 63.0 & feafetff > 53.5 ~ 0.04,
    feafetff > 63 ~ 0)) %>% 
  count(feprob) %>% 
  mutate(feprev=feprob*n) %>% 
  summarise(sumprov=sum(feprev))
  
# by urban/rural
fe_probffur <- nafetFAw %>%
  mutate(feprob = case_when(
    feafetff <= 15 ~ 1,
    feafetff <= 16.7 & feafetff > 15 ~ 0.96,
    feafetff <= 18.7 & feafetff > 16.7 ~ 0.93,
    feafetff <= 21.4 & feafetff > 18.7 ~ 0.85,
    feafetff <= 23.6 & feafetff > 21.4 ~ 0.75,
    feafetff <= 25.7 & feafetff > 23.6 ~ 0.65,
    feafetff <= 27.8 & feafetff > 25.7 ~ 0.55,
    feafetff <= 30.2 & feafetff > 27.8 ~ 0.45,
    feafetff <= 33.2 & feafetff > 30.2 ~ 0.35,
    feafetff <= 37.3 & feafetff > 33.2 ~ 0.25,
    feafetff <= 45.0 & feafetff > 37.3 ~ 0.15,
    feafetff <= 53.5 & feafetff > 45.0 ~ 0.08,
    feafetff <= 63.0 & feafetff > 53.5 ~ 0.04,
    feafetff > 63 ~ 0)) %>% 
  group_by(clustertype) %>% 
  count(feprob) %>% 
  mutate(feprev=feprob*n) %>% 
  summarise(sumprov=sum(feprev)) 

# by strata
fe_probffdo <- nafetFAw %>%
  mutate(feprob = case_when(
    feafetff <= 15 ~ 1,
    feafetff <= 16.7 & feafetff > 15 ~ 0.96,
    feafetff <= 18.7 & feafetff > 16.7 ~ 0.93,
    feafetff <= 21.4 & feafetff > 18.7 ~ 0.85,
    feafetff <= 23.6 & feafetff > 21.4 ~ 0.75,
    feafetff <= 25.7 & feafetff > 23.6 ~ 0.65,
    feafetff <= 27.8 & feafetff > 25.7 ~ 0.55,
    feafetff <= 30.2 & feafetff > 27.8 ~ 0.45,
    feafetff <= 33.2 & feafetff > 30.2 ~ 0.35,
    feafetff <= 37.3 & feafetff > 33.2 ~ 0.25,
    feafetff <= 45.0 & feafetff > 37.3 ~ 0.15,
    feafetff <= 53.5 & feafetff > 45.0 ~ 0.08,
    feafetff <= 63.0 & feafetff > 53.5 ~ 0.04,
    feafetff > 63 ~ 0)) %>% 
  group_by(domain) %>% 
  count(feprob) %>% 
  mutate(feprev=feprob*n) %>% 
  summarise(sumprov=sum(feprev)) 
```