7  Data in rtemis

This chapter reads the three datasets used in the book from the files in data/, inspects them, and derives each step of the functions read_heart_failure(), read_diabetes() and read_parkinsons() defined in the shared setup file. The heart failure data receive the most attention, since their outcome must be redefined and one variable excluded (Section 6.3, Section 6.4).

7.1 Reading and checking data

read() reads delimited text, Parquet, Excel, and other formats into a data.table, and cleans column names into syntactic R names by default. The book uses data.table objects throughout. A data.table is a data.frame with additional syntax: dt[i, j] selects rows i and computes j with column names in scope, and := adds, modifies or removes columns in place.

hf_raw <- read("data/heart_failure_clinical_records_dataset.csv")

check_data() reports the number of cases and features, the type of each column, and common issues: constant columns, duplicate rows and missing values.

check_data(hf_raw)

  hf_raw: A data.table with 299 rows and 13 columns.

  Data types
  * 3 numeric features
  * 10 integer features
  * 0 factors
  * 0 character features
  * 0 date features

  Issues
  * 0 constant features
  * 0 duplicate cases
  * 0 missing values

  Recommendations
  * Everything looks good

The file contains 299 patients and 13 variables, with no duplicates and no missing values. Every variable is stored as a number, including the binary variables, which are stored as integer codes 0 and 1.

7.2 The heart failure data

The data describe 299 patients aged 40 years or older, admitted with heart failure due to left ventricular systolic dysfunction, NYHA class III or IV, to the Institute of Cardiology and Allied Hospital, Faisalabad, Pakistan, between April and December 2015 (Ahmad et al. 2017). Chicco and Jurman (2020) document the variables:

Table 7.1: Heart failure variables

Variable	Description	Unit or coding
age	Age	Years
anaemia	Decrease of red blood cells or hemoglobin	0, 1
creatinine_phosphokinase	Level of the CPK enzyme in the blood	mcg/L
diabetes	Diabetes	0, 1
ejection_fraction	Percentage of blood leaving the heart at each contraction	%
high_blood_pressure	Hypertension	0, 1
platelets	Platelets in the blood	Platelets per µL
serum_creatinine	Level of creatinine in the blood	mg/dL
serum_sodium	Level of sodium in the blood	mEq/L
sex	Woman or man	0, 1
smoking	Smoking	0, 1
time	Follow-up period	Days
DEATH_EVENT	Death during the follow-up period	0, 1

The sex coding is not stated in the table; Ahmad et al. (2017) report 105 women and 194 men, and sex equals 1 for 194 patients, so 1 denotes a man. Chicco and Jurman (2020) give the platelet unit as kiloplatelets per microliter; the values in the file, from 25,100 to 850,000, are platelets per microliter.

7.2.1 Follow-up time

time is the number of days from baseline to death or to the end of follow-up. Its distribution differs sharply between patients who died and patients who did not:

draw_boxplot(
  split(hf_raw[["time"]], ifelse(hf_raw[["DEATH_EVENT"]] == 1L, "Died", "Censored")),
  ylab = "Follow-up time (days)"
)

The median follow-up is 44.5 days among patients who died and 172 days among those censored. Follow-up of censored patients ranges from 12 to 285 days. Short follow-up is therefore largely a consequence of death, and time is not known at baseline. It is a proxy of the outcome (Section 6.4).

Its effect on assessment can be measured directly. The following fits a logistic regression of the published outcome DEATH_EVENT on all variables, with and without time, and estimates the AUC by 10-fold cross-validation:

Since := modifies a data.table in place, copy() creates an independent copy before modification, leaving hf_raw unchanged:

hf_published <- copy(hf_raw)
hf_published[, DEATH_EVENT := factor(DEATH_EVENT, levels = c(0L, 1L), labels = c("No", "Yes"))]
cv <- setup_KFold(n_resamples = 10L, seed = 2026L)
mod_with_time <- train(
  hf_published,
  hyperparameters = setup_GLM(),
  outer_resampling_config = cv,
  execution_config = setup_SerialExecution(seed = 2026L),
  verbosity = 0L
)
mod_without_time <- train(
  hf_published[, !"time"],
  hyperparameters = setup_GLM(),
  outer_resampling_config = cv,
  execution_config = setup_SerialExecution(seed = 2026L),
  verbosity = 0L
)
c(
  with_time = mod_with_time@metrics_test@mean_metrics[["auc"]],
  without_time = mod_without_time@metrics_test@mean_metrics[["auc"]]
)

   with_time without_time 
      0.8799       0.7793 

Including time raises the cross-validated AUC from 0.779 to 0.880. The higher figure describes a model that cannot be applied at baseline, since time is unknown until follow-up ends. Cross-validation does not detect this form of leakage: the leaking variable is available in every training and test set alike.

7.2.2 A fixed-horizon outcome

DEATH_EVENT records death at any point during follow-up, over periods that differ between patients. The outcome used in this book is death within 60 days of baseline, defined for each patient as:

• Yes, if the patient died on or before day 60;
• No, if the patient was followed for more than 60 days, whether or not the patient died later;
• unknown, if the patient was alive at the end of a follow-up of 60 days or less.

horizon <- 60L
died_by_horizon <- hf_raw[["DEATH_EVENT"]] == 1L & hf_raw[["time"]] <= horizon
followed_past_horizon <- hf_raw[["time"]] > horizon
status_unknown <- !died_by_horizon & !followed_past_horizon
c(
  yes = sum(died_by_horizon),
  no = sum(followed_past_horizon),
  unknown = sum(status_unknown)
)

    yes      no unknown 
     54     236       9 

The 9 patients with unknown status are excluded. The outcome is then a factor with levels "No" and "Yes", so that death is the positive class:

hf <- hf_raw[!status_unknown]
hf[, death := factor(
  ifelse(DEATH_EVENT == 1L & time <= horizon, "Yes", "No"),
  levels = c("No", "Yes")
)]

time and DEATH_EVENT are removed, since neither is known at baseline:

hf[, c("time", "DEATH_EVENT") := NULL]

7.2.3 Variable types and units

The binary variables are converted from integer codes to factors with labeled levels, and platelets are expressed in thousands per µL, the unit in which platelet counts are usually reported:

binary <- c("anaemia", "diabetes", "high_blood_pressure", "smoking")
hf[, (binary) := lapply(.SD, factor, levels = c(0L, 1L), labels = c("No", "Yes")),
  .SDcols = binary
]
hf[, sex := factor(sex, levels = c(0L, 1L), labels = c("Female", "Male"))]
hf[, platelets := platelets / 1000]

These are the steps of read_heart_failure(), which returns the same data:

all.equal(hf, read_heart_failure())

[1] TRUE

check_data(hf)

  hf: A data.table with 290 rows and 12 columns.

  Data types
  * 3 numeric features
  * 3 integer features
  * 6 factors, of which 0 are ordered
  * 0 character features
  * 0 date features

  Issues
  * 0 constant features
  * 0 duplicate cases
  * 0 missing values

  Recommendations
  * Everything looks good

7.2.4 Distributions

draw_histogram() shows the distribution of a variable, here by outcome group:

draw_histogram(
  hf[["ejection_fraction"]],
  group = hf[["death"]],
  xlab = "Ejection fraction (%)"
)

Serum creatinine and creatinine phosphokinase are strongly right-skewed:

draw_histogram(hf[["serum_creatinine"]], xlab = "Serum creatinine (mg/dL)")

draw_histogram(
  hf[["creatinine_phosphokinase"]],
  xlab = "Creatinine phosphokinase (mcg/L)"
)

The median serum creatinine is 1.1 mg/dL and the maximum 9.4 mg/dL; the median creatinine phosphokinase is 255 mcg/L and the maximum 7861 mcg/L. In models that are sensitive to the scale of the features, such as linear models and distance-based methods, a few extreme values can have a large influence on the fit (Chapter 12).

7.3 The diabetes data

The diabetes data contain ten baseline variables for 442 patients with diabetes and a quantitative measure of disease progression one year after baseline (Efron et al. 2004). The file is tab-delimited; read() detects the delimiter:

db_raw <- read("data/diabetes.tab.txt")
head(db_raw, 3)

     AGE   SEX   BMI    BP    S1    S2    S3    S4    S5    S6     Y
   <int> <int> <num> <num> <int> <num> <num> <num> <num> <int> <int>
1:    59     2  32.1   101   157  93.2    38     4 4.860    87   151
2:    48     1  21.6    87   183 103.2    70     3 3.892    69    75
3:    72     2  30.5    93   156  93.6    41     4 4.673    85   141

The source names the baseline variables age, sex, body mass index (BMI), average blood pressure (BP), and six blood serum measurements (S1 to S6). The columns are renamed to lower case, the outcome is named progression, and sex, coded 1 and 2 without a documented mapping, becomes a factor:

db <- copy(db_raw)
setnames(db, c(
  "age", "sex", "bmi", "bp", "s1", "s2", "s3", "s4", "s5", "s6",
  "progression"
))
db[, sex := factor(sex)]
all.equal(db, read_diabetes())

[1] TRUE

check_data(db)

  db: A data.table with 442 rows and 11 columns.

  Data types
  * 6 numeric features
  * 4 integer features
  * 1 factor, which is not ordered
  * 0 character features
  * 0 date features

  Issues
  * 0 constant features
  * 0 duplicate cases
  * 0 missing values

  Recommendations
  * Everything looks good

draw_histogram(db[["progression"]], xlab = "Disease progression at one year")

7.4 The Parkinson’s telemonitoring data

The Parkinson’s data contain 5875 voice recordings from 42 people with early-stage Parkinson’s disease, recorded at home during a six-month trial of a telemonitoring device (Tsanas et al. 2010). Each recording has 16 voice measures, derived from variation in the fundamental frequency (jitter), variation in amplitude (shimmer), noise-to-tonal ratios, and nonlinear dynamical measures (Little et al. 2009). read() cleans names such as subject# and Jitter(%) to subject and Jitter:

pk_raw <- read("data/parkinsons_updrs.data")
names(pk_raw)

 [1] "subject"       "age"           "sex"           "test_time"    
 [5] "motor_UPDRS"   "total_UPDRS"   "Jitter"        "Jitter_Abs"   
 [9] "Jitter_RAP"    "Jitter_PPQ5"   "Jitter_DDP"    "Shimmer"      
[13] "Shimmer_dB"    "Shimmer_APQ3"  "Shimmer_APQ5"  "Shimmer_APQ11"
[17] "Shimmer_DDA"   "NHR"           "HNR"           "RPDE"         
[21] "DFA"           "PPE"          

The patient identifier, subject, shows the repeated-measures structure:

recordings <- table(pk_raw[["subject"]])
length(recordings)

[1] 42

range(recordings)

[1] 101 168

Each patient contributes between 101 and 168 recordings. Recordings from one patient share the patient’s age, sex and disease severity, so they are not independent cases. Chapter 11 uses subject to keep each patient’s recordings together when dividing the data.

The dataset has two UPDRS outcomes: the motor score, which assesses motor function, and the total score, of which the motor score is one component.

draw_scatter(
  pk_raw[["motor_UPDRS"]], pk_raw[["total_UPDRS"]],
  xlab = "Motor UPDRS", ylab = "Total UPDRS"
)

The two scores have a correlation of 0.947. With total UPDRS as the outcome, the motor score is a component of the outcome and is not admissible as a feature. test_time, the number of days since recruitment, is also removed, since the question concerns prediction from the voice measures, age and sex. Sex is coded 0 for men and 1 for women in the source:

pk <- copy(pk_raw)
pk[, sex := factor(sex, levels = c(0L, 1L), labels = c("Male", "Female"))]
pk[, c("motor_UPDRS", "test_time") := NULL]
pk <- set_outcome(pk, "total_UPDRS")
all.equal(pk, read_parkinsons())

[1] TRUE

7.5 Summary

Table 7.2: Datasets used in the book

Dataset	Reader	Cases	Features	Outcome
Heart failure	read_heart_failure()	290 patients	11	death: death within 60 days (No, Yes)
Diabetes	read_diabetes()	442 patients	10	progression: progression at one year
Parkinson’s	read_parkinsons()	5875 recordings, 42 patients	18, plus subject	total_UPDRS

Ahmad, Tanvir, Assia Munir, Sajjad Haider Bhatti, Muhammad Aftab, and Muhammad Ali Raza. 2017. “Survival Analysis of Heart Failure Patients: A Case Study.” PLOS ONE 12 (7): e0181001. https://doi.org/10.1371/journal.pone.0181001.

Chicco, Davide, and Giuseppe Jurman. 2020. “Machine Learning Can Predict Survival of Patients with Heart Failure from Serum Creatinine and Ejection Fraction Alone.” BMC Medical Informatics and Decision Making 20 (1). https://doi.org/10.1186/s12911-020-1023-5.

Efron, Bradley, Trevor Hastie, Iain Johnstone, and Robert Tibshirani. 2004. “Least Angle Regression.” The Annals of Statistics 32 (2). https://doi.org/10.1214/009053604000000067.

Little, M. A., P. E. McSharry, E. J. Hunter, J. Spielman, and L. O. Ramig. 2009. “Suitability of Dysphonia Measurements for Telemonitoring of Parkinson’s Disease.” IEEE Transactions on Biomedical Engineering 56 (4): 1015–22. https://doi.org/10.1109/tbme.2008.2005954.

Tsanas, A., M. A. Little, P. E. McSharry, and L. O. Ramig. 2010. “Accurate Telemonitoring of Parkinson’s Disease Progression by Noninvasive Speech Tests.” IEEE Transactions on Biomedical Engineering 57 (4): 884–93. https://doi.org/10.1109/tbme.2009.2036000.