Read Chapter 3 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.
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, 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?
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 > 100, ] # patient 1630L
cardiac$triglyceride[cardiac$triglyceride == 0]
cardiac$hdlchol[cardiac$hdlchol == 0]
# 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 > 100] <- 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 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?
# 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:
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 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?
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 met 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?
# 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 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?
# 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 for adding a column). Compare the mean and the median before
and after. What has happened to the gap between them?
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.
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 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
Data 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?
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 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 TRUEs 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.
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 Data 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.
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.
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.)
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 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.
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