R/hei_2010_PerPerson_pratio.R

In yikeshu0611/nhanesR: Tools for "NHANES" Database

hei_2010_PerPerson_pratio <- function(seqn=NULL,
                                      years=years,
                                      dietary='tot',
                                      seed=NULL){
    if (dietary=='iff') join <- c('seqn','iline') else join <- 'seqn'

    # day=1
    fped <- fped_read(years = years,day = 1,dietary = dietary,version = 2010)
    tsv <- nhs_tsv(sprintf('drx%s|dr1%s',dietary,dietary),years = years,cat = FALSE)
    dt <- nhs_read(tsv,'wtdr4yr','wtdrd1',
                   "drxtkcal,drxikcal,dr1tkcal,dr2tkcal,dr1ikcal,dr2ikcal:kcal",
                   "drxtsfat,drxisfat,dr1tsfat,dr2tsfat,dr1isfat,dr2isfat:sfat",
                   "drxtalco,drxialco,dr1talco,dr2talco,dr1ialco,dr2ialco:alco",
                   "drdtsodi,drdisodi,dr1tsodi,dr1isodi,dr2tsodi,dr2isodi:sodi",
                   "drxtmfat,drximfat,dr1tmfat,dr2tmfat,dr1imfat,dr2imfat:mfat",
                   "drxtpfat,drxipfat,dr1tpfat,dr2tpfat,dr1ipfat,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
    maxalcgr <- 13*(dt$kcal/1000)
    dt$exalccal <- ifelse(dt$alco <= maxalcgr,0,7*( dt$alco - maxalcgr))
    indat <- dplyr::inner_join(dt,fped,join)
    indat$emptycal10 <- indat$addsugc + indat$solfatc + indat$exalccal
    indat1 <- indat

    # day=2
    fped <- fped_read(years = years,day = 2,dietary = dietary,version = 2010)
    tsv <- nhs_tsv(sprintf('drx%s|dr2%s',dietary,dietary),years = years,cat = FALSE)
    dt <- nhs_read(tsv,
                   "drxtkcal,drxikcal,dr1tkcal,dr2tkcal,dr1ikcal,dr2ikcal:kcal",
                   "drxtsfat,drxisfat,dr1tsfat,dr2tsfat,dr1isfat,dr2isfat:sfat",
                   "drxtalco,drxialco,dr1talco,dr2talco,dr1ialco,dr2ialco:alco",
                   "drdtsodi,drdisodi,dr1tsodi,dr1isodi,dr2tsodi,dr2isodi:sodi",
                   "drxtmfat,drximfat,dr1tmfat,dr2tmfat,dr1imfat,dr2imfat:mfat",
                   "drxtpfat,drxipfat,dr1tpfat,dr2tpfat,dr1ipfat,dr2ipfat:pfat",
                   codebook = FALSE,varLabel = FALSE,cat = FALSE,Year = FALSE)

    colnames(dt) <- rename_line(colnames(dt))
    colnames(dt) <- rename_fdcd(colnames(dt))
    dt <- drop_col(dt,'fdcd')

    dt$monopoly <- dt$mfat + dt$pfat
    maxalcgr <- 13*(dt$kcal/1000)
    dt$exalccal <- ifelse(dt$alco <= maxalcgr,0,7*( dt$alco - maxalcgr))
    indat <- dplyr::inner_join(dt,fped,join)
    indat$emptycal10 <- indat$addsugc + indat$solfatc + indat$exalccal
    indat2 <- indat


    indat <- dplyr::inner_join(indat1,indat2,join)

    indat <- wt_dr_day1(indat,wtname = 'wtdr',cat=FALSE)

    if (!is.null(seqn)) indat <- indat[indat$seqn %in% seqn,]

    choice <- colnames(indat2)[-c(1:length(join))]
    for (i in choice) {
        which <- which(colnames(indat) %in% paste0(i,c('.x','.y')))
        indat$last <- rowSums(indat[,which],na.rm = TRUE)
        indat <- indat[,-which]
        colnames(indat)[ncol(indat)] <- i
    }
    colnames(indat) |> do::increase()
    ck <- indat$Year %in% c('1999-2000','2001-2002')
    indat[ck,choice] <- indat[ck,choice]/2 # these four years is the same

    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", "v_total", "v_drkgr", "v_legumes",
                "f_total", "f_whole", "g_whole", "d_total", "pf_total",
                "seaplant", "monopoly", "sfat", "sodi", "g_refined",
                "emptycal10")
    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+v_total+v_drkgr+v_legumes+f_total+f_whole+g_whole+d_total+pf_total+seaplant+monopoly+sfat+sodi+g_refined+emptycal10,
                          design=dg)
    sigma <- glm$cov.unscaled

    dg <- survey::svydesign(id = ~sdmvpsu, weights = ~wtdr, strata = ~sdmvstra,
                            nest = TRUE, data = indat)
    mean <- survey::svymean(~kcal+v_total+v_drkgr+v_legumes+f_total+f_whole+g_whole+d_total+pf_total+seaplant+monopoly+sfat+sodi+g_refined+emptycal10,dg)
    mu <- as.data.frame(mean)[,1]
    if (!is.null(seed)) set.seed(seed)
    indat <- mvtnorm::rmvnorm(n = 10000,mean = mu,sigma = sigma) |>
        as.data.frame()
    colnames(indat) <- colnames(sigma)



    # leg2010a
    afterleg <- leg2010a(indat = indat,
                         kcal = 'kcal',
                         allmeat ='pf_total',
                         seaplant = 'seaplant',
                         v_total = 'v_total',
                         v_drkgr = 'v_drkgr',
                         legumes = 'v_legumes')
    aftermac <- hei_2010(
        indat=afterleg,
        kcal = 'kcal',
        lv_total = 'legume_added_v_total',
        lbeangrn = 'legume_added_beangrn',
        f_total = 'f_total',
        wholefrt = 'f_whole',
        g_whl = 'g_whole',
        d_total = 'd_total',
        lallmeat = 'legume_added_allmeat',
        lseaplant = 'legume_added_seaplant',
        monopoly = 'monopoly',
        sfat = 'sfat',
        sodi = 'sodi',
        g_nwhl = 'g_refined',
        emptycal10 = 'emptycal10',
        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))
    )
}