---
title: "A Shifting Landscape for Laboratory Science"
author: "John Maindonald, MSI, Australian National University"
date: "10 May 2016"
output: html_document
subtitle: "Reproducibility and Related Issues"
bibliography: osc1.bib
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, fig.retina=NULL)
```

##

There are increasing calls, in Nature and elsewhere, 
for reform of major aspects of scientific processes --- 
to address concerns about reproducibility, to improve 
efficiency, and to leverage current technology more effectively. 
In some major areas, an unacceptably high proportion 
of published work is not reproducible.  Reasons for this, and
initiatives that are designed to address the problem, will be 
noted. A component of reproducibility is ensuring that data 
analysis details can be readily reproduced.  An ideal is to 
have as a primary record a document that combines text and code, 
and which can be processed to produce a final paper that has all
tables, graphs and analysis results. Individual memories, bottom
drawers, and stacks of files, are less than ideal for this
purpose. Data should be accessed from a properly maintained
database.  The benefits are wide-ranging --- for avoiding manual
errors, for communication of details of work within and outside 
the organization, and for preserving data and code for any 
needed revision or updating or future use.

## Creaking Scientific Processes

* 'Obvious' nonsense: the hurricanes saga  
    + The _Growth in Time of Debt_ saga  
* Scientific processes need to leverage new technology  
    + Not just the lab equipment, but processes 
* The bloggers are coming, already here!    
* Complexity is everywhere.  Manage it!
    + Transfer complexity to the outside world!  
* Demonstrated failures in reproducibility
    + Amgen: 6 'successes' from 53 'landmark' cancer studies.   
    NB: Issues with the studies that failed . . .  
    + Bayer: Drug studies also come off badly   
* Reproducible reporting

## US Hurricane Deaths vs $ Damage

> Jung, Shavitta, Viswanathana, and Hilbe (2013).
> "Female hurricanes are deadlier than male hurricanes".
> http://www.pnas.org/cgi/doi/10.1073/pnas.1402786111.  
> PNAS, May 2014

A 'female' name caused less concern than a 'male' name?

Damage data (in 2013 $) are from ICAT: "ICAT . . . was founded in 1998 to provide production, underwriting, and risk management services for insurance companies."  

Deaths were related to damage, in 2013 dollars, for a comparable
hurricane in 2013 (but my graph will use 2014 data)

## Explanatory Variables

- Damage variables (in either case, in 2013 or 2014 $)
    + Base Damage, i.e., Damage at the time
    + ICAT damage estimate, if storm were in 2013 or 2014  
- Gender; factor with levels `f`, `m`  
  _(Jung et al used a femaleness panel rating, on a 1 to 11 scale.)_
  
[Also available was Barometric pressure at (first) US landfall]
  
__Question__: Which makes more sense?

- Use as explanatory variable: Base damage or ICAT damage?
    + Jung et al: ICAT damage; next slide: Base damage
- With deaths on log scale,  
   measure damage on untransformed or log scale?

## 
```{r cf-graphs1, echo=FALSE, fig.width=3.75, fig.height=2.75, results='hide', out.width="780px"}
dset <- DAAG::hurricNamed
parset1 <- modifyList(lattice::simpleTheme(pch=16, lty=2,
                                col.points=c('red','blue'),
                                col.line=c('red','blue'),
                                alpha.line=0.2), 
                      list(fontsize=list(text=10, points=5),
                           clip=list(panel='off', strip='on')))
library(lattice)
atxlab <- c(1, 10, 100, 1000, 10000, 100000)
atlab <- c(1,10,100, 1000)
omitdset <- dset[c(13,84),]
gph1 <- xyplot(log(deaths+1) ~ BaseDam2014, data=dset, groups=mf,
       type=c('p'), auto.key=list(space="right"), xlim=c(0,120000),
       par.settings=parset1,
       xlab="Damage at time (US$000 in 2014)", ylab="Deaths",
       scales=list(x=list(at=c(0, 25000, 50000, 75000, 100000)),
         y=list(at=log(atlab+1), labels=paste(atlab)), tck=0.5),
       main=list("Regression on `BaseDam2014`!!", cex=0.8))
hD1.lm <- lm(log(deaths+1) ~ mf + BaseDam2014,
             data=dset[-c(13,84),])
b <- coef(hD1.lm)
gph1+latticeExtra::layer(panel.abline(c(b[1],b[3]), 
                                      pch=16, col='red', 
                                      lty=2, lwd=2, alpha=0.4)) +
  latticeExtra::layer(panel.abline(c(b[1]+b[2],b[3]), 
                                   pch=16, col='blue', 
                                   lty=2, lwd=2, alpha=0.4)) +
  latticeExtra::layer(ltext(BaseDam2014, log(deaths+1), cex=0.75, labels=Year, pos=4, col='gray40'), ltext(BaseDam2014, log(deaths+1), cex=0.75, labels=Name, pos=2, col='gray40'), data=omitdset)
```

## 
```{r cf-graphs2, echo=FALSE, fig.width=3.75, fig.height=2.75, results='hide', out.width="780px"}
parset2 <- modifyList(lattice::simpleTheme(pch=16, lty=2,
                                col.points=c('red','blue'),
                                col.line=c('red','blue'),
                                alpha.line=0.2), 
                      list(fontsize=list(text=10, points=5), clip=list(panel='off', strip='on')))
atxlab <- c(1, 100, 10000, 1000000)
gph2 <- xyplot(log(deaths+1) ~ log(BaseDam2014), data=dset, groups=mf,
       type=c('p'), auto.key=list(space="right"),
       par.settings=parset2, xlim=log(c(0.5,110000)),
       xlab="Damage at time (US$000 in 2014)", ylab="Deaths",
       scales=list(y=list(at=log(atlab+1), labels=paste(atlab)),
                   x=list(at=log(atxlab), labels=paste(atxlab)),
                   tck=0.5),
       main=list("Regress on `log(BaseDam2014)`", cex=0.8))
hD2.lm <- lm(log(deaths+1) ~ mf + log(BaseDam2014), 
             data=dset[-c(13,84),])
b <- coef(hD2.lm)
gph2+latticeExtra::layer(panel.abline(c(b[1],b[3]), 
                                      pch=16, col='red', 
                                      lty=2, lwd=2, alpha=0.4)) +
  latticeExtra::layer(panel.abline(c(b[1]+b[2],b[3]), 
                                   pch=16, col='blue', 
                                   lty=2, lwd=2, alpha=0.4))+
  latticeExtra::layer(ltext(log(BaseDam2014), log(deaths+1), cex=0.75, labels=Year, pos=4, col='gray40'), ltext(log(BaseDam2014), log(deaths+1), cex=0.75, labels=Name, pos=2, col='gray40'), data=omitdset)
```

## The _Growth in Time of Debt_ Saga

The <font color="purple">[_Herndon et al_]</font> debunking of Reinhart and Rogoff: "_Growth in Time of Debt_" (2010) created a huge stir.
   
"... [Identified errors] include spreadsheet [mistakes] ..., omission of available data, weighting, and transcription, ... [affecting] countries in the high public debt category.<br>   
... [The result is] a false image that high public debt ratios inevitably entail sharp declines in GDP growth. ... there is a wide range of GDP growth ... at every level of public debt among the 20 advanced economies that RR survey."  

<font size="5">http://nymag.com/daily/intelligencer/2013/04/grad-student-who-shook-global-austerity-movement.html</font>
  
[_Herndon et al_]:http://nymag.com/daily/intelligencer/2013/04/grad-student-who-shook-global-austerity-movement.html  

##
```{r debt2gdp, echo=FALSE, fig.width=5.5, fig.height=4.25, results='hide', message=FALSE, warning=FALSE, out.width="750px"}
library(ggplot2, quietly=TRUE)
library(mgcv)
load('RR.RData')
RR.gam <- gam(dRGDP ~ s(debtgdp, bs="cs"),data=RR)
res <- residuals(RR.gam)
m <- ggplot(RR, aes(x=debtgdp,y=dRGDP))
m <- m + geom_vline(xintercept=90,color='lightgray',size=1.0)
m <- m + geom_point(color='darkgray', size=0.75) + ylab("Real GDP Growth") + xlab("Public Debt/GDP Ratio") + scale_x_continuous(breaks=seq(0,240,30)) + theme_bw()
m <- m + geom_smooth(method=gam, color='black', size=0.75, formula= y ~ s(x, bs = "cs"))
m <- m + geom_segment(aes(x = 90, y = 2.2, xend = 245, yend = 2.2), linetype="dashed", size=0.25, alpha=0.25, color="sienna") + geom_segment(aes(x = 90, y = -0.1, xend = 245, yend = -0.1), linetype="dashed", size=0.25, alpha=0.25, color="sienna")
keep <- (res>0 & abs(res)>11.15)|(res<0 & abs(res)>10.15) |
  (RR[,'debtgdp'] > 180 & res<0)
RRkeep <- RR[keep,]
hj <- rep(-0.085, sum(keep))
onleft <- with(RRkeep,(Country %in% c("NZ","Ireland"))|
                 (debtgdp > 180))
hj[onleft] <- 1.075
txlab <- with(RRkeep, paste(iso2c, Year))
m1 <- m + geom_text(data = RRkeep, aes(x=debtgdp,y=dRGDP, label = txlab), size=2, hjust = hj)
tagged <- as.character(unique(RRkeep$Country))
tag2 <- RRkeep[match(tagged, RRkeep[, 'Country']), 'iso2c']
tag2 <- as.character(tag2)
leg <- cbind(tag2, tagged)[tag2!=tagged,]
colnames(leg) <- NULL
library(gridExtra)
print(m1+annotation_custom(tableGrob(leg,
                           theme=ttheme_default(base_size=6)),
                           xmin=180, xmax=250, 
                           ymin=14, ymax=25.5))
```

## Herndon et al's Conclusions

"... RR’s findings have served as an intellectual bulwark in support of austerity politics. The fact that RR’s findings are wrong should therefore lead us to reassess the austerity agenda itself in both Europe and the United States."  

"For econometricians a lesson from the problems in RR is the advantages of reproducible code relative to working spreadsheets. ...   
  
[O]ur simplified version of the spreadsheet and R code that reproduces RR and corrected results ... [is] on our website."  
  
<font color="blue">Herndon (28yr old Grad Student at UMass) is now famous!</font>

## Science With or Versus the Bloggers 

* If the media finds work interesting, beware!
    + Deficiencies may get unwelcome scrutiny  
* Formalize processes for post-publication review?
    + Make Publish and Perish a reality!
* The blogosphere can be an early warning mechanism
    + It can identify failures of scientific process
    + Would it catch a new Andrew Wakefield event?
* Scientific processes can/must use the blogosphere?
    + Sure, there are dangers, but . . .
    + Good things are happening:   
    Stack OverFlow; CrossValidated  

## Framework Issues --- c.f. Disinfestation

* What can insect physiologists tell us about insect death?
* What are the effects of temperature?  
    + On insect death?
    + On fruit damage?
* What is an appropriate measure of dose?
    + Does the concentration effect change with time?
    + How far can we trust the C-T product?  
* What is the form of the dose-response?
    + Much of the literature has it wrong
* Results must be reproducible across years

## Move complexity to the outside world

NB: Daniel Levitin: "The Organized Mind: Thinking  
    Straight in an Age of Information Overload"

* Data --- put complexity into well-managed databases
    + GeoNet, with access via QuakeSearch is a good model
    + Use for both internal and external access
* Laboratory processes
    + Monitor automatically wherever possible, . . .
* Complexity in working through the details of analyses
    + Marry together text and code in a "markup" version
    that can be processed to give the paper or report  
    + Serious users of results will often want the    
    filled-out details of the "markup" version.
    + This aids communication both within and outside.

## Reproducibility --- Evidence of Failures
* Ioannidis (2005)^1^ --- '... Most Published ... Findings Are False' 
    + This paper put reproducibility issues "on the map"
* Direct evidence --- results do not reproduce
    + Examples shortly, best evidence is in psychology
    + Most worrying evidence is in cancer studies
* Warning signals, from examination of papers
    + 'lack of +ve & -ve controls', faulty stats,   
  'inappropriate use of key reagents', 'failure to repeat', ...
* Results may be unreplicible
    + Key information may be omitted or wrong 

<div style="font-size:smaller">
1: @r19_ioannidis_2005, 'Why Most Published Research Findings Are False' </div>

## Direct Evidence

* Amgen: Reproduced 6 only of 53 'landmark' cancer studies.$^1$
    + @r2_begley_2013 notes issues with the studies that failed
* Bayer: Main results from 19 of 65 'seminal' drug studies
    + NB, journal impact factor was not a good predictor!$^2$  
* fMRI studies: 57 of 134 papers (42%) had >= 1 case lacking
check on separate test image. 
Another 14%, unclear ...$^4$  
* The Reproducibility:Psychology Project (~40% replicated)
    + Summary of results: 28 Aug 2015 issue of Science

<div style="font-size:smaller">
1: @r23_begley_ellis_2012, 'Raise standards ...'; NB also @r2_begley_2013  
2: @r9_prinz_schlange_asadullah_2011, 'Believe it or not ... drug targets'  
3: @r8_kriegeskorte_simmons_bellgowan_baker_2009, '... dangers of double dipping'  
4: @r1_osc_EstRep2015, 'Estimating the reproducibility of psychological science'
</div>

<div class="notes">
Begley's 6 red flags  

- Blinding?
- Were basic expts repeated?
- Are all results given?  
  (Cross-check images)
- +ve & -ve controls?
- Were reagents validated?
- Flawed statistical analysis
</div>

## Collins & Tabak^1^ --- Factors include ... 

<div style="color:blue">
<ul style="list-style-type:none">
<li> poor training . . . in experimental design</li>  
<li> making provocative statements rather than presenting technical details</li> 
<li> Crucial experimental design elements that are too frequently ignored include blinding, randomization, replication, sample-size calculation and the effect of sex differences</li>  
<li> some scientists reputedly use a 'secret sauce' to make their experiments work — and withhold details . . . or describe them only vaguely . . .</li>
</ul>
</div>

Also: Deviations from stated protocol; errors in data; selective use of data;
selection effects

<div style="font-size:smaller">
1. @r15_collins_tabak_2014, '... NIH plans to enhance reproducibility'
</div>

## What are the issues?
* Faulty design/reporting/execution (Begley, Collins & Tabak) 
    + Repeating the same mistakes will not help
    + Some results are in principle unreplicable
* Selection effects, where mostly there is no effect
    + More enlightened use of p-values will help
    + More crucially, the total process is not transparent
    
We do not know, certainly not with any certainty, where the balance 
lies between these two different types of issue.  In this area, science lacks a scientific understanding of its own processes.

## Aside: Fisher on _p_-values

<p style="color:blue">If one in twenty does not seem high enough 
odds, we may, if we prefer it, draw the line at one in fifty 
(the 2 per cent point), or one in a hundred (the 1 per cent point). 
... <strong>A scientific fact should be regarded as experimentally 
established only if a properly designed experiment rarely 
fails to give this level (0.05 or 0.02) of significance.</strong>^1^
</p>

<p style="color:blue"> ... we may say that a phenomenon is 
experimentally demonstrable when we know how to conduct 
an experiment which will rarely fail to give us 
a statistically significant result.$^2$</p>  
<br>

<div style="font-size:smaller">
1: @r17_fisher_1926  
2: @r18_fisher_1937
</div>

## Science Community Responses

* The most effective responses have come from psychology.
* The next slide will summarize one major study.
* Studies in other areas are under way
    + Wet lab studies are harder to do
    + They interfere with regular work programs
    + Greater defensiveness than among psychologists.

##
```{r repPsych, echo=FALSE, message=FALSE, fig.width=6.5, fig.height=6.5, out.width="600px"}
df <- read.csv('repPsych.csv')
ord <- order(df$T_r.O)
df <- df[ord,]

library(mgcv, quietly=TRUE)
library(latticeExtra, quietly=TRUE)
gph <- xyplot(T_r.R ~ T_r.O, data=df, groups=Discipline,
     xlab = list("Effect size r (original)", cex=1.25), alpha=0.75,
     ylab = list("Effect size r (replication)", cex=1.25),
     xlim = c(-.2, 1), ylim = c(-.3, 1),
     xaxs="i", yaxs="i", aspect = 1, 
     auto.key=list(text=c(levels(df$Discipline),levels(df$Journal)),
                   points=TRUE, space="right"
                   ),
     scales=list(cex=1.0, tck=0.5),
     par.settings=simpleTheme(pch=c(16,16,c(3,5,4)), cex=1.5,
                              col=c(trellis.par.get()$superpose.symbol$col[1:3], rep("black",2))),
     main="Psychology: Open Science Collaboration Results"
)
gph2 <- xyplot(T_r.R ~ T_r.O, data=df, groups=Journal,
     xlim = c(-.2, 1), ylim = c(-.3, 1),
     xaxs="i", yaxs="i", 
     par.settings=simpleTheme(pch=c(3,5,4), col="white", lwd=2, cex=1.25)
)
# Fitted curve
RO.gam <- gam(T_r.R ~ Discipline+s(T_r.O, by=Discipline), 
              family=scat, data=df)
df[, c('fit','se')] <- predict(RO.gam, se=TRUE)
df <- within(df, {lower=fit-se
                  upper=fit+se})
gph+as.layer(gph2)+
  layer(data=subset(df, Discipline=="Social"), 
          col2 <- trellis.par.get()$superpose.symbol$col[2],
          llines(T_r.O, fit, col=col2),
          llines(T_r.O, lower, col=col2, lty=2),
          llines(T_r.O, upper, col=col2, lty=2),
        panel.abline(h=0, col="gray"),
        panel.abline(v=0, col="gray"),
        panel.abline(a=0, b=1, col="gray"))
```

## Is the criticism overblown?

* _The scientific process does finally identify the chaff_  
    _This contrasts with, e.g., alternative medicine._
    
Sure, but the process is far too protracted & tortuous.    
Too many methodological failures go undetected.  
Rewards systems encourage work that is poor quality.

* _One can never achieve 100% certainty_  

Sure, but we can do a lot better than at present.
Science should not be ignorant of its own processes.

* _Present processes are pretty much OK_ (implied, not said)

We have the technology needed to do a vastly better job, but are not using it!

## Journal & refereeing failures

* Referees & readers do not have the information needed
    + Know exactly what was done; check data, code
* Referees are at or beyond the limits of their expertise
    + Statistical analysis is an especial difficulty
* The system is not making good use of modern technology
    + plus, it interacts with rewards systems in malign ways
* Savvy critics are a huge untapped resource that is wasted
    + Other experts, in the area or in relevant areas
* Put papers out for comment pre-publication.
    
__Open Science__'s response: make all processes transparent 
    
## Commentary in _Science_ (June 26 2015)
### 1. Self-correction in Science at work^1^

* Publish replications (PPS now has a section for this)
* Highlight & reward completeness of information 
* Encourage publishing well (not often), ...
* Create a culture that is willing to admit mistakes, ...
* School scientists in research ethics

## Communication and Storage

* Spoken word 
* Writing
* Printing (as it relates to reporting & publishing)
    + For scientific work, a fairly complete record
    + Nowadays, typically, a very incomplete record
* Computer-based systems and the world wide web
    + Rethink (NOT adapt) paper-based systems
    + Old excuses for the deficiencies of paper  
    based sytems no longer wash.
    
Note the importance of organizational memory.

## Promiscuous Publication?

<p style="color:blue">
"Now we are witnessing the transition to yet another scholarly communication system — one that will harness the technology of the Web to vastly improve dissemination.  ... The Web opens the workshop windows to disseminate scholarship as it happens, erasing the artificial distinction between process and product. ...</p>

<p style="color:blue">
Today's publication silos will be replaced by a set of decentralized, interoperable services that are built on a core infrastructure of open data and evolving standards — like the Web itself .... This 'decoupled journal' publishes promiscuously, then subjects products to rigorous review through the aggregated judgements of expert communities, supporting both rapid, fine-grained filtering and consistent, meaningful evaluation."
</p>

<div style="font-size:smaller">
Jason Priem:
Nature 495, 437–440 (28 March 2013) doi:10.1038/495437a
</div>

## An Open Source Model for Science

* Open Source Malaria --- think “Linux for Malaria Research"^1^  
  This follows a successful Schistosomiasis project.
  
* The Validation Science Exchange's Reproducibility Initiative^2^  
* Cancer Studies --- 50 "most impactful" from 2010-2012^3^
<br>

<br>

<div style="font-size:smaller">
1: @r32_todd_others_2015 (Matt Todd & others), @r33_osbr_2015  
2: @r34_iorns_others_2015 (Iorns & others)  
3: @r25_errington_iorns_gunn_tan_lomax_nosek_2014; @r5_kaiser_2015, in June 26 2015 _Science_    
</div>

## Hard to Reproduce Reporting

Steps in preparing a report or paper include:

* Microsoft Word or LaTeX or . . . 
* Code files: one for each table, graph, . . .
* Data files
* Requires work to adapt code to new data, or . . .
* Manual steps required are prone to induce errors
* Difficult to retrace steps a year or more later.

## Reproducible Reporting

* One file for text, code and data or data access
* To make revisions, just revise and re-process
* Easy for a later worker to check code & re-use
* Allows a consultant to know just what was done
* Ideal setup: 
    + Code in the file accesses packaged code
    + Check help files to know what code does
    + Data is likewise from a 'packaged' source  
    (Such sources are called 'databases')
* Both R and Python can be used, others?

## Reporting with R Markdown

* Very simple to use with RStudio
    + 10 minutes tuition can get one started
* For finer control, can use HTML markup
* Can output to HTML or PDF or Word or . . .

## Slides

* Slides for this talk (pdf + R Markdown sources) will be posted at:
http://maths-people.anu.edu.au/~johnm/stats-issues/
* The technology use for the R Markdown sources for this talk 
can also be used for reports and papers
    + Simple text markup: use R Markdown
    + Sophisticated: Use Sweave = LaTeX with markup
    + Further possibility: HTML with markup
    + At least one CRI is using Sweave for reports

## References


