Obesity criteria: NHANES data

1. NHANES data (4 cycles)

1.1 Data management

load("C:/Thach/Research projects/BMI_NHANES/w - create_merged_nhanes_with_dxa.RData")

Libraries

suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(dplyr))  
suppressPackageStartupMessages(library(expss))

## Warning: package 'expss' was built under R version 4.3.2

## Warning: package 'maditr' was built under R version 4.3.2

suppressPackageStartupMessages(library(table1)) 
suppressPackageStartupMessages(library(compareGroups)) 
suppressPackageStartupMessages(library(ggplot2)) 
suppressPackageStartupMessages(library(gridExtra))
suppressPackageStartupMessages(library(grid))
suppressPackageStartupMessages(library(ggthemes))

## Warning: package 'ggthemes' was built under R version 4.3.2

Merge datasets

demo = subset(demographics_clean, select = c(SEQN, SDDSRVYR, RIDRETH3))
whr = subset(response_clean, select = c(SEQN, SDDSRVYR, BMXWAIST))
  temp1 = merge(demo, whr, by = c("SEQN", "SDDSRVYR"))
dm = subset(questionnaire_clean, select = c(SEQN, SDDSRVYR, DIQ010, DIQ160))
  temp2 = merge(temp1, dm, by = c("SEQN", "SDDSRVYR"))
ob = merge(df_merged_nhanes, temp2, by = c("SEQN", "SDDSRVYR"))
ob = subset(ob, RIDAGEYR >= 20 & RIDAGEYR <= 59)
head(ob)

##      SEQN SDDSRVYR SEQN_new RIDRETH1 RIAGENDR RIDAGEYR BMXHT BMXWT BMXBMI
## 1  100001       10 C-100001        3        1       39 176.4  96.0   30.9
## 5  100008       10 C-100008        5        2       56 156.5  73.3   29.9
## 7  100012       10 C-100012        3        1       51 170.2 130.7   45.1
## 8  100013       10 C-100013        5        2       43 153.6  59.0   25.0
## 10 100015       10 C-100015        3        1       27 171.8  69.8   23.6
## 11 100016       10 C-100016        5        1       36 167.4  82.8   29.5
##    VNALCOHOLPERYEAR  LBXCOT DXDTOBMD DXDTOBMC DXDTOLE DXDTOLI DXDTOFAT DXDTOPF
## 1                NA 220.000    1.075  2536.13 62147.1 64683.3  30550.0    32.1
## 5                NA   0.016    1.024  2028.52 42335.1 44363.6  29061.5    39.6
## 7                NA   0.011    1.084  2641.44 80905.5 83546.9  46193.1    35.6
## 8                NA   0.011    1.096  1948.22 33436.1 35384.3  24176.7    40.6
## 10               NA 320.000    1.242  2662.50 51456.9 54119.4  16255.1    23.1
## 11               NA   0.035    1.162  2528.12 52889.6 55417.7  28026.7    33.6
##    RIDRETH3 BMXWAIST DIQ010 DIQ160
## 1         3    102.5      2      2
## 5         6    101.1      2      2
## 7         3    140.9      1     NA
## 8         6     86.2      2      2
## 10        3     87.0      2      2
## 11        6    100.6      2      2

dim(ob)

## [1] 14346    21

Subgroup definition:

ob$sex <- factor(ob$RIAGENDR,
                 levels = c(1, 2),
                 labels = c("Male", "Female"))
var_lab(ob$sex) = "Sex"

ob = ob %>%
  mutate(age.grp = case_when(RIDAGEYR >= 20 & RIDAGEYR <= 29 ~ "20-29",
                             RIDAGEYR >= 30 & RIDAGEYR <= 39 ~ "30-39",
                             RIDAGEYR >= 40 & RIDAGEYR <= 49 ~ "40-49",
                             RIDAGEYR >= 50 & RIDAGEYR <= 59 ~ "50-59"))
var_lab(ob$age.grp) = "Age group"

ob = ob %>% 
  mutate(race = case_when(RIDRETH3 == 1 | RIDRETH3 == 2 ~ "Hispanic",
                          RIDRETH3 == 3 ~ "Non-Hispanic White",
                          RIDRETH3 == 4 ~ "Non-Hispanic Black",
                          RIDRETH3 == 6 ~ "Non-Hispanic Asian",
                          RIDRETH3 == 7 ~ "Other Race, incl. Multi-Racial"))
ob$race = factor(ob$race, levels = c("Hispanic", "Non-Hispanic Asian", "Non-Hispanic White", "Non-Hispanic Black", "OtherRace, incl. Multi-Racial"))
var_lab(ob$race) = "Race and ethnicity"

var_lab(ob$SDDSRVYR) = "NHANES cycles"
val_lab(ob$SDDSRVYR) = num_lab("
                              7 2011-2012
                              8 2013-2014
                              9 2015-2016
                              10 2017-2018")

Obesity definition:

Obesity (Rubio et al. Lancet Diabetes Endocrinol 2025): (i) BMI (30 kg/m2 or 27.5 kg/m2 [Asian]); (ii) BMI and (percentage fat [25% in M, 35% in W] or WHR [102cm in M, 88cm in W /90, 80 for Asian])

# Obesity defined by only BMI
ob$ob.bmi = ifelse((ob$BMXBMI >= 30 & ob$RIDRETH3 !=6) | (ob$BMXBMI >= 27.5 & ob$RIDRETH3 == 6), 1, 0)
ob$ob.bmi <- factor(ob$ob.bmi, levels = c(0, 1), labels = c("Non-obese_BMI", "Obese_BMI"))
var_lab(ob$ob.bmi) = "Obesity_only BMI"
  table(ob$ob.bmi)

## 
## Non-obese_BMI     Obese_BMI 
##          8392          5815

# Temporal:  
ob$perc.fat = ifelse((ob$DXDTOPF > 25 & ob$RIAGENDR == 1) | (ob$DXDTOPF > 35 & ob$RIAGENDR == 2), 1, 0)  
  table(ob$perc.fat)

## 
##    0    1 
## 3350 7513

ob$waist = ifelse((ob$BMXWAIST > 102 & ob$RIAGENDR == 1 & ob$RIDRETH3 != 6) | (ob$BMXWAIST > 88 & ob$RIAGENDR == 2 & ob$RIDRETH3 != 6) | 
                  (ob$BMXWAIST > 90 & ob$RIAGENDR == 1 & ob$RIDRETH3 == 6) | (ob$BMXWAIST > 80 & ob$RIAGENDR == 2 & ob$RIDRETH3 == 6), 1, 0)
  table(ob$waist)

## 
##    0    1 
## 5838 7836

ob$excess.fat = ifelse(ob$perc.fat == 1 | ob$waist == 1, 1, 0)
ob$excess.fat2 = ifelse((ob$perc.fat == 1 | (ob$waist == 1 | is.na(ob$waist))) | 
                          ((ob$perc.fat == 1 | is.na(ob$perc.fat)) | ob$waist == 1), 1, 0)
ob$excess.fat2[(ob$perc.fat == 0 & is.na(ob$waist)) | 
                 (is.na(ob$perc.fat) & ob$waist == 0)] = NA
table(ob$excess.fat, ob$excess.fat2)

##    
##        0    1
##   0 2941    0
##   1    0 9793

  ### Exactly the same
  
# Obesity defined by BMI + (percentage fat or waist)
ob$ob.bmi_fat = ifelse((ob$ob.bmi == "Obese_BMI" & (ob$perc.fat == 1 | ob$waist == 1)) | 
                         ob$BMXBMI >= 40, 1, 0)
ob$ob.bmi_fat <- factor(ob$ob.bmi_fat, levels = c(0, 1), labels = c("Non-obese_BMI.Fat", "Obese_BMI.Fat"))
var_lab(ob$ob.bmi_fat) = "Obesity_BMI.Fat"
  table(ob$ob.bmi_fat)   

## 
## Non-obese_BMI.Fat     Obese_BMI.Fat 
##              8454              5591

# Check missing data
missing = subset(ob, is.na(BMXBMI) | is.na(perc.fat) | is.na(waist), select = c(SEQN, SDDSRVYR, BMXBMI, perc.fat, waist, ob.bmi, ob.bmi_fat))
missing_non.obese = subset(missing, ob.bmi == "Obese_BMI" & (is.na(perc.fat) | is.na(waist)))
### Correct data management (12/7/2025)

1.2 Characteristics by NHANS cycles

table1(~ sex + RIDAGEYR + age.grp + race + BMXWT + BMXHT + BMXBMI + BMXWAIST + DXDTOBMD + DXDTOFAT + DXDTOPF + ob.bmi + as.factor(ob$excess.fat) + ob.bmi_fat + as.factor(perc.fat) + as.factor(waist) | SDDSRVYR, render.continuous = c(. = "Mean (SD)", . = "Median [Q1, Q3]", . = "Median [min, max]"), data = ob)

## Warning in table1.formula(~sex + RIDAGEYR + age.grp + race + BMXWT + BMXHT + :
## Terms to the right of '|' in formula 'x' define table columns and are expected
## to be factors with meaningful labels.

	2011-2012 (N=3632)	2013-2014 (N=3803)	2015-2016 (N=3665)	2017-2018 (N=3246)	Overall (N=14346)
Sex
Male	1781 (49.0%)	1821 (47.9%)	1727 (47.1%)	1519 (46.8%)	6848 (47.7%)
Female	1851 (51.0%)	1982 (52.1%)	1938 (52.9%)	1727 (53.2%)	7498 (52.3%)
RIDAGEYR
Mean (SD)	38.9 (11.5)	39.4 (11.4)	39.3 (11.4)	40.0 (11.6)	39.4 (11.5)
Median [Q1, Q3]	39.0 [29.0, 49.0]	40.0 [30.0, 49.0]	39.0 [29.0, 49.0]	40.0 [30.0, 50.0]	39.0 [30.0, 49.0]
Median [min, max]	39.0 [20.0, 59.0]	40.0 [20.0, 59.0]	39.0 [20.0, 59.0]	40.0 [20.0, 59.0]	39.0 [20.0, 59.0]
Age group
20-29	954 (26.3%)	919 (24.2%)	931 (25.4%)	775 (23.9%)	3579 (24.9%)
30-39	924 (25.4%)	961 (25.3%)	933 (25.5%)	813 (25.0%)	3631 (25.3%)
40-49	875 (24.1%)	1007 (26.5%)	913 (24.9%)	778 (24.0%)	3573 (24.9%)
50-59	879 (24.2%)	916 (24.1%)	888 (24.2%)	880 (27.1%)	3563 (24.8%)
Race and ethnicity
Hispanic	748 (20.6%)	883 (23.2%)	1108 (30.2%)	799 (24.6%)	3538 (24.7%)
Non-Hispanic Asian	582 (16.0%)	478 (12.6%)	510 (13.9%)	544 (16.8%)	2114 (14.7%)
Non-Hispanic White	1241 (34.2%)	1525 (40.1%)	1089 (29.7%)	971 (29.9%)	4826 (33.6%)
Non-Hispanic Black	943 (26.0%)	772 (20.3%)	811 (22.1%)	744 (22.9%)	3270 (22.8%)
OtherRace, incl. Multi-Racial	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
Missing	118 (3.2%)	145 (3.8%)	147 (4.0%)	188 (5.8%)	598 (4.2%)
BMXWT
Mean (SD)	81.7 (22.0)	82.7 (23.1)	82.7 (22.7)	84.2 (24.3)	82.8 (23.0)
Median [Q1, Q3]	78.4 [65.8, 93.1]	79.0 [66.6, 94.1]	79.2 [66.7, 94.7]	80.0 [67.3, 97.2]	79.1 [66.6, 94.6]
Median [min, max]	78.4 [34.7, 216]	79.0 [32.8, 223]	79.2 [39.0, 199]	80.0 [37.5, 243]	79.1 [32.8, 243]
Missing	39 (1.1%)	29 (0.8%)	29 (0.8%)	31 (1.0%)	128 (0.9%)
BMXHT
Mean (SD)	168 (9.96)	168 (9.99)	167 (10.0)	167 (10.1)	168 (10.0)
Median [Q1, Q3]	168 [161, 175]	168 [161, 175]	167 [160, 174]	167 [160, 175]	167 [160, 175]
Median [min, max]	168 [135, 205]	168 [141, 203]	167 [137, 203]	167 [138, 196]	167 [135, 205]
Missing	36 (1.0%)	28 (0.7%)	29 (0.8%)	31 (1.0%)	124 (0.9%)
BMXBMI
Mean (SD)	28.8 (7.12)	29.2 (7.52)	29.6 (7.34)	30.0 (7.85)	29.4 (7.46)
Median [Q1, Q3]	27.6 [23.7, 32.2]	27.8 [23.8, 32.8]	28.4 [24.3, 33.5]	28.7 [24.5, 33.9]	28.1 [24.0, 33.1]
Median [min, max]	27.6 [13.6, 80.6]	27.8 [14.1, 82.9]	28.4 [14.5, 64.6]	28.7 [14.8, 86.2]	28.1 [13.6, 86.2]
Missing	41 (1.1%)	32 (0.8%)	31 (0.8%)	35 (1.1%)	139 (1.0%)
BMXWAIST
Mean (SD)	97.0 (16.9)	98.0 (17.2)	98.9 (17.2)	99.3 (18.0)	98.3 (17.3)
Median [Q1, Q3]	95.3 [84.3, 107]	95.8 [85.6, 108]	97.0 [86.4, 109]	97.4 [86.2, 110]	96.4 [85.6, 108]
Median [min, max]	95.3 [61.9, 176]	95.8 [55.5, 178]	97.0 [58.7, 172]	97.4 [57.9, 170]	96.4 [55.5, 178]
Missing	175 (4.8%)	150 (3.9%)	188 (5.1%)	159 (4.9%)	672 (4.7%)
DXDTOBMD
Mean (SD)	1.12 (0.109)	1.11 (0.111)	1.10 (0.113)	1.12 (0.107)	1.11 (0.110)
Median [Q1, Q3]	1.11 [1.05, 1.18]	1.11 [1.04, 1.18]	1.10 [1.03, 1.17]	1.11 [1.05, 1.18]	1.11 [1.04, 1.18]
Median [min, max]	1.11 [0.762, 1.88]	1.11 [0.729, 1.55]	1.10 [0.733, 1.52]	1.11 [0.697, 1.73]	1.11 [0.697, 1.88]
Missing	914 (25.2%)	679 (17.9%)	865 (23.6%)	989 (30.5%)	3447 (24.0%)
DXDTOFAT
Mean (SD)	27000 (11900)	27600 (11800)	27800 (12100)	27700 (11600)	27500 (11900)
Median [Q1, Q3]	24800 [18600, 32600]	25300 [19200, 33500]	25500 [19200, 34200]	25700 [19300, 34200]	25300 [19100, 33600]
Median [min, max]	24800 [5910, 102000]	25300 [4900, 89300]	25500 [7230, 77500]	25700 [6670, 86900]	25300 [4900, 102000]
Missing	917 (25.2%)	689 (18.1%)	880 (24.0%)	997 (30.7%)	3483 (24.3%)
DXDTOPF
Mean (SD)	32.4 (8.69)	33.1 (8.37)	33.1 (8.65)	33.5 (8.56)	33.0 (8.57)
Median [Q1, Q3]	31.9 [25.9, 39.3]	32.8 [26.8, 39.7]	33.2 [26.6, 40.2]	33.4 [27.3, 40.5]	32.8 [26.6, 39.9]
Median [min, max]	31.9 [12.6, 55.4]	32.8 [12.8, 54.7]	33.2 [11.7, 54.7]	33.4 [12.1, 56.1]	32.8 [11.7, 56.1]
Missing	917 (25.2%)	689 (18.1%)	880 (24.0%)	997 (30.7%)	3483 (24.3%)
Obesity_only BMI
Non-obese_BMI	2262 (62.3%)	2288 (60.2%)	2096 (57.2%)	1746 (53.8%)	8392 (58.5%)
Obese_BMI	1329 (36.6%)	1483 (39.0%)	1538 (42.0%)	1465 (45.1%)	5815 (40.5%)
Missing	41 (1.1%)	32 (0.8%)	31 (0.8%)	35 (1.1%)	139 (1.0%)
as.factor(ob$excess.fat)
0	813 (22.4%)	826 (21.7%)	726 (19.8%)	576 (17.7%)	2941 (20.5%)
1	2364 (65.1%)	2642 (69.5%)	2516 (68.6%)	2271 (70.0%)	9793 (68.3%)
Missing	455 (12.5%)	335 (8.8%)	423 (11.5%)	399 (12.3%)	1612 (11.2%)
Obesity_BMI.Fat
Non-obese_BMI.Fat	2271 (62.5%)	2305 (60.6%)	2118 (57.8%)	1760 (54.2%)	8454 (58.9%)
Obese_BMI.Fat	1279 (35.2%)	1434 (37.7%)	1471 (40.1%)	1407 (43.3%)	5591 (39.0%)
Missing	82 (2.3%)	64 (1.7%)	76 (2.1%)	79 (2.4%)	301 (2.1%)
as.factor(perc.fat)
0	911 (25.1%)	950 (25.0%)	835 (22.8%)	654 (20.1%)	3350 (23.4%)
1	1804 (49.7%)	2164 (56.9%)	1950 (53.2%)	1595 (49.1%)	7513 (52.4%)
Missing	917 (25.2%)	689 (18.1%)	880 (24.0%)	997 (30.7%)	3483 (24.3%)
as.factor(waist)
0	1594 (43.9%)	1590 (41.8%)	1451 (39.6%)	1203 (37.1%)	5838 (40.7%)
1	1863 (51.3%)	2063 (54.2%)	2026 (55.3%)	1884 (58.0%)	7836 (54.6%)
Missing	175 (4.8%)	150 (3.9%)	188 (5.1%)	159 (4.9%)	672 (4.7%)

table(ob$excess.fat, ob$ob.bmi)

##    
##     Non-obese_BMI Obese_BMI
##   0          2878        62
##   1          4240      5520

table(ob$excess.fat[ob$SDDSRVYR == 7], ob$ob.bmi[ob$SDDSRVYR == 7])

##    
##     Non-obese_BMI Obese_BMI
##   0           803         9
##   1          1088      1268

table(ob$excess.fat[ob$SDDSRVYR == 8], ob$ob.bmi[ob$SDDSRVYR == 8])

##    
##     Non-obese_BMI Obese_BMI
##   0           809        17
##   1          1221      1410

table(ob$excess.fat[ob$SDDSRVYR == 9], ob$ob.bmi[ob$SDDSRVYR == 9])

##    
##     Non-obese_BMI Obese_BMI
##   0           704        22
##   1          1057      1453

table(ob$excess.fat[ob$SDDSRVYR == 10], ob$ob.bmi[ob$SDDSRVYR == 10])

##    
##     Non-obese_BMI Obese_BMI
##   0           562        14
##   1           874      1389

1.3 Prevalence of obesity by BMI only

Calculate %obesity (#obesity/total) by sex, age group, ethnicity and wave

# Overall
ob %>%
  group_by(SDDSRVYR) %>%
  summarize(
    n = n(),
    ob_bmi = sum(ob.bmi == "Obese_BMI", na.rm = TRUE),
    prev.ob_bmi = ob_bmi / n,
    se.ob_bmi = sqrt(prev.ob_bmi * (1 - prev.ob_bmi) / n),
    ob_bmi.fat = sum(ob.bmi_fat == "Obese_BMI.Fat", na.rm = TRUE),
    prev.ob_bmi.fat = ob_bmi.fat / n,
    se.ob_bmi,fat = sqrt(prev.ob_bmi.fat * (1 - prev.ob_bmi.fat) / n),
    .groups = "drop"
  )

## # A tibble: 4 × 8
##   SDDSRVYR     n ob_bmi prev.ob_bmi se.ob_bmi ob_bmi.fat prev.ob_bmi.fat     fat
##   <labell> <int>  <int>       <dbl>     <dbl>      <int>           <dbl>   <dbl>
## 1  7        3632   1329       0.366   0.00799       1279           0.352 0.00793
## 2  8        3803   1483       0.390   0.00791       1434           0.377 0.00786
## 3  9        3665   1538       0.420   0.00815       1471           0.401 0.00810
## 4 10        3246   1465       0.451   0.00873       1407           0.433 0.00870

# By Sex
ob %>%
  group_by(SDDSRVYR, sex) %>%
  summarize(
    n = n(),
    ob_bmi = sum(ob.bmi == "Obese_BMI", na.rm = TRUE),
    prev.ob_bmi = ob_bmi / n,
    se.ob_bmi = sqrt(prev.ob_bmi * (1 - prev.ob_bmi) / n),
    ob_bmi.fat = sum(ob.bmi_fat == "Obese_BMI.Fat", na.rm = TRUE),
    prev.ob_bmi.fat = ob_bmi.fat / n,
    se.ob_bmi,fat = sqrt(prev.ob_bmi.fat * (1 - prev.ob_bmi.fat) / n),
    .groups = "drop"
  )

## # A tibble: 8 × 9
##   SDDSRVYR   sex       n ob_bmi prev.ob_bmi se.ob_bmi ob_bmi.fat prev.ob_bmi.fat
##   <labelled> <fct> <int>  <int>       <dbl>     <dbl>      <int>           <dbl>
## 1  7         Male   1781    592       0.332    0.0112        561           0.315
## 2  7         Fema…  1851    737       0.398    0.0114        718           0.388
## 3  8         Male   1821    634       0.348    0.0112        601           0.330
## 4  8         Fema…  1982    849       0.428    0.0111        833           0.420
## 5  9         Male   1727    667       0.386    0.0117        622           0.360
## 6  9         Fema…  1938    871       0.449    0.0113        849           0.438
## 7 10         Male   1519    671       0.442    0.0127        638           0.420
## 8 10         Fema…  1727    794       0.460    0.0120        769           0.445
## # ℹ 1 more variable: fat <dbl>

# By Age group
ob %>%
  group_by(SDDSRVYR, age.grp) %>%
  summarize(
    n = n(),
    ob_bmi = sum(ob.bmi == "Obese_BMI", na.rm = TRUE),
    prev.ob_bmi = ob_bmi / n,
    se.ob_bmi = sqrt(prev.ob_bmi * (1 - prev.ob_bmi) / n),
    ob_bmi.fat = sum(ob.bmi_fat == "Obese_BMI.Fat", na.rm = TRUE),
    prev.ob_bmi.fat = ob_bmi.fat / n,
    se.ob_bmi,fat = sqrt(prev.ob_bmi.fat * (1 - prev.ob_bmi.fat) / n),
    .groups = "drop"
  )

## # A tibble: 16 × 9
##    SDDSRVYR   age.grp        n ob_bmi prev.ob_bmi se.ob_bmi ob_bmi.fat
##    <labelled> <labelled> <int>  <int>       <dbl>     <dbl>      <int>
##  1  7         20-29        954    287       0.301    0.0148        275
##  2  7         30-39        924    316       0.342    0.0156        301
##  3  7         40-49        875    359       0.410    0.0166        346
##  4  7         50-59        879    367       0.418    0.0166        357
##  5  8         20-29        919    277       0.301    0.0151        262
##  6  8         30-39        961    402       0.418    0.0159        390
##  7  8         40-49       1007    397       0.394    0.0154        384
##  8  8         50-59        916    407       0.444    0.0164        398
##  9  9         20-29        931    316       0.339    0.0155        295
## 10  9         30-39        933    395       0.423    0.0162        378
## 11  9         40-49        913    431       0.472    0.0165        411
## 12  9         50-59        888    396       0.446    0.0167        387
## 13 10         20-29        775    301       0.388    0.0175        289
## 14 10         30-39        813    367       0.451    0.0175        345
## 15 10         40-49        778    387       0.497    0.0179        376
## 16 10         50-59        880    410       0.466    0.0168        397
## # ℹ 2 more variables: prev.ob_bmi.fat <dbl>, fat <dbl>

# By Race
ob %>%
  group_by(SDDSRVYR, race) %>%
  summarize(
    n = n(),
    ob_bmi = sum(ob.bmi == "Obese_BMI", na.rm = TRUE),
    prev.ob_bmi = ob_bmi / n,
    se.ob_bmi = sqrt(prev.ob_bmi * (1 - prev.ob_bmi) / n),
    ob_bmi.fat = sum(ob.bmi_fat == "Obese_BMI.Fat", na.rm = TRUE),
    prev.ob_bmi.fat = ob_bmi.fat / n,
    se.ob_bmi.fat = sqrt(prev.ob_bmi.fat * (1 - prev.ob_bmi.fat) / n),
    .groups = "drop"
  )

## # A tibble: 20 × 9
##    SDDSRVYR  race      n ob_bmi prev.ob_bmi se.ob_bmi ob_bmi.fat prev.ob_bmi.fat
##    <labelle> <fct> <int>  <int>       <dbl>     <dbl>      <int>           <dbl>
##  1  7        Hisp…   748    310       0.414    0.0180        300           0.401
##  2  7        Non-…   582    119       0.204    0.0167        111           0.191
##  3  7        Non-…  1241    428       0.345    0.0135        418           0.337
##  4  7        Non-…   943    435       0.461    0.0162        413           0.438
##  5  7        <NA>    118     37       0.314    0.0427         37           0.314
##  6  8        Hisp…   883    385       0.436    0.0167        372           0.421
##  7  8        Non-…   478    113       0.236    0.0194        108           0.226
##  8  8        Non-…  1525    559       0.367    0.0123        549           0.36 
##  9  8        Non-…   772    367       0.475    0.0180        349           0.452
## 10  8        <NA>    145     59       0.407    0.0408         56           0.386
## 11  9        Hisp…  1108    531       0.479    0.0150        512           0.462
## 12  9        Non-…   510    129       0.253    0.0192        118           0.231
## 13  9        Non-…  1089    416       0.382    0.0147        402           0.369
## 14  9        Non-…   811    389       0.480    0.0175        369           0.455
## 15  9        <NA>    147     73       0.497    0.0412         70           0.476
## 16 10        Hisp…   799    364       0.456    0.0176        345           0.432
## 17 10        Non-…   544    190       0.349    0.0204        171           0.314
## 18 10        Non-…   971    429       0.442    0.0159        421           0.434
## 19 10        Non-…   744    385       0.517    0.0183        376           0.505
## 20 10        <NA>    188     97       0.516    0.0364         94           0.5  
## # ℹ 1 more variable: se.ob_bmi.fat <dbl>

1.4 Prevalence of confirmed excess adiposity among obesity by BMI only

By sex, age group, ethnicity and wave

ob.bmi = subset(ob, ob.bmi == "Obese_BMI")

# Overall
ob.bmi %>%
  group_by(SDDSRVYR) %>%
  summarize(
    n = n(),
    excess_fat = sum(excess.fat == 1, na.rm = TRUE),
    prev.excess_fat = excess_fat / n,
    se.excess_fat = sqrt(prev.excess_fat * (1 - prev.excess_fat) / n),
    .groups = "drop"
  )

## # A tibble: 4 × 5
##   SDDSRVYR       n excess_fat prev.excess_fat se.excess_fat
##   <labelled> <int>      <int>           <dbl>         <dbl>
## 1  7          1329       1268           0.954       0.00574
## 2  8          1483       1410           0.951       0.00562
## 3  9          1538       1453           0.945       0.00583
## 4 10          1465       1389           0.948       0.00579

# By Sex
ob.bmi %>%
  group_by(SDDSRVYR, sex) %>%
  summarize(
    n = n(),
    excess_fat = sum(excess.fat == 1, na.rm = TRUE),
    prev.excess_fat = excess_fat / n,
    se.excess_fat = sqrt(prev.excess_fat * (1 - prev.excess_fat) / n),
    .groups = "drop"
  )

## # A tibble: 8 × 6
##   SDDSRVYR   sex        n excess_fat prev.excess_fat se.excess_fat
##   <labelled> <fct>  <int>      <int>           <dbl>         <dbl>
## 1  7         Male     592        557           0.941       0.00969
## 2  7         Female   737        711           0.965       0.00680
## 3  8         Male     634        594           0.937       0.00966
## 4  8         Female   849        816           0.961       0.00663
## 5  9         Male     667        619           0.928       0.0100 
## 6  9         Female   871        834           0.958       0.00683
## 7 10         Male     671        629           0.937       0.00935
## 8 10         Female   794        760           0.957       0.00718

# By Age group
ob.bmi %>%
  group_by(SDDSRVYR, age.grp) %>%
  summarize(
    n = n(),
    excess_fat = sum(excess.fat == 1, na.rm = TRUE),
    prev.excess_fat = excess_fat / n,
    se.excess_fat = sqrt(prev.excess_fat * (1 - prev.excess_fat) / n),
    .groups = "drop"
  )

## # A tibble: 16 × 6
##    SDDSRVYR   age.grp        n excess_fat prev.excess_fat se.excess_fat
##    <labelled> <labelled> <int>      <int>           <dbl>         <dbl>
##  1  7         20-29        287        272           0.948       0.0131 
##  2  7         30-39        316        299           0.946       0.0127 
##  3  7         40-49        359        341           0.950       0.0115 
##  4  7         50-59        367        356           0.970       0.00890
##  5  8         20-29        277        258           0.931       0.0152 
##  6  8         30-39        402        379           0.943       0.0116 
##  7  8         40-49        397        380           0.957       0.0102 
##  8  8         50-59        407        393           0.966       0.00903
##  9  9         20-29        316        294           0.930       0.0143 
## 10  9         30-39        395        372           0.942       0.0118 
## 11  9         40-49        431        403           0.935       0.0119 
## 12  9         50-59        396        384           0.970       0.00861
## 13 10         20-29        301        285           0.947       0.0129 
## 14 10         30-39        367        338           0.921       0.0141 
## 15 10         40-49        387        374           0.966       0.00916
## 16 10         50-59        410        392           0.956       0.0101

# By Race
ob.bmi %>%
  group_by(SDDSRVYR, race) %>%
  summarize(
    n = n(),
    excess_fat = sum(excess.fat == 1, na.rm = TRUE),
    prev.excess_fat = excess_fat / n,
    se.excess_fat = sqrt(prev.excess_fat * (1 - prev.excess_fat) / n),
    .groups = "drop"
  )

## # A tibble: 20 × 6
##    SDDSRVYR   race                   n excess_fat prev.excess_fat se.excess_fat
##    <labelled> <fct>              <int>      <int>           <dbl>         <dbl>
##  1  7         Hispanic             310        298           0.961       0.0110 
##  2  7         Non-Hispanic Asian   119        111           0.933       0.0230 
##  3  7         Non-Hispanic White   428        415           0.970       0.00830
##  4  7         Non-Hispanic Black   435        408           0.938       0.0116 
##  5  7         <NA>                  37         36           0.973       0.0267 
##  6  8         Hispanic             385        366           0.951       0.0110 
##  7  8         Non-Hispanic Asian   113        108           0.956       0.0193 
##  8  8         Non-Hispanic White   559        545           0.975       0.00661
##  9  8         Non-Hispanic Black   367        338           0.921       0.0141 
## 10  8         <NA>                  59         53           0.898       0.0393 
## 11  9         Hispanic             531        508           0.957       0.00883
## 12  9         Non-Hispanic Asian   129        117           0.907       0.0256 
## 13  9         Non-Hispanic White   416        400           0.962       0.00943
## 14  9         Non-Hispanic Black   389        359           0.923       0.0135 
## 15  9         <NA>                  73         69           0.945       0.0266 
## 16 10         Hispanic             364        343           0.942       0.0122 
## 17 10         Non-Hispanic Asian   190        170           0.895       0.0223 
## 18 10         Non-Hispanic White   429        414           0.965       0.00887
## 19 10         Non-Hispanic Black   385        369           0.958       0.0102 
## 20 10         <NA>                  97         93           0.959       0.0202

1.4b Sensitivity analysis for a subgroup of participants with all data available

ob.all = subset(ob, !is.na(ob.bmi) & !is.na(ob.bmi_fat))
table1(~ RIDAGEYR + sex + age.grp + race + ob.bmi + ob.bmi_fat | SDDSRVYR, data = ob.all)

## Warning in table1.formula(~RIDAGEYR + sex + age.grp + race + ob.bmi +
## ob.bmi_fat | : Terms to the right of '|' in formula 'x' define table columns
## and are expected to be factors with meaningful labels.

	2011-2012 (N=3550)	2013-2014 (N=3739)	2015-2016 (N=3589)	2017-2018 (N=3167)	Overall (N=14045)
RIDAGEYR
Mean (SD)	38.9 (11.6)	39.4 (11.4)	39.3 (11.4)	40.0 (11.6)	39.4 (11.5)
Median [Min, Max]	39.0 [20.0, 59.0]	40.0 [20.0, 59.0]	39.0 [20.0, 59.0]	40.0 [20.0, 59.0]	39.0 [20.0, 59.0]
Sex
Male	1745 (49.2%)	1790 (47.9%)	1688 (47.0%)	1479 (46.7%)	6702 (47.7%)
Female	1805 (50.8%)	1949 (52.1%)	1901 (53.0%)	1688 (53.3%)	7343 (52.3%)
Age group
20-29	939 (26.5%)	906 (24.2%)	911 (25.4%)	759 (24.0%)	3515 (25.0%)
30-39	900 (25.4%)	951 (25.4%)	909 (25.3%)	790 (24.9%)	3550 (25.3%)
40-49	851 (24.0%)	985 (26.3%)	892 (24.9%)	758 (23.9%)	3486 (24.8%)
50-59	860 (24.2%)	897 (24.0%)	877 (24.4%)	860 (27.2%)	3494 (24.9%)
Race and ethnicity
Hispanic	730 (20.6%)	864 (23.1%)	1087 (30.3%)	770 (24.3%)	3451 (24.6%)
Non-Hispanic Asian	567 (16.0%)	469 (12.5%)	497 (13.8%)	520 (16.4%)	2053 (14.6%)
Non-Hispanic White	1221 (34.4%)	1511 (40.4%)	1069 (29.8%)	962 (30.4%)	4763 (33.9%)
Non-Hispanic Black	914 (25.7%)	752 (20.1%)	791 (22.0%)	731 (23.1%)	3188 (22.7%)
OtherRace, incl. Multi-Racial	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
Missing	118 (3.3%)	143 (3.8%)	145 (4.0%)	184 (5.8%)	590 (4.2%)
Obesity_only BMI
Non-obese_BMI	2262 (63.7%)	2288 (61.2%)	2096 (58.4%)	1746 (55.1%)	8392 (59.8%)
Obese_BMI	1288 (36.3%)	1451 (38.8%)	1493 (41.6%)	1421 (44.9%)	5653 (40.2%)
Obesity_BMI.Fat
Non-obese_BMI.Fat	2271 (64.0%)	2305 (61.6%)	2118 (59.0%)	1760 (55.6%)	8454 (60.2%)
Obese_BMI.Fat	1279 (36.0%)	1434 (38.4%)	1471 (41.0%)	1407 (44.4%)	5591 (39.8%)

1.4b.1 Prevalence of obesity by BMI only

# Overall
ob.all %>%
  group_by(SDDSRVYR) %>%
  summarize(
    n = n(),
    ob_bmi = sum(ob.bmi == "Obese_BMI", na.rm = TRUE),
    prev.ob_bmi = ob_bmi / n,
    se.ob_bmi = sqrt(prev.ob_bmi * (1 - prev.ob_bmi) / n),
    ob_bmi.fat = sum(ob.bmi_fat == "Obese_BMI.Fat", na.rm = TRUE),
    prev.ob_bmi.fat = ob_bmi.fat / n,
    se.ob_bmi,fat = sqrt(prev.ob_bmi.fat * (1 - prev.ob_bmi.fat) / n),
    .groups = "drop"
  )

## # A tibble: 4 × 8
##   SDDSRVYR     n ob_bmi prev.ob_bmi se.ob_bmi ob_bmi.fat prev.ob_bmi.fat     fat
##   <labell> <int>  <int>       <dbl>     <dbl>      <int>           <dbl>   <dbl>
## 1  7        3550   1288       0.363   0.00807       1279           0.360 0.00806
## 2  8        3739   1451       0.388   0.00797       1434           0.384 0.00795
## 3  9        3589   1493       0.416   0.00823       1471           0.410 0.00821
## 4 10        3167   1421       0.449   0.00884       1407           0.444 0.00883

# By Sex
ob.all %>%
  group_by(SDDSRVYR, sex) %>%
  summarize(
    n = n(),
    ob_bmi = sum(ob.bmi == "Obese_BMI", na.rm = TRUE),
    prev.ob_bmi = ob_bmi / n,
    se.ob_bmi = sqrt(prev.ob_bmi * (1 - prev.ob_bmi) / n),
    ob_bmi.fat = sum(ob.bmi_fat == "Obese_BMI.Fat", na.rm = TRUE),
    prev.ob_bmi.fat = ob_bmi.fat / n,
    se.ob_bmi,fat = sqrt(prev.ob_bmi.fat * (1 - prev.ob_bmi.fat) / n),
    .groups = "drop"
  )

## # A tibble: 8 × 9
##   SDDSRVYR   sex       n ob_bmi prev.ob_bmi se.ob_bmi ob_bmi.fat prev.ob_bmi.fat
##   <labelled> <fct> <int>  <int>       <dbl>     <dbl>      <int>           <dbl>
## 1  7         Male   1745    570       0.327    0.0112        561           0.321
## 2  7         Fema…  1805    718       0.398    0.0115        718           0.398
## 3  8         Male   1790    618       0.345    0.0112        601           0.336
## 4  8         Fema…  1949    833       0.427    0.0112        833           0.427
## 5  9         Male   1688    644       0.382    0.0118        622           0.368
## 6  9         Fema…  1901    849       0.447    0.0114        849           0.447
## 7 10         Male   1479    651       0.440    0.0129        638           0.431
## 8 10         Fema…  1688    770       0.456    0.0121        769           0.456
## # ℹ 1 more variable: fat <dbl>

# By Age group
ob.all %>%
  group_by(SDDSRVYR, age.grp) %>%
  summarize(
    n = n(),
    ob_bmi = sum(ob.bmi == "Obese_BMI", na.rm = TRUE),
    prev.ob_bmi = ob_bmi / n,
    se.ob_bmi = sqrt(prev.ob_bmi * (1 - prev.ob_bmi) / n),
    ob_bmi.fat = sum(ob.bmi_fat == "Obese_BMI.Fat", na.rm = TRUE),
    prev.ob_bmi.fat = ob_bmi.fat / n,
    se.ob_bmi,fat = sqrt(prev.ob_bmi.fat * (1 - prev.ob_bmi.fat) / n),
    .groups = "drop"
  )

## # A tibble: 16 × 9
##    SDDSRVYR   age.grp        n ob_bmi prev.ob_bmi se.ob_bmi ob_bmi.fat
##    <labelled> <labelled> <int>  <int>       <dbl>     <dbl>      <int>
##  1  7         20-29        939    278       0.296    0.0149        275
##  2  7         30-39        900    303       0.337    0.0158        301
##  3  7         40-49        851    350       0.411    0.0169        346
##  4  7         50-59        860    357       0.415    0.0168        357
##  5  8         20-29        906    268       0.296    0.0152        262
##  6  8         30-39        951    396       0.416    0.0160        390
##  7  8         40-49        985    387       0.393    0.0156        384
##  8  8         50-59        897    400       0.446    0.0166        398
##  9  9         20-29        911    303       0.333    0.0156        295
## 10  9         30-39        909    383       0.421    0.0164        378
## 11  9         40-49        892    418       0.469    0.0167        411
## 12  9         50-59        877    389       0.444    0.0168        387
## 13 10         20-29        759    292       0.385    0.0177        289
## 14 10         30-39        790    352       0.446    0.0177        345
## 15 10         40-49        758    378       0.499    0.0182        376
## 16 10         50-59        860    399       0.464    0.0170        397
## # ℹ 2 more variables: prev.ob_bmi.fat <dbl>, fat <dbl>

# By Race
ob.all %>%
  group_by(SDDSRVYR, race) %>%
  summarize(
    n = n(),
    ob_bmi = sum(ob.bmi == "Obese_BMI", na.rm = TRUE),
    prev.ob_bmi = ob_bmi / n,
    se.ob_bmi = sqrt(prev.ob_bmi * (1 - prev.ob_bmi) / n),
    ob_bmi.fat = sum(ob.bmi_fat == "Obese_BMI.Fat", na.rm = TRUE),
    prev.ob_bmi.fat = ob_bmi.fat / n,
    se.ob_bmi.fat = sqrt(prev.ob_bmi.fat * (1 - prev.ob_bmi.fat) / n),
    .groups = "drop"
  )

## # A tibble: 20 × 9
##    SDDSRVYR  race      n ob_bmi prev.ob_bmi se.ob_bmi ob_bmi.fat prev.ob_bmi.fat
##    <labelle> <fct> <int>  <int>       <dbl>     <dbl>      <int>           <dbl>
##  1  7        Hisp…   730    302       0.414    0.0182        300           0.411
##  2  7        Non-…   567    112       0.198    0.0167        111           0.196
##  3  7        Non-…  1221    419       0.343    0.0136        418           0.342
##  4  7        Non-…   914    418       0.457    0.0165        413           0.452
##  5  7        <NA>    118     37       0.314    0.0427         37           0.314
##  6  8        Hisp…   864    377       0.436    0.0169        372           0.431
##  7  8        Non-…   469    108       0.230    0.0194        108           0.230
##  8  8        Non-…  1511    554       0.367    0.0124        549           0.363
##  9  8        Non-…   752    354       0.471    0.0182        349           0.464
## 10  8        <NA>    143     58       0.406    0.0411         56           0.392
## 11  9        Hisp…  1087    516       0.475    0.0151        512           0.471
## 12  9        Non-…   497    119       0.239    0.0191        118           0.237
## 13  9        Non-…  1069    408       0.382    0.0149        402           0.376
## 14  9        Non-…   791    379       0.479    0.0178        369           0.466
## 15  9        <NA>    145     71       0.490    0.0415         70           0.483
## 16 10        Hisp…   770    350       0.455    0.0179        345           0.448
## 17 10        Non-…   520    172       0.331    0.0206        171           0.329
## 18 10        Non-…   962    424       0.441    0.0160        421           0.438
## 19 10        Non-…   731    381       0.521    0.0185        376           0.514
## 20 10        <NA>    184     94       0.511    0.0369         94           0.511
## # ℹ 1 more variable: se.ob_bmi.fat <dbl>

1.4b.2 Prevalence of confirmed excess adiposity among obesity by BMI only

By sex, age group, ethnicity and wave

ob.all.bmi = subset(ob.all, ob.bmi == "Obese_BMI")

# Overall
ob.all.bmi %>%
  group_by(SDDSRVYR) %>%
  summarize(
    n = n(),
    excess_fat = sum(excess.fat == 1, na.rm = TRUE),
    prev.excess_fat = excess_fat / n,
    se.excess_fat = sqrt(prev.excess_fat * (1 - prev.excess_fat) / n),
    .groups = "drop"
  )

## # A tibble: 4 × 5
##   SDDSRVYR       n excess_fat prev.excess_fat se.excess_fat
##   <labelled> <int>      <int>           <dbl>         <dbl>
## 1  7          1288       1268           0.984       0.00345
## 2  8          1451       1410           0.972       0.00435
## 3  9          1493       1453           0.973       0.00418
## 4 10          1421       1389           0.977       0.00394

# By Sex
ob.all.bmi %>%
  group_by(SDDSRVYR, sex) %>%
  summarize(
    n = n(),
    excess_fat = sum(excess.fat == 1, na.rm = TRUE),
    prev.excess_fat = excess_fat / n,
    se.excess_fat = sqrt(prev.excess_fat * (1 - prev.excess_fat) / n),
    .groups = "drop"
  )

## # A tibble: 8 × 6
##   SDDSRVYR   sex        n excess_fat prev.excess_fat se.excess_fat
##   <labelled> <fct>  <int>      <int>           <dbl>         <dbl>
## 1  7         Male     570        557           0.977       0.00625
## 2  7         Female   718        711           0.990       0.00367
## 3  8         Male     618        594           0.961       0.00777
## 4  8         Female   833        816           0.980       0.00490
## 5  9         Male     644        619           0.961       0.00761
## 6  9         Female   849        834           0.982       0.00452
## 7 10         Male     651        629           0.966       0.00708
## 8 10         Female   770        760           0.987       0.00408

# By Age group
ob.all.bmi %>%
  group_by(SDDSRVYR, age.grp) %>%
  summarize(
    n = n(),
    excess_fat = sum(excess.fat == 1, na.rm = TRUE),
    prev.excess_fat = excess_fat / n,
    se.excess_fat = sqrt(prev.excess_fat * (1 - prev.excess_fat) / n),
    .groups = "drop"
  )

## # A tibble: 16 × 6
##    SDDSRVYR   age.grp        n excess_fat prev.excess_fat se.excess_fat
##    <labelled> <labelled> <int>      <int>           <dbl>         <dbl>
##  1  7         20-29        278        272           0.978       0.00872
##  2  7         30-39        303        299           0.987       0.00656
##  3  7         40-49        350        341           0.974       0.00846
##  4  7         50-59        357        356           0.997       0.00280
##  5  8         20-29        268        258           0.963       0.0116 
##  6  8         30-39        396        379           0.957       0.0102 
##  7  8         40-49        387        380           0.982       0.00677
##  8  8         50-59        400        393           0.982       0.00656
##  9  9         20-29        303        294           0.970       0.00975
## 10  9         30-39        383        372           0.971       0.00853
## 11  9         40-49        418        403           0.964       0.00910
## 12  9         50-59        389        384           0.987       0.00571
## 13 10         20-29        292        285           0.976       0.00895
## 14 10         30-39        352        338           0.960       0.0104 
## 15 10         40-49        378        374           0.989       0.00526
## 16 10         50-59        399        392           0.982       0.00657

# By Race
ob.all.bmi %>%
  group_by(SDDSRVYR, race) %>%
  summarize(
    n = n(),
    excess_fat = sum(excess.fat == 1, na.rm = TRUE),
    prev.excess_fat = excess_fat / n,
    se.excess_fat = sqrt(prev.excess_fat * (1 - prev.excess_fat) / n),
    .groups = "drop"
  )

## # A tibble: 20 × 6
##    SDDSRVYR   race                   n excess_fat prev.excess_fat se.excess_fat
##    <labelled> <fct>              <int>      <int>           <dbl>         <dbl>
##  1  7         Hispanic             302        298           0.987       0.00658
##  2  7         Non-Hispanic Asian   112        111           0.991       0.00889
##  3  7         Non-Hispanic White   419        415           0.990       0.00475
##  4  7         Non-Hispanic Black   418        408           0.976       0.00747
##  5  7         <NA>                  37         36           0.973       0.0267 
##  6  8         Hispanic             377        366           0.971       0.00867
##  7  8         Non-Hispanic Asian   108        108           1           0      
##  8  8         Non-Hispanic White   554        545           0.984       0.00537
##  9  8         Non-Hispanic Black   354        338           0.955       0.0110 
## 10  8         <NA>                  58         53           0.914       0.0369 
## 11  9         Hispanic             516        508           0.984       0.00544
## 12  9         Non-Hispanic Asian   119        117           0.983       0.0118 
## 13  9         Non-Hispanic White   408        400           0.980       0.00686
## 14  9         Non-Hispanic Black   379        359           0.947       0.0115 
## 15  9         <NA>                  71         69           0.972       0.0196 
## 16 10         Hispanic             350        343           0.98        0.00748
## 17 10         Non-Hispanic Asian   172        170           0.988       0.00817
## 18 10         Non-Hispanic White   424        414           0.976       0.00737
## 19 10         Non-Hispanic Black   381        369           0.969       0.00895
## 20 10         <NA>                  94         93           0.989       0.0106

1.5 Changes over time

1.5.1 Changes of BMI

p1 = ggplot(data = ob, aes(y = BMXBMI, x = as.factor(SDDSRVYR), fill = as.factor(SDDSRVYR))) + geom_boxplot() + theme_bw() + theme(legend.position="none") + labs(title = "Changes of BMI", x = "NHANES Cycles", y = "BMI (kg/m2)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))

p2 = ggplot(data = ob, aes(y = DXDTOPF, x = as.factor(SDDSRVYR), fill = as.factor(SDDSRVYR))) + geom_boxplot() + theme_bw() + theme(legend.position="none") + labs(title = "Changes of percentage of fat", x = "NHANES Cycles", y = "Percentage of fat (%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))

p3 = ggplot(data = ob, aes(y = BMXWAIST, x = as.factor(SDDSRVYR), fill = as.factor(SDDSRVYR))) + geom_boxplot() + theme_bw() + theme(legend.position="none") + labs(title = "Changes of waist circumference", x = "NHANES Cycles", y = "Waist circumference (cm)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))

grid.arrange(p1, p2, p3, nrow = 1)

## Warning: Removed 139 rows containing non-finite values (`stat_boxplot()`).

## Warning: Removed 3483 rows containing non-finite values (`stat_boxplot()`).

## Warning: Removed 672 rows containing non-finite values (`stat_boxplot()`).

1.5.2 Relationship between BMI and pcfat by NHANES cycles

p1 = ggplot(data = ob, aes(x = BMXBMI, y = DXDTOPF, col = as.factor(SDDSRVYR))) + geom_point() + geom_smooth() + theme_bw() + theme(legend.position="none") + labs(title = "BMI-Fat relationship", x = "BMI (kg/m2)", y = "Percentage of fat (%)")

p2 = ggplot(data = ob, aes(x = BMXBMI, y = BMXWAIST, col = as.factor(SDDSRVYR))) + geom_point() + geom_smooth() + theme_bw() + theme(legend.position="none") + labs(title = "BMI-Waist relationship", x = "BMI (kg/m2)", y = "Waist circumference (cm)")

p3 = ggplot(data = ob, aes(x = BMXWAIST, y = DXDTOPF, col = as.factor(SDDSRVYR))) + geom_point() + geom_smooth() + theme_bw() + theme(legend.position="none") + labs(title = "Waist circ-Fat relationship", x = "Waist circumference (cm)", y = "Percentage of fat (%)")

grid.arrange(p1, p2, p3, nrow = 1)

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

## Warning: Removed 3518 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 3518 rows containing missing values (`geom_point()`).

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

## Warning: Removed 691 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 691 rows containing missing values (`geom_point()`).

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

## Warning: Removed 3593 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 3593 rows containing missing values (`geom_point()`).

2. VOS data

1.1 Data management

temp = read.csv("C:/Thach/Research projects/BMI_NHANES/Alldata.csv", header = TRUE)
vos = subset(temp, select = c(id, sex, height, weight, bmi, wc, yob, ethnic, wbbmd, wbtotfat, pcfat))
head(vos)

##     id sex height weight      bmi   wc  yob ethnic    wbbmd wbtotfat    pcfat
## 1 2303   M  167.0   64.0 22.94812 68.0 1989      1 1.180955 19444.57 30.03026
## 2 2433   M  161.5   50.0 19.17013 23.6 1957      1 1.119966 12574.17 24.78205
## 3 3172   F  145.1   43.8 20.80364 43.5 1970      1 1.009092 16444.02 37.14034
## 4 3830   M  166.8   72.5 26.05829 50.6 1969      1 1.142285 24241.46 33.71228
## 5 2094   F  156.0   40.0 16.43655 53.0 1989      1 0.851204 14871.17 36.75267
## 6 2087   M  170.0   66.0 22.83737 54.4 1970      1 1.114502 23273.87 35.18074

vos$age = 2016 - vos$yob
vos = subset(vos, age>= 20 & age <= 59)

Subgroup definition:

vos = vos %>%
  mutate(age.grp = case_when(age >= 20 & age <= 29 ~ "20-29",
                             age >= 30 & age <= 39 ~ "30-39",
                             age >= 40 & age <= 49 ~ "40-49",
                             age >= 50 & age <= 59 ~ "50-59"))
var_lab(vos$age.grp) = "Age group"

vos = vos %>%
  mutate(race = case_when(ethnic == 1 ~ "Viet",
                          ethnic == 2 ~ "Chinese",
                          ethnic == 3 ~ "Other"))

Obesity definition:

Obesity (Rubio et al. Lancet Diabetes Endocrinol 2025): (i) BMI (27.5 kg/m2 [Asian]); (ii) BMI and (percentage fat [25% in M, 35% in W] or WHR [90 cm in M, 80 in W for Asian])

# Obesity defined by only BMI
vos$ob.bmi = ifelse(vos$bmi >= 27.5, 1, 0)
vos$ob.bmi <- factor(vos$ob.bmi, levels = c(0, 1), labels = c("Non-obese_BMI", "Obese_BMI"))
var_lab(vos$ob.bmi) = "Obesity_only BMI"
  table(vos$ob.bmi)

## 
## Non-obese_BMI     Obese_BMI 
##          3035           256

# Temporal:  
vos$perc.fat = ifelse((vos$pcfat > 25 & vos$sex == "M") | (vos$pcfat > 35 & vos$sex == "F"), 1, 0)  
  table(vos$perc.fat)

## 
##    0    1 
##  390 2905

  #### pcfat in VOS 3% higher than NHANES_2018 (???): 41.2% (W), 30.1% (M) versus 38.9% (W), 27.5% (M)
  #### though BMI much lower: 22.5 (W), 23.4 (M) versus 30.3 (W), 29.7 (M)
  
vos$waist = ifelse((vos$wc > 90 & vos$sex == "M") | (vos$wc > 80 & vos$sex == "F"), 1, 0)
  table(vos$waist)

## 
##    0    1 
## 2247 1041

vos$excess.fat = ifelse(vos$perc.fat == 1 | vos$waist == 1, 1, 0)
  
# Obesity defined by BMI + (percentage fat or waist)
vos$ob.bmi_fat = ifelse((vos$ob.bmi == "Obese_BMI" & (vos$perc.fat == 1 | vos$waist == 1)) | 
                         vos$bmi >= 40, 1, 0)
vos$ob.bmi_fat <- factor(vos$ob.bmi_fat, levels = c(0, 1), labels = c("Non-obese_BMI.Fat", "Obese_BMI.Fat"))
var_lab(vos$ob.bmi_fat) = "Obesity_BMI.Fat"
  table(vos$ob.bmi_fat)   

## 
## Non-obese_BMI.Fat     Obese_BMI.Fat 
##              3035               256

2.2 Characteristics

table1(~ age + age.grp + race + weight + height + bmi + wc + wbbmd + wbtotfat + pcfat + ob.bmi + ob.bmi_fat + as.factor(perc.fat) + as.factor(waist) + as.factor(excess.fat) | sex, render.continuous = c(. = "Mean (SD)", . = "Median [Q1, Q3]", . = "Median [min, max]"), data = vos) 

	F (N=2106)	M (N=1189)	Overall (N=3295)
age
Mean (SD)	43.8 (11.4)	41.3 (12.1)	42.9 (11.7)
Median [Q1, Q3]	46.0 [36.0, 53.0]	43.0 [30.0, 52.0]	46.0 [33.0, 53.0]
Median [min, max]	46.0 [20.0, 59.0]	43.0 [20.0, 59.0]	46.0 [20.0, 59.0]
Age group
20-29	347 (16.5%)	281 (23.6%)	628 (19.1%)
30-39	317 (15.1%)	228 (19.2%)	545 (16.5%)
40-49	609 (28.9%)	284 (23.9%)	893 (27.1%)
50-59	833 (39.6%)	396 (33.3%)	1229 (37.3%)
race
Chinese	68 (3.2%)	47 (4.0%)	115 (3.5%)
Viet	2035 (96.6%)	1142 (96.0%)	3177 (96.4%)
Missing	3 (0.1%)	0 (0%)	3 (0.1%)
weight
Mean (SD)	53.1 (7.95)	63.5 (10.1)	56.9 (10.1)
Median [Q1, Q3]	52.0 [48.0, 57.1]	63.0 [57.0, 70.0]	55.5 [50.0, 62.5]
Median [min, max]	52.0 [30.0, 90.0]	63.0 [37.0, 121]	55.5 [30.0, 121]
Missing	3 (0.1%)	1 (0.1%)	4 (0.1%)
height
Mean (SD)	154 (5.24)	165 (6.03)	158 (7.66)
Median [Q1, Q3]	154 [150, 157]	165 [161, 168]	157 [152, 163]
Median [min, max]	154 [131, 173]	165 [140, 197]	157 [131, 197]
Missing	3 (0.1%)	1 (0.1%)	4 (0.1%)
bmi
Mean (SD)	22.5 (3.18)	23.4 (3.29)	22.8 (3.25)
Median [Q1, Q3]	22.2 [20.3, 24.3]	23.2 [21.2, 25.3]	22.5 [20.6, 24.7]
Median [min, max]	22.2 [14.0, 39.3]	23.2 [14.2, 38.4]	22.5 [14.0, 39.3]
Missing	3 (0.1%)	1 (0.1%)	4 (0.1%)
wc
Mean (SD)	77.8 (8.77)	83.2 (9.22)	79.8 (9.30)
Median [Q1, Q3]	77.0 [71.8, 83.4]	83.5 [77.4, 89.2]	79.5 [73.2, 86.0]
Median [min, max]	77.0 [43.5, 108]	83.5 [23.6, 125]	79.5 [23.6, 125]
Missing	4 (0.2%)	3 (0.3%)	7 (0.2%)
wbbmd
Mean (SD)	1.01 (0.0949)	1.08 (0.0898)	1.03 (0.0982)
Median [Q1, Q3]	1.01 [0.947, 1.07]	1.07 [1.01, 1.13]	1.04 [0.971, 1.10]
Median [min, max]	1.01 [0.654, 1.35]	1.07 [0.767, 1.45]	1.04 [0.654, 1.45]
wbtotfat
Mean (SD)	22400 (5250)	19600 (5900)	21400 (5650)
Median [Q1, Q3]	22000 [18900, 25300]	19500 [15700, 23300]	21100 [17700, 24600]
Median [min, max]	22000 [8140, 48300]	19500 [6680, 54800]	21100 [6680, 54800]
pcfat
Mean (SD)	41.2 (4.58)	30.1 (5.20)	37.2 (7.21)
Median [Q1, Q3]	41.3 [38.3, 44.3]	30.6 [26.8, 33.6]	38.2 [32.2, 42.6]
Median [min, max]	41.3 [23.1, 54.9]	30.6 [14.1, 47.2]	38.2 [14.1, 54.9]
Obesity_only BMI
Non-obese_BMI	1971 (93.6%)	1064 (89.5%)	3035 (92.1%)
Obese_BMI	132 (6.3%)	124 (10.4%)	256 (7.8%)
Missing	3 (0.1%)	1 (0.1%)	4 (0.1%)
Obesity_BMI.Fat
Non-obese_BMI.Fat	1971 (93.6%)	1064 (89.5%)	3035 (92.1%)
Obese_BMI.Fat	132 (6.3%)	124 (10.4%)	256 (7.8%)
Missing	3 (0.1%)	1 (0.1%)	4 (0.1%)
as.factor(perc.fat)
0	185 (8.8%)	205 (17.2%)	390 (11.8%)
1	1921 (91.2%)	984 (82.8%)	2905 (88.2%)
as.factor(waist)
0	1316 (62.5%)	931 (78.3%)	2247 (68.2%)
1	786 (37.3%)	255 (21.4%)	1041 (31.6%)
Missing	4 (0.2%)	3 (0.3%)	7 (0.2%)
as.factor(excess.fat)
0	177 (8.4%)	203 (17.1%)	380 (11.5%)
1	1929 (91.6%)	985 (82.8%)	2914 (88.4%)
Missing	0 (0%)	1 (0.1%)	1 (0.0%)

Commparison with NHANES 2017-2018

Missing data excluded

# Whole population:
ob.all.10 = subset(ob.all, SDDSRVYR == 10)
  table(ob.all.10$race)

## 
##                      Hispanic            Non-Hispanic Asian 
##                           770                           520 
##            Non-Hispanic White            Non-Hispanic Black 
##                           962                           731 
## OtherRace, incl. Multi-Racial 
##                             0

table1(~ RIDAGEYR + age.grp + race + BMXWT + BMXHT + BMXBMI + BMXWAIST + DXDTOBMD + DXDTOFAT + DXDTOPF + ob.bmi + ob.bmi_fat + as.factor(perc.fat) + as.factor(waist) + as.factor(excess.fat) | sex, render.continuous = c(. = "Mean (SD)", . = "Median [Q1, Q3]", . = "Median [min, max]"), data = ob.all.10)

	Male (N=1479)	Female (N=1688)	Overall (N=3167)
RIDAGEYR
Mean (SD)	39.8 (11.7)	40.2 (11.5)	40.0 (11.6)
Median [Q1, Q3]	40.0 [30.0, 50.0]	40.0 [30.0, 51.0]	40.0 [30.0, 50.0]
Median [min, max]	40.0 [20.0, 59.0]	40.0 [20.0, 59.0]	40.0 [20.0, 59.0]
Age group
20-29	368 (24.9%)	391 (23.2%)	759 (24.0%)
30-39	355 (24.0%)	435 (25.8%)	790 (24.9%)
40-49	360 (24.3%)	398 (23.6%)	758 (23.9%)
50-59	396 (26.8%)	464 (27.5%)	860 (27.2%)
Race and ethnicity
Hispanic	351 (23.7%)	419 (24.8%)	770 (24.3%)
Non-Hispanic Asian	241 (16.3%)	279 (16.5%)	520 (16.4%)
Non-Hispanic White	457 (30.9%)	505 (29.9%)	962 (30.4%)
Non-Hispanic Black	329 (22.2%)	402 (23.8%)	731 (23.1%)
OtherRace, incl. Multi-Racial	0 (0%)	0 (0%)	0 (0%)
Missing	101 (6.8%)	83 (4.9%)	184 (5.8%)
BMXWT
Mean (SD)	90.6 (23.4)	78.6 (23.8)	84.2 (24.4)
Median [Q1, Q3]	86.2 [74.1, 102]	73.7 [61.1, 90.8]	79.7 [67.0, 97.2]
Median [min, max]	86.2 [43.7, 243]	73.7 [37.5, 204]	79.7 [37.5, 243]
BMXHT
Mean (SD)	175 (7.71)	161 (7.06)	167 (10.1)
Median [Q1, Q3]	175 [169, 180]	161 [156, 165]	167 [160, 175]
Median [min, max]	175 [148, 196]	161 [138, 189]	167 [138, 196]
BMXBMI
Mean (SD)	29.6 (6.97)	30.3 (8.60)	30.0 (7.88)
Median [Q1, Q3]	28.6 [25.0, 33.1]	28.6 [24.0, 34.7]	28.6 [24.5, 33.9]
Median [min, max]	28.6 [15.5, 86.2]	28.6 [14.8, 74.8]	28.6 [14.8, 86.2]
BMXWAIST
Mean (SD)	101 (17.2)	97.7 (18.6)	99.4 (18.0)
Median [Q1, Q3]	99.3 [89.5, 111]	95.4 [83.7, 109]	97.4 [86.2, 110]
Median [min, max]	99.3 [62.3, 164]	95.4 [57.9, 170]	97.4 [57.9, 170]
Missing	32 (2.2%)	59 (3.5%)	91 (2.9%)
DXDTOBMD
Mean (SD)	1.16 (0.105)	1.08 (0.0959)	1.12 (0.107)
Median [Q1, Q3]	1.15 [1.09, 1.22]	1.08 [1.02, 1.14]	1.11 [1.05, 1.18]
Median [min, max]	1.15 [0.876, 1.73]	1.08 [0.697, 1.40]	1.11 [0.697, 1.73]
Missing	409 (27.7%)	511 (30.3%)	920 (29.0%)
DXDTOFAT
Mean (SD)	24800 (10200)	30300 (12200)	27700 (11600)
Median [Q1, Q3]	23700 [17700, 30100]	28600 [20900, 37400]	25700 [19300, 34200]
Median [min, max]	23700 [6670, 72500]	28600 [8830, 86900]	25700 [6670, 86900]
Missing	413 (27.9%)	515 (30.5%)	928 (29.3%)
DXDTOPF
Mean (SD)	27.5 (6.17)	38.9 (6.57)	33.5 (8.54)
Median [Q1, Q3]	27.9 [23.5, 31.7]	39.9 [35.0, 43.4]	33.4 [27.3, 40.5]
Median [min, max]	27.9 [12.1, 48.6]	39.9 [17.7, 56.1]	33.4 [12.1, 56.1]
Missing	413 (27.9%)	515 (30.5%)	928 (29.3%)
Obesity_only BMI
Non-obese_BMI	828 (56.0%)	918 (54.4%)	1746 (55.1%)
Obese_BMI	651 (44.0%)	770 (45.6%)	1421 (44.9%)
Obesity_BMI.Fat
Non-obese_BMI.Fat	841 (56.9%)	919 (54.4%)	1760 (55.6%)
Obese_BMI.Fat	638 (43.1%)	769 (45.6%)	1407 (44.4%)
as.factor(perc.fat)
0	354 (23.9%)	298 (17.7%)	652 (20.6%)
1	712 (48.1%)	875 (51.8%)	1587 (50.1%)
Missing	413 (27.9%)	515 (30.5%)	928 (29.3%)
as.factor(waist)
0	716 (48.4%)	479 (28.4%)	1195 (37.7%)
1	731 (49.4%)	1150 (68.1%)	1881 (59.4%)
Missing	32 (2.2%)	59 (3.5%)	91 (2.9%)
as.factor(excess.fat)
0	327 (22.1%)	249 (14.8%)	576 (18.2%)
1	973 (65.8%)	1290 (76.4%)	2263 (71.5%)
Missing	179 (12.1%)	149 (8.8%)	328 (10.4%)

# Asians:
ob.all.10.asian = subset(ob.all, SDDSRVYR == 10 & race == "Non-Hispanic Asian")
  table(ob.all.10.asian$race)

## 
##                      Hispanic            Non-Hispanic Asian 
##                             0                           520 
##            Non-Hispanic White            Non-Hispanic Black 
##                             0                             0 
## OtherRace, incl. Multi-Racial 
##                             0

table1(~ RIDAGEYR + age.grp + race + BMXWT + BMXHT + BMXBMI + BMXWAIST + DXDTOBMD + DXDTOFAT + DXDTOPF + ob.bmi + ob.bmi_fat + as.factor(perc.fat) + as.factor(waist) + as.factor(excess.fat) | sex, render.continuous = c(. = "Mean (SD)", . = "Median [Q1, Q3]", . = "Median [min, max]"), data = ob.all.10.asian)

	Male (N=241)	Female (N=279)	Overall (N=520)
RIDAGEYR
Mean (SD)	40.7 (11.1)	41.2 (11.3)	41.0 (11.2)
Median [Q1, Q3]	41.0 [31.0, 50.0]	42.0 [31.0, 52.0]	41.0 [31.0, 51.0]
Median [min, max]	41.0 [20.0, 59.0]	42.0 [20.0, 59.0]	41.0 [20.0, 59.0]
Age group
20-29	52 (21.6%)	53 (19.0%)	105 (20.2%)
30-39	56 (23.2%)	72 (25.8%)	128 (24.6%)
40-49	64 (26.6%)	67 (24.0%)	131 (25.2%)
50-59	69 (28.6%)	87 (31.2%)	156 (30.0%)
Race and ethnicity
Hispanic	0 (0%)	0 (0%)	0 (0%)
Non-Hispanic Asian	241 (100%)	279 (100%)	520 (100%)
Non-Hispanic White	0 (0%)	0 (0%)	0 (0%)
Non-Hispanic Black	0 (0%)	0 (0%)	0 (0%)
OtherRace, incl. Multi-Racial	0 (0%)	0 (0%)	0 (0%)
BMXWT
Mean (SD)	79.4 (16.4)	62.0 (11.9)	70.1 (16.6)
Median [Q1, Q3]	77.4 [69.6, 85.9]	59.7 [53.3, 69.1]	68.8 [58.1, 78.9]
Median [min, max]	77.4 [43.7, 160]	59.7 [39.6, 109]	68.8 [39.6, 160]
BMXHT
Mean (SD)	171 (6.77)	157 (5.71)	164 (9.31)
Median [Q1, Q3]	170 [167, 177]	157 [154, 161]	163 [156, 170]
Median [min, max]	170 [152, 186]	157 [140, 177]	163 [140, 186]
BMXBMI
Mean (SD)	27.1 (5.11)	25.1 (4.80)	26.1 (5.04)
Median [Q1, Q3]	26.3 [24.1, 29.0]	24.3 [21.5, 28.0]	25.5 [22.6, 28.5]
Median [min, max]	26.3 [17.0, 60.6]	24.3 [16.5, 43.2]	25.5 [16.5, 60.6]
BMXWAIST
Mean (SD)	95.0 (12.6)	85.7 (11.5)	90.1 (12.9)
Median [Q1, Q3]	94.2 [87.0, 100]	84.8 [76.5, 93.5]	89.3 [81.1, 97.2]
Median [min, max]	94.2 [65.3, 152]	84.8 [65.0, 140]	89.3 [65.0, 152]
Missing	5 (2.1%)	13 (4.7%)	18 (3.5%)
DXDTOBMD
Mean (SD)	1.11 (0.0849)	1.05 (0.0883)	1.08 (0.0908)
Median [Q1, Q3]	1.10 [1.04, 1.16]	1.05 [1.00, 1.10]	1.08 [1.02, 1.13]
Median [min, max]	1.10 [0.918, 1.34]	1.05 [0.804, 1.35]	1.08 [0.804, 1.35]
Missing	38 (15.8%)	68 (24.4%)	106 (20.4%)
DXDTOFAT
Mean (SD)	22700 (7020)	23100 (7150)	22900 (7080)
Median [Q1, Q3]	22100 [18000, 26300]	22200 [18100, 28200]	22100 [18000, 27200]
Median [min, max]	22100 [9640, 52800]	22200 [9220, 48700]	22100 [9220, 52800]
Missing	37 (15.4%)	67 (24.0%)	104 (20.0%)
DXDTOPF
Mean (SD)	28.2 (4.58)	37.0 (5.84)	32.7 (6.84)
Median [Q1, Q3]	28.1 [25.0, 31.4]	37.3 [33.6, 41.3]	31.9 [27.3, 37.6]
Median [min, max]	28.1 [16.4, 41.3]	37.3 [19.0, 49.5]	31.9 [16.4, 49.5]
Missing	37 (15.4%)	67 (24.0%)	104 (20.0%)
Obesity_only BMI
Non-obese_BMI	149 (61.8%)	199 (71.3%)	348 (66.9%)
Obese_BMI	92 (38.2%)	80 (28.7%)	172 (33.1%)
Obesity_BMI.Fat
Non-obese_BMI.Fat	150 (62.2%)	199 (71.3%)	349 (67.1%)
Obese_BMI.Fat	91 (37.8%)	80 (28.7%)	171 (32.9%)
as.factor(perc.fat)
0	52 (21.6%)	75 (26.9%)	127 (24.4%)
1	152 (63.1%)	137 (49.1%)	289 (55.6%)
Missing	37 (15.4%)	67 (24.0%)	104 (20.0%)
as.factor(waist)
0	86 (35.7%)	89 (31.9%)	175 (33.7%)
1	150 (62.2%)	177 (63.4%)	327 (62.9%)
Missing	5 (2.1%)	13 (4.7%)	18 (3.5%)
as.factor(excess.fat)
0	40 (16.6%)	53 (19.0%)	93 (17.9%)
1	184 (76.3%)	198 (71.0%)	382 (73.5%)
Missing	17 (7.1%)	28 (10.0%)	45 (8.7%)

Original data (missing data INCLUDED)

# Whole population:
ob.10 = subset(ob, SDDSRVYR == 10)
  table(ob.10$race)

## 
##                      Hispanic            Non-Hispanic Asian 
##                           799                           544 
##            Non-Hispanic White            Non-Hispanic Black 
##                           971                           744 
## OtherRace, incl. Multi-Racial 
##                             0

table1(~ RIDAGEYR + age.grp + race + BMXWT + BMXHT + BMXBMI + BMXWAIST + DXDTOBMD + DXDTOFAT + DXDTOPF + ob.bmi + ob.bmi_fat + as.factor(perc.fat) + as.factor(waist) + as.factor(excess.fat) | sex, render.continuous = c(. = "Mean (SD)", . = "Median [Q1, Q3]", . = "Median [min, max]"), data = ob.10)

	Male (N=1519)	Female (N=1727)	Overall (N=3246)
RIDAGEYR
Mean (SD)	39.8 (11.7)	40.1 (11.5)	40.0 (11.6)
Median [Q1, Q3]	40.0 [30.0, 50.0]	40.0 [30.0, 51.0]	40.0 [30.0, 50.0]
Median [min, max]	40.0 [20.0, 59.0]	40.0 [20.0, 59.0]	40.0 [20.0, 59.0]
Age group
20-29	377 (24.8%)	398 (23.0%)	775 (23.9%)
30-39	365 (24.0%)	448 (25.9%)	813 (25.0%)
40-49	368 (24.2%)	410 (23.7%)	778 (24.0%)
50-59	409 (26.9%)	471 (27.3%)	880 (27.1%)
Race and ethnicity
Hispanic	363 (23.9%)	436 (25.2%)	799 (24.6%)
Non-Hispanic Asian	251 (16.5%)	293 (17.0%)	544 (16.8%)
Non-Hispanic White	464 (30.5%)	507 (29.4%)	971 (29.9%)
Non-Hispanic Black	337 (22.2%)	407 (23.6%)	744 (22.9%)
OtherRace, incl. Multi-Racial	0 (0%)	0 (0%)	0 (0%)
Missing	104 (6.8%)	84 (4.9%)	188 (5.8%)
BMXWT
Mean (SD)	90.6 (23.3)	78.6 (23.7)	84.2 (24.3)
Median [Q1, Q3]	86.3 [74.3, 103]	74.0 [61.4, 90.8]	80.0 [67.3, 97.2]
Median [min, max]	86.3 [43.7, 243]	74.0 [37.5, 204]	80.0 [37.5, 243]
Missing	17 (1.1%)	14 (0.8%)	31 (1.0%)
BMXHT
Mean (SD)	175 (7.71)	161 (7.06)	167 (10.1)
Median [Q1, Q3]	175 [169, 180]	161 [156, 165]	167 [160, 175]
Median [min, max]	175 [148, 196]	161 [138, 189]	167 [138, 196]
Missing	17 (1.1%)	14 (0.8%)	31 (1.0%)
BMXBMI
Mean (SD)	29.7 (6.94)	30.3 (8.55)	30.0 (7.85)
Median [Q1, Q3]	28.7 [25.0, 33.1]	28.7 [24.1, 34.8]	28.7 [24.5, 33.9]
Median [min, max]	28.7 [15.5, 86.2]	28.7 [14.8, 74.8]	28.7 [14.8, 86.2]
Missing	20 (1.3%)	15 (0.9%)	35 (1.1%)
BMXWAIST
Mean (SD)	101 (17.1)	97.7 (18.6)	99.3 (18.0)
Median [Q1, Q3]	99.3 [89.5, 111]	95.4 [83.7, 109]	97.4 [86.2, 110]
Median [min, max]	99.3 [62.3, 164]	95.4 [57.9, 170]	97.4 [57.9, 170]
Missing	63 (4.1%)	96 (5.6%)	159 (4.9%)
DXDTOBMD
Mean (SD)	1.16 (0.104)	1.08 (0.0959)	1.12 (0.107)
Median [Q1, Q3]	1.15 [1.09, 1.22]	1.08 [1.02, 1.14]	1.11 [1.05, 1.18]
Median [min, max]	1.15 [0.876, 1.73]	1.08 [0.697, 1.40]	1.11 [0.697, 1.73]
Missing	443 (29.2%)	546 (31.6%)	989 (30.5%)
DXDTOFAT
Mean (SD)	24800 (10300)	30300 (12200)	27700 (11600)
Median [Q1, Q3]	23700 [17700, 30100]	28600 [20900, 37400]	25700 [19300, 34200]
Median [min, max]	23700 [6670, 72500]	28600 [8830, 86900]	25700 [6670, 86900]
Missing	447 (29.4%)	550 (31.8%)	997 (30.7%)
DXDTOPF
Mean (SD)	27.5 (6.19)	38.9 (6.58)	33.5 (8.56)
Median [Q1, Q3]	27.9 [23.5, 31.7]	39.9 [35.0, 43.4]	33.4 [27.3, 40.5]
Median [min, max]	27.9 [12.1, 48.6]	39.9 [17.7, 56.1]	33.4 [12.1, 56.1]
Missing	447 (29.4%)	550 (31.8%)	997 (30.7%)
Obesity_only BMI
Non-obese_BMI	828 (54.5%)	918 (53.2%)	1746 (53.8%)
Obese_BMI	671 (44.2%)	794 (46.0%)	1465 (45.1%)
Missing	20 (1.3%)	15 (0.9%)	35 (1.1%)
Obesity_BMI.Fat
Non-obese_BMI.Fat	841 (55.4%)	919 (53.2%)	1760 (54.2%)
Obese_BMI.Fat	638 (42.0%)	769 (44.5%)	1407 (43.3%)
Missing	40 (2.6%)	39 (2.3%)	79 (2.4%)
as.factor(perc.fat)
0	356 (23.4%)	298 (17.3%)	654 (20.1%)
1	716 (47.1%)	879 (50.9%)	1595 (49.1%)
Missing	447 (29.4%)	550 (31.8%)	997 (30.7%)
as.factor(waist)
0	722 (47.5%)	481 (27.9%)	1203 (37.1%)
1	734 (48.3%)	1150 (66.6%)	1884 (58.0%)
Missing	63 (4.1%)	96 (5.6%)	159 (4.9%)
as.factor(excess.fat)
0	327 (21.5%)	249 (14.4%)	576 (17.7%)
1	977 (64.3%)	1294 (74.9%)	2271 (70.0%)
Missing	215 (14.2%)	184 (10.7%)	399 (12.3%)

# Asians:
ob.10.asian = subset(ob, SDDSRVYR == 10 & race == "Non-Hispanic Asian")
  table(ob.10.asian$race)

## 
##                      Hispanic            Non-Hispanic Asian 
##                             0                           544 
##            Non-Hispanic White            Non-Hispanic Black 
##                             0                             0 
## OtherRace, incl. Multi-Racial 
##                             0

table1(~ RIDAGEYR + age.grp + race + BMXWT + BMXHT + BMXBMI + BMXWAIST + DXDTOBMD + DXDTOFAT + DXDTOPF + ob.bmi + ob.bmi_fat + as.factor(perc.fat) + as.factor(waist) + as.factor(excess.fat) | sex, render.continuous = c(. = "Mean (SD)", . = "Median [Q1, Q3]", . = "Median [min, max]"), data = ob.10.asian)

	Male (N=251)	Female (N=293)	Overall (N=544)
RIDAGEYR
Mean (SD)	40.7 (11.0)	41.0 (11.3)	40.8 (11.2)
Median [Q1, Q3]	41.0 [31.5, 50.0]	41.0 [31.0, 52.0]	41.0 [31.0, 51.0]
Median [min, max]	41.0 [20.0, 59.0]	41.0 [20.0, 59.0]	41.0 [20.0, 59.0]
Age group
20-29	53 (21.1%)	55 (18.8%)	108 (19.9%)
30-39	59 (23.5%)	81 (27.6%)	140 (25.7%)
40-49	68 (27.1%)	68 (23.2%)	136 (25.0%)
50-59	71 (28.3%)	89 (30.4%)	160 (29.4%)
Race and ethnicity
Hispanic	0 (0%)	0 (0%)	0 (0%)
Non-Hispanic Asian	251 (100%)	293 (100%)	544 (100%)
Non-Hispanic White	0 (0%)	0 (0%)	0 (0%)
Non-Hispanic Black	0 (0%)	0 (0%)	0 (0%)
OtherRace, incl. Multi-Racial	0 (0%)	0 (0%)	0 (0%)
BMXWT
Mean (SD)	79.7 (16.6)	62.6 (12.0)	70.5 (16.6)
Median [Q1, Q3]	77.8 [69.7, 86.2]	60.6 [53.8, 70.5]	69.2 [58.4, 79.2]
Median [min, max]	77.8 [43.7, 160]	60.6 [39.6, 109]	69.2 [39.6, 160]
Missing	3 (1.2%)	1 (0.3%)	4 (0.7%)
BMXHT
Mean (SD)	171 (6.74)	157 (5.71)	163 (9.32)
Median [Q1, Q3]	170 [167, 177]	157 [154, 161]	163 [156, 170]
Median [min, max]	170 [152, 186]	157 [140, 177]	163 [140, 186]
Missing	5 (2.0%)	1 (0.3%)	6 (1.1%)
BMXBMI
Mean (SD)	27.2 (5.15)	25.4 (4.88)	26.2 (5.08)
Median [Q1, Q3]	26.4 [24.1, 29.1]	24.5 [21.8, 28.4]	25.6 [22.8, 28.7]
Median [min, max]	26.4 [17.0, 60.6]	24.5 [16.5, 43.2]	25.6 [16.5, 60.6]
Missing	5 (2.0%)	1 (0.3%)	6 (1.1%)
BMXWAIST
Mean (SD)	95.0 (12.5)	85.7 (11.5)	90.1 (12.9)
Median [Q1, Q3]	94.4 [87.0, 100]	84.8 [76.5, 93.5]	89.4 [81.2, 97.2]
Median [min, max]	94.4 [65.3, 152]	84.8 [65.0, 140]	89.4 [65.0, 152]
Missing	13 (5.2%)	27 (9.2%)	40 (7.4%)
DXDTOBMD
Mean (SD)	1.11 (0.0847)	1.05 (0.0883)	1.08 (0.0907)
Median [Q1, Q3]	1.10 [1.04, 1.16]	1.05 [1.00, 1.10]	1.08 [1.02, 1.13]
Median [min, max]	1.10 [0.918, 1.34]	1.05 [0.804, 1.35]	1.08 [0.804, 1.35]
Missing	46 (18.3%)	82 (28.0%)	128 (23.5%)
DXDTOFAT
Mean (SD)	22700 (6980)	23100 (7150)	22900 (7060)
Median [Q1, Q3]	22200 [18000, 26200]	22200 [18100, 28200]	22200 [18000, 27100]
Median [min, max]	22200 [9640, 52800]	22200 [9220, 48700]	22200 [9220, 52800]
Missing	45 (17.9%)	81 (27.6%)	126 (23.2%)
DXDTOPF
Mean (SD)	28.3 (4.61)	37.0 (5.84)	32.7 (6.83)
Median [Q1, Q3]	28.1 [25.0, 31.5]	37.3 [33.6, 41.3]	31.9 [27.4, 37.6]
Median [min, max]	28.1 [16.4, 41.3]	37.3 [19.0, 49.5]	31.9 [16.4, 49.5]
Missing	45 (17.9%)	81 (27.6%)	126 (23.2%)
Obesity_only BMI
Non-obese_BMI	149 (59.4%)	199 (67.9%)	348 (64.0%)
Obese_BMI	97 (38.6%)	93 (31.7%)	190 (34.9%)
Missing	5 (2.0%)	1 (0.3%)	6 (1.1%)
Obesity_BMI.Fat
Non-obese_BMI.Fat	150 (59.8%)	199 (67.9%)	349 (64.2%)
Obese_BMI.Fat	91 (36.3%)	80 (27.3%)	171 (31.4%)
Missing	10 (4.0%)	14 (4.8%)	24 (4.4%)
as.factor(perc.fat)
0	52 (20.7%)	75 (25.6%)	127 (23.3%)
1	154 (61.4%)	137 (46.8%)	291 (53.5%)
Missing	45 (17.9%)	81 (27.6%)	126 (23.2%)
as.factor(waist)
0	86 (34.3%)	89 (30.4%)	175 (32.2%)
1	152 (60.6%)	177 (60.4%)	329 (60.5%)
Missing	13 (5.2%)	27 (9.2%)	40 (7.4%)
as.factor(excess.fat)
0	40 (15.9%)	53 (18.1%)	93 (17.1%)
1	186 (74.1%)	198 (67.6%)	384 (70.6%)
Missing	25 (10.0%)	42 (14.3%)	67 (12.3%)

2.3 Prevalence of obesity by BMI only

# Overall
vos %>%
  summarize(
    n = n(),
    ob_bmi = sum(ob.bmi == "Obese_BMI", na.rm = TRUE),
    prev.ob_bmi = ob_bmi / n,
    se.ob_bmi = sqrt(prev.ob_bmi * (1 - prev.ob_bmi) / n),
    ob_bmi.fat = sum(ob.bmi_fat == "Obese_BMI.Fat", na.rm = TRUE),
    prev.ob_bmi.fat = ob_bmi.fat / n,
    se.ob_bmi,fat = sqrt(prev.ob_bmi.fat * (1 - prev.ob_bmi.fat) / n),
    .groups = "drop"
  )

##      n ob_bmi prev.ob_bmi   se.ob_bmi ob_bmi.fat prev.ob_bmi.fat         fat
## 1 3295    256  0.07769347 0.004663395        256      0.07769347 0.004663395

# By Sex
vos %>%
  group_by(sex) %>%
  summarize(
    n = n(),
    ob_bmi = sum(ob.bmi == "Obese_BMI", na.rm = TRUE),
    prev.ob_bmi = ob_bmi / n,
    se.ob_bmi = sqrt(prev.ob_bmi * (1 - prev.ob_bmi) / n),
    ob_bmi.fat = sum(ob.bmi_fat == "Obese_BMI.Fat", na.rm = TRUE),
    prev.ob_bmi.fat = ob_bmi.fat / n,
    se.ob_bmi,fat = sqrt(prev.ob_bmi.fat * (1 - prev.ob_bmi.fat) / n),
    .groups = "drop"
  )

## # A tibble: 2 × 8
##   sex       n ob_bmi prev.ob_bmi se.ob_bmi ob_bmi.fat prev.ob_bmi.fat     fat
##   <chr> <int>  <int>       <dbl>     <dbl>      <int>           <dbl>   <dbl>
## 1 F      2106    132      0.0627   0.00528        132          0.0627 0.00528
## 2 M      1189    124      0.104    0.00886        124          0.104  0.00886

# By Age group
vos %>%
  group_by(age.grp) %>%
  summarize(
    n = n(),
    ob_bmi = sum(ob.bmi == "Obese_BMI", na.rm = TRUE),
    prev.ob_bmi = ob_bmi / n,
    se.ob_bmi = sqrt(prev.ob_bmi * (1 - prev.ob_bmi) / n),
    ob_bmi.fat = sum(ob.bmi_fat == "Obese_BMI.Fat", na.rm = TRUE),
    prev.ob_bmi.fat = ob_bmi.fat / n,
    se.ob_bmi,fat = sqrt(prev.ob_bmi.fat * (1 - prev.ob_bmi.fat) / n),
    .groups = "drop"
  )

## # A tibble: 4 × 8
##   age.grp      n ob_bmi prev.ob_bmi se.ob_bmi ob_bmi.fat prev.ob_bmi.fat     fat
##   <labell> <int>  <int>       <dbl>     <dbl>      <int>           <dbl>   <dbl>
## 1 20-29      628     30      0.0478   0.00851         30          0.0478 0.00851
## 2 30-39      545     38      0.0697   0.0109          38          0.0697 0.0109 
## 3 40-49      893     79      0.0885   0.00950         79          0.0885 0.00950
## 4 50-59     1229    109      0.0887   0.00811        109          0.0887 0.00811

# By Race
vos %>%
  group_by(race) %>%
  summarize(
    n = n(),
    ob_bmi = sum(ob.bmi == "Obese_BMI", na.rm = TRUE),
    prev.ob_bmi = ob_bmi / n,
    se.ob_bmi = sqrt(prev.ob_bmi * (1 - prev.ob_bmi) / n),
    ob_bmi.fat = sum(ob.bmi_fat == "Obese_BMI.Fat", na.rm = TRUE),
    prev.ob_bmi.fat = ob_bmi.fat / n,
    se.ob_bmi.fat = sqrt(prev.ob_bmi.fat * (1 - prev.ob_bmi.fat) / n),
    .groups = "drop"
  )

## # A tibble: 3 × 8
##   race        n ob_bmi prev.ob_bmi se.ob_bmi ob_bmi.fat prev.ob_bmi.fat
##   <chr>   <int>  <int>       <dbl>     <dbl>      <int>           <dbl>
## 1 Chinese   115      8      0.0696   0.0237           8          0.0696
## 2 Viet     3177    248      0.0781   0.00476        248          0.0781
## 3 <NA>        3      0      0        0                0          0     
## # ℹ 1 more variable: se.ob_bmi.fat <dbl>

2.4 Prevalence of confirmed excess adiposity among obesity by BMI only

vos.bmi = subset(vos, ob.bmi == "Obese_BMI")

# Overall
vos.bmi %>%
  summarize(
    n = n(),
    excess_fat = sum(excess.fat == 1, na.rm = TRUE),
    prev.excess_fat = excess_fat / n,
    se.excess_fat = sqrt(prev.excess_fat * (1 - prev.excess_fat) / n),
    .groups = "drop"
  )

##     n excess_fat prev.excess_fat se.excess_fat
## 1 256        256               1             0

# By Sex
vos.bmi %>%
  group_by(sex) %>%
  summarize(
    n = n(),
    excess_fat = sum(excess.fat == 1, na.rm = TRUE),
    prev.excess_fat = excess_fat / n,
    se.excess_fat = sqrt(prev.excess_fat * (1 - prev.excess_fat) / n),
    .groups = "drop"
  )

## # A tibble: 2 × 5
##   sex       n excess_fat prev.excess_fat se.excess_fat
##   <chr> <int>      <int>           <dbl>         <dbl>
## 1 F       132        132               1             0
## 2 M       124        124               1             0

# By Age group
vos.bmi %>%
  group_by(age.grp) %>%
  summarize(
    n = n(),
    excess_fat = sum(excess.fat == 1, na.rm = TRUE),
    prev.excess_fat = excess_fat / n,
    se.excess_fat = sqrt(prev.excess_fat * (1 - prev.excess_fat) / n),
    .groups = "drop"
  )

## # A tibble: 4 × 5
##   age.grp        n excess_fat prev.excess_fat se.excess_fat
##   <labelled> <int>      <int>           <dbl>         <dbl>
## 1 20-29         30         30               1             0
## 2 30-39         38         38               1             0
## 3 40-49         79         79               1             0
## 4 50-59        109        109               1             0

# By Race
vos.bmi %>%
  group_by(race) %>%
  summarize(
    n = n(),
    excess_fat = sum(excess.fat == 1, na.rm = TRUE),
    prev.excess_fat = excess_fat / n,
    se.excess_fat = sqrt(prev.excess_fat * (1 - prev.excess_fat) / n),
    .groups = "drop"
  )

## # A tibble: 2 × 5
##   race        n excess_fat prev.excess_fat se.excess_fat
##   <chr>   <int>      <int>           <dbl>         <dbl>
## 1 Chinese     8          8               1             0
## 2 Viet      248        248               1             0

2.5 Relationship between BMI and pcfat by NHANES cycles

p1 = ggplot(data = vos, aes(x = bmi, y = pcfat, col = sex)) + geom_point() + geom_smooth() + theme_bw() + theme(legend.position="none") + labs(title = "BMI-Fat relationship", x = "BMI (kg/m2)", y = "Percentage of fat (%)")

p2 = ggplot(data = vos, aes(x = bmi, y = wc, col = sex)) + geom_point() + geom_smooth() + theme_bw() + theme(legend.position="none") + labs(title = "BMI-Waist relationship", x = "BMI (kg/m2)", y = "Waist circumference (cm)")

p3 = ggplot(data = vos, aes(x = wc, y = pcfat, col = sex)) + geom_point() + geom_smooth() + theme_bw() + theme(legend.position="none") + labs(title = "Waist circ-Fat relationship", x = "Waist circumference (cm)", y = "Percentage of fat (%)")

grid.arrange(p1, p2, p3, nrow = 1)

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 4 rows containing missing values (`geom_point()`).

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

## Warning: Removed 7 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 7 rows containing missing values (`geom_point()`).

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

## Warning: Removed 7 rows containing non-finite values (`stat_smooth()`).
## Removed 7 rows containing missing values (`geom_point()`).

Estimated percent body fat for BMI threshold of 27.5 kg/m2

pcfat = -18.9 - 10.9Men + 0.044Age + 3.473BMI - 0.051BMI^2

# Aged 50:

## Men
-18.9 - 10.9*1 + 0.044*50 + 3.473*27.5 - 0.051*27.5^2

## [1] 29.33875

### BMI= 27.5 kg/m2 equivalent to Percent body fat of 29.3%

## Women
-18.9 - 10.9*0 + 0.044*50 + 3.473*27.5 - 0.051*27.5^2

## [1] 40.23875

### BMI= 27.5 kg/m2 equivalent to Percent body fat of 40.2%