<jupyter_start><jupyter_text>Statistical analysis on NHANES data collected between 1999 and 2006 Before starting the process, let's have a look to the **software environment** :<jupyter_code>sessionInfo()<jupyter_output><empty_output><jupyter_text>We can see that the `R version` in the Jupyter notebook is `3.4.1` and the OS is `Ubuntu 18.04 LTS`. The `Jupyter notebook version` not displayed in the output is `5.6.0`. These elements might be important if a user tries to re-create the exact software-environment. Also for each library that will be loaded, the versions will be specified for the same reason.  The present notebook is divided in six parts : 1. Connect to the MRAN repository and load the packages  2. Dataset provenance 3. Load the data from the local `nhanes_extracted` file 4. Data management 5. Descriptive statistics 6. Hypothesis testing  Assuming you are currently using either Jupyter or mybinder and want to reproduce the workflow, you need to erase all the cells outputs and then execute the code again. In order to so, click on Kernel > Restart & clear output. Then you can execute the code by selecting a cell and clicking on the run button (or by using the appropriate shortcut mentioned in the Help section). You can interact with the code and edit it all along the notebook._More information on how to use Jupyter at http://jupyter.org/documentation_  1) Connect to the M.R.A.N repository and load the packages  C.R.A.N (Comprehensive R Archive Network) is a package repository of the programming language R. Packages tend to get many updates over time and M.R.A.N (Microsoft R Application Network) enables to retrieve packages with 'frozen' versions corresponding to a specific date (in our illustration the date is set to 2018-09-14). MRAN basically consists of daily snapshots of CRAN. Extracting versions of different packages set on a specific date facilitates consistency between packages and therefore reproducibility. More information on : https://cran.r-project.org/ and https://mran.microsoft.com/ .The line of code below sets the repository to MRAN instead of the default CRAN. (_When executing the notebook on mybinder it is facultative to run the next cell since the `runtime.txt` file in the Github repository already specifies it_)<jupyter_code>options(repos = c(CRAN = "https://mran.revolutionanalytics.com/snapshot/2018-09-14"))<jupyter_output><empty_output><jupyter_text>We can check with the line of code below that the repository used to install packages is correctly set to MRAN. You should get a URL that contains 'mran' and the date mentioned above.<jupyter_code>getOption("repos")<jupyter_output><empty_output><jupyter_text>If running the notebook in Mybinder,there is no need to execute the `install.packages()` function since the installation of every packages is already done by the `install.R` file from the Github repository (except the `Rcupcake` package which is requiring the `install_github()`function). However, we did mention the line of code necessary for the installation of packages after the `` symbol in case of using a local Jupyter notebook.  dplyr package<jupyter_code>#install.packages("dplyr")
library(dplyr)
packageVersion("dplyr")<jupyter_output>Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union<jupyter_text>prettyR package<jupyter_code>#install.packages("prettyR")
library(prettyR)
packageVersion("prettyR")<jupyter_output><empty_output><jupyter_text>2) Dataset provenance The `NHANES_extracted.csv` which is published in the same repository, is NHANES data collected between 1999 and 2006 with only the seven variables of interest for this work (female, male, age, BMI, waist circumference and sugar-free diet). It was extracted from the NHANES dataset (`MainTable.csv file`) archived in dryad digital repository at http://dx.doi.org/10.5061/dryad.d5h62. This large dataset containing data collected between 1999 and 2006, was reassembled by the PIC-SURE team from CDC NHANES files (https://www.cdc.gov/nchs/nhanes/index.htm) to form what is called the 'NHANES unified dataset'.This second section describes how we extracted the dataset (`NHANES_extracted.csv`) which we will be working on. **The extraction process itself is voluntary displayed in 'markdown cell' meaning that it is not executable because we do not provide the whole dataset in the Github repository. However you are free to reproduce the whole steps after downloading the starting point dataset if needed.**_For more information on NHANES unified dataset, see the article 'A database of human exposomes and phenomes from the US National Health and Nutrition Examination Survey' on https://www.nature.com/articles/sdata201696_  a) Load the whole dataset in the object `nhanes_dataset` nhanes_dataset <- read.table("MainTable.csv", sep=",", header = TRUE)  b) Create a smaller dataset with only the variables of interest : `nhanes_extracted` nhanes_extracted % select('SEQN', 'female', 'male', 'RIDAGEYR', 'BMXBMI', 'BMXWAIST', 'DRDSDT4')  c) Save the dataset in `NHANES_extracted.csv` file write.csv(nhanes_extracted, file = "NHANES_extracted.csv")  3) Load the data from the local `nhanes_extracted` file Make sure the `nhanes_extracted.csv` file is in the directory you are working in, and execute the following line of code :<jupyter_code>nhanes_extracted <- read.csv("NHANES_extracted.csv")<jupyter_output><empty_output><jupyter_text>Here's a view of the six first rows of `nhanes_extracted` :<jupyter_code>head(nhanes_extracted)<jupyter_output><empty_output><jupyter_text>`X` and `SEQN` being redundant, we can delete the `X` column :<jupyter_code>nhanes_extracted <- nhanes_extracted[-1]<jupyter_output><empty_output><jupyter_text>4) Data management  a) Basic quick overview of variables structure<jupyter_code>glimpse(nhanes_extracted)<jupyter_output>Observations: 41,474
Variables: 7
$ SEQN     <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,...
$ female   <int> 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1...
$ male     <int> 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0...
$ RIDAGEYR <int> 2, 77, 10, 1, 49, 19, 59, 13, 11, 43, 15, 37, 70, 81, 38, ...
$ BMXBMI   <dbl> 14.90, 24.90, 17.63, NA, 29.10, 22.56, 29.39, 15.51, 18.48...
$ BMXWAIST <dbl> 45.7, 98.0, 64.7, NA, 99.9, 81.6, 90.7, 64.1, 64.6, 108.0,...
$ DRDSDT4  <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...<jupyter_text>We can notice that there are two variables for the sex variable (`female` and `male`) and some variables names are not intuitive. Also the variable adressing the sugar-fee diet `DRDSDT4` (yes or no) is considered as integer instead of categorical. We will then modify these in the next steps.  b) Modify variables names to make it more understandable _Variables correspondence are findable in the dryad digital repository mentioned in section 2 and in cdc.gov website (https://www.cdc.gov/nchs/nhanes/index.htm)._<jupyter_code>variables <- c('index', 'female', 'male', 'age', 'body_mass_index', 'waist_circumference', 'sugarfree_diet' )

names(nhanes_extracted) <- variables

head(nhanes_extracted)<jupyter_output><empty_output><jupyter_text>c) Transform the variable `sugarfree_diet` into a categorical variable and label it<jupyter_code>nhanes_extracted$sugarfree_diet <- factor(nhanes_extracted$sugarfree_diet, levels=c('0','1'), labels=c('No','Yes'))<jupyter_output><empty_output><jupyter_text>d)  Create one categorical variable `sex` First let's check that `female` and `male` are mutually exclusive as expected:<jupyter_code>table(nhanes_extracted$female, nhanes_extracted$male, deparse.level = 2, useNA = 'always')<jupyter_output><empty_output><jupyter_text>Now let's creat a single categorical variable `sex` :<jupyter_code>nhanes_extracted <- nhanes_extracted %>% mutate(sex = ifelse(female == 1, "female", "male")) %>% select(-c('female', 'male'))
nhanes_extracted$sex <- as.factor(nhanes_extracted$sex)<jupyter_output><empty_output><jupyter_text>5) Descriptive statistics Here's a view of the six first rows of `nhanes_extracted` now that we have modified the variables:<jupyter_code>head(nhanes_extracted)<jupyter_output><empty_output><jupyter_text>a) Description of variables The following meaning of the variables are provided according to the CDC NHANES website (https://www.cdc.gov/nchs/nhanes/index.htm):- `index` : index number of patients- `sex` : female or male- `age` : age of individuals in years- `body_mass_index` : Body mass index of individual in kg/m^2- `waist_circumference` : waist circumference of indviduals in cm- `sugarfree_diet` : yes (the individual has a sugar-free diet) / no (the individual does not have a sugar free diet) Below you can see a basic statistics summary of our dataset :<jupyter_code>describe(nhanes_extracted[, -1], num.desc = c("mean", "sd", "median", "min", "max", "valid.n"))<jupyter_output>Description of nhanes_extracted[, -1]<jupyter_text>b) Gender distribution<jupyter_code>gender_distribution <- table(nhanes_extracted$sex)

gender_distribution
gender_proportion<- round(100*gender_distribution/sum(gender_distribution), 1)

gender_proportion<jupyter_output><empty_output><jupyter_text>There are 21210 females (51.1%) and 20264 males (48.9%) in this sample.<jupyter_code>gender_percent <- paste(gender_proportion,"%",sep="")

pie(gender_distribution, labels = gender_percent, main = "Gender distribution", col = c("pink", "blue"))
legend("topright", c("female", "male"), fill = c("pink", "blue"))<jupyter_output><empty_output><jupyter_text>c) Age according to gender distribution<jupyter_code>nhanes_extracted %>% group_by(sex)%>% summarise(first_quantile_age = quantile(age, probs = 0.25, na.rm = TRUE), median_age = median(age, na.rm = TRUE), third_quantile_age = quantile(age, probs = 0.75, na.rm = TRUE))

boxplot(age~sex, data = nhanes_extracted, main = "Boxplot of age distribution by gender", xlab = "Gender", ylab = "Age in years", col = c("pink", "blue"))<jupyter_output><empty_output><jupyter_text>The output above shows similar distributions of age in males and females in this sample.  d) Body mass index (BMI) according to gender<jupyter_code>nhanes_extracted %>% group_by(sex)%>% summarise(first_quantile_bmi = quantile(body_mass_index, probs = 0.25, na.rm = TRUE), median_bmi = median(body_mass_index, na.rm = TRUE), third_quantile_bmi = quantile(body_mass_index, probs = 0.75, na.rm = TRUE))

boxplot(body_mass_index~sex, data = nhanes_extracted, main = "Boxplot of BMI distibution by gender", xlab = "Gender", ylab = "BMI in kg/m^2", col = c("pink", "blue"))<jupyter_output><empty_output><jupyter_text>BMI distribution of females is slightly more widespread than males distribution.  6) Hypothesis testing  We are here exploring if there is a relationship between having a sugarfree diet and the waist circumference in this NHANES sample. In this study, measurements are performed on individuals over 2 years old. For better clarity, let's subset a new dataset comprising only those individuals.<jupyter_code>nhanes_over_two <- nhanes_extracted  %>% filter(age >=2)
#quick check of the number of observations and variables
dim(nhanes_over_two)<jupyter_output><empty_output><jupyter_text>We have now generated our new dataset with 38 124 individuals (over 2 years old) and the same 6 variables from the previous dataset.You can see below a basic statistics summary :<jupyter_code>describe(nhanes_over_two[, -1], num.desc = c("mean", "sd", "median", "min", "max", "valid.n"))<jupyter_output>Description of nhanes_over_two[, -1]<jupyter_text>We can see that there are missing values in every variables except from `age` and `sex`.**Missing values in `sugarfree_diet` and `waist_circumference` represent more than 50% and about 10% of the observations respectively.**  a) What are the means of waist circumference in both groups ?<jupyter_code>nhanes_over_two %>% group_by(sugarfree_diet)%>% summarise(count = n(), mean = mean(waist_circumference, na.rm = TRUE))<jupyter_output><empty_output><jupyter_text>**3716** observations with missing values in `waist_circumference` were exluded (`na.rm = TRUE`) in order to compute means (_run the code below to get this number_).<jupyter_code>nhanes_over_two %>% filter(is.na(nhanes_over_two$waist_circumference))%>% nrow()<jupyter_output><empty_output><jupyter_text>For the purpose of this work, we will assume that all the missing data are **MCAR** (missing completely at random), therefore we will not impute any missing data.  b) Test the hypothesis ( wilcoxon test ) Let's test the association between having a sugar-free diet and the waist circumference with the Wilcoxon test :<jupyter_code>wilcox.test(waist_circumference~sugarfree_diet, conf.int=TRUE, data = nhanes_over_two)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09988358556461001, 'fraction_numerical': 0.023282887077997673, 'mean_word_length': 5.905680600214362, 'pattern_counts': {'":': 0, '<': 112, '<?xml version=': 0, '>': 113, 'https://': 7, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '373532', 'license': 'Apache-2.0', 'n_tokens_mistral': 4367, 'n_tokens_neox': 3942, 'n_words': 1492, 'path': 'Statistical_analysis_on_NHANES_data_(1999-2006).ipynb', 'repo_name': 'hms-dbmi/Reproducible_and_FAIR_statistical_analysis', 'size': '12885'}