Streamlining the analysis of survey data using functional programming in R

1. Simple example

##-------- STEP 0: Data preparation

## Ensuring categorical grouping variables are factors with descriptive labels

data <- mtcars %>% 
          mutate(
            fullsample = factor("All cars"), 
            cyl = factor(cyl, levels = c(4, 6, 8), labels = c("4 cylinders", "6 cylinders", "8 cylinders")), 
            am = factor(am, levels = c(0, 1), labels = c("Automatic transmission", "Manual transmission"))
          )

##-------- STEP 1-2

## Named vector of all the numeric variables
## we want the sample means for:
indicators <- c("mpg" = "Miles per gallon", 
               "hp" = "Gross horsepower", 
               "disp" = "Displacement (cu. in.)", 
               "wt" = "Weight (1,000 lbs)")

# Named vector of all the categorical variables
# containing population sub-groups we want sample means for:
groups <- c("fullsample" = "All cars", 
            "cyl" = "Number of cylinders", 
            "am" = "Type of transmission")

The code section below creates a function called svy_summary which requires 5 inputs: data (your dataset), i (a string - the column name of an indicator in data), g (a string - the column name of a grouping variable in data), iref (a lookup vector with the descriptive names of the indicator), gref (a lookup vector with the descriptive names of the grouping variable). In the code for this function you will see some functions and operators that you may not have seen before, specifically sym, as_string and !!. These are all part of the rlang library which provides a set of tools for tidy evaluation - to help you write functions using tidyverse pipelines and grammars. The svy_summary function returns a data frame with the estimates for the mean and the 95% confidence interval for the indicators and groups specified, as well as descriptive information. The function assumes that there is a variable called “fullsample” to accomodate the case where you want to estimate the mean and confidence interval of an indicator for the population as a whole.

#----------STEP 3:
# Creating function to compute mean & standard deviation (no sampling weights or complex survey design) 
#given an indicator and grouping variable

# INPUTS:
# data --> dnenotes data frame you are using
# i --> indicator(s)
# g --> population sub-group(s)
# iref --> named vector which provides labels for all of your indicators (variables)
# gref --> named vector which provides labels for all of the variables that contain your sub-groups

svy_summary <- function(data, i, g, iref, gref) { 
  
  i <- sym(i)
  g <- sym(g)

  # Computing mean and 95% confidence interval, 
  # adding a name for the indicator, 
  # and group names for group variable
  
  if (as_string(g) %in% c("fullsample")) { 
    results <- data %>%
      summarize(
        mean = mean(!!i, na.rm=TRUE), 
        sd = sd(!!i, na.rm = TRUE), 
        se = sd/sqrt(length(!!i))
      ) %>% 
      mutate(
        group_cat_val = gref[as_string(g)]
      )
  }
  else {
    results <- data %>%
      group_by(!!g) %>%
      summarize(
        mean = mean(!!i, na.rm=TRUE), 
        sd = sd(!!i, na.rm = TRUE), 
        se = sd/sqrt(length(!!i))
      ) %>% 
      mutate(
        group_cat_val = !!g
      )
  }
  
  # Computing lower & upper 95%, adding descriptive info
  results <- results %>%  mutate(
        lci = mean - qnorm(0.975)*se, 
        uci = mean + qnorm(0.975)*se, 
        indicator = as_string(i), 
        indicator_name = iref[as_string(i)], 
        group = as_string(g), 
        group_name = gref[as_string(g)]
      ) %>%
      dplyr::select(
        indicator, indicator_name, group, group_name, group_cat_val, mean, lci, uci
      )
  
  return(results)
  
}

#---------- STEP 4:
# Apply the function from step 3 to each combination of indicator and grouping variable

combinations <- expand.grid(names(indicators), names(groups), stringsAsFactors = FALSE)
# Indicators
is <- combinations[[1]]
# Groups
gs <- combinations[[2]]

results <- dplyr::bind_rows(
              map2(is, gs, svy_summary, data = data, iref = indicators, gref = groups)
            )
 

results %>% 
  kable() %>% 
  kable_styling(bootstrap_options = "responsive", font_size = 10) %>% 
  scroll_box(height = "250px")

2. Example with sampling weights

If you have data from a more complex sampling procedure, you can use the set of functions in the srvyr package to first define the survey design using as_survey_design.

#-------- STEP 0: Data preparation

# Adding two columns, one for survey weights, another for the primary sampling unit

data <- data %>% 
          mutate(
            probweights = runif(n(), min = 0, max = 1), 
            psu = rep(1:3, length.out = n()), 
            strata = 1
          )

#---------- STEP 3:
# Creating function to compute mean & standard deviation 
# (for a data set with sampling weights or complex survey design) 
# given an indicator and grouping variable

# INPUTS:
# same as before, except now you need: 
# w = survey weights
# psu = primary sampling units
# strata = survey strata (if applicable)

svy_summary_weights <- function(data, i, g, psu, strata, w, iref, gref) { 
  
  i <- sym(i)
  g <- sym(g)
  w <- sym(w)
  psu <- sym(psu)
  strata <- sym(strata)
  
  # Defining survey design
  svydesign <- data %>%
    as_survey_design(ids = !!psu, weight = !!w, strata = !!strata, variables = c(!!i, !!g))
  
  # Computing mean and 95% confidence interval, adding a name fot the indicator, and group namees for group variable
  
  if (as_string(g) %in% c("fullsample")) { 
    results <- svydesign %>%
      summarise(
        mean = survey_mean(!!i, vartype = "ci", na.rm=TRUE)
      ) %>% 
      mutate(
        group_cat_val = gref[as_string(g)]
      )
  }
  
  else {
    results <- svydesign %>%
      group_by(!!g) %>%
      summarise(
        mean = survey_mean(!!i, vartype = "ci", na.rm=TRUE)
      ) %>% 
      mutate(
        group_cat_val = !!g
      ) 
  }
  
  results <- results %>% 
    rename(lci = mean_low, uci = mean_upp) %>% 
      mutate(
        indicator = as_string(i), 
        indicator_name = iref[as_string(i)], 
        group = as_string(g), 
        group_name = gref[as_string(g)]
      ) %>%
      dplyr::select(
        indicator, indicator_name, group, group_name, group_cat_val, mean, lci, uci
      )
  
}

#---------- STEP 4:
# Apply the function from step 3 to each combination of indicator and grouping variable

combinations <- expand.grid(names(indicators), names(groups), stringsAsFactors = FALSE)
# Indicators
is <- combinations[[1]]
# Groups
gs <- combinations[[2]]

results <- dplyr::bind_rows(
              map2(is, 
                   gs, 
                   svy_summary_weights, 
                   data = data, 
                   iref = indicators, 
                   gref = groups,
                   psu = "psu", 
                   strata = "strata", 
                   w = "probweights")
            )
 

results %>% 
  kable() %>% 
  kable_styling(bootstrap_options = "responsive", font_size = 10) %>% 
  scroll_box(height = "250px")

You can easily create a more publication friendly of this table (with a subset of the results):

makeci <- function(m, l, u) { 
  paste0(round(m, 1), " (", round(l, 1), " - ", round(u, 1), ")")
}

results %>% 
  mutate(meanci = makeci(mean, lci, uci)) %>% 
  filter(indicator_name %in% c("Miles per gallon", "Gross horsepower")) %>% 
  arrange(indicator_name, group_name, group_cat_val) %>%
  select(indicator_name, group_name, group_cat_val, meanci) %>% 
  rename(`Performance Indicator` = indicator_name, Feature = group_name, Type = group_cat_val,`Mean (95% Confidence Interval)` = meanci) %>% 
  kable() %>% 
  kable_styling(bootstrap_options = "responsive", font_size = 10) %>% 
  collapse_rows(columns = 1:2, valign = "top")