Tutorial 4.1: Exploring the NHANES dataset

In this tutorial, we will learn how to import, tidy, wrangle and visualize data. We will work with one specific dataset;

1 The NHANES dataset

The National Health and Nutrition Examination Survey (NHANES) contains data that has been collected since 1960. For this tutorial, we will make use of the data that were collected between 2009 and 2012, for 10.000 U.S. civilians. The dataset contains a large number of physical, demographic, nutritional and life-style-related parameters.

However, before we can actually start working with data, we will need to learn how to import the required datasets into our Rstudio environment.

2 Data Import with the readr R package

library(readr)

Let’s try reading in some data. We will begin by reading in the NHANES.csv dataset.

NHANES <- read_csv("https://raw.githubusercontent.com/GTPB/PSLS20/master/data/NHANES.csv")

## Parsed with column specification:
## cols(
##   .default = col_double(),
##   SurveyYr = col_character(),
##   Gender = col_character(),
##   AgeDecade = col_character(),
##   Race1 = col_character(),
##   Race3 = col_logical(),
##   Education = col_character(),
##   MaritalStatus = col_character(),
##   HHIncome = col_character(),
##   HomeOwn = col_character(),
##   Work = col_character(),
##   BMICatUnder20yrs = col_logical(),
##   BMI_WHO = col_character(),
##   Testosterone = col_logical(),
##   Diabetes = col_character(),
##   HealthGen = col_character(),
##   LittleInterest = col_character(),
##   Depressed = col_character(),
##   SleepTrouble = col_character(),
##   PhysActive = col_character(),
##   TVHrsDay = col_logical()
##   # ... with 12 more columns
## )

## See spec(...) for full column specifications.

## Warning: 20122 parsing failures.
##  row              col           expected     actual                                                                   file
## 5001 Race3            1/0/T/F/TRUE/FALSE Asian      'https://raw.githubusercontent.com/GTPB/PSLS20/master/data/NHANES.csv'
## 5001 BMICatUnder20yrs 1/0/T/F/TRUE/FALSE NormWeight 'https://raw.githubusercontent.com/GTPB/PSLS20/master/data/NHANES.csv'
## 5001 Testosterone     1/0/T/F/TRUE/FALSE 274.95     'https://raw.githubusercontent.com/GTPB/PSLS20/master/data/NHANES.csv'
## 5001 TVHrsDay         1/0/T/F/TRUE/FALSE 2_hr       'https://raw.githubusercontent.com/GTPB/PSLS20/master/data/NHANES.csv'
## 5001 CompHrsDay       1/0/T/F/TRUE/FALSE 3_hr       'https://raw.githubusercontent.com/GTPB/PSLS20/master/data/NHANES.csv'
## .... ................ .................. .......... ......................................................................
## See problems(...) for more details.

head(NHANES) ## take a look at the first rows of the dataset

## # A tibble: 6 x 76
##      ID SurveyYr Gender   Age AgeDecade AgeMonths Race1 Race3 Education
##   <dbl> <chr>    <chr>  <dbl> <chr>         <dbl> <chr> <lgl> <chr>    
## 1 51624 2009_10  male      34 30-39           409 White NA    High Sch…
## 2 51624 2009_10  male      34 30-39           409 White NA    High Sch…
## 3 51624 2009_10  male      34 30-39           409 White NA    High Sch…
## 4 51625 2009_10  male       4 0-9              49 Other NA    <NA>     
## 5 51630 2009_10  female    49 40-49           596 White NA    Some Col…
## 6 51638 2009_10  male       9 0-9             115 White NA    <NA>     
## # … with 67 more variables: MaritalStatus <chr>, HHIncome <chr>,
## #   HHIncomeMid <dbl>, Poverty <dbl>, HomeRooms <dbl>, HomeOwn <chr>,
## #   Work <chr>, Weight <dbl>, Length <dbl>, HeadCirc <dbl>, Height <dbl>,
## #   BMI <dbl>, BMICatUnder20yrs <lgl>, BMI_WHO <chr>, Pulse <dbl>,
## #   BPSysAve <dbl>, BPDiaAve <dbl>, BPSys1 <dbl>, BPDia1 <dbl>, BPSys2 <dbl>,
## #   BPDia2 <dbl>, BPSys3 <dbl>, BPDia3 <dbl>, Testosterone <lgl>,
## #   DirectChol <dbl>, TotChol <dbl>, UrineVol1 <dbl>, UrineFlow1 <dbl>,
## #   UrineVol2 <dbl>, UrineFlow2 <dbl>, Diabetes <chr>, DiabetesAge <dbl>,
## #   HealthGen <chr>, DaysPhysHlthBad <dbl>, DaysMentHlthBad <dbl>,
## #   LittleInterest <chr>, Depressed <chr>, nPregnancies <dbl>, nBabies <dbl>,
## #   Age1stBaby <dbl>, SleepHrsNight <dbl>, SleepTrouble <chr>,
## #   PhysActive <chr>, PhysActiveDays <dbl>, TVHrsDay <lgl>, CompHrsDay <lgl>,
## #   TVHrsDayChild <dbl>, CompHrsDayChild <dbl>, Alcohol12PlusYr <chr>,
## #   AlcoholDay <dbl>, AlcoholYear <dbl>, SmokeNow <chr>, Smoke100 <chr>,
## #   Smoke100n <chr>, SmokeAge <dbl>, Marijuana <chr>, AgeFirstMarij <dbl>,
## #   RegularMarij <chr>, AgeRegMarij <dbl>, HardDrugs <chr>, SexEver <chr>,
## #   SexAge <dbl>, SexNumPartnLife <dbl>, SexNumPartYear <dbl>, SameSex <chr>,
## #   SexOrientation <chr>, PregnantNow <chr>

tail(NHANES) ## take a look at the last rows of the dataset

## # A tibble: 6 x 76
##      ID SurveyYr Gender   Age AgeDecade AgeMonths Race1 Race3 Education
##   <dbl> <chr>    <chr>  <dbl> <chr>         <dbl> <chr> <lgl> <chr>    
## 1 71909 2011_12  male      28 20-29            NA Mexi… NA    9 - 11th…
## 2 71909 2011_12  male      28 20-29            NA Mexi… NA    9 - 11th…
## 3 71910 2011_12  female     0 0-9               5 White NA    <NA>     
## 4 71911 2011_12  male      27 20-29            NA Mexi… NA    College …
## 5 71915 2011_12  male      60 60-69            NA White NA    College …
## 6 71915 2011_12  male      60 60-69            NA White NA    College …
## # … with 67 more variables: MaritalStatus <chr>, HHIncome <chr>,
## #   HHIncomeMid <dbl>, Poverty <dbl>, HomeRooms <dbl>, HomeOwn <chr>,
## #   Work <chr>, Weight <dbl>, Length <dbl>, HeadCirc <dbl>, Height <dbl>,
## #   BMI <dbl>, BMICatUnder20yrs <lgl>, BMI_WHO <chr>, Pulse <dbl>,
## #   BPSysAve <dbl>, BPDiaAve <dbl>, BPSys1 <dbl>, BPDia1 <dbl>, BPSys2 <dbl>,
## #   BPDia2 <dbl>, BPSys3 <dbl>, BPDia3 <dbl>, Testosterone <lgl>,
## #   DirectChol <dbl>, TotChol <dbl>, UrineVol1 <dbl>, UrineFlow1 <dbl>,
## #   UrineVol2 <dbl>, UrineFlow2 <dbl>, Diabetes <chr>, DiabetesAge <dbl>,
## #   HealthGen <chr>, DaysPhysHlthBad <dbl>, DaysMentHlthBad <dbl>,
## #   LittleInterest <chr>, Depressed <chr>, nPregnancies <dbl>, nBabies <dbl>,
## #   Age1stBaby <dbl>, SleepHrsNight <dbl>, SleepTrouble <chr>,
## #   PhysActive <chr>, PhysActiveDays <dbl>, TVHrsDay <lgl>, CompHrsDay <lgl>,
## #   TVHrsDayChild <dbl>, CompHrsDayChild <dbl>, Alcohol12PlusYr <chr>,
## #   AlcoholDay <dbl>, AlcoholYear <dbl>, SmokeNow <chr>, Smoke100 <chr>,
## #   Smoke100n <chr>, SmokeAge <dbl>, Marijuana <chr>, AgeFirstMarij <dbl>,
## #   RegularMarij <chr>, AgeRegMarij <dbl>, HardDrugs <chr>, SexEver <chr>,
## #   SexAge <dbl>, SexNumPartnLife <dbl>, SexNumPartYear <dbl>, SameSex <chr>,
## #   SexOrientation <chr>, PregnantNow <chr>

#knitr::kable(NHANES[c(1,4,5,6,7,8),c(1,3,4,7,17,20,21,25)],format = "markdown")
NHANES[c(1,4,5,6,7,8),c(1,3,4,7,17,20,21,25)] ## take a look at a subset of the dataset

## # A tibble: 6 x 8
##      ID Gender   Age Race1 Weight Height   BMI BPSysAve
##   <dbl> <chr>  <dbl> <chr>  <dbl>  <dbl> <dbl>    <dbl>
## 1 51624 male      34 White   87.4   165.  32.2      113
## 2 51625 male       4 Other   17     105.  15.3       NA
## 3 51630 female    49 White   86.7   168.  30.6      112
## 4 51638 male       9 White   29.8   133.  16.8       86
## 5 51646 male       8 White   35.2   131.  20.6      107
## 6 51647 female    45 White   75.7   167.  27.2      118

2.1 Take a glimpse() at your data

The glimpse function allows us to (obviously) take a first, informative glimpse at our data. The function is part of the dplyr, which we will explore in much more detail below!

dplyr::glimpse(NHANES[,1:10])

## Observations: 10,000
## Variables: 10
## $ ID            <dbl> 51624, 51624, 51624, 51625, 51630, 51638, 51646, 51647,…
## $ SurveyYr      <chr> "2009_10", "2009_10", "2009_10", "2009_10", "2009_10", …
## $ Gender        <chr> "male", "male", "male", "male", "female", "male", "male…
## $ Age           <dbl> 34, 34, 34, 4, 49, 9, 8, 45, 45, 45, 66, 58, 54, 10, 58…
## $ AgeDecade     <chr> "30-39", "30-39", "30-39", "0-9", "40-49", "0-9", "0-9"…
## $ AgeMonths     <dbl> 409, 409, 409, 49, 596, 115, 101, 541, 541, 541, 795, 7…
## $ Race1         <chr> "White", "White", "White", "Other", "White", "White", "…
## $ Race3         <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Education     <chr> "High School", "High School", "High School", NA, "Some …
## $ MaritalStatus <chr> "Married", "Married", "Married", NA, "LivePartner", NA,…

# or glimpse(NHANES) to see all the variables in the dataset

3 Data Tidying

library(tidyverse)

## ── Attaching packages ────────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.2.1     ✓ dplyr   0.8.3
## ✓ tibble  2.1.3     ✓ stringr 1.4.0
## ✓ tidyr   1.0.0     ✓ forcats 0.4.0
## ✓ purrr   0.3.3

## ── Conflicts ───────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

** Important ** If you are not familiar yet with the concepts of tidy data, have a look at the preliminary_tidyverse.Rmd file!

If we consider our NHANES dataframe, we see it is already in a tidy format, as;

• Each variable forms a column.
• Each observation forms a row.
• Each type of observational unit forms a table.

Each row contains all of the information on a single subject (US civilian) in the study.

In the next tutorial, we will work with a dataset on the effects of a certain drug, captopril, on the systolic and diastolic blood pressure of patients. This will not be a tidy dataset. As such, the details of tidying data with tidyverse will be described there

4 Data wrangling with dplyr

library(dplyr)

The dplyr package provides us with a large set of functions for handling our data. We refer to the preliminary_tidyverse.Rmd file for a more detailed description of these functions.

The most important dplyr functions to remember are:

dplyr verbs	Description
select()	select columns
filter()	filter rows
arrange()	re-order or arrange rows
mutate()	create new columns
summarize()	summarize values
group_by()	allows for group operations in the “split-apply-combine” concept

In addition, dplyr uses the pipe operator %>% to use the output of one function as an input to the next function. Instead of nesting functions (reading from the inside to the outside), the idea of of piping is to read the functions from left to right.

4.1 Select and filter

Here we will demonstrate the dplyr functionalities with a couple of small examples.

Let’s say we want to investigate how many subjects are 1. men that 2. are older than 18 years old, 3. are taller than 150 cm, and 4. have weight of less than 80kg

NHANES %>% 
  select(c("Gender","Age","Weight","Height")) %>% ## select the columns of interest
  filter(Gender == "male", 
         Age > 18,
         Height > 150,
         Weight < 80) %>% ## filter observations (rows) based on the required values
  nrow()

## [1] 1286

Filtering the dataset based on thee three filtering criteria retain 3.601 subjects. Note that the select step was not strictly necessary, but since we could answer the question based on only these 4 columns, there is no need to retain the the rest of the data (if this is the only question).

4.2 Arrange

Next question; Within the Race1 category of Hispanics, select men that are not married and display the first five by descending height.

Hint: use the desc() function inside of arrange() to order rows in a descending order. Use the base R function head to display the first five.

NHANES %>%
  select(c("Race1", "Gender","MaritalStatus","Height")) %>%
  filter(MaritalStatus != "Married", Race1 == "Hispanic", Gender == "male") %>% 
  arrange(desc(Height)) %>% 
  head(n=5)

## # A tibble: 5 x 4
##   Race1    Gender MaritalStatus Height
##   <chr>    <chr>  <chr>          <dbl>
## 1 Hispanic male   NeverMarried    197.
## 2 Hispanic male   NeverMarried    190.
## 3 Hispanic male   NeverMarried    182.
## 4 Hispanic male   NeverMarried    181.
## 5 Hispanic male   NeverMarried    180.

We have successfully combined three of the main dplyr functionalities. Let’s try to explore another one!

4.3 Mutate

Assume that we don’t trust the BMI values in our dataset and we decide to calculate them ourselves. BMI is typically calculated by taking a person’s weight (in kg) and dividing it by the its height (in m) squared. Based on this rule, generate a new column, BMI_self, for the subset of the NHANES dataset that we obtained from the previous question. To create new columns, we will use the mutate() function in dplyr.

NHANES %>%
  select(c("Race1", "Gender","MaritalStatus","Height","Weight","BMI")) %>%
  filter(MaritalStatus != "Married", Race1 == "Hispanic", Gender == "male") %>% 
  mutate(BMI_self = Weight/(Height/100)**2)

## # A tibble: 107 x 7
##    Race1    Gender MaritalStatus Height Weight   BMI BMI_self
##    <chr>    <chr>  <chr>          <dbl>  <dbl> <dbl>    <dbl>
##  1 Hispanic male   Divorced        177.  114.   36.3     36.3
##  2 Hispanic male   Divorced        177.  114.   36.3     36.3
##  3 Hispanic male   NeverMarried    190.  125.   34.7     34.7
##  4 Hispanic male   NeverMarried    176.  161.   52.1     52.1
##  5 Hispanic male   LivePartner     164.   88.2  32.8     32.8
##  6 Hispanic male   LivePartner     164.   88.2  32.8     32.8
##  7 Hispanic male   NeverMarried    162.   78.5  30.0     30.0
##  8 Hispanic male   Divorced        173.   95.4  31.8     31.8
##  9 Hispanic male   LivePartner     179.   92.8  28.9     28.9
## 10 Hispanic male   LivePartner     179.   92.8  28.9     28.9
## # … with 97 more rows

Note that now we need to include the Weight and BMI columns in the select statement, because we need that input for the mutate function. Good news: It turns out that the BMI column in the original dataset was computed correctly after all! But now we are interested in the mean value of the BMI_Self column. This requires a fifth dplyr function: summarize().

4.4 summarize

Given the filtering of the question above, compute the mean value of the BMI_self column.

NHANES %>%
  select(c("Race1", "Gender","MaritalStatus","Height","Weight","BMI")) %>%
  filter(MaritalStatus != "Married", Race1 == "Hispanic", Gender == "male") %>% 
  mutate(BMI_self = Weight/(Height/100)**2) %>%
  summarize(avg_BMI_self = mean(BMI_self, na.rm = TRUE))

## # A tibble: 1 x 1
##   avg_BMI_self
##          <dbl>
## 1         28.8

## the addition argument na.rm = TRUE makes sure to remove missing
## values for the purpose of calculating the mean

For this particular subset of the data, we find an average BMI value of 28.59 kg/m**2.

There are many other summary statistics you could consider such sd(), min(), median(), max(), sum(), n() (returns the length of vector), first() (returns first value in vector), last() (returns last value in vector) and n_distinct() (number of distinct values in vector). We will elaborate on these functions later.

Note that choosing the most informative summary statistic is very important! This can be shown with the following example;

** change location to github location ** ![/Users/jg/Desktop/PhD/Teaching/Lisbon_2020/psls20/data/Figure_partners.png){width=60%}

In this study, men and woman were asked about their “ideal number of partners desired over 30 years”. While almost all subjects desired a number between 0 and 50, three male subjects selected a number above 100. These three outliers in the data can have a large impact on the data analysis, especially when we work with summary statistics that are sensitive to these outliers.

When we look at the mean, for instance, we see that on average woman desire 2.8 partners, while men desire 64.3 partners on average, suggesting a large discrepancy between male and female desires.

However, if look at a more robust summary statistic such as the median, we see that the result is 1 for both men and woman. It is clear that the mean value was completely distorted by the three outliers in the data.

Another example of a more robust summary statistic is the geometric mean.

4.5 Group

We have already combined 5 very important functions. Here we will add a final one: group_by().

The group_by() verb is and incredibly powerful function in dplyr. It allows us, for example, to calculate summary statistics for different groups of observations.

If we take our example from above, let’s say we want to split the data frame by some variable (e.g. MaritalStatus), apply a function (mean) to a column (e.g. BMI_self) of the individual
data frames and then combine the output back into a summary data frame.

NHANES %>% 
  select(c("Race1", "Gender","MaritalStatus","Height","Weight","BMI")) %>%
  filter(MaritalStatus != "Married", Race1 == "Hispanic", Gender == "male") %>%
  mutate(BMI_self = Weight/(Height/100)**2) %>%
  group_by(MaritalStatus) %>% ##  group the subjects by their marital status
  summarize(avg_BMI_self = mean(BMI_self, na.rm = TRUE)) ## calculate the mean BMI_self value of each group

## # A tibble: 5 x 2
##   MaritalStatus avg_BMI_self
##   <chr>                <dbl>
## 1 Divorced              28.4
## 2 LivePartner           29.0
## 3 NeverMarried          28.1
## 4 Separated             31.4
## 5 Widowed               30.6

We have successfully combined six of the most important dplyr functionalities!

Now that we have all the required functions for importing, tidying and wrangling data in place, we will learn how to visualize our data with the ggplot2 package.

5 Data Visualization

library(ggplot2)

As you might have have already seen, there are many functions available in base R that can create plots (e.g. plot(), boxplot()). Others include: hist(), qqplot(), etc. These functions are great because they come with a basic installation of R and can be quite powerful when you need a quick visualization of something when you are exploring data.

We are choosing to introduce ggplot2 because, in our opinion, it’s one of the most simple ways for beginners to create relatively complicated plots that are intuitive and aesthetically pleasing.

5.1 Univariate statistics

In univariate statistics, we focus on a single variable of interest. Here, we will show different ways to visualize the BMI variable from the NHANES dataset. Importantly, different types of visualizations will provide us with different types of information!

Here we will visualize the same data using a dotplot, a histogram and a boxplot.

5.1.1 Dotplot

NHANES %>%
  filter(!is.na(BMI)) %>%
  head(NHANES, n=100) %>% ## to make the visualization more clear, we take only the first 100 subjects
  ggplot(aes(x=BMI)) +
  geom_dotplot(method = "histodot") 

## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

In stead of making a visualization “on the fly” as we did in the code chunk above, we may also store the visualization in an object. This allows us to: - Call upon the object again later - Add extra layers to the plot.

## store the plot from above in the "p" object
p <- NHANES %>%
  filter(!is.na(BMI)) %>%
  head(NHANES, n=100) %>%
  ggplot(aes(x=BMI)) +
  geom_dotplot(method = "histodot") 

## Add additional info to the plot
p <- p + xlab("BMI (kg/m2)") + 
  theme_bw()

## Display the plot
p

## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

Note that if we are sure we do not want to store the plot for later retrieval, we shouldn’t: it would simple be a waste of storage space!

In the dotplot, the BMI value for each subject (in this case, the first 100 subjects of the NHANES study), is plotted on the x-axis. The y-axis allows us to count how many subjects have a BMI value that falls within a certain counting bin. For instance, in the BMI interval [29.5, 30.[ we observe 4 subjects.

Note, however, that the dotplot does not provide us with information on the mean/median values of the data, which are typically features of high interest.

5.1.2 Histogram

NHANES %>%
  filter(!is.na(BMI)) %>%
  head(NHANES, n=100) %>%
  ggplot(aes(x=BMI)) +
  geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

The information that can be obtained from this histogram is similar to that of the dotplot. Here, we can read immediately (i.e. without counting) that the BMI interval [29.5, 30.0[ contains 4 subjects. Again, the histogram does not provide us with information on the mean/median values of the data, but only on its distribution.

5.1.3 Boxplot

set.seed(2) ## to make the horizontal position of the jitter non-random

give.n <- function(x){
  return(c(y = 50, label = length(x))) 
}

NHANES %>%
  filter(!is.na(BMI)) %>%
  head(NHANES, n=100) %>%
  ggplot(aes(x="", y=BMI)) +
  geom_boxplot(outlier.shape=NA) + 
  geom_jitter(width=0.2) +
  stat_summary(geom="text", fun.y=quantile,
               aes(label=sprintf("%1.1f", ..y..)),
               position=position_nudge(x=0.5), size=4.5) +
  annotate("text", x = c(1.5,1.5,1.5,1.5,1.5,1.05,1.05), y = c(12.5,19,24.5,28.5,45,18,35), label = c("(minimum)","(25%)","(median)","(75%)","(maximum)","whisker","whisker"),size=3) +
  stat_summary(fun.data = give.n, geom = "text") 

Arguably, the boxplot is the most informative default visualization strategy. First, it provides us with similar insights to the shape of the distribution as the histogram. Second, It clearly displays several useful summary statistics such as the median, the interquartile range, whiskers and outliers. Third, with the geom_jitter functionality we can also project each individual value of the dataset on the plot.

The only disadvantage of the boxplot as compared to the dotplot or histogram is that we cannot see (from the figure) how many subjects have a BMI value within a certain interval. However, this feature usually is not of primary interest.

5.2 Bivariate statistics

In bivariate statistics, the goal is to study two variables, including the relationship between both variables.

In terms of visualizations, the scatterplot is the baseline method for displaying two variables.

5.2.1 Create scatter plots using geom_point()

For the NHANES dataset, we can for instance look at the relationship between a person’s height and weight values.

NHANES %>%
  ggplot(aes(x = Height, y = Weight)) +
  geom_point() + 
  xlab("Height (cm)") + 
  ylab("Weight (kg)")

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

We used the xlab() and ylab() functions in ggplot2 to specify the x-axis and y-axis labels. Note that NA values were automatically remove by the geom_point function.

ggplot2 also has a very broad panel of aesthetic features for improving your plot. One very basic feature is that we can give colors to the ggplot object. For instance, we can give different colors to the dots in the previous scatterplot based on a subject’s gender.

NHANES %>%
  ggplot(aes(x = Height, y = Weight, color = Gender)) +
  geom_point() + 
  xlab("Height (cm)") + 
  ylab("Weight (kg)")

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

5.3 Combining dplyr and ggplot2

Note that the previous functions from the dplyr package can be easily combined with ggplot through the concept of pipes %>%.

For instance, we could make the same scatterplot as above, but only for white, married adults.

In addition, we here also show some convenient ggplot features; 1. Setting the colors manually 2. Set to have a white background 3. Set a (main) title for the plot 4. Pick a different shape for the dots in the scatterplot 5. Manually set the limits of the x- and y-axes

Note that this a only the tip of the iceberg of the the ggplot functionalities!

NHANES %>%
  filter(Age >= 18, Race1 == "White", MaritalStatus == "Married") %>% ## select the required data
  ggplot(aes(x = Height, y = Weight, color = Gender)) +
  geom_point(shape=17, size = 1) + ## set to a different shape (triangle) and size
  ggtitle("Height versus Weight") + ## for the main title
  xlab("Height (cm)") + 
  ylab("Weight (kg)") +
  xlim(0,220) + ## set limit of x-axis
  ylim(0,220) + ## set limit of y-axis
  scale_color_manual(values = c("red","blue")) + ## manually set colors
  theme_bw() ## set white background

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

Next to scatterplots, we have a large number of other types of plots at our disposal. Importantly, some plots will be more informative than others, depending on the research question. Therefore, choosing the plot that is most informative is crucial. We show this with a more elaborate example later (Data_exploration_captopril.Rmd).

5.4 Final example

5.4.1 Goal

Set up a reference interval for the systolic blood pressure in the NHANES dataset.

5.4.2 Background

The captopril dataset, which we will explore and analyse later, holds information on 15 patients that have increased blood pressure values.

Before we may conduct such an experiment, we first need to know which values should be considered increased and which ones should be considered normal. To find an interval for values that are normal, we can set up a reference interval. To set up this interval, we will use the NHANES dataset. The BPSysAve column holds data on the systolic blood pressure. To select healthy subjects, we will need to subset the data.

5.4.3 Analysis

First, we will plot the data for all subjects for which we have all the required data and that are between 40 and 65 years old.

Note, that we calculate the standard deviation on the average blood pressure measurements. If you have the same number of “technical repeats” on the same person you are allowed to average them first. Note, however, that if there are some technical measurements missing, averaging is no longer allowed because then the averaged values will be heteroscedastic: some averages are based on more technical repeats than others so they will have a different standard deviation.

## histogram of BPSysAve for subjects wit age between 40 and 65 years old
NHANES %>%
  filter(!is.na(Race1), !is.na(Smoke100n), !is.na(BMI_WHO), !is.na(Age), !is.na(HardDrugs), !is.na(HealthGen), !is.na(Gender), !is.na(AlcoholYear), !is.na(BPSys1), !is.na(BPSys2), !is.na(BPSys3), !is.na(SleepTrouble)) %>% ## retains 4660 subjects with the required data
  filter(between(Age,40,65)) %>% ## filter the subjects on age, retains 2522 subjects
  distinct(ID,.keep_all=TRUE) %>% ## removes duplicated IDs, retains 1467 subjects
  ggplot(aes(x=BPSysAve)) +
  geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

An important requirement of calculating reference intervals is that the data are normally distributed. this is clearly not the case; the data has a long right tail, which is quite common in biological data (in this easier to have extreme values on the right-hand side than on the left-hand size, which is additionally bounded by zero for most biological variables).

By selecting only healthy subjects, we expect the data to be distributed more normally. We define healthy as being a non-smoker, without a history of diabetes, hard drugs or sleeping trouble, with a general health that is not considered poor and that has a BMI between 18.5 and 29.9.

## histogram of BPSysAve for HEALTHY male subjects with age between 40 and 65 years old
NHANES %>%
  filter(!is.na(Race1), !is.na(Smoke100n), !is.na(BMI_WHO), !is.na(Age), !is.na(HardDrugs), !is.na(HealthGen), !is.na(Gender), !is.na(AlcoholYear), !is.na(BPSys1), !is.na(BPSys2), !is.na(BPSys3), !is.na(SleepTrouble)) %>% ## retains 4660 subjects with the required data
  filter(between(Age,40,65)) %>% ## filter the subjects on age, retains 2522 subjects
  distinct(ID,.keep_all=TRUE) %>% ## removes duplicated IDs, retains 1467 subjects
  filter(Smoke100n == "Non-Smoker", Diabetes == "No", HardDrugs == "No", HealthGen != "Poor", SleepTrouble == "No", BMI_WHO%in%c("18.5_to_24.9","25.0_to_29.9")) %>% ## filter to have only healthy subjects, based on multiple health criteria. Retains 275 subjects.
  ggplot(aes(x=BPSysAve)) +
  geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

We will evaluate if the distribution is normally distributed using a QQ-plot. In this plot we compare the quantiles in the data with the quantiles of the normal distribution. If the data are normally distributed then the data in the QQ-plot approximately follows a straight line. Deviations from the straight line indicate deviations of normality.

NHANES %>%
  filter(!is.na(Race1), !is.na(Smoke100n), !is.na(BMI_WHO), !is.na(Age), !is.na(HardDrugs), !is.na(HealthGen), !is.na(Gender), !is.na(AlcoholYear), !is.na(BPSys1)&!is.na(BPSys2)&!is.na(BPSys3), !is.na(SleepTrouble)) %>%
  filter(between(Age,40,65)) %>% 
  distinct(ID,.keep_all=TRUE) %>% 
  filter(Smoke100n == "Non-Smoker", Diabetes == "No", HardDrugs == "No", HealthGen != "Poor", SleepTrouble == "No", BMI_WHO%in%c("18.5_to_24.9","25.0_to_29.9")) %>% 
  ggplot(aes(sample=BPSysAve))+
  geom_qq()+
  geom_qq_line()

We observe that the distribution is slightly skewed too the right. We will therefore log transform the bloodpressure measurements. We use the log2 transformation: a difference of 1 on the log scale corresponds to doubling the blood pressure on the original scale: \(log2(B)-log2(A)=log2(B/A)\) \(2^{1}=2\).

NHANES %>%
  filter(!is.na(Race1), !is.na(Smoke100n), !is.na(BMI_WHO), !is.na(Age), !is.na(HardDrugs), !is.na(HealthGen), !is.na(Gender), !is.na(AlcoholYear), !is.na(BPSys1)&!is.na(BPSys2)&!is.na(BPSys3), !is.na(SleepTrouble)) %>%
  filter(between(Age,40,65)) %>% 
  distinct(ID,.keep_all=TRUE) %>% 
  filter(Smoke100n == "Non-Smoker", Diabetes == "No", HardDrugs == "No", HealthGen != "Poor", SleepTrouble == "No", BMI_WHO%in%c("18.5_to_24.9","25.0_to_29.9")) %>% 
  ggplot(aes(sample=BPSysAve %>% log2))+
  geom_qq()+
  geom_qq_line()

We see that the data is normally distributed upon log transformation. We can now use the normal distribution to construct a 95% interval at log scale that we will backtransform to the original scale.

## to get the 95% reference interval for the healthy group;

summary_NHANES <- NHANES %>%
  filter(!is.na(Race1), !is.na(Smoke100n), !is.na(BMI_WHO), !is.na(Age), !is.na(HardDrugs), !is.na(HealthGen), !is.na(Gender), !is.na(AlcoholYear), !is.na(BPSys1), !is.na(BPSys2), !is.na(BPSys3), !is.na(SleepTrouble)) %>% ## retains 4660 subjects with the required data
  filter(between(Age,40,65)) %>% ## filter the subjects on age, retains 2522 subjects
  distinct(ID,.keep_all=TRUE) %>% ## removes duplicated IDs, retains 1467 subjects
  filter(Smoke100n == "Non-Smoker", Diabetes == "No", HardDrugs == "No", HealthGen != "Poor", SleepTrouble == "No", BMI_WHO%in%c("18.5_to_24.9","25.0_to_29.9")) %>% 
  mutate(BPSysAveLog=BPSysAve %>% log2) %>%
  summarize_at("BPSysAveLog",
               list(mean=~mean(.,na.rm=TRUE),
                    sd=~sd(.,na.rm=TRUE),
                    n=function(x) x%>%is.na%>%`!`%>%sum)) %>%
  mutate(se=sd/sqrt(n))

paste("Mean value (log2):", summary_NHANES$mean)

## [1] "Mean value (log2): 6.88522582539534"

paste("Geometric mean value (log2):", 2^summary_NHANES$mean)

## [1] "Geometric mean value (log2): 118.211437989542"

paste("Standard deviation (log2):", summary_NHANES$sd)

## [1] "Standard deviation (log2): 0.168454138831745"

paste0("Reference interval (log2): [", summary_NHANES$mean - 2*summary_NHANES$sd, ";", summary_NHANES$mean + 2*summary_NHANES$sd, "]")

## [1] "Reference interval (log2): [6.54831754773185;7.22213410305883]"

paste0("Reference interval (log2): [", 2^(summary_NHANES$mean - 2*summary_NHANES$sd), ";", 2^(summary_NHANES$mean + 2*summary_NHANES$sd), "]")

## [1] "Reference interval (log2): [93.5922747749331;149.306597207509]"

The log2 systolic blood pressure for healthy subjects is approximately normally distributed. We have calculated the mean and standard deviation of blood pressure values in this healthy subset at the log2 scale, which allows us to set up the (95%) reference interval, for what we can consider to be normal blood pressure values. Any of the later patients who’s values do not fall within this interval, we can consider abnormal.

The log2 mean blood pressure value of the healthy subset is 6.89, with a standard deviation of 0.17. When we replace the population average and population standard deviation with these estimates, we obtain a 95% reference interval of [6.55;7.22] on log scale and [93.59,149.31] mmHg on the original scale.

Note, that in the literature a value 140 mmHg for the systolic blood pressure is typically considered to be the upper limit of normality.

---

Home