ecoinfoPlots.Rmd

---
title: "EVOSTC data collection paper - Data Visualizations"
author: "Jessica Couture"
date: "10/31/2017"
output:
  pdf_document:
    toc: yes
    toc_depth: '3'
  html_document:
    toc: yes
    toc_depth: 3
---

#Background

The data and plots produced here contribute to a publication ("Funder imposed data publication requirements seldom inspire data sharing") to report results of a two year data archiving effort by a small group of researchers and students at the National Center for Ecological Analysis and Synthesis at UC Santa Barbara. The [Exxon Valdez Oil Spill Trustee Council](http://www.evostc.state.ak.us) (EVOSTC) was formed following the Exon Valdez oil spill in Alaska in 1989. Since then, the EVOSTC has funded hundreds of projects and in 2012 we began a project to recover and archive the data collected in these EVOSTC funded projects.

For this paper we ask 5 main questions about the data collected from the Exxon Valdez Oil Spill Trustee Council funded projects:

1. Twenty-five years after the EVOS, for how many projects funded by EVOSTC can we collect data?
2. Are there certain **research fields** that are more likely to make data available than others?
3. Are there certain **sectors** that are more likely to make data available than others?
4. Is the availability of data correlated to how old the data are? (temporal relationships?)
5.Why did people refuse to share their data?

***

```{r data prep,echo=FALSE,message=FALSE}

library(knitr)
library(dplyr)
library(readr)
library(ggplot2)
library(reshape2)
library(grid)
library(gridExtra)
library(gtable)
library(cowplot)
library(party)

rslt=read.csv('data/evosArchResults.csv',header=T,stringsAsFactors=F)
rslt[rslt$Status=='published','Status']<-'Published'

rslt$succ=ifelse(rslt$Status %in% c('Published','SentData','Revised'),0,1) # success=0, not successful=1....why?
rslt$statSucc=ifelse(rslt$Status %in% c('Published','SentData','Revised'),1,2) # 1 = success, 2 = not successful

agDes=read.csv('data/ecoInfoAgencyDesig.csv',header=T,stringsAsFactors=T,col.names=c('agency','agGrp','agSubGrp'))
agDes$agGrp=gsub('nonProf','private',agDes$agGrp)

rslt2=merge(rslt,agDes,all.x=T)
rslt2=rslt2 %>%
  filter(!end>2010) 

rslt2$dataType<-gsub("biological-",'',rslt2$dataType)
```

# MS Plots:

## 1. Project Status Reporting
> Twenty-five years after the EVOS, for how many projects funded by EVOSTC can we collect data?</h4>

__Figure 1:__ Final status of all projects requested. Black are bars are projects for which data were successfully acquired, grey bars represent projects for which no data were acquired.

_NOTES: For this first plot, we decided to use the original plot generated from redmine that simply displays the number of datasets in each staus at the end to the 2 year data collection effort. Still having issues getting the arrow in the x-axis title without the period or other character after it._

```{r status plot, fig.width=10,echo=FALSE}

overall=table(rslt2$Status)
overallLong=melt(overall)
colnames(overallLong)=c('Status','nDatasets')
overall2=overallLong %>%
  mutate(prop=nDatasets/sum(nDatasets))

possible=overallLong %>%
  filter(!Status=='Unrecoverable') %>%
  mutate(prop=nDatasets/sum(nDatasets))

statusOrd=c('Unrecoverable','Emailed','Replied','SentData','Revised','Published')
#rslt2$Status=as.factor(rslt2$Status,levels=statusOrd,ordered=T)
statGrey=c('Unrecoverable'='gray67',
             'Emailed'='gray67',
             'Replied'='gray67',
             'SentData'='black',
             'Revised'='black',
             'Published'='black')

stSuccCols=c("2"='gray67',"1"="black")

statPlot=ggplot()+
  geom_bar(data=rslt2,
                 aes(x=factor(Status,levels=statusOrd,ordered=T),fill=factor(statSucc)))+
  scale_fill_manual(name="Data obtained?",
                    breaks=c("1","2"),
                    values=stSuccCols,
                    labels=c("Data acquired", "No data acquired"))+
  theme(legend.title=element_blank())+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())+
  xlab("Recovery stage")+
  ylab('Number of projects\n')+
  theme(axis.text=element_text(size=14),
        title=element_text(size=16,face="bold"),
        legend.text = element_text(size = 14, face = "bold"))+
  theme(axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black"))

statPlot

```

***

## 2. Research field
> Are there certain **research fields** that are more likely to make data available than others?
 
__Figure 2:__ Percent success of projects by research field. The top plot is total number of projects funded for each field, the bottom plot shows the percent success for the given field.

_NOTES: We stack nProjects/percent graphs to display variation in nProjects between groups as well as variation in relative responses for each group._

```{r research field success plots,fig.height=9, echo=FALSE}

catYN=as.data.frame(table(rslt2$dataType,rslt2$succ))
colnames(catYN)=c('category','succ','nDatasets')

fldYN=catYN %>%
  group_by(category) %>%
  mutate(grpTots=sum(nDatasets)) %>%
  mutate(prop=nDatasets/grpTots)

fldYN$category=factor(fldYN$category,levels=c("oil","physical","habitat",'plankton','benthicInverts','fish','birds','mammals',"social","modeling"))

##### all fields

fieldNDs=ggplot(fldYN[fldYN$succ==0,],aes(x=category,y=grpTots))+
  geom_bar(stat = "identity",fill="grey")+
  theme(axis.line=element_line('black'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())+
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank())+
  ylim(0,100)+
  ylab('Number of projects\n')+
  theme(axis.text=element_text(size=14),
        title=element_text(size=16,face="bold"))+
  theme(axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black"))
  
fieldPerc=ggplot(fldYN[fldYN$succ==0,],aes(x=category,y=(prop*100)))+
  geom_bar(stat = "identity",fill="black")+
  scale_x_discrete(labels=c("oil","physical","habitat",'plankton','benthicInverts','fish','birds','mammals',"social","modeling"))+
  theme(axis.line=element_line('black'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.text.x = element_text(angle=45,size = 14,vjust = 1,hjust = 1))+
  xlab('Research field')+
  ylab('% recovered\n')+
  ylim(0,63)+
  theme(title=element_text(size=18,face="bold"),
        axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black"))


plot_grid(fieldNDs,fieldPerc, ncol = 1, rel_heights=c(1,1.5),align="v")

```

***

## 3. Classification tree
> How do the important characteristics influence the output (success)? 

__Figure 3:__ Classification tree of which data variables and characterizations predict successful and non-successful data acquisition. Data field was the only variable that could be used to predict data availability. 


```{r class tree bio, message=F, warning=FALSE}
rslt2$dataType<-factor(rslt2$dataType)
partreeBio<-ctree(succ~agSubGrp+end+dataType,data=rslt2)

plot(partreeBio)
```

***

## 4. Reasons for not sharing

> Why don't people share data?

__Figure 4:__ Reasons given for not providing data. Communication loss and lack of contact information were the main reasons data were not obtained.

_NOTES: Here we flip the axes to separate this analysis from the others. These data are not a sub-grouping of our results like the others, but rather a deeper analysis of only the unrecovered datasets. Also, unlike the other data above, we did not run stats here._

Data details:
Here I used all of the data that were not collected and sorted them into categories based on the redmine notes we took:

* All datasets labeled "emailed" were grouped with a other datasets that were labeled "unrecoverable" because no contact information could be found.
* All datasets labeled "replied" we put into their own category called "communication lost".
* The "data lost" category includes all datasets in which someone confirmed that the data no longer existed either due to damage, non-persistent formats (not including printed/non-digital), etc.
* The other are self-explanatory, but represent confirmed responses as to one of these reasons.
* I tried to separate out the datasets that were not recovered but "should be on the CD" to sort these into their own category, but we acked sufficient notes to confidently add this grouping.

```{r Unrecoverable reasons, echo=F}
reas=rslt2 %>%
  filter(Status=='Unrecoverable') %>%
  filter(!reason=="no data generated") %>%
  filter(!reason=="data uvailable") %>%
  select(Status,reason)
reasTab=as.data.frame(table(reas$reason))
colnames(reasTab)=c('reason','nDatasets')

## add in status = 'emailed' to 'no contact info' & status = "replied" to unwilling to share? ...discuss with group

reasTab[reasTab$reason=='no contact info',"nDatasets"] <- (reasTab[reasTab$reason=='no contact info',"nDatasets"] + length(rslt2[rslt2$Status=='Emailed','Status'])) 

levels(reasTab$reason)=c(levels(reasTab$reason),'communication lost')

reasTab[6,]<-c('communication lost', length(rslt2[rslt2$Status=='Replied','Status'])+1)# +1 for one of the 'data unavailable' datasets, add the other to "data lost"
reasTab[1,2]<-as.numeric(as.character(reasTab[reasTab$reason=="data lost","nDatasets"]))+1

reasTab<-arrange(reasTab,as.numeric(nDatasets))
reasTab$reason<-factor(reasTab$reason,levels = reasTab$reason,ordered = T)

reasTot=sum(as.numeric(reasTab$nDatasets))
reasTab$perc=as.numeric(reasTab$nDatasets)/reasTot

whyNot<-ggplot(reasTab,aes(x=reason,y=as.numeric(perc)*100))+
  geom_bar(stat = "identity")+
  coord_flip()+
  #scale_x_discrete(breaks = seq(0,60,by=10), labels=seq(0,60,by=10)) +
  theme(axis.line=element_line('black'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())+
  theme_light()+
  ylim(0,60)+
  xlab('')+
  ylab('\n% of unrecovered projects')+
  #ggtitle('Why are data not recovered?')+
  theme(axis.text=element_text(size=14),
        title=element_text(size=16,face="bold"))+
  theme(axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black")) ## *** add each n as text above bar
whyNot
```

***

\pagebreak

# Supplementary figures

## S1 Figure. Awardee agency sector
> Are there certain **sectors** that are more likely to make data available than others?
(based on PI's affiliation)

__SuppFigure 1:__ Percent success of projects by agency sector. The top plot is total number of projects funded for each sector, the bottom plot shows the percent success for the given sector.

_NOTES: We stack nProjects/percent graphs to display variation in nProjects between groups as well as variation in relative responses for each group. These data show interesting opposing trends in nProjects and %recovery_

```{r, subSector stacked, fig.height=10,echo=FALSE}
secSubYN=as.data.frame(table(rslt2$agSubGrp,rslt2$succ))
colnames(secSubYN)=c('subSector','succ','nDatasets')

secSubYN2=secSubYN %>%
  group_by(subSector) %>%
  mutate(grpTots=sum(nDatasets)) %>%
  mutate(prop=nDatasets/grpTots)
secSubYN2$subSector=factor(secSubYN2$subSector,levels=c("gov_fed","gov_state","academia","nonProf","private","tribe"),ordered = T)

subSecNds=ggplot(secSubYN2,aes(x=subSector,y=nDatasets))+
  geom_bar(stat = "identity",fill="grey")+
  theme(axis.line=element_line('black'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())+
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank())+
  ylab('Number of projects\n')+
  theme(axis.text=element_text(size=14),
        title=element_text(size=16,face="bold"))+
  theme(axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black"))

subSecPerc=ggplot(secSubYN2[secSubYN2$succ==0,],aes(x=subSector,y=as.numeric(prop*100)))+
  geom_bar(stat = "identity",fill="black")+
  scale_x_discrete(breaks=levels(secSubYN2$subSector),
                   labels=c("Gov - fed","Gov - State","academia","non-profit","for-profit","AK native"))+
  theme(axis.line=element_line('black'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.text.x = element_text(angle=45,size = 14,vjust = 1,hjust = 1))+
  xlab('\nSector')+
  ylab('% recovered\n')+
  ylim(0,63)+
  #ggtitle('Percent success by research sector type')+
  theme(axis.text=element_text(size=14),
        title=element_text(size=16,face="bold"))+
  theme(axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black"))
  # annotate("text", x = 1, y=secSubYN2$prop[[2]]*100+3, label = as.character(secSubYN2[secSubYN2$subSector=="gov_fed" & secSubYN2$succ==0,"grpTots"]))+
  # annotate("text", x = 2, y=secSubYN2$prop[[3]]*100+3, label = as.character(secSubYN2[secSubYN2$subSector=="gov_state" & secSubYN2$succ==0,"grpTots"]))+
  # annotate("text", x = 3, y=secSubYN2$prop[[1]]*100+3, label = as.character(secSubYN2[secSubYN2$subSector=="academia" & secSubYN2$succ==0,"grpTots"]))+
  # annotate("text", x = 4, y=secSubYN2$prop[[4]]*100+3, label = as.character(secSubYN2[secSubYN2$subSector=="nonProf" & secSubYN2$succ==0,"grpTots"]))+
  # annotate("text", x = 5, y=secSubYN2$prop[[5]]*100+3, label = as.character(secSubYN2[secSubYN2$subSector=="private" & secSubYN2$succ==0,"grpTots"]))+
  # annotate("text", x = 6, y=secSubYN2$prop[[6]]*100+3, label = as.character(secSubYN2[secSubYN2$subSector=="tribe" & secSubYN2$succ==0,"grpTots"]))

# grid.arrange(subSecNds, subSecPerc, ncol = 1, heights=c(1,1.75))
plot_grid(subSecNds, subSecPerc, ncol = 1, rel_heights=c(1,1.5),align="v")

```

***

## S2 Figure. Temporal trends?
> Is the availability of data correlated to how old the data are?

__SuppFigure 2:__ Percent success of projects by age of data. Age is calculated based on number of years between the last year of EVOSTC funding and start of the archiving project (2012). The top plot is total number of projects that ended each year, the bottom plot shows the percent success for the given year.

_NOTES: We stack nProjects/percent graphs to display variation in nProjects between groups as well as variation in relative responses for each group._

_The x-axis might not be initially intuitive but we are representing years since a project ended in opposite-chronological order to show any effect of increasing age of a dataset. We chose this design to be able to compare to [Michener et al. 1997](http://onlinelibrary.wiley.com/doi/10.1890/1051-0761(1997)007[0330:NMFTES]2.0.CO;2/abstract;jsessionid=57D1E6F52B947ABD71ABEA56185DA7D3.f01t04) - fig 1._

```{r temporal - stacked plots,fig.height=10, echo=F, message=F, warning=F}

rslt2$succ<-ifelse(rslt2$Status %in% c('Unrecoverable','Emailed','Replied'),1,0)

endYr=rslt2 %>%
  filter(!end<1989) %>%
  filter(!end>2010) %>%
  select(succ,end)

endYrs=as.data.frame(table(endYr$end,endYr$succ))
colnames(endYrs)=c('year','succ','freq')
endYrs<- mutate(endYrs,since=(2012-as.numeric(paste(endYrs$year))))

yrsPropEnd<-endYrs %>%
  group_by(as.character(year)) %>%
  mutate(yrSum=sum(freq)) %>%
  mutate(prop=freq/yrSum)

yrsNds=ggplot(yrsPropEnd[yrsPropEnd$succ==0,],aes(x=since,y=yrSum)) +
  geom_bar( stat="identity",fill="grey")+
  theme(axis.line=element_line('black'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.title.y=element_text(color="black"),
        axis.title.x=element_blank())+
  xlim(min(endYrs$since),max(endYrs$since))+
  #ylim(0,100)+ 
  #xlab('Years since funding award')+
  ylab('Number\nof projects')+
  theme(legend.position="none",
        axis.text=element_text(size=14),
        title=element_text(size=16,face="bold"),
        axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black"))


yrsPerc=ggplot(yrsPropEnd[yrsPropEnd$succ==0,],aes(x=since,y=as.numeric(prop*100)))+
  geom_bar(stat = "identity",fill="black")+
  xlim(min(endYrs$since),max(endYrs$since))+
  theme(axis.line=element_line('black'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())+
  xlab('Age of data (years)')+
  ylab('% recovered\n')+
  ylim(0,63)+
  theme(axis.text=element_text(size=14),
        title=element_text(size=16,face="bold"))+
  theme(axis.line.x = element_line(color="black"),
        axis.line.y = element_line(color="black"))

# grid.arrange(yrsNds, yrsPerc, ncol = 1, heights=c(1,1.5))
plot_grid(yrsNds, yrsPerc, ncol = 1, rel_heights=c(1,1.5),align="v")
```

***

## S3 Figure. Temporal trends in reasons?
> Are there temporal trends in the reasons given for not sharing data?

__SuppFigure 3:__ Counts of reasons for not sharing data based on the age of the data. 

_To test if there were temporal trends in reasons given for not sharing data, we plotted the frequency of each reason for not sharing based on the age of the data._

```{r reasons - temp, echo=F, message=F, warning=F}
reas2a=rslt2 %>%
  #filter(!reason=="no data generated") %>%
  #filter(!reason=="data uvailable") %>%
  filter(Status %in% c('Unrecoverable',"Emailed","Replied")) %>%
  select(Status,reason,end,agSubGrp,dataType)
reas2a[reas2a$Status=="Emailed","reason"]<-"no contact info"
reas2a[reas2a$Status=="Replied","reason"]<-"lost communication"

reas2<-filter(reas2a,!reason %in% c("no data generated","data uvailable"))

reasYr=as.data.frame(with(reas2,table(reason,end)))

DL = ggplot(data=reasYr[reasYr$reason=="data lost",],aes(as.numeric(as.character((end))),Freq)) +
  geom_line() +
  ylim(0,15) +
  xlim(1989,2010) +
  xlab("") +
  ylab("")+
  ggtitle("data lost")

NCI = ggplot(data=reasYr[reasYr$reason=="no contact info",],aes(as.numeric(as.character((end))),Freq)) +
  geom_line() +
  ylim(0,15) +
  xlim(1989,2010) +
  xlab("") +
  ylab("")+
  ggtitle("no contact info")

ND = ggplot(data=reasYr[reasYr$reason=="non digital",],aes(as.numeric(as.character((end))),Freq)) +
  geom_line() +
  ylim(0,15) +
  xlim(1989,2010) +
  xlab("") +
  ylab("")+
  ggtitle("non digital")

RF = ggplot(data=reasYr[reasYr$reason=="requested funding",],aes(as.numeric(as.character((end))),Freq)) +
  geom_line() +
  ylim(0,15) +
  xlim(1989,2010) +
  xlab("") +
  ylab("") +
  ggtitle("requested funding")

US = ggplot(data=reasYr[reasYr$reason=="unwilling to share",],aes(as.numeric(as.character((end))),Freq)) +
  geom_line() +
  ylim(0,15) +
  xlim(1989,2010) +
  ylab("")+
  xlab("") +
  ggtitle("unwilling to share")

LC = ggplot(data=reasYr[reasYr$reason=="lost communication",],aes(as.numeric(as.character((end))),Freq)) +
  geom_line() +
  ylim(0,15) +
  xlab("") +
  ylab("")+
  ggtitle("lost communication")

yAxLab <- ggdraw() + draw_label("Frequency", fontface='bold',angle = 90,hjust=0)
xAxLab <- ggdraw() + draw_label("Year", fontface='bold',vjust=0)

pYr<-plot_grid(DL,NCI,ND,RF,US,LC,labels="AUTO", ncol = 2,align="v")
xPlot<-plot_grid(pYr,xAxLab,ncol=1,rel_heights = c(1,0.05))
plot_grid(yAxLab,xPlot,ncol=2,rel_widths = c(0.05,1))
```