August 2016 – The Lab-R-torian

The Problem

As you know in Clinical Chemistry, we are not always writing a major paper but sometimes just preparing a short-report to answer a technical question that we've encounted at work. For shorter papers, journals often have more stringent rules about how many figures you can submit and even sometimes forbid multipanelled figures. In these situations, we might want to cram a little more into your figure than we might otherwise. In a recent submission, I wanted to produce a difference plot of immunoassay results before and after storage but I also wanted to show the distribution of the results using a histogram–but this would have counted as two separate figures.

However, thanks to some fine work by Dean Attali of UBC Department of Statistics where he works with R-legend Jenny Bryan, it is quite easy to add marginal histograms to a Bland Altman (or any other scatter) plot using the ggExtra package.

How To

Let's make some fake data for a Bland Altman plot. Let's pretend that we are measuring the same quantity by immunoassay at baseline and after 1 year of storage at -80 degrees. We'll add some heteroscedastic error and create some apparent degradation of about 20%:

set.seed(10) #make predictable random data
baseline <- rlnorm(100, 0, 1)
post <- 0.8*baseline + rnorm(100, 0, 0.10*baseline)
plot(baseline,post)
abline(lm(post ~ baseline))
abline(0, 1, col="red", lty = 2)

set.seed(10) #make predictable random data

baseline <- rlnorm(100, 0, 1)

post <- 0.8*baseline + rnorm(100, 0, 0.10*baseline)

plot(baseline,post)

abline(lm(post ~ baseline))

abline(0, 1, col="red", lty = 2)

plot of chunk unnamed-chunk-1

Or if we plot this in the ggplot() paradigm

library(ggplot2)
my.data <- data.frame(baseline, post)
ggplot(my.data, aes(x=baseline, y=post)) +
    theme_bw() + 
    geom_point(shape=1) +    # Use hollow circles
    geom_smooth(method=lm) +  # Add linear regression line 
    geom_abline(slope = 1, intercept = 0, linetype = 2, colour = "red")

library(ggplot2)

my.data <- data.frame(baseline, post)

ggplot(my.data, aes(x=baseline, y=post)) +

theme_bw() +

geom_point(shape=1) + # Use hollow circles

geom_smooth(method=lm) + # Add linear regression line

geom_abline(slope = 1, intercept = 0, linetype = 2, colour = "red")

plot of chunk unnamed-chunk-2

Now we will prepare the difference data:

diff <- (post - baseline)
diffp <- (post - baseline)/baseline*100
sd.diff <- sd(diff)
sd.diffp <- sd(diffp)
my.data <- data.frame(baseline, post, diff, diffp)

diff <- (post - baseline)

diffp <- (post - baseline)/baseline*100

sd.diff <- sd(diff)

sd.diffp <- sd(diffp)

my.data <- data.frame(baseline, post, diff, diffp)

In standard Bland Altman plots, one plots the difference between methods against the average of the methods, but in this case, the x-axis should be the baseline result, because that is the closest thing we have to the truth.

library(ggExtra)
diffplot <- ggplot(my.data, aes(baseline, diff)) + 
  geom_point(size=2, colour = rgb(0,0,0, alpha = 0.5)) + 
  theme_bw() + 
  #when the +/- 2SD lines will fall outside the default plot limits 
  #they need to be pre-stated explicitly to make the histogram line up properly. 
  #Thanks to commenter for noticing this.
  ylim(mean(my.data$diff) - 3*sd.diff, mean(my.data$diff) + 3*sd.diff) +
  geom_hline(yintercept = 0, linetype = 3) +
  geom_hline(yintercept = mean(my.data$diff)) +
  geom_hline(yintercept = mean(my.data$diff) + 2*sd.diff, linetype = 2) +
  geom_hline(yintercept = mean(my.data$diff) - 2*sd.diff, linetype = 2) +
  ylab("Difference pre and post Storage (mg/L)") +
  xlab("Baseline Concentration (mg/L)")

#And now for the magic - we'll use 25 bins
ggMarginal(diffplot, type="histogram", bins = 25)

library(ggExtra)

diffplot <- ggplot(my.data, aes(baseline, diff)) +

geom_point(size=2, colour = rgb(0,0,0, alpha = 0.5)) +

theme_bw() +

#when the +/- 2SD lines will fall outside the default plot limits

#they need to be pre-stated explicitly to make the histogram line up properly.

#Thanks to commenter for noticing this.

ylim(mean(my.data$diff) - 3*sd.diff, mean(my.data$diff) + 3*sd.diff) +

geom_hline(yintercept = 0, linetype = 3) +

geom_hline(yintercept = mean(my.data$diff)) +

geom_hline(yintercept = mean(my.data$diff) + 2*sd.diff, linetype = 2) +

geom_hline(yintercept = mean(my.data$diff) - 2*sd.diff, linetype = 2) +

ylab("Difference pre and post Storage (mg/L)") +

xlab("Baseline Concentration (mg/L)")

#And now for the magic - we'll use 25 bins

ggMarginal(diffplot, type="histogram", bins = 25)

plot of chunk unnamed-chunk-4

So that is the difference plot for the absolute difference. We can also obviously do the percent difference.

diffplotp <- ggplot(my.data, aes(baseline, diffp)) + 
  geom_point(size=2, colour = rgb(0,0,0, alpha = 0.5)) + 
  theme_bw() + 
  geom_hline(yintercept = 0, linetype = 3) +
  geom_hline(yintercept = mean(my.data$diffp)) +
  geom_hline(yintercept = mean(my.data$diffp) + 2*sd.diffp, linetype = 2) +
  geom_hline(yintercept = mean(my.data$diffp) - 2*sd.diffp, linetype = 2) +
  ylab("Difference pre and post Storage (%)") +
  xlab("Baseline Concentration (mg/L)")


ggMarginal(diffplotp, type="histogram", bins = 25)

diffplotp <- ggplot(my.data, aes(baseline, diffp)) +

geom_point(size=2, colour = rgb(0,0,0, alpha = 0.5)) +

theme_bw() +

geom_hline(yintercept = 0, linetype = 3) +

geom_hline(yintercept = mean(my.data$diffp)) +

geom_hline(yintercept = mean(my.data$diffp) + 2*sd.diffp, linetype = 2) +

geom_hline(yintercept = mean(my.data$diffp) - 2*sd.diffp, linetype = 2) +

ylab("Difference pre and post Storage (%)") +

xlab("Baseline Concentration (mg/L)")

ggMarginal(diffplotp, type="histogram", bins = 25)

plot of chunk unnamed-chunk-5

Kickin' it Old School

You can also do this in a non-ggplot() paradigm using base plotting utilities as described in this R-bloggers post.

Conclusion

And that, friends, is a way of squishing in a histogram of your sample concentrations into your difference plot which allows you to graphically display your sampling distribution and justify whether you would use parametric or non-parametric statistics to assess the extent of loss of immunoreactivity from storage.

And speaking of scatterplots

“…then the Lord your God will restore your fortunes and have compassion on you and gather you again from all the nations where he scattered you.”
Deut 30:3

Background

Back in 2011 I was not aware of any tool in R for Passing Bablok (PB) regression, a form of robust regression described in a series of three papers in Clinical Chemistry and Laboratory Medicine (then J Clin Chem and Biochem) available here, here and here. For reasons that are not entirely clear to me, this regression methodology is favoured by clinical chemists but seems largely ignored by other disciplines. However since reviewers clinical chemistry journals will demand the use of PB regression, it seemed expeditious to me to code it in R. This is what spawned a small project for a piece of software to do PB (and Deming and ordinary least squares) regression using a self-contained executable that could be downloaded, unzipped on a Windows Desktop and just ran. You can download here and instructions for installation and use are here and here respectively. The calculations are all done in R, the GUI is built with Python and Py-Qt4 and the executable with cx_freeze. I made it run without an installer because hospital IT often refuse to install software that has not been officially vetted and purchased. The tool was a lot more popular than I anticipated now having about 2000 downloads. In any case, maintenance, upgrades, bug fixing and dealing with operating system updates that break things (like OSX El Capitan's security policies) are no-fun so a Shiny based solution to the same problem makes a lot of sense.

Update

Since 2011, statisticians at Roche Diagnostics programmed the mcr package for PB and Deming regression. Additionally, there is also the MethComp package and the deming package from the Mayo Clinic which both offer PB regression.

Shiny App

Enter Burak Bahar, a like-minded Clinical Pathologist who is currently doing a fellowship at Yale. He liked my cp-R program but he saw the need for a web-based equivalent.

Burak and his wife Ayse, also a physician, have coded a Shiny App for doing Deming, PB and least squares regression in R which is capable of producing publication quality figures and provides all the regression statistics you would need for method-validation or publication. It can also produce a regression report in PDF, Word or HTML. The dynamic duo of the Bahar-MDs deserve all credit here as my only contribution related to suggestions related to usability. This project was presented at the 2016 American Association of Clinical Chemistry meeting in Philadelphia.

The app URL is bahar.shinyapps.io/method_compare. Go to the data tab on the left and then cut and paste your data from an spreadsheet program. Shortcuts CTRL-C (copy) and CTRL-V (paste) work natively in the table. The table is pre-populated with some random data for demonstration purposes. Once your data is pasted in, click on the Plots tab and choose the Bland-Altman or Scatter Plot.

Example

Here is an image generated with the Bahar Shiny app using method comparison data obtained from St. Paul's Hospital Laboratory in migrating from Siemens Immulite 2000 XPi to Roche Cobas e601 for Calcitonin determination. Don't worry, we did more than 33 comparison–I am just showing the low end.

Try adjusting some of the plot parameters. The figures will update in real time. Thanks to Burak and Ayse Bahar for your work!

(Dan's) Parting Thought

There are straight lines that matter a lot more than regression.

I will make justice the measuring line and righteousness the plumb line
(Isa 28:17)

The Lab-R-torian

Month: August 2016

Make Bland Altman Plots with Marginal Histograms using ggExtra

The Problem

How To

Kickin' it Old School

Conclusion

A Shiny App for Passing Bablok and Deming Regression

Background

Update

Shiny App

Example