exercise_4_solution.Rmd

---
title: 'Exercises'
output: 
  html_document:
    toc: false
    code_folding: hide
---


```{r setup, echo=FALSE, purl = FALSE}
knitr::opts_chunk$set(echo=TRUE, message = FALSE, warning = FALSE, eval = FALSE, cache = FALSE)
options(width = 80)
SOLUTIONS <- TRUE
```

\  

## Exercise 4: Cleaning, summarising and linking data

\  

Read [Chapter 3](https://intro2r.com/data_r.html) to help you complete the questions in this exercise.

\  

This exercise picks up exactly where Exercise 3 left off and uses the same `cardiac` dataframe. If you are starting a fresh R session, either continue with the R script you wrote for Exercise 3 or run the two lines below before you go any further. They are the import from **Exercise 3, Q5** and the two factors from **Exercise 3, Q7**.

```{r Q0, echo=TRUE}
cardiac <- read.table('data/cardiacdata.txt', header = TRUE, sep = "\t", stringsAsFactors = TRUE)

cardiac$Fsex <- factor(cardiac$sex, levels = c(1, 2),
                       labels = c("Female", "Male"))

cardiac$Fsmoking <- factor(cardiac$smoking, levels = c(1, 2, 3),
                           labels = c("Current", "Ex", "Never"))
```

\  

1\. In Exercise 3 you learned how to extract rows that meet a condition. You can use that skill to check whether your data are believable. This matters a lot: a value that is simply wrong will not stop your code running, will not produce a warning, and will quietly change every result you calculate from it. Run `summary()` on your dataframe again and look at the minimum and maximum of each numeric variable. Three variables contain values that are not merely unusual but impossible. Two are impossibly low and one is impossibly high (hint: a body mass index below about 15 or above about 60 is essentially unheard of in a living adult, and nobody can survive with no cholesterol in their blood at all, so a blood fat concentration of exactly 0 is not a low reading, it is an impossible one). Find them, work out which patient the impossibly high one belongs to, and then set all of the offending values to `NA` so they cannot contaminate anything you calculate later (you will need the square bracket `[ ]` notation from [Section 3.4.2](https://intro2r.com/wrangling-data-frames.html#logical-indexes), this time on the left hand side of an assignment). Re-run `summary()` afterwards to check it worked. Why is `NA` the right answer here, rather than deleting the whole patient record or guessing what the value should have been?

```{r Q1, echo=SOLUTIONS}
summary(cardiac)

# hdlchol has a minimum of 0.000, triglyceride has a minimum of 0.000, and
# bmi has a maximum of 514.60. A blood concentration of zero is not a low
# reading, it is impossible; you cannot have no cholesterol in your blood and
# still be alive to take part in a study. And a body mass index of 514 is not
# a very large person. It is almost certainly 51.46 with the decimal point in
# the wrong place, although we have no way of confirming that.

# which patient has the impossible bmi?
cardiac[cardiac$bmi > 60, ]     # patient 1630L

# the two blood fats. Each of these prints the NA values already in the
# variable alongside the zeros you are looking for, for the same reason as
# Exercise 3, Q11. hdlchol and triglyceride have 2 missing values each and R
# won't guess whether a missing value is a zero, so read the 0s and ignore the
# NAs.
cardiac$triglyceride[cardiac$triglyceride == 0]   # one 0, and 2 NAs
cardiac$hdlchol[cardiac$hdlchol == 0]             # two 0s, and 2 NAs

# set the impossible values to NA. Putting the condition inside [ ] on the left
# of the arrow means 'the elements that match this condition become NA', and
# every other value is left alone.
cardiac$bmi[cardiac$bmi > 60] <- NA                  # 1 patient
cardiac$triglyceride[cardiac$triglyceride == 0] <- NA # 1 patient
cardiac$hdlchol[cardiac$hdlchol == 0] <- NA           # 2 patients

summary(cardiac)   # check: the minima and maximum are now sensible

#      tchol           hdlchol       triglyceride        bmi       
#  Min.   : 4.120   Min.   :0.720   Min.   :0.520   Min.   :17.57  
#  Median : 6.830   Median :1.360   Median :1.380   Median :25.20  
#  Mean   : 6.978   Mean   :1.421   Mean   :1.556   Mean   :25.72  
#  Max.   :11.660   Max.   :3.000   Max.   :4.670   Max.   :44.44  
#  NA's   :2        NA's   :4       NA's   :3       NA's   :1      

# Notice what that one bmi value was doing to the mean: 28.72 before, 25.72
# after. The median barely moved. 

# Why NA and not something else? Deleting the whole record throws away all the
# other measurements for that patient, which are potentially perfectly good.
# Guessing the value - even a very reasonable guess like 51.46 - means
# inventing data, and
# nobody reading your results afterwards could tell which numbers you measured
# and which you made up. NA says exactly what you know: there should be a value
# here, and it isn't usable. Every R function has an na.rm
# argument to cope with it.

# Note that you have changed the dataframe in your R session, not the file on
# disk. data/cardiacdata.txt still contains the original values, which is
# exactly as it should be: your raw data stays raw, and your script is the
# record of what you changed.
```

\  

2\. Another useful way to manipulate your dataframes is to sort the rows based on the value of a variable, or on a combination of variables. Rather counter-intuitively you should use the `order()` function to sort your dataframes, not the `sort()` function (see [Section 3.4.3](https://intro2r.com/wrangling-data-frames.html#ordering-data-frames) of the Introduction to R book for an explanation). Ordering dataframes uses the same logic you practised in Q12 in Exercise 2. Sort all rows in the `cardiac` dataframe by ascending order of systolic pressure within each level of smoking status, and assign the result to a variable with a sensible name. The trick is to remember that `order()` will take more than one variable, and that the order you give them in matters. Now take a look at the bottom of your sorted dataframe. Where have the patients with a missing smoking status ended up, and why?

```{r Q2, echo=SOLUTIONS}
# systolic within smoking status. smoking comes first because it is the
# grouping you want, systolic second because it sorts within each group
cardiac_sorted <- cardiac[order(cardiac$Fsmoking, cardiac$systolic), ]

# order() puts NAs last by default, so the 7 patients with no smoking status
# are all at the bottom rather than mixed in with the ones you can interpret.
tail(cardiac_sorted)
```

\    

3\. Often, we would like to summarise variables by, for example, calculating a mean, median or counting the number of observations. To do this for a single variable it's fairly straightforward:

\  

```{r Q3a, echo=TRUE}
mean(cardiac$age)          # mean age
median(cardiac$systolic)   # median systolic blood pressure
length(cardiac$tchol)      # number of observations
```

\  

Perhaps more interestingly, you might want to summarise one variable conditional on the level of another categorical variable, and to do several variables at once. The `aggregate()` function does exactly this (see [Section 3.5](https://intro2r.com/summarising-data-frames.html#summarising-data-frames) of the Introduction to R book, or `?aggregate`). Use `aggregate()` to calculate the mean age, systolic pressure, diastolic pressure and total cholesterol for each smoking category.

The grouping variable goes inside the `by` argument wrapped in `list()`, because a list is how `aggregate()` accepts more than one of them at once. Use `Fsmoking`, the factor you created in **Exercise 3, Q7**, rather than the original `smoking` column, so that your output comes back labelled Current, Ex and Never instead of 1, 2 and 3.

Run it exactly as it comes, with no special handling of anything, and then look hard at the `tchol` column. Something has gone wrong. What, and why?

```{r Q3b, echo=SOLUTIONS}
aggregate(cardiac[, c(2, 4, 5, 6)],
          by = list(smoking = cardiac$Fsmoking), FUN = mean)

#   smoking      age systolic diastolic tchol
# 1 Current 63.10060 143.5400  76.40000    NA
# 2      Ex 66.40769 145.0385  77.48077  6.82
# 3   Never 64.82019 139.2778  77.50000    NA

# Two of the three cholesterol means have come back as NA. Remember those
# missing values you have been tripping over since Exercise 3? They have not
# gone away: tchol still has two, one belonging to a current smoker and one to a
# never smoker, and mean() returns NA if even a single value handed to it is
# missing. Two missing values out of 163 patients have wiped out two of the
# three means, the age and blood pressure columns look perfectly healthy either
# side of them, and nothing warned you.
```

\  

4\. Now fix it. `aggregate()` hands any extra arguments straight on to the function you asked it to use, so the `na.rm` argument you first met in **Exercise 2, Q13** and used again in **Exercise 3, Q11** works here too. Add it and check that all three cholesterol means appear.

Then calculate the same means for each combination of smoking category and sex. Two grouping variables go inside the same `list()`, and again use the factor versions, `Fsmoking` and `Fsex`, rather than the original `smoking` and `sex` columns. When you have your answer, count up the patients in it. Are they all there?

```{r Q4, echo=SOLUTIONS}
# na.rm = TRUE is passed straight through to mean()
aggregate(cardiac[, c(2, 4, 5, 6)],
          by = list(smoking = cardiac$Fsmoking), FUN = mean, na.rm = TRUE)

#   smoking      age systolic diastolic    tchol
# 1 Current 63.10060 143.5400  76.40000 6.927143
# 2      Ex 66.40769 145.0385  77.48077 6.820000
# 3   Never 64.82019 139.2778  77.50000 7.157170

# two grouping variables, both named
aggregate(cardiac[, c(2, 4, 5, 6)],
          by = list(smoking = cardiac$Fsmoking, sex = cardiac$Fsex),
          FUN = mean, na.rm = TRUE)

# no, they are not all there. The 7 patients with no smoking status are dropped
# from every one of these summaries, because aggregate() has no group to put
# them in, and it does not tell you it has left them out. na.rm = TRUE fixed the
# missing cholesterol values; it does nothing at all about missing groups.
```

\  

5\. Knowing how many observations are present for each category (or combinations of categories) is useful to determine whether you have an adequate sample size (for subsequent modelling for example). Use the `table()` function to determine the number of patients in each smoking category (see [Section 3.5](https://intro2r.com/summarising-data-frames.html#summarising-data-frames) again for more information). Next use the same function to display the number of patients for each combination of smoking category and sex. Use `Fsmoking` and `Fsex` here as well, so that your table is labelled rather than numbered. Does `table()` tell you about the patients whose smoking status is missing?

```{r Q5, echo=SOLUTIONS}
# using table
table(cardiac$Fsmoking)

# Current      Ex   Never 
#      50      52      54 

table(cardiac$Fsmoking, cardiac$Fsex)

#           Female Male
#   Current     26   24
#   Ex          15   37
#   Never       37   17

# by default table() silently drops the missing values - 50 + 52 + 54 = 156,
# not 163. Ask for them explicitly:
table(cardiac$Fsmoking, useNA = "ifany")
```

\  

6\. Not every variable has its values distributed symmetrically about the centre. Look back at `triglyceride` in your `summary()` output. The mean, 1.556, sits above the median, 1.380, and the rest of the summary is lopsided in the same direction: the third quartile is 0.50 above the median while the first is only 0.32 below it, and the maximum is 3.29 above while the minimum is 0.86 below. Everything is stretched out to the right, which is what a long right hand tail looks like in a table of numbers, a few patients with much higher values than everybody else. Transforming to a log scale pulls that tail in, and it is a common way of dealing with skewed variables.

Create a new variable in the `cardiac` dataframe called `log_triglyceride`, holding the base 10 logarithm of `triglyceride` (see `?log10`, and [Section 3.4.4](https://intro2r.com/wrangling-data-frames.html#adding-columns-and-rows) for adding a column). Compare the mean and the median before and after. What has happened to the gap between them?

```{r Q6, echo=SOLUTIONS}
cardiac$log_triglyceride <- log10(cardiac$triglyceride)

mean(cardiac$triglyceride, na.rm = TRUE)         # 1.556
median(cardiac$triglyceride, na.rm = TRUE)       # 1.380

mean(cardiac$log_triglyceride, na.rm = TRUE)     # 0.149
median(cardiac$log_triglyceride, na.rm = TRUE)   # 0.140

# On the raw scale the mean is about 13% above the median. On the log scale the
# two are almost identical, which is what you expect once the long tail has been
# pulled in.

# Worth noticing: this only works because you set that triglyceride of 0.00 to
# NA in Q1. log10(0) is -Inf, which is not a number you can do anything useful
# with, and it would have spread into everything you calculated next.
```

\  

### Data linkage

Real analyses almost never involve a single file. The measurements you want are usually spread across several sources and you have to bring them together first. When those sources hold records about the same people, combining them into one dataset is called data linkage, and it is an important technique in health data science. Hospital admissions linked to prescribing records, a birth cohort linked to school attainment, a patient list linked to the death register. Each linkage answers a question that neither source could answer on its own, and this is the term we will use for it throughout the course.

The idea is simple. You match records on something that identifies the same person in both sources, usually an identifier such as a patient number. However, if implemented carelessly, analyses can go quietly wrong, which is what the next three questions are about. In R the function that does the work is `merge()` (see [Section 3.4.5](https://intro2r.com/wrangling-data-frames.html#merging-data-frames) of the Introduction to R book), and the operation it performs is called a join.

\  

7\. The patients in this study were followed up ten years after their first examination, and those follow-up measurements are held in a separate file. Download *'cardiac_followup.txt'* from the **[<i class="fa fa-download"></i> Data](data.html)** link, save it to your `data` directory and import it with `read.table()` into a variable called `followup`. It holds the patient number, `patno`, along with the systolic blood pressure and body mass index measured at follow-up, `systolic10` and `bmi10`. How many rows does it have, and how does that compare with `cardiac`? The `nrow()` function will tell you, and it does exactly what the name suggests: give it a dataframe and it returns the number of rows. 

Now use the `merge()` function to link the two together on the patient number, and assign the result to `cardiac_fu`. Left to its own devices `merge()` performs what is called an inner join, which keeps only the patients who appear in both sources. How many rows does the linked dataframe have, and can you explain why?

```{r Q7, echo=SOLUTIONS}
followup <- read.table('data/cardiac_followup.txt', header = TRUE,
                       sep = "\t", stringsAsFactors = TRUE)

nrow(cardiac)     # 163 patients at the start of the study
nrow(followup)    # 108 patients with follow-up measurements

cardiac_fu <- merge(cardiac, followup, by = "patno")
nrow(cardiac_fu)  # 108

# By default merge() keeps only the patients who appear in BOTH dataframes,
# which is called an inner join. 55 of the original patients have no follow-up
# measurements, so they have quietly disappeared. Nothing warned you about this.
# ALWAYS check the number of rows before and after a linkage.
```

\  

8\. Losing 55 patients without being told is exactly the sort of thing that can ruin an analysis, and it is why you check row counts every single time you link two sources. Often what you actually want is a left join, which keeps every patient from the first dataframe whether or not they have a match in the second. Look at [Section 3.4.5](https://intro2r.com/wrangling-data-frames.html#merging-data-frames) and the help file for `merge()`, and find the argument that switches `merge()` from an inner join to a left join. Redo the linkage keeping all 163 patients, this time assigning the result to a new dataframe called `cardiac_all`, and then check how many of those patients have a missing follow-up systolic pressure.

That last part needs a way of counting missing values, which you have not met yet. `is.na()` takes a variable and returns `TRUE` or `FALSE` for every value in it, `TRUE` where the value is missing. `sum()` then counts those `TRUE`s for you, because R treats `TRUE` as 1 and `FALSE` as 0. So `sum(is.na(cardiac_all$systolic10))` reads as 'how many follow-up systolic values are missing'. You will need it again in Q9.

```{r Q8, echo=SOLUTIONS}
cardiac_all <- merge(cardiac, followup, by = "patno", all.x = TRUE)
nrow(cardiac_all)                      # 163 - everybody is kept

sum(is.na(cardiac_all$systolic10))     # 55 patients have no follow-up

# all.x = TRUE keeps every row of x (the first dataframe) and fills the missing
# columns with NA. This is a left join.
```

\  

9\. One last linkage, and this one has a twist in it. The patients were also followed up for hospital admissions over the same ten years. Download *'cardiac_admissions.txt'* from the **[<i class="fa fa-download"></i> Data](data.html)** link, save it to your `data` directory and import it into a variable called `admissions`. It holds one row for each patient who was admitted at least once, and the number of times they were admitted. Note these admission records are simulated but reflect how these types of data are typically formatted.

a) Merge the `admissions` dataframe and `cardiac_all`, keeping all 163 patients exactly as you did in Q8. Check you still have 163 rows, then count how many patients ended up with a missing `n_admissions`.

b) Now the twist. Up until now, an `NA` has meant 'there should be a value here and we do not have it (missing)'. Does it mean that here? Work out what a missing `n_admissions` actually tells you about that patient (hint: should it actually be a value, and if so, what value). Once you have figured out what `NA` should actually represent, replace these `NA` values (hint: you cannot identify missing values with a conditional statement like `cardiac_all$n_admissions == NA`. Use the `is.na()` function from Q8 instead.)

```{r Q9, echo=SOLUTIONS}
admissions <- read.table('data/cardiac_admissions.txt', header = TRUE,
                         sep = "\t", stringsAsFactors = TRUE)

# a)
cardiac_all <- merge(cardiac_all, admissions, by = "patno", all.x = TRUE)
nrow(cardiac_all)                        # 163 - everybody is still here

sum(is.na(cardiac_all$n_admissions))     # 53 patients have no admissions record

# b)
# Those 53 NAs do NOT mean 'we do not know'. Those patients are absent from the
# admissions file because they were never admitted, so the right number for them
# is 0, and we know that for certain. The NA is an artefact of how a left join
# fills gaps, not a statement about our knowledge.
cardiac_all$n_admissions[is.na(cardiac_all$n_admissions)] <- 0

# Always ask what a missing value means before you decide how to handle it. The
# answer is not always the same.
```

\  

10\. Ok, we have spent quite a bit of time (and energy) learning how to clean, summarise and link dataframes. The last thing we need to cover is how to export a dataframe from R to an external file (see [Section 3.6](https://intro2r.com/exporting-data.html#exporting-data) of the book for more details). Export your cleaned and linked dataframe `cardiac_all` to a file called 'cardiac_clean.txt' in the `output` directory you created in Exercise 1. To do this you will need to use the `write.table()` function. You want to include the variable names in the first row of the file, but you don't want to include the row names. Also, make sure the file is a tab delimited file. Once you have created your file, try to open it in Microsoft Excel (or open source equivalent). Finally, and this is the part people forget: add a comment block at the top of your R script listing every change you made to the data and why. Which values did you set to `NA`, in which variables, and which patients did they belong to? Someone reading your work in six months, including you, needs to be able to answer that without re-running anything.


```{r Q10, echo=SOLUTIONS}
write.table(cardiac_all, "output/cardiac_clean.txt", col.names = TRUE,
            row.names = FALSE, sep = "\t")

# Decision log - the sort of thing that belongs at the top of your script:
#
# Data: data/cardiacdata.txt, 163 patients, imported unchanged.
# 1. Fsex and Fsmoking created as factor versions of sex (Female, Male) and
#    smoking (Current, Ex, Never). The original coded columns were left
#    untouched.
# 2. hdlchol of 0.00 set to NA (2 patients). A HDL cholesterol of zero is
#    physiologically impossible and is almost certainly a missing value that
#    was recorded as 0.
# 3. triglyceride of 0.00 set to NA (1 patient). Same reasoning.
# 4. bmi of 514.60 set to NA (patient 1630L). Almost certainly 51.46 with the
#    decimal point misplaced, but as we cannot confirm that, NA rather than
#    a correction.
# 5. log_triglyceride added, the base 10 log of triglyceride, taken after the
#    cleaning in step 3.
# 6. Ten year follow-up measurements linked on from data/cardiac_followup.txt,
#    keeping all 163 patients; 55 have no follow-up data.
# 7. Admissions linked on from data/cardiac_admissions.txt (SIMULATED data),
#    keeping all 163 patients. n_admissions set to 0, not left as NA, for the 53
#    patients with no admission record, because a patient who was never admitted
#    has zero admissions rather than an unknown number.
```

\  

End of Exercise 4