This post shows, side by side, the Python and R code to import data from the US National Health and Nutrition Examination Survey (NHANES) and run a linear and a logistic regression. NHANES is a complex survey with sampling weights, so an estimate meant for the US population has to account for its design. I first fit the models the naive way, treating every row equally, then add the survey design so the standard errors and coefficients are correct for the population.
The R code needs three packages:
library(haven) # read SAS .xpt files
library(dplyr) # data wrangling
library(survey) # complex-survey models
Import the data
Each NHANES component is a SAS transport file (.xpt), which both languages read directly from the CDC for this post. I pull four files from the 2017 to 2018 cycle: demographics (which also holds the survey-design columns), body measures, blood pressure, and glycohemoglobin (HbA1c).
Python uses pandas.read_sas:
import pandas as pd
import numpy as np
import urllib.request
import io
base = "https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2017/DataFiles/"
def nhanes(f):
req = urllib.request.Request(base + f + ".xpt",
headers={"User-Agent": "Mozilla/5.0"})
raw = urllib.request.urlopen(req).read()
return pd.read_sas(io.BytesIO(raw), format="xport")
demo = nhanes("DEMO_J")
bmx = nhanes("BMX_J")
bpx = nhanes("BPX_J")
ghb = nhanes("GHB_J")
demo.shape
## (9254, 46)
R uses haven::read_xpt:
base <- "https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2017/DataFiles/"
nhanes <- function(f) read_xpt(paste0(base, f, ".xpt"))
demo <- nhanes("DEMO_J")
bmx <- nhanes("BMX_J")
bpx <- nhanes("BPX_J")
ghb <- nhanes("GHB_J")
dim(demo)
## [1] 9254 46
Merge and prepare the variables
Every file shares an id SEQN. I average the four blood-pressure readings, left-join the components onto demographics, rename the predictors, code sex as a label, and derive a diabetes outcome (HbA1c of 6.5% or higher). A missing HbA1c stays missing rather than counting as "no diabetes".
Python:
bpx = bpx.assign(SBP=bpx[["BPXSY1", "BPXSY2", "BPXSY3", "BPXSY4"]].mean(axis=1))
df = (demo[["SEQN", "RIAGENDR", "RIDAGEYR",
"WTMEC2YR", "SDMVPSU", "SDMVSTRA"]]
.merge(bmx[["SEQN", "BMXBMI"]], on="SEQN", how="left")
.merge(bpx[["SEQN", "SBP"]], on="SEQN", how="left")
.merge(ghb[["SEQN", "LBXGH"]], on="SEQN", how="left"))
df = df.rename(columns={"RIDAGEYR": "age", "BMXBMI": "bmi", "SBP": "sbp"})
df = df.assign(
sex = np.where(df["RIAGENDR"] == 1, "Male", "Female"),
diabetes = np.where(df["LBXGH"] >= 6.5, 1.0,
np.where(df["LBXGH"].notna(), 0.0, np.nan)),
)
adult = df[df["age"] >= 20]
adult[["age", "sex", "bmi", "sbp", "diabetes"]].head()
## age sex bmi sbp diabetes
## 2 66.0 Female 31.7 200.000000 0.0
## 5 66.0 Female 23.7 142.000000 0.0
## 6 75.0 Female 38.9 118.666667 0.0
## 8 56.0 Male 21.3 101.333333 0.0
## 10 67.0 Male 23.5 104.666667 0.0
R:
bpx <- bpx |>
mutate(
SBP = rowMeans(across(c(BPXSY1, BPXSY2, BPXSY3, BPXSY4)), na.rm = TRUE)
)
df <- demo |>
select(SEQN, RIAGENDR, RIDAGEYR, WTMEC2YR, SDMVPSU, SDMVSTRA) |>
left_join(select(bmx, SEQN, BMXBMI), by = "SEQN") |>
left_join(select(bpx, SEQN, SBP), by = "SEQN") |>
left_join(select(ghb, SEQN, LBXGH), by = "SEQN") |>
rename(age = RIDAGEYR, bmi = BMXBMI, sbp = SBP) |>
mutate(
sex = if_else(RIAGENDR == 1, "Male", "Female"),
diabetes = case_when(LBXGH >= 6.5 ~ 1, !is.na(LBXGH) ~ 0, TRUE ~ NA_real_)
)
adult <- filter(df, age >= 20)
head(adult[c("age", "sex", "bmi", "sbp", "diabetes")])
## # A tibble: 6 × 5
## age sex bmi sbp diabetes
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 66 Female 31.7 200 0
## 2 66 Female 23.7 142 0
## 3 75 Female 38.9 119. 0
## 4 56 Male 21.3 101. 0
## 5 67 Male 23.5 105. 0
## 6 54 Female 39.9 162 1
Linear regression
I model systolic blood pressure on age, BMI, and sex. Both languages use a formula with a categorical sex (Female is the reference).
Python (statsmodels, missing="drop" drops rows with missing data):
import statsmodels.formula.api as smf
lin = smf.ols("sbp ~ age + bmi + C(sex)", data=adult, missing="drop").fit()
lin.params.round(3)
## Intercept 87.734
## C(sex)[T.Male] 1.691
## age 0.528
## bmi 0.378
## dtype: float64
R (lm drops incomplete rows by default):
lin <- lm(sbp ~ age + bmi + sex, data = adult)
round(coef(lin), 3)
## (Intercept) age bmi sexMale
## 87.734 0.528 0.378 1.691
lin.summary() in Python and summary(lin) in R print the full table with standard errors and p-values.
Both languages report about 0.53 mmHg per year of age, 0.38 per BMI point, and 1.7 higher for men.
Logistic regression
Same formula for the binary diabetes outcome. The coefficients are on the log-odds scale, so I exponentiate them to read odds ratios.
Python (smf.logit):
logit = smf.logit("diabetes ~ age + bmi + C(sex)",
data=adult, missing="drop").fit(disp=0)
np.exp(logit.params).round(3)
## Intercept 0.001
## C(sex)[T.Male] 1.283
## age 1.052
## bmi 1.068
## dtype: float64
R (glm with family = binomial):
logit <- glm(diabetes ~ age + bmi + sex, data = adult, family = binomial)
round(exp(coef(logit)), 3)
## (Intercept) age bmi sexMale
## 0.001 1.052 1.068 1.283
Both languages report an odds ratio of about 1.05 per year of age, 1.07 per BMI point, and 1.28 for men versus women.
Add the survey weights
The fits above treat every row equally. NHANES is a complex survey: it oversamples some groups and samples in clusters, so any estimate meant for the US population needs the sampling weight (WTMEC2YR), the stratum (SDMVSTRA), and the cluster (SDMVPSU). In both languages I describe the design once, then reuse it, and I build it on the full data and restrict to adults afterward so the design stays intact for correct standard errors.
Python uses the svy package (built on polars). where= restricts to the adult domain, and family is "gaussian" for the linear model and "binomial" for the logistic one:
import polars as pl
import svy
design = svy.Design(stratum="SDMVSTRA", psu="SDMVPSU", wgt="WTMEC2YR")
sample = svy.Sample(pl.from_pandas(df), design)
lin_w = sample.glm.fit(y="sbp", x=["age", "bmi", svy.Cat("sex", ref="Female")],
family="gaussian", where=svy.col("age") >= 20)
{c.term: round(float(c.est), 3) for c in lin_w.coefs}
## {'_intercept_': 87.602, 'age': 0.458, 'bmi': 0.418, 'sex_Male': 3.1}
log_w = sample.glm.fit(y="diabetes", x=["age", "bmi", svy.Cat("sex", ref="Female")],
family="binomial", where=svy.col("age") >= 20)
{c.term: round(float(np.exp(c.est)), 3) for c in log_w.coefs}
## {'_intercept_': 0.0, 'age': 1.057, 'bmi': 1.08, 'sex_Male': 1.552}
R uses the survey package. svydesign sets up the design, subset restricts to adults, and svyglm fits the models (quasibinomial for logistic to avoid a harmless weights warning):
options(survey.lonely.psu = "adjust")
des <- svydesign(
ids = ~SDMVPSU,
strata = ~SDMVSTRA,
weights = ~WTMEC2YR,
nest = TRUE,
data = df
)
des_ad <- subset(des, age >= 20)
lin_w <- svyglm(sbp ~ age + bmi + sex, design = des_ad)
round(coef(lin_w), 3)
## (Intercept) age bmi sexMale
## 87.602 0.458 0.418 3.100
log_w <- svyglm(
diabetes ~ age + bmi + sex,
design = des_ad,
family = quasibinomial()
)
round(exp(coef(log_w)), 3)
## (Intercept) age bmi sexMale
## 0.000 1.057 1.080 1.552
The weighted estimates are population estimates with design-based standard errors, and the two toolkits agree: the male-versus-female blood-pressure gap is about 3.1 mmHg and the diabetes odds ratio for men is about 1.55, in both Python and R.
That’s it. I hope you find this tutorial useful, and if you have questions, leave a comment below.
