MATH2349 Data Wrangling

Assignment 2

Required packages

The following packages were called:

require(rmarkdown)
require(readr)
require(readxl)
require(tidyr)
require(dplyr)
require(knitr)
require(stringr)
require(foreign)
require(ggplot2)
require(lubridate)
require(outliers)
require(forecast)

Executive Summary

Three data sets have been selected for this assignment:

• Two meet Hadley Wickham’s “Tidy Data” principles [1] and will be merged to create a larger sample of individual observations; and
• One data set is messy and will be imported to create a data frame that meets the “Tidy Data” principles.

Preprocessing steps include:

• Checking data frame dimension and structure using dim() and str();
• Converting variable types (e.g. numeric to factor with ordered levels);
• Creating and converting variables to standardise data formats for merging;
  
• Sub-setting data frames using select() and filter();
  
• Renaming variables using rename(); and
  
• Merging data sets using rbind() and cbind().

Data sets were then scanned for missing values, special characters and other obvious errors, and data visualisation and z-score tests were used to identify, evaluate and manage outliers.

Within the merged data set, two variables were transformed to:

• Convert a non-linear relationship to a linear one; and
  
• Decrease a left-skew in an attempt to convert the sub-set to a normal distribution.

Data Description

Data set 1:

Sourced from Heart Disease UCI on Kaggle, the data contains 303 observations of 14 variables. “heart.csv” was downloaded as a ZIP file, uncompressed, and the CSV file saved locally, prior to importing as heart1 using read_csv().

Variables in the data set:

• age: Age (in years)
  
• sex: Sex (1=Male, 0=Female)
  
• cp: Classification of chest pain symptoms (4 categorical levels)
  
• trestbps: Resting blood pressure (mm Hg)
  
• chol: Total cholesterol level (mg/dL)
  
• fbs: Fasting blood sugar (if >120 mg/dL, 1=true, 0=false)
  
• restecg: Resting electrocardiographic results (3 categorical levels)
  
• thalach: Maximum heart rate achieved (bpm)
  
• exang: Presence of exercise induced angina in the person (1=yes, 0=no)
  
• oldpeak: Indicates whether ST depression was induced by exercise, relative to rest
  
• slope: Slope of the peak exercise ST segment (1=sloping up, 2=flat, 3=sloping down)
  
• ca: Number of major vessels (from 0 to 3) coloured by fluorosopy
  
• thal: Presence of thalassemia, a blood disorder (1=normal, 2=fixed defect, 3=reversable defect)
  
• target: Indicates whether the person has been diagnosed with heart disease (0=no, 1=yes)

heart1 <- read_csv("~/RStudio Projects/RStudio Data/heart.csv")

Parsed with column specification:
cols(
  age = col_double(),
  sex = col_double(),
  cp = col_double(),
  trestbps = col_double(),
  chol = col_double(),
  fbs = col_double(),
  restecg = col_double(),
  thalach = col_double(),
  exang = col_double(),
  oldpeak = col_double(),
  slope = col_double(),
  ca = col_double(),
  thal = col_double(),
  target = col_double()
)

Data Description cont.

Data set 2:

Sourced from Cardiovascular Disease on Kaggle. The data contains 70,000 observations of 13 variables. “cardio_train.csv” was downloaded and saved locally, prior to importing as heart2 using read_delim(). Variable “id” is not needed and was removed using a subset command.

Variables in the data set:

• id: Identification number of observation
  
• age: Age (in days)
  
• gender: Sex (1=women, 2=men)
  
• height: Height (cm)
  
• Weight: Weight (kg)
  
• ap_hi: Systolic blood pressure
  
• ap_lo: Diastolic blood pressure
  
• cholesterol: Cholesterol group classification (1=normal, 2=above normal, 3=well above normal)
  
• gluc: Glucose level (1=normal, 2=above normal, 3=well above normal)
  
• smoke: Whether subject smokes (0=no, 1=yes)
  
• alco: Whether subject consumes alcohol (0=no, 1=yes)
  
• active: Whether subject engages in physical activity (0=no, 1=yes)
  
• cardio: Presence or absence of cardiovascular disease (0=no, 1=yes)

heart2 <- read_delim("~/RStudio Projects/RStudio Data/cardio_train.csv", ";",
                     escape_double = FALSE, col_names = TRUE, trim_ws = TRUE)

Parsed with column specification:
cols(
  id = col_double(),
  age = col_double(),
  gender = col_double(),
  height = col_double(),
  weight = col_double(),
  ap_hi = col_double(),
  ap_lo = col_double(),
  cholesterol = col_double(),
  gluc = col_double(),
  smoke = col_double(),
  alco = col_double(),
  active = col_double(),
  cardio = col_double()
)

heart2 <- heart2 %>% subset(select = 2:13)

Data Description cont.

Data set 3:

Sourced from the ABS, the Australian Health Survey: Biomedical Results for Chronic Diseases, 2011–12 data file was downloaded and saved locally. Variables were imported individually and will be remerged as heart3 under Tidy and Manipulate below.

Selected variables for import from the data set:

• Total cholesterol: Cholesterol group based on reading (mmol/L)
  
• Males: Number of males per cholesterol group (’000)
  
• Females: Number of females per cholesterol group (’000)
  
• Persons: Number of persons per cholesterol group (’000)

abs_file <- "C:/Users/Wendy/Documents/RStudio Projects/Rstudio Data/43640do004_20112012.xls"

group <- as.data.frame(read_excel(path = abs_file, sheet = "Table_4_1",
                                  range = "A8:A15", col_names = "chol_group"))

males <- as.data.frame(read_excel(path = abs_file, sheet = "Table_4_1",
                                  range = "B8:B15", col_names = "count"))

females <- as.data.frame(read_excel(path = abs_file, sheet = "Table_4_1",
                                    range = "C8:C15", col_names = "count"))

persons <- as.data.frame(read_excel(path = abs_file, sheet = "Table_4_1",
                                    range = "D8:D15", col_names = "count"))

Understanding the Data

After importing, and before preprocessing the data, the dim() and str() functions are used to summarise the variables, understand the size, structures, and data types.

It is apparent that both data sets have imported with all variables classified as numeric. This is not surprising due to the use of numbers to represent ordinal variables.

heart1 %>% dim()

[1] 303  14

heart1 %>% str(give.attr = FALSE)

tibble [303 x 14] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ age     : num [1:303] 63 37 41 56 57 57 56 44 52 57 ...
 $ sex     : num [1:303] 1 1 0 1 0 1 0 1 1 1 ...
 $ cp      : num [1:303] 3 2 1 1 0 0 1 1 2 2 ...
 $ trestbps: num [1:303] 145 130 130 120 120 140 140 120 172 150 ...
 $ chol    : num [1:303] 233 250 204 236 354 192 294 263 199 168 ...
 $ fbs     : num [1:303] 1 0 0 0 0 0 0 0 1 0 ...
 $ restecg : num [1:303] 0 1 0 1 1 1 0 1 1 1 ...
 $ thalach : num [1:303] 150 187 172 178 163 148 153 173 162 174 ...
 $ exang   : num [1:303] 0 0 0 0 1 0 0 0 0 0 ...
 $ oldpeak : num [1:303] 2.3 3.5 1.4 0.8 0.6 0.4 1.3 0 0.5 1.6 ...
 $ slope   : num [1:303] 0 0 2 2 2 1 1 2 2 2 ...
 $ ca      : num [1:303] 0 0 0 0 0 0 0 0 0 0 ...
 $ thal    : num [1:303] 1 2 2 2 2 1 2 3 3 2 ...
 $ target  : num [1:303] 1 1 1 1 1 1 1 1 1 1 ...

heart1 %>% head()

heart2 %>% dim()

[1] 70000    12

heart2 %>% str()

tibble [70,000 x 12] (S3: tbl_df/tbl/data.frame)
 $ age        : num [1:70000] 18393 20228 18857 17623 17474 ...
 $ gender     : num [1:70000] 2 1 1 2 1 1 1 2 1 1 ...
 $ height     : num [1:70000] 168 156 165 169 156 151 157 178 158 164 ...
 $ weight     : num [1:70000] 62 85 64 82 56 67 93 95 71 68 ...
 $ ap_hi      : num [1:70000] 110 140 130 150 100 120 130 130 110 110 ...
 $ ap_lo      : num [1:70000] 80 90 70 100 60 80 80 90 70 60 ...
 $ cholesterol: num [1:70000] 1 3 3 1 1 2 3 3 1 1 ...
 $ gluc       : num [1:70000] 1 1 1 1 1 2 1 3 1 1 ...
 $ smoke      : num [1:70000] 0 0 0 0 0 0 0 0 0 0 ...
 $ alco       : num [1:70000] 0 0 0 0 0 0 0 0 0 0 ...
 $ active     : num [1:70000] 1 1 0 1 0 0 1 1 1 0 ...
 $ cardio     : num [1:70000] 0 1 1 1 0 0 0 1 0 0 ...

heart2 %>% head()

Understanding the Data cont.

Data conversions are needed before merging heart1 and heart2 in order to standardise data formats. The following actions have been taken:

• Convert heart1 ‘sex’ and heart2 ‘gender’ to a common variable ‘sex’ and values converted from numeric to factor with labels ‘male’ and ‘female’

heart1$sex <- heart1$sex %>% factor(levels = c("0", "1"), labels = c("female", "male"))

heart2 <- heart2 %>% dplyr::rename(sex = gender)
heart2$sex <- heart2$sex %>% factor(levels = c("1", "2"), labels = c("female", "male"))

• Convert variable heart2 ‘cholesterol’ from numeric to factor with ordered levels by changing cholesterol levels 1, 2, 3 to normal, above normal, and well above normal, aligned with Total Cholesterol levels = Desirable, Borderline High, and High (respectively) reproduced in the table below, based on a WebMD table [2]:

Total Cholesterol	Category	New values
Less than 200	Desirable	Normal
200 - 239	Borderline High	Above normal
240 and above	High	Well above normal

heart2$cholesterol <- heart2$cholesterol %>% factor(levels = c("1", "2", "3"),
         labels = c("normal", "above normal", "well above normal"), ordered = TRUE)

• Rename variable heart1 ‘target’ to ‘hdisease’ and change values for easier interpretation

heart1 <- heart1 %>% dplyr::rename(hdisease = target)

heart1$hdisease <- heart1$hdisease %>% factor(levels = c("0", "1"), labels = c("no", "yes"))

• Convert 5 additional variables from numeric to factors.

heart1$fbs <- heart1$fbs %>% as.factor()
heart1$cp <- heart1$cp %>% as.factor()
heart1$exang <- heart1$exang %>% as.factor()
heart1$slope <- heart1$slope %>% as.factor()
heart1$thal <- heart1$thal %>% as.factor()

Understanding the Data cont.

The data is reviewed to confirm the changes have been made:

heart1 %>% head()

heart2 %>% head()

Tidy and Manipulate Data

Finding an untidy data set was challenging. To meet the assignment criteria a third data set was selected and which is ‘untidy’. It was not scraped directly from the web. The data set does not comply with Hadley Wickham’s Tidy Data principle that “(e)ach variable forms a column”[1]. Specifically, it contains variables in both rows (cholesterol groups) and columns (males, females, persons). Steps taken to prepare this data set:

• New variables were created (country and sex) and appended to data frames

group$country <- "aus"
males$sex <- "male"
females$sex <- "female"
persons$sex <- "all"

• Data frame group merged with males, females and persons data frames using cbind() to create three dfs of dimension (8 x 4)

males <- cbind(group, males)
females <- cbind(group, females)
persons <- cbind(group, persons)

• The resulting 3 data frames males, females and persons are merged using rbind() to create heart3

heart3 <- rbind(males, females, persons)

• The variable ‘chol_group’ is converted to an ordered factor with 8 levels

heart3$chol_group <- heart3$chol_group %>% factor(ordered = TRUE)

Review the changes using the str() function:

heart3 %>% str()

'data.frame':   24 obs. of  4 variables:
 $ chol_group: Ord.factor w/ 8 levels "<4.0"<"≥4.0 to <4.5"<..: 1 2 3 4 5 6 7 8 1 2 ...
 $ country   : chr  "aus" "aus" "aus" "aus" ...
 $ count     : num  1193 1198 1585 1608 1189 ...
 $ sex       : chr  "male" "male" "male" "male" ...

Tidy and Manipulate Data cont.

To merge data from two countries using different cholesterol scores, it is necessary to create a new variable for chol_group to convert the score ranges for Australia (mmol/L) to the US score ranges (mg/dL). The mg/dL thresholds (200, 239, 240 mg/dL) are converted to mmol/L to determine where they fall within the ABS ranges. This is done with the following conversion formula [6]: \(chol (mmol/L) = \frac{chol (mg/dL)}{38.67}\)

norm_au <- 200 / 38.67  # Determine mmol/L value for max of 'normal'
abnorm_au <- 239 / 38.67  # Determine mmol/L value for max of 'above normal'
wanorm_au <- 240 /38.67  # Determine mmol/L value for min of 'well above normal'
norm_au

[1] 5.171968

abnorm_au

[1] 6.180502

wanorm_au

[1] 6.206362

It is clear that the ABS cholesterol groupings cannot be accurately mapped to cholesterol level groups in the WebMD table above as summarised below:

• ABS <4.0 = normal
• ABS =4.0 < 4.5 = normal
• ABS =4.5 < 5.0 = normal
• ABS =5.0 < 5.5 cannot be categorised, 200mg/dL = 5.17 mmol/L
• ABS =5.5 < 6.0 = above normal
• ABS =6.0 < 6.5 cannot be categorised, 240mg/dL = 6.21 mmol/L
• ABS =6.5 < 7.0 = well above normal
• ABS > 7.0 = well above normal

As a result, the decision is made not to approximate, therefore heart3 will not be preprocessed further.

Tidy and Manipulate Data cont.

Returning to the heart1 and heart2 data sets:

• Create a new factor variable heart1 ‘chol_group’ with ordered levels using cut() function

heart1$chol_group <- heart1$chol %>% 
  cut(breaks = c(0, 200, 239, 600), labels = c("normal","above normal","well above normal"),
      ordered = TRUE)
heart1 %>% head()

• Create a new variable heart2 ‘age_yr’ by calculating \({{age\space (days)} \over {365.25}} = {age\space (years)}\)

heart2$age_yr <- (heart2$age / 365.25) %>% round(0)
heart2 %>% head()

Tidy and Manipulate Data cont.

It is time to create subsets of heart1 and heart2, standardise column names. Verify the subsets are ready to be merged by reviewing the head(). The key variables are:

• Age (in years)
  
• Sex
  
• Cholesterol (by group)

heart1_sub <- heart1 %>% select(c(age, sex, chol_group))

heart2_sub <- heart2 %>% select(c(age_yr, sex, cholesterol))
heart2_sub <- heart2_sub %>% dplyr::rename(age = age_yr)
heart2_sub <- heart2_sub %>% dplyr::rename(chol_group = cholesterol)

heart1_sub %>% head()

heart2_sub %>% head()

Tidy and Manipulate Data cont.

• Merge the subsets using the rbind() function and retain the variable names of heart1

heart_merge <- rbind(heart1_sub, heart2_sub)

• Confirm the changes using dim() and str() functions

heart_merge %>% dim()

[1] 70303     3

heart_merge %>% str()

tibble [70,303 x 3] (S3: tbl_df/tbl/data.frame)
 $ age       : num [1:70303] 63 37 41 56 57 57 56 44 52 57 ...
 $ sex       : Factor w/ 2 levels "female","male": 2 2 1 2 1 2 1 2 2 2 ...
 $ chol_group: Ord.factor w/ 3 levels "normal"<"above normal"<..: 2 3 2 2 3 1 3 3 1 1 ...

Scan the Data

The data is checked for issues:

• The colSums() function was used to check for NA or missing values across all columns in the data frames;

colSums(is.na(heart_merge))

       age        sex chol_group 
         0          0          0 

• The summary() function was used to check for non-numeric values or special characters in numeric variable fields

heart_merge %>% summary()

      age            sex                    chol_group   
 Min.   :29.00   female:45626   normal           :52436  
 1st Qu.:48.00   male  :24677   above normal     : 9646  
 Median :54.00                  well above normal: 8221  
 Mean   :53.31                                           
 3rd Qu.:58.00                                           
 Max.   :77.00                                           

Scan the Data cont.

• To check for outliers, Tukey’s method [7] is used via a boxplot to visualise the data.

heart_merge %>% ggplot(mapping = aes(x = chol_group, y = age)) + 
  geom_boxplot(col = "black", outlier.colour = "red") +
  labs(title = "Boxplot of Age Ranges by Cholesterol Groups (heart_merge)", x = "Cholesterol Group", y = "Age")

The boxplot indicates possible outliers (age per cholesterol group) which may be skewing the results.

Assess outliers by calculating observations with an absolute z-score greater than 3. This method is performed with reference to MATH2349 Module 6 Notes.[7]

z.scores <- heart_merge$age %>% scores(type = "z")
z.scores %>% summary()

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-3.5894 -0.7838  0.1022  0.0000  0.6929  3.4985 

which(abs(z.scores) > 3)

[1]    73   130   145   239  6523 22647 30970 56209

Scan the Data cont.

Eight observations (out of 70,303) were identified as outliers with a Z-score > 3, so these are removed and a new data frame heart_clean is created for further analysis.

heart_clean <- heart_merge[ - which( abs(z.scores) >3 ), ]
heart_clean %>% dim()

[1] 70295     3

heart_clean %>% head()

Rerun the boxplot, after removing outliers with a z-score > 3. It is apparent that the ‘well above normal’ group still has observations which are outside the expected range, and relates to people aged less than 40 with ‘well above normal’ cholesterol levels.

heart_clean %>% ggplot(mapping = aes(x = chol_group, y = age)) +
  geom_boxplot(col = "black", outlier.colour = "red") +
  labs(title = "Boxplot of Age Ranges by Cholesterol Groups (heart_clean)",
       x = "Cholesterol Group", y = "Age")

Transform the Data

The first data transformation occurred in Tidy and Manipulate to create a new variable age_yr from heart2 age (days) through a mathematical calculation. This needed to occur prior to merging the data sets to standardise the data (i.e. to convert a non-linear relation into a linear one).

Noting that heart_clean chol_group “well above normal” still has outliers, a subset heart_wan is created and visualised with a histogram. To make the histogram clearer, an outline/fill has been added to clearly see each column [3]. The histogram shows a left-skew.

heart_wan <- heart_clean %>% filter(chol_group == "well above normal")
heart_wan %>% dim()

[1] 8219    3

ggplot(heart_wan) + geom_histogram(aes(x = age), bins = 15, color="black", fill="white") +
  labs(main = "Age Distribution for Cholesterol Group (well above normal)", x = "Age", y = "Count")

To reduce the skew, the Box-Cox transformation is applied to the subset and visualised by histogram using the function hist(). Note: The transformation would be applied to the complete data set, this is for demonstration purposes only. This method is performed with reference to MATH2349 Module 7 Notes. To reduce the number of \(\lambda\) values printed, the ‘max.print’ argument is added to limit to <100 values [4], and the histogram is manually adjusted [5].

options('max.print' = 60)
getOption('max.print')

[1] 60

boxcox_wan <- BoxCox(heart_wan$age, lambda = "auto")
boxcox_wan

 [1]  683.8387 1623.5640 1567.0813  967.2592 1151.2058 1199.6916 1681.0464 1681.0464
 [9]  923.7716 1859.4911 2519.2817 2111.4121 1457.1147 1457.1147 2111.4121 2111.4121
[17] 1299.6619 1151.2058 1403.6308  759.8178 1351.1466 1103.7197 1299.6619 2176.8913
[25] 1983.4524 1011.7464 1623.5640 2519.2817 1457.1147 1457.1147 1299.6619 1511.5982
[33] 1299.6619 1681.0464  839.7954 1011.7464 1799.0099 1351.1466  881.2837 2243.3702
[41] 2310.8486 1681.0464 1151.2058 1457.1147 1920.9720  923.7716 1403.6308  881.2837
[49] 1739.5284  881.2837 1249.1769 1249.1769 2046.9324 2046.9324 1511.5982 1057.2332
[57] 2046.9324 1739.5284  839.7954 1457.1147
 [ reached getOption("max.print") -- omitted 8159 entries ]
attr(,"lambda")
[1] 1.999924

mean_boxcox <- boxcox_wan %>% mean()
median_boxcox <- boxcox_wan %>% median()

hist_boxcox <- hist(boxcox_wan, breaks = seq(500, 3000, 100), col = "lightgrey",
                    border = "black", xlim = c(500, 3000), ylim = c(0, 1200),
                    main = "Histogram of 'well above normal' Cholesterol Group (Box-Cox)",
                    xlab = "Box-Cox Transformed Values") %>%
  abline(v=c(mean_boxcox, median_boxcox), lty=c(2,3), lwd=2) %>%
  legend(x = "topright", legend = c("Mean ","Median"), lty=c(2,3), lwd=2)

The resulting histogram indicates a left-skew remains in the normalised data. If this were a real study, additional tests would be required (for example, to test for homogeneity of variance) before conducting further statistical analyses.

For this assignment, this concludes the study.

References

[1] Wickham, H., “Tidy Data” (p4), Journal of Statistical Software, Volume 59, Issue 10, August 2014

[2] “Heart Disease and Lowering Cholesterol” (p2), WebMD

[3] “Adding Space between my geom_histogram bars-not barplot” stackoverflow

[4] “Increasing-max-print-does-not-always-work-in-r-and-neither-does-global-options-f” stackoverflow

[5] Davies, Tilman M., 2016 The Book of R: A First Course in Programming and Statistics (Chapter 14, pp294-296), No Starch Press

[6] Lipid Conversion Factors, National Center for Biotechnology Information

[7] Dr Anil Dolgun, 2020 MATH2349 Module Notes

Long, J.D., Teetor, P., 2019 R Cookbook, O’Reilly Media, Inc. 

Xie, Y., Dervieux, C., Riederer, E. 2020 R Markdown Cookbook