Demo_Manuscript/zz_regression_analyses.R

In PAHPLabResearch/FLASH: FLASH: Free Living Activity Study for Health

# Establish which pairs of variables to test ------------------------------

  var_pairs <-
    list(c("act24_kcal", "AG_kcal", "swa_kcal")) %>%
    rep(2) %>%
    expand.grid(
      KEEP.OUT.ATTRS = FALSE,
      stringsAsFactors = FALSE
    ) %>%
    .[.$Var1 != .$Var2, ] %>%
    stats::setNames(c("xvar", "yvar")) %>%
    as.list(.)

# Functions ---------------------------------------------------------------

  read_format <- function(file) {

    d <- readRDS(file)

    d %<>%
      names(.) %>%
      .[grepl("kcal", .)] %>%
      d[ ,.] %>%
      apply(1, anyNA) %>%
      {d[!., ]}

    if (nrow(d) < 900) {
      warning(
        "Less than 15 h of wear time for ",
        d$id[1], " -- skipping", call. = FALSE
      )
      return(NULL)
    }

    d

  }

  regression_wrapper <- function(file, var_pairs) {

    cat("\nProcessing", basename(file))
    d <- read_format(file)

    if (is.null(d)) return (NULL)

    {mapply(
      regression_test,
      xvar = var_pairs$xvar,
      yvar = var_pairs$yvar,
      MoreArgs = list(dataset = d),
      SIMPLIFY = FALSE
    )} %>%
    lapply(coef_summary) %>%
    c(make.row.names = FALSE) %>%
    do.call(rbind, .)

  }

  regression_test <- function(xvar, yvar, dataset) {

    dataset[ ,xvar] %>%
    mean(.) %>%
    {within(dataset, {center_val = .})} %>%
    lm(
      I(eval(parse(text = yvar)) - center_val) ~
      I(eval(parse(text = xvar)) - center_val),
      .
    ) %>%
    summary(.) %>%
    test_beta(c(0,1)) %>%
    c(id = dataset$id[1], xvar = xvar, yvar = yvar, .)

  }

  test_beta <- function(summary, H0 = 0, ...) {

    summary$coefficients[ ,"t value"] <-
      (summary$coefficients[ ,"Estimate"] - H0) /
      summary$coefficients[ ,"Std. Error"]

    summary$coefficients[ ,"Pr(>|t|)"] <-
      2 * pt(
        abs(summary$coefficients[ ,"t value"]),
        summary$df[2],
        lower.tail = FALSE
      )

    summary

  }

  coef_summary <- function(x) {

    data.frame(
      id = x$id,
      xvar = x$xvar,
      yvar = x$yvar,
      intercept = x$coefficients[1,"Estimate"],
      slope = x$coefficients[2,"Estimate"],
      r2_adjusted = x$adj.r.squared,
      stringsAsFactors = FALSE,
      row.names = NULL
    )

  }

  d_sum <- function(d, description) {

    paste(d$xvar, d$yvar) %>%
    split(d, .) %>%
    lapply(function(x) {
      varsum(x) %>%
      data.frame(
        xvar = x$xvar[1],
        yvar = x$yvar[1],
        .,
        stringsAsFactors = FALSE,
        row.names = NULL
      )
    }) %>%
    c(make.row.names = FALSE) %>%
    do.call(rbind, .) %>%
    data.frame(
      description = description,
      .,
      stringsAsFactors = FALSE,
      row.names = NULL
    )

  }

  varsum <- function(d) {

    mapply(
      function(d, variable) PAutilities::mean_sd(
        d[ ,variable], give_df = FALSE,
        digits = 2, nsmall = 2
      ),
      variable = list("intercept", "slope", "r2_adjusted"),
      MoreArgs = list(d = d),
      SIMPLIFY = FALSE
    ) %>%
    c(stringsAsFactors = FALSE) %>%
    do.call(data.frame, .) %>%
    stats::setNames(c("intercept", "slope", "r2_adjusted"))

  }