Tutorial 4.3: Exploring the captopril dataset

In this tutorial, you will learn how to explore paired data.

1 The captopril dataset

The captopril dataset stems from a small experiment with 15 patients with hypertension. For each patient systolic and diasystolic blood pressure measurements where taken before and after administering captopril.

2 Import the data

library(tidyverse)

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

## Parsed with column specification:
## cols(
##   id = col_double(),
##   SBPb = col_double(),
##   DBPb = col_double(),
##   SBPa = col_double(),
##   DBPa = col_double()
## )

head(captopril)

## # A tibble: 6 x 5
##      id  SBPb  DBPb  SBPa  DBPa
##   <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     1   210   130   201   125
## 2     2   169   122   165   121
## 3     3   187   124   166   121
## 4     4   160   104   157   106
## 5     5   167   112   147   101
## 6     6   176   101   145    85

3 Data visualization

One straightforward possibility to visualize the four types of blood pressure values is by adopting the gather function from tidyverse. It will reshape the dataframe, such that we have have a single variable type, which points at one of the four blood pressure types, and bp, which points at the actual value for each type for each patient.

captopril %>% 
  gather(type,bp,-id)

## # A tibble: 60 x 3
##       id type     bp
##    <dbl> <chr> <dbl>
##  1     1 SBPb    210
##  2     2 SBPb    169
##  3     3 SBPb    187
##  4     4 SBPb    160
##  5     5 SBPb    167
##  6     6 SBPb    176
##  7     7 SBPb    185
##  8     8 SBPb    206
##  9     9 SBPb    173
## 10    10 SBPb    146
## # … with 50 more rows

3.1 Barplot

A Barplot is a plot that you will commonly find in papers.

captopril %>% 
  gather(type,bp,-id) %>% 
  group_by(type) %>%
     summarize_at("bp",
               list(mean=~mean(.,na.rm=TRUE),
                    sd=~sd(.,na.rm=TRUE),
                    n=function(x) x%>%is.na%>%`!`%>%sum)) %>%
  mutate(se=sd/sqrt(n)) %>%
  ggplot(aes(x=type,y=mean,fill=type)) + 
  scale_fill_brewer(palette="RdGy") +
  theme_bw() +
  geom_bar(stat="identity") + 
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),width=.2) +
  ggtitle("Barplot of different blood pressure measures") +
  ylab("blood pressure (mmHg)")

A barplot, however, is not very informative. The height of the bars only provides us with information of the mean blood pressure. However, we don’t see the actual underlying values, so we for instance don’t have any information on the spread of the data. It is usually more informative to represent to underlying values as raw as possible. Boxplots are ideal for this!

Note that it is possible to add the raw data on the barplot, but we still would not see any measures of the spread, such as the interquartile range.

3.2 Boxplot

With the boxplot, we get a lot of useful information. We immediately see multiple features of the spread of the data, such as the median, the 25% and 75% quantiles and outliers. Since we only have 15 raw values (patients), we can easily add them to the plot without getting messy.

captopril %>% 
  gather(type,bp,-id) %>% 
  ggplot(aes(x=type,y=bp,fill=type)) + 
  scale_fill_brewer(palette="RdGy") +
  theme_bw() +
  geom_boxplot(outlier.shape=NA) + 
  geom_jitter(width = 0.2) +
  ggtitle("Boxplot of different blood pressure measures") +
  ylab("blood pressure (mmHg)") + stat_summary(fun.y=mean, geom="point", shape=5, size=3, color="black", fill="black")

In terms of interpretation, we can see that the median systolic and diastolic blood pressure values are lower after treatment with captopril than before. If we want to have visual inference on the mean values (cfr. class on t-test), we can add them to the plot stat_summary function.

4 Paired data

An important feature of this dataset is that it contains paired data; for each patient, we have blood pressure values (systolic and diastolic) before and after treatment with captopril.

We first plot the systolic blood pressure measurements. We make a plot with the data points for the two systolic blood pressure measurements in the patient id.

captopril %>% 
  gather(type,bp,-id) %>% 
  filter(type%in%c("SBPa","SBPb")) %>%
  ggplot(aes(x=id,y=bp,color=type)) +
  geom_point() +
  scale_color_manual(values = c("darkgreen","red"))

We see that for all patients, the systolic blood pressure is lower after captopril treatment than before. Note that we could not see this from the boxplot, directly.

Another good way to visualise the paired data is by using a line plot that is grouped by subject. Note, that we make a factor of the type and that we choose order of the levels so that the level “b” before becomes before the level “a”. Otherwise the after measurements will appear first in the plot.

captopril %>% 
  gather(type,bp,-id) %>%
  mutate(type=factor(type,levels=c("SBPb","SBPa","DBPb","DBPa"))) %>%
  filter(type%in%c("SBPa","SBPb")) %>%
  ggplot(aes(x=type,y=bp)) +
  geom_line(aes(group = id)) +
  geom_point()

captopril %>% 
  gather(type,bp,-id) %>%
  mutate(type=factor(type,levels=c("SBPb","SBPa","DBPb","DBPa"))) %>%
  filter(type%in%c("DBPa","DBPb")) %>%
  ggplot(aes(x=type,y=bp)) +
  geom_line(aes(group = id)) +
  geom_point()

Analogously, we may immediately difference the after measurement from the before measurement because we have paired data.

captopril %>% 
  mutate(bp_diff = SBPa-SBPb) %>%
  ggplot(aes(x="",y=bp_diff)) +
  geom_boxplot(outlier.shape=NA) +
  geom_jitter(width = 0.2) + 
  ggtitle("Boxplot of the difference in blood pressure") +
  ylab("blood pressure difference (mmHg)") + stat_summary(fun.y=mean, geom="point", shape=5, size=3, color="black", fill="black") +
  theme_bw()+
  ylim(-40,10) +
  geom_hline(yintercept = 0, color="red") # adds horizontal line to plot

The mean difference is lower than zero. It seems that on average the captopril treatment has lowered the blood pressure by about 20 mmHg.

5 Assumptions: QQ-plot

In the tutorial on hypothesis testing we will assess if this difference is statistically significant. We will introduce a statistical test that relies on the assumption that the differences are normally distributed. We can check this with a QQ-plot:

captopril %>% 
  mutate(bp_diff = SBPa-SBPb) %>%
  select(bp_diff) %>%
  ggplot(aes(sample=bp_diff)) +
  geom_qq() +
  geom_qq_line()

The plot shows no large deviations of the straight line indicating that the data are approximately normally distributed.

6 Descriptive statistics

Descriptive statistics for all blood pressure measurements are found below.

captopril %>% 
  gather(type,bp,-id) %>% 
  group_by(type) %>%
     summarize_at("bp",
               list(mean=~mean(.,na.rm=TRUE),
                    sd=~sd(.,na.rm=TRUE),
                    n=function(x) x%>%is.na%>%`!`%>%sum))%>%
  mutate(se=sd/sqrt(n))

## # A tibble: 4 x 5
##   type   mean    sd     n    se
##   <chr> <dbl> <dbl> <int> <dbl>
## 1 DBPa   103.  12.6    15  3.24
## 2 DBPb   112.  10.5    15  2.70
## 3 SBPa   158   20.0    15  5.16
## 4 SBPb   177.  20.6    15  5.31

It is also very useful to calculate the descriptive statistics for the systolic blood pressure difference because we can assess the effect size of administering captopril for every patient.

captopril %>% 
  mutate(bp_diff = SBPa-SBPb) %>%
     summarize_at("bp_diff",
               list(mean=~mean(.,na.rm=TRUE),
                    sd=~sd(.,na.rm=TRUE),
                    n=function(x) x%>%is.na%>%`!`%>%sum)) %>%
  mutate(se=sd/sqrt(n))

## # A tibble: 1 x 4
##    mean    sd     n    se
##   <dbl> <dbl> <int> <dbl>
## 1 -18.9  9.03    15  2.33

7 Wrap-up

Next, time if you explore paired data, you can immediately adopt the line plot and the boxplot of the differences from Section 3.3. so as to acknowledge the pairing and perform the descriptive statistics in section 3.4. It might also be useful to plot boxplots of the raw blood pressure measurements to get an idea on their varability.

It is not useful to report the barplot.

7.1 Advantage of paired design

Every patient acts as their own control: we can estimate the effect of the treatment for each patient

captopril %>% 
  ggplot(aes(SBPb,SBPa)) +
  geom_point() + 
  ggtitle(paste("correlation = ",round(cor(captopril$SBPb,captopril$SBPa),3)))

captopril %>% 
  mutate(bp_diff = SBPa-SBPb) %>%
  select(SBPb,SBPa,bp_diff) %>%
  summarize_all(sd)

## # A tibble: 1 x 3
##    SBPb  SBPa bp_diff
##   <dbl> <dbl>   <dbl>
## 1  20.6  20.0    9.03

Calculate variance based on fact that data are paired.

captopril %>% summarise(sd_diff=sqrt(sd(SBPb)^2 + sd(SBPa)^2-2*sd(SBPa)*sd(SBPb)*cor(SBPa,SBPb)))

## # A tibble: 1 x 1
##   sd_diff
##     <dbl>
## 1    9.03

Note, that the standard deviations on the blood pressure measurements are more than twice as large as the standard deviations on the difference. This is because each patient acts as their own control. As such, we can correct for their baseline blood pressure when estimating the effect of the treatment.

---

Home