---
title: "estimation of micronutrient intakes under no fortification (parts of Tables 4, 5 and 6)"
author: "R Goto"
date: "2024-07-01"
output: html_document
---
## 5.1. Settings
```{r eval=FALSE}
library(tidyverse)
library(summarytools)
library(survey)
library(srvyr)
```

## 5.2. Dataset
- open data of cleaned 60 food consumption per day per AFE and weighting factors
```{r eval=FALSE}
f60afec <- read_csv("f60afec.csv")  # cleaned data of 60 food consumption per day per AFE
TZFCT3 <- read_csv("TFNC_NCT_NTPS20_v.3.0.0.csv") # organised food composition table (FCT) for 60 food items

dataA <- read_csv("hh_sec_a.csv")
dataAw <- dataA %>% 
  select(y4_hhid, domain, strataid, clusterid, clustertype, y4_weights)  # organise dataA with weighting factors

```

## 5.3. Calculate estimated micronutrient intakes
```{r eval=FALSE}
# organise FCT
names(TZFCT3)
TZFCT3 <- TZFCT3 %>% 
  rename("itemcode"="item_id",
         "itemname"="item_desc",
         "VAmcg"="VITA_RAEmcg",
         "FOLmcg"="FOLmcg_standardised",
         "VB12mcg"="VITB12mcg")

TZFCT3 <- TZFCT3 %>% 
  select(itemcode, ENERCkcal, FEmg, ZNmg, VAmcg, FOLmcg, VB12mcg)

write_csv(TZFCT3, 'TZFCT3.csv')

# merge TZFCT with estimated food consumption in g
nf60afec <- left_join(f60afec, TZFCT3, by = "itemcode")

# calculate estimated micronutrient intakes per day per AFE in HHs considering that FCT values are calculated nutrient amount per 100g

options(scipen = 10, digits=2)

nf60afec <- nf60afec %>% 
  mutate(feafe = g_afec*FEmg*0.01,
         znafe = g_afec*ZNmg*0.01,
         vaafe = g_afec*VAmcg*0.01,
         folafe = g_afec*FOLmcg*0.01,
         vb12afe = g_afec*VB12mcg*0.01,
         kcalafe = g_afec*ENERCkcal*0.01)

# calculate sum of nutrients from 60 food items by HHs
nafet <- nf60afec %>% 
  group_by(y4_hhid) %>% 
  summarise(feafet=sum(feafe, na.rm=TRUE),
            znafet=sum(znafe, na.rm=TRUE),
            vaafet=sum(vaafe, na.rm=TRUE),
            folafet=sum(folafe, na.rm=TRUE),
            vb12afet=sum(vb12afe, na.rm=TRUE),
            kcalafet=sum(kcalafe, na.rm=TRUE))

write_csv(nafet, "nafet.csv")

# check the energy distribution
install.packages("cowplot")
library(cowplot)
kcal1 <- nafet %>% 
  ggplot() +
  geom_histogram(aes(kcalafet))
kcal2 <- nafet%>% 
  ggplot() +
  geom_boxplot(aes(kcalafet))
cowplot::plot_grid(kcal1, kcal2, nrow = 2, labels = "AUTO")

png("kcal.png")

```

# get estimated micronutrient intakes by 60 food items in median and the interquartile range (IQR)
# in national (Table 4)
```{r eval=FALSE}
# weighting
nmA <- left_join(nafet, dataAw, by = 'y4_hhid')   # merge with weighting factor

nmA_w <- nmA %>%
  as_survey_design(ids = clusterid, 
                   weights = y4_weights, 
                   strata = strataid,
                   nest=TRUE)

# get estimated micronutrient intakes by 60 food items in median and the interquartile range (IQR)
# in national (Table 4, Supplement Table 3) - included kcal
options(scipen = 4, digits=4)
nmA_wmd <- nmA_w %>%
  summarise(femd = survey_median(feafet, na.rm=TRUE),
            feqt = survey_quantile(feafet, c(0.25, 0.75), na.rm=TRUE),            ,
            znmd = survey_median(znafet, na.rm=TRUE),
            znqt = survey_quantile(znafet, c(0.25, 0.75), na.rm=TRUE),
            vamd = survey_median(vaafet, na.rm=TRUE),
            vaqt = survey_quantile(vaafet, c(0.25, 0.75), na.rm=TRUE),
            folmd = survey_median(folafet, na.rm=TRUE),
            folqt = survey_quantile(folafet, c(0.25, 0.75), na.rm=TRUE),
            vb12md = survey_median(vb12afet, na.rm=TRUE),
            vb12qt = survey_quantile(vb12afet, c(0.25, 0.75), na.rm=TRUE),
            kcalmd = survey_median(kcalafet, na.rm=TRUE),
            kcalqt = survey_quantile(kcalafet, c(0.25, 0.75, 0.90, 0.95), na.rm=TRUE)
            )

# by urban/rural
nmA_wmdur <- nmA_w %>%
  group_by(clustertype) %>% 
  summarise(femd = survey_median(feafet, na.rm=TRUE),
            feqt = survey_quantile(feafet, c(0.25, 0.75), na.rm=TRUE),
            znmd = survey_median(znafet, na.rm=TRUE),
            znqt = survey_quantile(znafet, c(0.25, 0.75), na.rm=TRUE),
            vamd = survey_median(vaafet, na.rm=TRUE),
            vaqt = survey_quantile(vaafet, c(0.25, 0.75), na.rm=TRUE),
            folmd = survey_median(folafet, na.rm=TRUE),
            folqt = survey_quantile(folafet, c(0.25, 0.75), na.rm=TRUE),
            vb12md = survey_median(vb12afet, na.rm=TRUE),
            vb12qt = survey_quantile(vb12afet, c(0.25, 0.75), na.rm=TRUE)
            )

# by strata
nmA_wmddo <- nmA_w %>%
  group_by(domain) %>% 
  summarise(femd = survey_median(feafet, na.rm=TRUE),
            feqt = survey_quantile(feafet, c(0.25, 0.75), na.rm=TRUE),
            znmd = survey_median(znafet, na.rm=TRUE),
            znqt = survey_quantile(znafet, c(0.25, 0.75), na.rm=TRUE),
            vamd = survey_median(vaafet, na.rm=TRUE),
            vaqt = survey_quantile(vaafet, c(0.25, 0.75), na.rm=TRUE),
            folmd = survey_median(folafet, na.rm=TRUE),
            folqt = survey_quantile(folafet, c(0.25, 0.75), na.rm=TRUE),
            vb12md = survey_median(vb12afet, na.rm=TRUE),
            vb12qt = survey_quantile(vb12afet, c(0.25, 0.75), na.rm=TRUE)
            )

```


### 5.5. calculate the prevalance of micronutrient inadequecy in no fortification scenario (Table 5 - no fortification)
```{r eval=FALSE}
nafet <- nafet %>% 
  mutate(znadq=if_else(znafet<10.2, "inadequate", "adequate"),
         vaadq=if_else(vaafet<490, "inadequate", "adequate"),
         foladq=if_else(folafet<250, "inadequate", "adequate"),
         vb12adq=if_else(vb12afet<2.0, "inadequate", "adequate"))

# weighting 
nafet_w <- left_join(nafet, dataAw, by = 'y4_hhid')   # merge with weighting factor

nafet_w <- nafet_w %>%
  as_survey_design(ids = clusterid, 
                   weights = y4_weights, 
                   strata = strataid,
                   nest=TRUE)

# proportion at the zinc inadequacy
options(digits=3)
# nationally
znadqn <- nafet_w %>%
  group_by(znadq) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by urban/rural
znadqur <- nafet_w %>%
  group_by(clustertype, znadq) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by strata
znadqdo <- nafet_w %>%
  group_by(domain, znadq) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))

# proportion at the VA inadequacy
# nationally
vaadqn <- nafet_w %>%
  group_by(vaadq) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by urban/rural
vaadqur <- nafet_w %>%
  group_by(clustertype, vaadq) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by strata
vaadqdo <- nafet_w %>%
  group_by(domain, vaadq) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))

# proportion at the folate inadequacy
# nationally
foladqn <- nafet_w %>%
  group_by(foladq) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by urban/rural
foladqur <- nafet_w %>%
  group_by(clustertype, foladq) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by strata
foladqdo <- nafet_w %>%
  group_by(domain, foladq) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))

# proportion at the VB12 inadequacy
# nationally
vb12adqn <- nafet_w %>%
  group_by(vb12adq) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by urban/rural
vb12adqur <- nafet_w %>%
  group_by(clustertype, vb12adq) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by strata
vb12adqdo <- nafet_w %>%
  group_by(domain, vb12adq) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))

```


### 5.6. calculate the prevalance of micronutrient excess intakes in no fortification scenario (Table 6)
- H-UL for iron 45mg, zinc used both EFSA 25mg and IOM 40mg, VA 3000mcg, folate 1700mcg (converted from folic acid 1000mcg in Allen et al. 2020), VB12 - no UL
```{r eval=FALSE}
nafet <- nafet %>% 
  mutate(feex=if_else(feafet<45, "noexcess", "excess"),
         znex25=if_else(znafet<25, "noexcess", "excess"),
         znex40=if_else(znafet<40, "noexcess", "excess"),
         vaex=if_else(vaafet<3000, "noexcess", "excess"),
         folex=if_else(folafet<1700, "noexcess", "excess"))

# weighting 
nafet_w <- left_join(nafet, dataAw, by = 'y4_hhid')   # merge with weighting factor

nafet_w <- nafet_w %>%
  as_survey_design(ids = clusterid, 
                   weights = y4_weights, 
                   strata = strataid,
                   nest=TRUE)

# proportion at the iron excess intake
options(digits=3)
# nationally
feex <- nafet_w %>%
  group_by(feex) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by urban/rural
feexur <- nafet_w %>%
  group_by(clustertype, feex) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by strata
feexdo <- nafet_w %>%
  group_by(domain, feex) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))

# proportion at the zinc excess intake (using 25mg and 40mg)
options(digits=3)
# nationally
znex25 <- nafet_w %>%
  group_by(znex25) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by urban/rural
znex25ur <- nafet_w %>%
  group_by(clustertype, znex25) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by strata
znex25do <- nafet_w %>%
  group_by(domain, znex25) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))

# 40mg
# nationally
znex40 <- nafet_w %>%
  group_by(znex40) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by urban/rural
znex40ur <- nafet_w %>%
  group_by(clustertype, znex40) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by strata
znex40do <- nafet_w %>%
  group_by(domain, znex40) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))


# proportion at the VA excess intake
# nationally
vaexn <- nafet_w %>%
  group_by(vaex) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by urban/rural
vaexur <- nafet_w %>%
  group_by(clustertype, vaex) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by strata
vaexdo <- nafet_w %>%
  group_by(domain, vaex) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))

# proportion at the folate excess intake (used 1700mcg converted from folic acid 1000mcg in Allen et al. 2020)
# nationally
folexn <- nafet_w %>%
  group_by(folex) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by urban/rural
folexur <- nafet_w %>%
  group_by(clustertype, folex) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))
# by strata
folexdo <- nafet_w %>%
  group_by(domain, folex) %>% 
  summarise(prop=survey_prop( , vartype = c("ci"),  level = 0.95))

```