R/hei_2015_PerDay_pratio.R

In yikeshu0611/nhanesR: Tools for "NHANES" Database

hei_2015_PerDay_pratio <- function(seqn=NULL,
                                   years=2009,
                                   day='1',
                                   dietary='tot',
                                   seed=NULL,
                                   varLabel=FALSE){
    if (dietary=='iff') join <- c('seqn','line') else join <- 'seqn'

    fped <- fped_read(years = years,day = day,dietary = dietary,cat = FALSE,version = 2015)

    tsv <- nhs_tsv(ifelse(day=='1',
                          sprintf('drx%s|dr1%s',dietary,dietary),
                          sprintf('drx%s|dr2%s',dietary,dietary)),years = years,cat = FALSE)
    dt <- nhs_read(tsv,'wtdr4yr','wtdrd1',
                   "drxtkcal,drxikcal,dr1tkcal,dr1ikcal,dr2tkcal,dr2ikcal:kcal",
                   "drxtsfat,drxisfat,dr1tsfat,dr1isfat,dr2tsfat,dr2isfat:sfat",
                   "drdtsodi,drdisodi,dr1tsodi,dr1isodi,dr2tsodi,dr2isodi:sodi",
                   "drxtmfat,drximfat,dr1tmfat,dr1imfat,dr2tmfat,dr2imfat:mfat",
                   "drxtpfat,drxipfat,dr1tpfat,dr1ipfat,dr2tpfat,dr2ipfat:pfat",
                   codebook = FALSE,varLabel = FALSE,cat = FALSE)
    colnames(dt) <- rename_line(colnames(dt))
    colnames(dt) <- rename_fdcd(colnames(dt))
    dt <- drop_col(dt,'fdcd')

    dt$monopoly <- dt$mfat + dt$pfat

    indat <- dplyr::inner_join(dt,fped,join)
    if (!is.null(seqn)) indat <- indat[indat$seqn %in% seqn,]
    indat <- wt_dr_day1(indat,wtname = 'wtdr',cat=FALSE)

    demo <- nhs_read(nhs_tsv('demo',years = years,cat=FALSE),'sdmvstra','sdmvpsu',cat = FALSE,Year = FALSE)

    indat <- dplyr::inner_join(indat,demo,'seqn')

    choice <- c('kcal', 'vtotalleg', 'vdrkgrleg', 'f_total', 'f_whole',
                'g_whole', 'd_total','pfallprotleg', 'pfseaplantleg',
                'monopoly', 'sfat', 'sodi', 'g_refined', 'add_sugars')
    for (i in choice) attributes(indat[,i]) <- NULL
    indat1 <- reshape2::melt(indat[,c(join,choice)],id.vars=join)
    one <- dplyr::inner_join(indat[,c(join,choice,'sdmvstra','sdmvpsu','wtdr')],indat1,join)
    head(one)
    for (i in choice) one[,i] <- ifelse(one$variable==i,1,0)

    dg <- survey::svydesign(id = ~sdmvpsu, weights = ~wtdr, strata = ~sdmvstra,
                            nest = TRUE,data = one)
    glm <- survey::svyglm(value~-1+kcal+vtotalleg+vdrkgrleg+f_total+f_whole+g_whole+d_total+pfallprotleg+pfseaplantleg+monopoly+sfat+sodi+g_refined+add_sugars,
                  design=dg)
    sigma <- glm$cov.unscaled

    dg <- survey::svydesign(id = ~sdmvpsu, weights = ~wtdr, strata = ~sdmvstra,
                            nest = TRUE, data = indat)
    mean <- survey::svymean(~kcal+vtotalleg+vdrkgrleg+f_total+f_whole+g_whole+d_total+pfallprotleg+pfseaplantleg+monopoly+sfat+sodi+g_refined+add_sugars,dg)
    mu <- as.data.frame(mean)[,1]
    if (!is.null(seed)) set.seed(seed)
    sim_data <- mvtnorm::rmvnorm(n = 10000,mean = mu,sigma = sigma) |>
        as.data.frame()
    colnames(sim_data) <- colnames(sigma)
    aftermac <- hei_2015(
        indat = sim_data,
        kcal = 'kcal',
        vtotalleg = 'vtotalleg',
        vdrkgrleg = 'vdrkgrleg',
        f_total = 'f_total',
        fwholefrt = 'f_whole',
        g_whole = 'g_whole',
        d_total = 'd_total',
        pfallprotleg = 'pfallprotleg',
        pfseaplantleg = 'pfseaplantleg',
        monopoly = 'monopoly',
        satfat = 'sfat',
        sodium = 'sodi',
        g_refined = 'g_refined',
        add_sugars = 'add_sugars',
        varLabel = FALSE,
        energy=FALSE,component=TRUE,density=FALSE,join=NULL
    )
    head(aftermac)
    data.frame(
        min = sapply(aftermac, min),
        max = sapply(aftermac, max),
        mean = sapply(aftermac, mean),
        sd = sapply(aftermac, sd),
        lowci = sapply(aftermac, function(i) quantile(i,0.025)),
        upci = sapply(aftermac, function(i) quantile(i,0.975))
    )
}