NHANES Diabetes EWAS: Reproduction

This document replicates the main findings from Table 1 of the following publication: Patel, C.J., Bhattacharya, J. & Butte, A.J. (2010) An Environment-Wide Association Study (EWAS) on Type 2 Diabetes Mellitus. PLoS ONE. [Online] 5 (5). Available from: doi:10.1371/journal.pone.0010746

Data from https://github.com/chiragjp/nhanes_scidata

For a general overview of weighted survey regression see: https://github.com/chiragjp/nhanes_scidata/blob/master/User_Guide.Rmd

Normalization of relevant variables

“Most chemical exposure data arising from mass spectrometry or absorption measurements occurred within a very small range and had a right skew; thus, we log transformed these variables. Further, we applied a z-score transformation (adjusting each observation to the mean and scaling by the standard deviation) in order to compare odds ratios from the many regressions.”

load('nh_99-06.Rdata')

# Log transform and then normalise (later called through regression formula)
logNorm = function(x) {
  x = log(x)
  return((x - mean(x, na.rm=T)) / sd(x, na.rm=T))
}

Survey-weighted adjusted Log regression

Investigate NHANES 1999/2000, 2001/2002, 2003/2004 and 2005/2006 cohorts together.

Perform survey-weighted logistic regression for each of the following exposures from the above paper against diabetes (as defined by fasting blood sugar >= 126mg/dl): ‘cis-beta-carotene’, ‘trans-beta-carotene’, ‘gamma-tocopherol’, ‘Heptachlor Epoxide’ and ‘PCB170’

Adjust for BMI, age, sex, socio-economic-status and ethnicity. The numbers are nearly the same as in Table 1 from the above paper.

suppressMessages(library(survey))

nhanesLogRegDiab = function(exposure_var, adjust_for, data) {
  dsn = svydesign(id=~SDMVPSU, strata=~SDMVSTRA, weights=~WTMEC2YR, nest=T, data=data)
  
  form = paste("I(LBXGLU >=126) ~ logNorm(", exposure_var, ") + ", paste(adjust_for, collapse="+"))
  
  mod = svyglm(as.formula(form), design=dsn, family=quasibinomial())
  
  return(mod)
}

exposures_of_interest = c('LBXCBC', 'LBXBEC', 'LBXGTC', 'LBXHPE', 'LBX170')
names = list('LBXCBC'='cis-beta-carotene', 'LBXBEC'='trans-beta-carotene', 'LBXGTC'='gamma-tocopherol', 'LBXHPE'='Heptachlor Epoxide', 'LBX170'='PCB170')

#subset for one cohort (1: 99/00, 2: 01/02, 3: 03/04, 4: 05/06)
#dat = subset(MainTable, WTMEC2YR > 0 & SDDSRVYR == 3)

# take all cohorts together
dat = subset(MainTable, WTMEC2YR > 0)
findings = data.frame(Cohort = c("2001-2006*", "2001-2006*", "1999-2006*", "1999-2004*", "1999-2004*"), row.names = exposures_of_interest)

# Replicate key numbers of Table 1
for (exposure in exposures_of_interest) {
  adjust_for = c("BMXBMI", "RIDAGEYR", "female", "SES_LEVEL", "black", "white", "mexican", "other_hispanic", "other_eth", "SDDSRVYR") # adjusting for cohort (SDDSRVYR) only matters when taking all cohorts, when subsetting for one its coef will be 1 anyway
  
  # Identify subset of samples with all necessary variables
  samples = dat[complete.cases(dat[, c(exposure, adjust_for)]),]
  
  # Run model
  mod = nhanesLogRegDiab(exposure, adjust_for, samples)
  
  # Save findings
  n_diab    = sum(samples$LBXGLU >= 126, na.rm=T)
  n_no_diab = sum(samples$LBXGLU <  126, na.rm=T)
  findings[exposure, "N T2D, No T2D"] = paste0(n_diab, ", ", n_no_diab)
  
  p = coef(summary(mod))[2,4]
  findings[exposure, "P"] = signif(p,2)
  
  or = exp(coef(mod)[2])
  or_ci_lo = exp(coef(mod)[2] - 1.96 * coef(summary(mod))[2,2])
  or_ci_hi = exp(coef(mod)[2] + 1.96 * coef(summary(mod))[2,2])
  
  findings[exposure, "OR (95% CI)"] = paste0(signif(or,2), " (",
                                             signif(or_ci_lo,2), "-",
                                             signif(or_ci_hi,2), ")")
}

rownames(findings) = names[exposures_of_interest]
findings

LS0tDQp0aXRsZTogIk5IQU5FUyBEaWFiZXRlcyBFV0FTOiBSZXByb2R1Y3Rpb24iDQphdXRob3I6ICJUaW0gV29lbGZsZSINCmRhdGU6ICIxMi8wMS8yMDE5Ig0Kb3V0cHV0Og0KICBodG1sX25vdGVib29rOg0KICAgIGNvZGVfZm9sZGluZzogaGlkZQ0KLS0tDQoNClRoaXMgZG9jdW1lbnQgcmVwbGljYXRlcyB0aGUgbWFpbiBmaW5kaW5ncyBmcm9tIFtUYWJsZSAxXShodHRwczovL3d3dy5uY2JpLm5sbS5uaWguZ292L3BtYy9hcnRpY2xlcy9QTUMyODczOTc4L3RhYmxlL3BvbmUtMDAxMDc0Ni10MDAxLz9yZXBvcnQ9b2JqZWN0b25seSkgb2YgdGhlIGZvbGxvd2luZyBwdWJsaWNhdGlvbjogW1BhdGVsLCBDLkouLCBCaGF0dGFjaGFyeWEsIEouICYgQnV0dGUsIEEuSi4gKDIwMTApICpBbiBFbnZpcm9ubWVudC1XaWRlIEFzc29jaWF0aW9uIFN0dWR5IChFV0FTKSBvbiBUeXBlIDIgRGlhYmV0ZXMgTWVsbGl0dXMuKiBQTG9TIE9ORS4gW09ubGluZV0gNSAoNSkuIEF2YWlsYWJsZSBmcm9tOiBkb2k6MTAuMTM3MS9qb3VybmFsLnBvbmUuMDAxMDc0Nl0oaHR0cHM6Ly9kb2kub3JnLzEwLjEzNzEvam91cm5hbC5wb25lLjAwMTA3NDYpDQoNCkRhdGEgZnJvbSBodHRwczovL2dpdGh1Yi5jb20vY2hpcmFnanAvbmhhbmVzX3NjaWRhdGENCg0KRm9yIGEgZ2VuZXJhbCBvdmVydmlldyBvZiB3ZWlnaHRlZCBzdXJ2ZXkgcmVncmVzc2lvbiBzZWU6ICBodHRwczovL2dpdGh1Yi5jb20vY2hpcmFnanAvbmhhbmVzX3NjaWRhdGEvYmxvYi9tYXN0ZXIvVXNlcl9HdWlkZS5SbWQNCg0KIyMgTm9ybWFsaXphdGlvbiBvZiByZWxldmFudCB2YXJpYWJsZXMNCg0KIk1vc3QgY2hlbWljYWwgZXhwb3N1cmUgZGF0YSBhcmlzaW5nIGZyb20gbWFzcyBzcGVjdHJvbWV0cnkgb3IgYWJzb3JwdGlvbiBtZWFzdXJlbWVudHMgb2NjdXJyZWQgd2l0aGluIGEgdmVyeSBzbWFsbCByYW5nZSBhbmQgaGFkIGEgcmlnaHQgc2tldzsgdGh1cywgd2UgbG9nIHRyYW5zZm9ybWVkIHRoZXNlIHZhcmlhYmxlcy4gRnVydGhlciwgd2UgYXBwbGllZCBhIHotc2NvcmUgdHJhbnNmb3JtYXRpb24gKGFkanVzdGluZyBlYWNoIG9ic2VydmF0aW9uIHRvIHRoZSBtZWFuIGFuZCBzY2FsaW5nIGJ5IHRoZSBzdGFuZGFyZCBkZXZpYXRpb24pIGluIG9yZGVyIHRvIGNvbXBhcmUgb2RkcyByYXRpb3MgZnJvbSB0aGUgbWFueSByZWdyZXNzaW9ucy4iDQoNCmBgYHtyfQ0KbG9hZCgnbmhfOTktMDYuUmRhdGEnKQ0KDQojIExvZyB0cmFuc2Zvcm0gYW5kIHRoZW4gbm9ybWFsaXNlIChsYXRlciBjYWxsZWQgdGhyb3VnaCByZWdyZXNzaW9uIGZvcm11bGEpDQpsb2dOb3JtID0gZnVuY3Rpb24oeCkgew0KICB4ID0gbG9nKHgpDQogIHJldHVybigoeCAtIG1lYW4oeCwgbmEucm09VCkpIC8gc2QoeCwgbmEucm09VCkpDQp9DQpgYGANCg0KIyMgU3VydmV5LXdlaWdodGVkIGFkanVzdGVkIExvZyByZWdyZXNzaW9uDQoNCkludmVzdGlnYXRlIE5IQU5FUyBbMTk5OS8yMDAwXShodHRwczovL3d3d24uY2RjLmdvdi9uY2hzL25oYW5lcy9Db250aW51b3VzTmhhbmVzL0RlZmF1bHQuYXNweD9CZWdpblllYXI9MTk5OSksIFsyMDAxLzIwMDJdKGh0dHBzOi8vd3d3bi5jZGMuZ292L25jaHMvbmhhbmVzL0NvbnRpbnVvdXNOaGFuZXMvRGVmYXVsdC5hc3B4P0JlZ2luWWVhcj0yMDAxKSwgWzIwMDMvMjAwNF0oaHR0cHM6Ly93d3duLmNkYy5nb3YvbmNocy9uaGFuZXMvQ29udGludW91c05oYW5lcy9EZWZhdWx0LmFzcHg/QmVnaW5ZZWFyPTIwMDMpIGFuZCBbMjAwNS8yMDA2XShodHRwczovL3d3d24uY2RjLmdvdi9uY2hzL25oYW5lcy9Db250aW51b3VzTmhhbmVzL0RlZmF1bHQuYXNweD9CZWdpblllYXI9MjAwNSkgY29ob3J0cyB0b2dldGhlci4NCg0KUGVyZm9ybSBzdXJ2ZXktd2VpZ2h0ZWQgbG9naXN0aWMgcmVncmVzc2lvbiBmb3IgZWFjaCBvZiB0aGUgZm9sbG93aW5nIGV4cG9zdXJlcyBmcm9tIHRoZSBhYm92ZSBwYXBlciBhZ2FpbnN0IGRpYWJldGVzIChhcyBkZWZpbmVkIGJ5IGZhc3RpbmcgYmxvb2Qgc3VnYXIgPj0gMTI2bWcvZGwpOiAnY2lzLWJldGEtY2Fyb3RlbmUnLCAndHJhbnMtYmV0YS1jYXJvdGVuZScsICdnYW1tYS10b2NvcGhlcm9sJywgJ0hlcHRhY2hsb3IgRXBveGlkZScgYW5kICdQQ0IxNzAnDQoNCkFkanVzdCBmb3IgQk1JLCBhZ2UsIHNleCwgc29jaW8tZWNvbm9taWMtc3RhdHVzIGFuZCBldGhuaWNpdHkuIFRoZSBudW1iZXJzIGFyZSBuZWFybHkgdGhlIHNhbWUgYXMgaW4gVGFibGUgMSBmcm9tIHRoZSBhYm92ZSBwYXBlci4NCg0KYGBge3J9DQpzdXBwcmVzc01lc3NhZ2VzKGxpYnJhcnkoc3VydmV5KSkNCg0KbmhhbmVzTG9nUmVnRGlhYiA9IGZ1bmN0aW9uKGV4cG9zdXJlX3ZhciwgYWRqdXN0X2ZvciwgZGF0YSkgew0KICBkc24gPSBzdnlkZXNpZ24oaWQ9flNETVZQU1UsIHN0cmF0YT1+U0RNVlNUUkEsIHdlaWdodHM9fldUTUVDMllSLCBuZXN0PVQsIGRhdGE9ZGF0YSkNCiAgDQogIGZvcm0gPSBwYXN0ZSgiSShMQlhHTFUgPj0xMjYpIH4gbG9nTm9ybSgiLCBleHBvc3VyZV92YXIsICIpICsgIiwgcGFzdGUoYWRqdXN0X2ZvciwgY29sbGFwc2U9IisiKSkNCiAgDQogIG1vZCA9IHN2eWdsbShhcy5mb3JtdWxhKGZvcm0pLCBkZXNpZ249ZHNuLCBmYW1pbHk9cXVhc2liaW5vbWlhbCgpKQ0KICANCiAgcmV0dXJuKG1vZCkNCn0NCg0KZXhwb3N1cmVzX29mX2ludGVyZXN0ID0gYygnTEJYQ0JDJywgJ0xCWEJFQycsICdMQlhHVEMnLCAnTEJYSFBFJywgJ0xCWDE3MCcpDQpuYW1lcyA9IGxpc3QoJ0xCWENCQyc9J2Npcy1iZXRhLWNhcm90ZW5lJywgJ0xCWEJFQyc9J3RyYW5zLWJldGEtY2Fyb3RlbmUnLCAnTEJYR1RDJz0nZ2FtbWEtdG9jb3BoZXJvbCcsICdMQlhIUEUnPSdIZXB0YWNobG9yIEVwb3hpZGUnLCAnTEJYMTcwJz0nUENCMTcwJykNCg0KI3N1YnNldCBmb3Igb25lIGNvaG9ydCAoMTogOTkvMDAsIDI6IDAxLzAyLCAzOiAwMy8wNCwgNDogMDUvMDYpDQojZGF0ID0gc3Vic2V0KE1haW5UYWJsZSwgV1RNRUMyWVIgPiAwICYgU0REU1JWWVIgPT0gMykNCg0KIyB0YWtlIGFsbCBjb2hvcnRzIHRvZ2V0aGVyDQpkYXQgPSBzdWJzZXQoTWFpblRhYmxlLCBXVE1FQzJZUiA+IDApDQpmaW5kaW5ncyA9IGRhdGEuZnJhbWUoQ29ob3J0ID0gYygiMjAwMS0yMDA2KiIsICIyMDAxLTIwMDYqIiwgIjE5OTktMjAwNioiLCAiMTk5OS0yMDA0KiIsICIxOTk5LTIwMDQqIiksIHJvdy5uYW1lcyA9IGV4cG9zdXJlc19vZl9pbnRlcmVzdCkNCg0KIyBSZXBsaWNhdGUga2V5IG51bWJlcnMgb2YgVGFibGUgMQ0KZm9yIChleHBvc3VyZSBpbiBleHBvc3VyZXNfb2ZfaW50ZXJlc3QpIHsNCiAgYWRqdXN0X2ZvciA9IGMoIkJNWEJNSSIsICJSSURBR0VZUiIsICJmZW1hbGUiLCAiU0VTX0xFVkVMIiwgImJsYWNrIiwgIndoaXRlIiwgIm1leGljYW4iLCAib3RoZXJfaGlzcGFuaWMiLCAib3RoZXJfZXRoIiwgIlNERFNSVllSIikgIyBhZGp1c3RpbmcgZm9yIGNvaG9ydCAoU0REU1JWWVIpIG9ubHkgbWF0dGVycyB3aGVuIHRha2luZyBhbGwgY29ob3J0cywgd2hlbiBzdWJzZXR0aW5nIGZvciBvbmUgaXRzIGNvZWYgd2lsbCBiZSAxIGFueXdheQ0KICANCiAgIyBJZGVudGlmeSBzdWJzZXQgb2Ygc2FtcGxlcyB3aXRoIGFsbCBuZWNlc3NhcnkgdmFyaWFibGVzDQogIHNhbXBsZXMgPSBkYXRbY29tcGxldGUuY2FzZXMoZGF0WywgYyhleHBvc3VyZSwgYWRqdXN0X2ZvcildKSxdDQogIA0KICAjIFJ1biBtb2RlbA0KICBtb2QgPSBuaGFuZXNMb2dSZWdEaWFiKGV4cG9zdXJlLCBhZGp1c3RfZm9yLCBzYW1wbGVzKQ0KICANCiAgIyBTYXZlIGZpbmRpbmdzDQogIG5fZGlhYiAgICA9IHN1bShzYW1wbGVzJExCWEdMVSA+PSAxMjYsIG5hLnJtPVQpDQogIG5fbm9fZGlhYiA9IHN1bShzYW1wbGVzJExCWEdMVSA8ICAxMjYsIG5hLnJtPVQpDQogIGZpbmRpbmdzW2V4cG9zdXJlLCAiTiBUMkQsIE5vIFQyRCJdID0gcGFzdGUwKG5fZGlhYiwgIiwgIiwgbl9ub19kaWFiKQ0KICANCiAgcCA9IGNvZWYoc3VtbWFyeShtb2QpKVsyLDRdDQogIGZpbmRpbmdzW2V4cG9zdXJlLCAiUCJdID0gc2lnbmlmKHAsMikNCiAgDQogIG9yID0gZXhwKGNvZWYobW9kKVsyXSkNCiAgb3JfY2lfbG8gPSBleHAoY29lZihtb2QpWzJdIC0gMS45NiAqIGNvZWYoc3VtbWFyeShtb2QpKVsyLDJdKQ0KICBvcl9jaV9oaSA9IGV4cChjb2VmKG1vZClbMl0gKyAxLjk2ICogY29lZihzdW1tYXJ5KG1vZCkpWzIsMl0pDQogIA0KICBmaW5kaW5nc1tleHBvc3VyZSwgIk9SICg5NSUgQ0kpIl0gPSBwYXN0ZTAoc2lnbmlmKG9yLDIpLCAiICgiLA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgc2lnbmlmKG9yX2NpX2xvLDIpLCAiLSIsDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzaWduaWYob3JfY2lfaGksMiksICIpIikNCn0NCg0Kcm93bmFtZXMoZmluZGluZ3MpID0gbmFtZXNbZXhwb3N1cmVzX29mX2ludGVyZXN0XQ0KZmluZGluZ3MNCmBgYA0K

NHANES Diabetes EWAS: Reproduction

Tim Woelfle

12/01/2019

Normalization of relevant variables

Survey-weighted adjusted Log regression