---
title: "PRAMS Reproducible Report for GA in `r params$year`"
author: "Melinda Higgins"
date: "2025-06-23"
output:
  pdf_document: default
  word_document: default
  html_document: default
subtitle: "Example Report from TIDAL Course - Summer 2025"
params:
  year:
    label: Year
    value: 2020
    input: select
    choices:
    - 2017
    - 2018
    - 2019
    - 2020
    - 2021
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
# set up knitr options
# - set echo = FALSE to hide code in final report
# - set error = TRUE to knit document 
#   even if there are coding errors
# - set warning and message to FALSE to hide
#   these in the final report
knitr::opts_chunk$set(echo = FALSE,
                      error = TRUE,
                      warning = FALSE,
                      message = FALSE)

# load any packages needed here
library(dplyr)
library(survey)
library(gtsummary)
library(gt)

# Load dataset(s) needed
# NOTE: Change directory to location of your PRAMS data
load(file = "PRAMS/prams.RData")

# select GA and the year of interest
# use dplyr to filter out just GA
prams_ga_year <- prams %>%
  filter(STATE == "GA") %>%
  filter(NEST_YR == params$year)
```

## Pregnancy Risk Assessment Monitoring System (PRAMS)

[PRAMS](https://www.cdc.gov/prams/about/index.html) is "an ongoing, site-specific, and population-based surveillance system. PRAMS is designed to identify groups of women and infants at high risk for health problems, to monitor changes in health status, and to measure progress towards goals in improving the health of mothers and infants."

* You can request access to the data from the  [CDC's PRAMS Portal](https://www.cdc.gov/prams/php/data-research/index.html).
* Here is the link to the PRAMS codebook for [Phase 8 Automated Research File Codebook](https://www.cdc.gov/prams/php/data-research/codebooks-phase-8.html).

For this example report we will be looking at these variables from the PRAMS 8 dataset:

* `BF5EVER` - Breastfeed ever
* `VITAMIN` - Multivitamin -- # of times /wk
* `MAT_AGE_PU` - Maternal age grouped

Specifically, we will try to replicate the statistics for the [Selected 2016-2022 Maternal and Child Health (MCH) Indicators](https://www.cdc.gov/prams/php/data-research/mch-indicators-by-site.html) for GA from any of the years from 2017-2021 which had metrics reported for GA.

```{r}
# create any variables needed BEFORE 
# applying the complex survey sampling design weights

# create a factor variable
# and add labels for BF5EVER
prams_ga_year$BF5EVER.f <- factor(
  prams_ga_year$BF5EVER,
  levels = c(1, 2),
  labels = c("NO", "YES")
)

# create a factor variable
# and add labels for VITAMINS
prams_ga_year$VITAMIN.f <- factor(
  prams_ga_year$VITAMIN,
  levels = c(1, 2, 3, 4),
  labels = c("1 = DIDNT TAKE VITAMIN",
             "2 = 1-3 TIMES/WEEK",
             "3 = 4-6 TIMES/WEEK",
             "4 = EVERY DAY/WEEK")
)

# for the metrics, compute a
# variable for anyone who
# took vitamins 4+ times a week
prams_ga_year$VITAMIN_4plus <- 
  ifelse(prams_ga_year$VITAMIN > 2, 1, 0)

# add labels, make a factor
prams_ga_year$VITAMIN_4plus.f <- factor(
  prams_ga_year$VITAMIN_4plus,
  levels = c(0, 1),
  labels = c("3x/week or less",
             "4x/week or more")
)

# set up survey design object
# create the survey design file for GA
# for year specified
prams_ga_year.svy <- 
  svydesign(ids = ~0, strata = ~SUD_NEST, 
            fpc = ~TOTCNT, weights = ~WTANAL, 
            data = prams_ga_year)
```

```{r}
# MCH data pulled out for easy comparison
breastfed_mch <- data.frame(
  year = c(2017, 2018, 2019, 2020, 2021),
  n_denominator = c(919, 735, 766, 677, 757),
  n_numerator = c(749, 605, 615, 563, 590),
  wtpct = c(85.2, 83.0, 79.8, 81.2, 79.8),
  wtpct_95cilb = c(81.6, 79.2, 75.8, 77.1, 75.9),
  wtpct_95ciup = c(88.2, 86.3, 83.2, 84.7, 83.3)
)

vitamin_mch <- data.frame(
  year = c(2017, 2018, 2019, 2020, 2021),
  n_denominator = c(953, 766, 783, 690, 781),
  n_numerator = c(298, 280, 278, 247, 328),
  wtpct = c(31.1, 35.9, 38.3, 35.9, 41.1),
  wtpct_95cilb = c(26.1, 31.6, 34.0, 31.5, 36.7),
  wtpct_95ciup = c(34.5, 40.4, 42.7, 40.7, 45.6)
)
```

\newpage

### MCH Indicators for GA between 2017-2021 for "Ever Breastfed".

Here are the Maternal and Child Health (MCH) Indicators for GA between 2017-2021 for "Ever Breastfed".

```{r}
knitr::kable(
  breastfed_mch,
  col.names = c("Year",
                "Unweighted Total Sample",
                "Unweighted Frequency",
                "Weighted %",
                "Lower 95% Confidence Interval",
                "Upper 95% Confidence Interval"),
  caption = "Ever Breastfed"
)
```

### Estimates for "Ever Breastfed" in GA in `r params$year`.

Here are the estimates computed using the `survey` package and the `gtsummary` package to make a nicely formatted table.

```{r}
# write custom function for confidence intervals
confidence_intervals <- function(data, variable, by, ...) {
  ## extract the confidence intervals and multiply to get percentages
  props <- svyciprop(as.formula(paste0( "~" , variable)),
              data, na.rm = TRUE)
  
  ## extract the confidence intervals 
  as.numeric(confint(props) * 100) %>% ## make numeric and multiply for percentage
    round(., digits = 1) %>%           ## round to one digit
    c(.) %>%                           ## extract the numbers from matrix
    paste0(., collapse = "-")          ## combine to single character
}

tbl_svysummary(
  data = prams_ga_year.svy,
  include = c(BF5EVER.f),
  statistic = list(everything() ~ c("{n} ({p}%)"))
  ) %>%
  add_n() %>%  
  add_stat(fns = everything() ~ confidence_intervals) %>%
  modify_header(
    list(
      n ~ "**Weighted total (N)**",
      stat_0 ~ "**Weighted Count**",
      add_stat_1 ~ "**95%CI**"
    )) %>%
  as_gt() %>%
  tab_options(latex.tbl.pos = "h")
```

\newpage

### MCH Indicators for GA between 2017-2021 for "Multivitamin use 4+ times a week during the month before pregnancy".

Here are the Maternal and Child Health (MCH) Indicators for GA between 2017-2021 for "Multivitamin use 4+ times a week during the month before pregnancy".

```{r}
knitr::kable(
  vitamin_mch,
  col.names = c("Year",
                "Unweighted Total Sample",
                "Unweighted Frequency",
                "Weighted %",
                "Lower 95% Confidence Interval",
                "Upper 95% Confidence Interval"),
  caption = "Multivitamin use 4+ times a week during the month before pregnancy"
)
```

### Estimates for "Multivitamin use 4+ times a week during the month before pregnancy" in GA in `r params$year`.

Here are the estimates computed using the `survey` package and the `gtsummary` package to make a nicely formatted table.

```{r}
tbl_svysummary(
  data = prams_ga_year.svy,
  include = c(VITAMIN_4plus.f),
  statistic = list(everything() ~ c("{n} ({p}%)"))
  ) %>%
  add_n() %>%  
  add_stat(fns = everything() ~ confidence_intervals) %>%
  modify_header(
    list(
      n ~ "**Weighted total (N)**",
      stat_0 ~ "**Weighted Count**",
      add_stat_1 ~ "**95%CI**"
    ))  %>%
  as_gt() %>%
  tab_options(latex.tbl.pos = "h")
```

\newpage

### Logistic Regression - Breastfed Ever by Maternal Age and Multivitamin Use

Look at predicting "Ever Breastfed" by multivitamin use and maternal age.

#### Unweighted Model - not adjusted for complex survey sampling weights

```{r}
glm1 <-glm(BF5EVER.f ~ MAT_AGE_PU + VITAMIN_4plus.f,
           data = prams_ga_year.svy,
           family = "binomial")
gtsummary::tbl_regression(glm1,
                          exponentiate = TRUE)  %>%
  as_gt() %>%
  tab_options(latex.tbl.pos = "h")
```

#### Weighted Model - adjusted for complex survey sampling weights

```{r}
wtglm1 <- svyglm(BF5EVER.f ~ MAT_AGE_PU + VITAMIN_4plus.f, 
                 design = prams_ga_year.svy,
                 family=quasibinomial())
sm1 <- summary(wtglm1)
or1 <- exp(coef(wtglm1))

# pull out odds ratios and other stats
df1 <- data.frame(
  beta = sm1[["coefficients"]][,1],
  se = sm1[["coefficients"]][,2],
  or = or1,
  pval = sm1[["coefficients"]][,4]
)

row.names(df1) <- c(
  "Intercept",
  "Maternal Age",
  "Multivitamin 4x/wk or more"
)

knitr::kable(
  df1,
  col.names = c("Beta",
                "SE Beta",
                "Odds Ratio",
                "p-value"),
  caption = "Weighted Model for Ever Breastfed by Age and Multivitamin Use"
)
```

