Revolutions

NYT on the importance of reproducible research

Yesterday's New York Times includes a great article on the failure of some genetic tests for cancer detection, and the flaws in the research that led to them. The article features quotes from Keith Baggerly of MD Anderson Cancer Center, and includes a photo of him and colleague Kevin Coombes in front of a page of R code:

The article highlights the importance of reproducible research: unless others can have access to the data and code that backs up the statistical conclusions in the paper, many such errors are likely to continue to go undiscovered. I saw Keith give an amazing talk about reproducible research and data forensics at the BioConductor conference a couple of years ago (I'm pretty sure that's one of the slides from the talk in the image above):

The talk concerned a published article where the results seemed odd. Baggerly tried in vain to get hold of the source data to reproduce the analysis, but the article authors didn't cooperate. So in an amazing feat of data forensics, he managed to recreate the data by matching public sources to measurements from the printed graphs, and figured out that there were gross data errors in the article: labels transposed, data duplicated, that kind of thing. The conclusions were completely bunk, but the journal refused to print a correction, despite the fact that it meant actual patients were being trialled on inappropriate drugs.

I'm glad to see this important issue is getting some wider media attention -- check out the Times article at the link below for the story.

New York Times: How Bright Promise in Cancer Testing Fell Apart

Five things Biologists should know about Statistics

In a thoughtful blog post, Bioinformatician Ewan Birney (Head of Nucleotide Data at the European Bioinformatics Institute) talks about the importance of Statistics to biologists:

Biology is really about stats. Indeed, the foundation of much of frequentist statistics - RA Fisher and colleagues - were totally motivated by biological problems.

He also cites the "Five statistical things I wished I had been taught 20 years ago". In order, they are:

1. Non-parametric statistics
2. R
3. The problem of multiple testing
4. The relationship between P-value, effect size, and sample size, and
5. Linear models and PCA

Read Ewan's full post at the link below for his reasons why every biologist should learn about these five statistical things.

Bioinformatician at Large: Five statistical things I wished I had been taught 20 years ago (via @paulblaser)

 

R in the Bioinformatics Knowledgeblog

The Knowledge Blog progect is a new, light-weight way of publishing scientific, academic and technical knowledge on the web, across several scientific disciplines. One such discipline is bioinformatics, and the Bioinformatics Knowledgeblog contains useful scientific reference material for bioinformatics, including several resources for R users. There you'll find an R Guide for Complete Beginners, a tutorial on analyzing microarray data with Bioconductor, a guide to making volcano plots for microarray data, and much more. Check out all the content at the link below.

The BioInformatics Knowledgeblog: About (via Simon Cockell)

Parallel Computing with R for Life Sciences

I hadn't heard of the CloudAsia 2010 conference before, but from the programme the workshop Master Class on HPC Application For Life Sciences looked like it was interesting. One workshop session in particular caught my eye: Practical Parallel Computing in R by Xie Chao and Tan Tin Wee (from the National University of Singapore). The workshop notes (PDF) provide a practical and concise overview of the parallel programming options in R, including several examples using our own foreach package under the heading "One Ring to Rule Them All":

The foreach and iterators packages created by REvolution Computing [now Revolution Analytics -- ed.] provide us a convenient framework for parallel computing in R. With these two packages, you only need to write one version of your code for all parallel backends.

The document gives examples of using several parallel backends with foreach, including doMC, doSNOW and doMPI. 

By the way, If you're planning to try out foreach yourself, don't forget that a new parallel backend for Windows is now available: doSMP, included with all editions of Revolution R. doSMP is an open-source R package, we're haven't yet submitted it to CRAN. But if you want to get in early and try it with R 2.11, Tal Galili has compiled it and made a Windows binary compatible with R 2.11 at his R Statistics blog (thanks, Tal!).

CloudAsia 2010: Master Class on HPC Applications For Life Sciences

Clinical Reporting with R

One of the main goals of analyzing clinical data is to produce a report. (What, you thought it was to make the world a better place?) The R Project has, of course, all the tools you need to perform the statistical analysis, calculate the tables of results, and present conclusions graphically. But how can you assemble all of that into a report that someone can, you know, read?

You could go the cut-and-paste route: write the text in Word, export the data from R to format the tables in Excel, dress up the saved charts in Photoshop. But that's a complex, manual process, and manual processes can introduce errors. Worse yet, if the data ever changes, you've got to go through the whole process again to update the report. That means no interim reports, and conversely, a big barrier to correcting the data and the report after it's published.

Vanderbilt Biostatistics professor Frank Harrell has a different solution: the rreport package for R. (See an overview slide deck here.) It's designed to produce statistical reports for clinical trials, and is especially useful for producing interim reports for data monitoring committees (DMCs). You can use it to create a complete report document, fully automating the process of generating tables from your R analyses like this:

and integrating R graphics into the document, like this:

The rreport system is an example of literate programming: the tables and reports are interwoven into the narrative text in a single source document in the open-source LaTeX typesetting language, and the entire report can be redone, with all results recalculated and all charts regenerated from the source data, in a single step. Using LaTeX does take some getting used to -- it's not a WYSIWYG environment like Word -- but does make for very attractive reports and the ability to easily typeset the Greek-laden mathematical equations so prevalent in clinical trial reports. You can download rreport from the Vanderbilt CVS repository at the link below.

Department of Biostatistics, Vanderbilt University: rreport package

Using R to analyze milk yields

Gregor Gorjanc has written a package for R to fit lactation curves to milk yield data. The package isn't yet on CRAN (it's not yet quite complete), but it's available for download from Gregor's website. Animal science isn't my specialty, so I'd never heard of Guo-Swalve or the four other lactation functions the package can fit. But I did think these graphs of the curves fitted to goat-milk data were pretty interesting (click for large version):

I'm a city boy, so for me, milk comes out of a carton and it always tastes the same. Never knew there was so much variation in milk fat, protein, and lactose over time.

Gregor Gorjanc: Fitting lactation curves/functions in R

Visualizing Hierarchical Cluster Models

The blog Getting Genetics Done features tips for analyzing genomic data, such as this recent post on hierarchical clustering with R using the hclust function. One of the interesting wrinkles with genomic data is that there will generally be many clusters, in which case the traditional tree representation is likely to be an unreadable "fuzzball":

The post gives some useful tips on how to export the tree data from R and use specialized software to create alternate visualizations, like this circular cladogram:

Getting Genetics Done: Hierarchical Clustering in R

Because it's Friday: Unfunny Statistics Jokes

Reader DE points me to this joke, given by Oxford professor Peter Donnelly in a 2006 TED talk:

Q: How do you tell the introverted statistician from an extroverted statistician?
A: The extroverted statistician is the one staring at the other person’s shoes.

I note however from his bio that Peter Donnelly is listed as a mathematician (despite his work in statistical genetics). I recall telling this same joke while I was working department shared by both mathematicians and statisticians, except the subject of the joke was mathematicians. Everyone has their stereotypes, I guess.

Despite the bad joke, the talk itself is well worth a look. It's an engaging look at common fallacies in statistics and probability, in particular the problem of inverting conditional probabilities (aka the Prosecutor's Fallacy). The example was one I hadn't heard of before, where Sally Clark, a mother in the UK, was convicted of murder after losing two children to cot death, and a paediatrician citied in expert testimony that the chances of this happening by accident was 1 in 73 million. Not only does this estimate assume independence of the two deaths, it's irrelevant: it tells us nothing about the relative likelihood of murder versus tragic accident. Sally Clark was later exonerated after two appeals, but adding further to the tragedy died suddenly in 2007.

On a lighter note, there was also an interesting stats puzzle I hadn't seen before. Imagine tossing a coin repeatedly until a certain sequence of three tosses is seen, say: HTH or HTT. On average, do you think it takes more or fewer tosses to see HTH or HTT, or is the expected count the same? Guess first, and then check the video below for the answer.

The impatient can skip ahead to 5:20 for the answer. (For the record, I guessed wrong.) This kind of problem has important applications in genetics. Finding a specific sequence of ACTG in a genome of millions of base-pairs is a harder problem than tossing coins, but the R package Biostrings solves it very efficiently. 

TED: Peter Donnelly shows how stats fool juries

BioConductor 2009: Parallel Programming for Gene Sequences

I had a really great time at the BioC2009 conference in Seattle this week. The conference is dedicated to the analysis of gene expression data using the the BioConductor packages for R. It's an area I haven't really looked at for quite a while, but boy, has everything changed. Three years ago the hot topic was microarray analysis; now, that's old hat. 

The big thing this year is sequence analysis: scientists can now look inside a cell and identify which segments of DNA are being expressed (and at what rate). This provides a direct link between experimental factors like drugs, diseases, and risk factors and the actual behavior of specific parts of the body in terms of the proteins our cells create via genes in our DNA. Not only does this create a vast volume of data about known genes, but has spawned a whole new data-identification problem of finding previously unknown (but observed, via expression) genes in the reference genome. Suddenly, the ability to find a short sequence (a gene) in a VERY long string (the genome) using fuzzy-matching techniques (genes aren't always the exact same ACTG sequence) is an important computational problem.

I think that may be why there was such interest in my Parallel Programming in R lab session. It was a 2-hour lecture and practical session. Much of the material I've covered here before, such as using iterators and foreach to parallelize loops in R on multi-core systems or on clusters. It does seem to be the case that a lot of typical calculations in BioConductor are embarrassingly parallel, and it's easy to convert a for loop or an lapply loop to foreach and get some significant speedups. (Others including Michael Knudsen and William Webber have found this to be the case, too.) There was also some new material on using video cards (GPUs) for high-performance multi-threaded calculations, which I'll probably write up as a blog post some day. You can download the slides (2Mb PDF) and the R script file I used for the session.

All in all, it was a great event (despite the record heat in Seattle) and it's clear that BioConductor has a thriving community about it. Robert Gentleman announced at the conference he's stepping down as the head of the project to take up a new role in the Bay Area, but with Martin Morgan taking the reins of the team in Seattle and with the development community worldwide, the BioConductor will continue to be at the leading edge of gene expression data analysis for the foreseeable future.

Swine Flu: Epidemiology scores a point over politics

Vice-President Biden got into a bit of trouble today, saying he wouldn't advise his family to travel by plane or by subway to limit the chances of their being infected with swine flu. Some, like the United States Travel Association trade group, weren't happy. But the science may be on Biden's side, says Ezra Klein:

Biden may be right to induce a bit of panic. The United States Travel Association won't think so. It's their profit stream on the line, after all. But epidemiologists are probably quietly relieved by the Vice President's comments.

The danger of swine flu lies in whether or not it becomes a full-blown epidemic, and that in turn depends on the base reproduction rate: the number of secondary infections generated by a primary infection. As shown in this chart (which unless I am very much mistaken was generated with R), anything we can do to reduce the reproduction rate R₀ can dramatically reduce the infection attack rate.

So score one for science.

Ezra Klein: Is Joe Biden right about swine flu?