data-raw/gaussian_parameters.R

In dewittpe/pedbp: Pediatric Blood Pressure

# Build a data set for the mean and standard deviations for Gaussian
# approximations of the percentiles

# pull in the code for estimating the distribution means and sd
set.seed(42)
source("R/est_norm.R")
load("./data/gemelli1990.rda")
load("./data/lo2013.rda")
load("./data/nhlbi_bp_norms.rda")
load("./data/flynn2017.rda")

data.table::setDT(nhlbi_bp_norms)
data.table::setDT(gemelli1990)
data.table::setDT(lo2013)
data.table::setDT(flynn2017)

sbp <- nhlbi_bp_norms[, as.list(est_norm(sbp, bp_percentile/100)$par), by = .(male, age, height_percentile)]
dbp <- nhlbi_bp_norms[, as.list(est_norm(dbp, bp_percentile/100)$par), by = .(male, age, height_percentile)]
data.table::setnames(sbp, old = c("mean", "sd"), new = c("sbp_mean", "sbp_sd"))
data.table::setnames(dbp, old = c("mean", "sd"), new = c("dbp_mean", "dbp_sd"))
cdc_bp <- merge(sbp, dbp)

sbp <- flynn2017[, as.list(est_norm(sbp, bp_percentile/100)$par), by = .(male, age, height_percentile)]
dbp <- flynn2017[, as.list(est_norm(dbp, bp_percentile/100)$par), by = .(male, age, height_percentile)]
data.table::setnames(sbp, old = c("mean", "sd"), new = c("sbp_mean", "sbp_sd"))
data.table::setnames(dbp, old = c("mean", "sd"), new = c("dbp_mean", "dbp_sd"))
flynn2017_bp <- merge(sbp, dbp)

bp_parameters <-
  data.table::rbindlist(list(gemelli1990 = gemelli1990, nhlbi = cdc_bp, lo2013 = lo2013, flynn2017 = flynn2017_bp), idcol = "source", use.names = TRUE, fill = TRUE)

bp_parameters <- as.data.frame(bp_parameters)

################################################################################
##                                 Save Data                                  ##
save(bp_parameters, file = "./data/bp_parameters.rda")

################################################################################
##                              Save cpp version                              ##
cat("// Generated by data-raw/gaussian_parameters.R",
    "// Do not edit by hand",
    "// [[Rcpp::depends(RcppArmadillo)]]",
    "#include <RcppArmadillo.h>",
    "#include <Rcpp.h>",
    "#ifndef pedbp_bp_params_H",
    "#define pedbp_bp_params_H",
    sep = "\n",
    file = "src/bp_params.h",
    append = FALSE)

cat("// Generated by data-raw/gaussian_parameters.R",
    "// Do not edit by hand",
    "// [[Rcpp::depends(RcppArmadillo)]]",
    "#include <RcppArmadillo.h>",
    "#include <Rcpp.h>",
    "#include \"bp_params.h\"",
    sep = "\n",
    file = "src/bp_params.cpp",
    append = FALSE)

for (src in unique(bp_parameters$source)) {
  for (sex in 0:1) {
    d <- subset(bp_parameters, bp_parameters$source == src & bp_parameters$male == sex)
    d <- as.matrix(d[, c("age", "sbp_mean", "sbp_sd", "dbp_mean", "dbp_sd", "height_percentile")])
    if (all(is.na(d[, "height_percentile"]))) {
      d[, "height_percentile"] <- 101 # if there is no percentile, then everyone is in this percentile
    }
    d <- cbind(d, source = ifelse(src == "gemelli1990", 1, ifelse(src == "lo2013", 2, ifelse(src == "nhlbi", 3, ifelse(src == "flynn2017", 4, stop("unknown src"))))))
    d <- apply(d, MARGIN = 1, function(x) paste("{", paste(x, collapse = ", "), "}"))
    d <- paste(d, collapse = ", \n")
    nm <- paste0(src, "_", ifelse(sex == 0, "female", "male"))
    cat(paste0("arma::mat ", nm , "();\n"), file = "src/bp_params.h", append = TRUE)
    cat("",
        paste0("arma::mat ", nm, "() {"),
        "\tarma::mat X = {",
        d,
        "\t};",
        "\treturn X;",
        "}",
        sep = "\n",
        file = "src/bp_params.cpp",
        append = TRUE
        )
  }
}

cat("#endif", sep = "\n", file = "src/bp_params.h", append = TRUE)

################################################################################
##                                End of File                                 ##
################################################################################