analysis-datasharing.Rmd

---
title: "Data Sharing Policies - Analysis"
author: "Analysis by Jessica Minnier, OHSU; github:jminnier"
date: "February 2, 2017"
output:
  html_document:
    code_folding: hide
    number_sections: yes
    theme: cerulean
    toc: yes
    toc_float: yes
---

```{r setup, include=FALSE}
list.of.packages <- c("broom","devtools","dplyr","ggplot2",
                      "tibble","tidyr","reshape2","stringr","knitr","pander","readxl")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
# install.packages(list.of.packages) # run this if some packages are too old and cause errors
# from http://stackoverflow.com/questions/4090169/elegant-way-to-check-for-missing-packages-and-install-them


library(ggplot2)
library(reshape2)
library(knitr)
library(pander)
library(readxl)
library(tidyverse)
library(tableone)
library(magrittr)
library(irr)
library(janitor)

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, tidy=TRUE, message = FALSE,fig.path = "figs/",fig.keep=TRUE)

```

# Data Processing

Here we read in the data and do some name and format manipulation.

Note, the abbreviation DSM denotes "Data Sharing Mark" and OAM denotes "Open Access Mark."

```{r data_process}
# read journal data
exceldat <- read_excel("../Data-Sharing-Policies_2017-01-25.xlsx", sheet = "Data",
                    na = "N/A")
jdata <- exceldat[-1,] # remove line of info
jdata <- jdata[-which(is.na(jdata$Journal)),] # remove NAs

# add labels for DSM
dsm_labels = c(
  "Required as condition of publication, barring exceptions",
  "Required but no explicit statement regarding effect on publication/editorial decisions",
  "Explicitly encouraged/addressed, but not required.",
  "Mentioned indirectly",
  "Only protein, proteomic, and/or genomic data sharing are addressed.",
  "No mention"
)

# Recommended Sharing Method
rsm_labels = c(
  "Public Online Repository", "Journal Hosted",
  "By Reader Request to Authors","Multiple methods equally recommended",
  "Unspecified"
)

journal_host_labels = c(
  "Journal will host regardless of size",
"Journal has data hosting file/s size limit","Unspecified"
)

# Fix typos:
# One line has an end line character when should be 'b'
jdata$`Copyright/Licensing Mark`[jdata$`Copyright/Licensing Mark`=="b\n"] = "b"
jdata$`Reproducibility Noted Mark`[jdata$`Reproducibility Noted Mark`=="a\n"] = "a"
# Space after an a
jdata$`Archival/Retention Mark`[jdata$`Archival/Retention Mark`=="a "] = "a"
jdata$`Protein, Proteomic, Genomic, or Microaray (Sequence or Structural) Data Sharing Addressed/Required with Deposit to Specific Data Banks`[jdata$`Protein, Proteomic, Genomic, or Microaray (Sequence or Structural) Data Sharing Addressed/Required with Deposit to Specific Data Banks`%in%c("a\n","a ")] = "a"
jdata$`Protein, Proteomic, Genomic, or Microaray (Sequence or Structural) Data Sharing Addressed/Required with Deposit to Specific Data Banks`[jdata$`Protein, Proteomic, Genomic, or Microaray (Sequence or Structural) Data Sharing Addressed/Required with Deposit to Specific Data Banks`%in%c("b\n","b ")] = "b"

# easier shorter names
jdata <- jdata%>%mutate(
  citable_2013 = `{2013} Citable Articles`,
  citable_2014 = `2014 Citable Items`,
  citable_2015 = `2015 Citable Articles`,
  dsm=as.numeric(`Data Sharing Mark`), 
  dsm_fac = factor(dsm,levels = 1:6,
                   labels=dsm_labels),
  dsm2 = 1*(dsm<3),
  dsm2_fac = factor(dsm2,levels=0:1,
                    labels=c("Not Required","Required")),
  dsm2_fac_flip = factor(dsm2,levels=c(1,0), # reverse levels for for plotting only
                    labels=c("Required","Not Required")),
  oam=as.numeric(`Open Access Mark`),
  oam_fac = factor(oam,levels=0:1,
                   labels=c("Subscription","Open Access")),
  if_2013=`2013 Impact Factor`,
  if_2014=`2014 Impact Factor`,
  if_2015 = `2015 Impact Factor`,
  total_cites_2013 = `2013 Total Cites`,
  total_cites_2014 = `2014 Total Cites`,
  total_cites_2015 = `2015 Total Cites`,
  omic_data_sharing_fac = factor(`Protein, Proteomic, Genomic, or Microaray (Sequence or Structural) Data Sharing Addressed/Required with Deposit to Specific Data Banks`,
                                 levels=c("a","b"),
                                 labels=c("Yes","No")),
  rec_pref_sharing_fac = factor(`Recommended/Preferred Sharing Mark`,levels=c("A","B","C","D","E"),
                   labels=rsm_labels),
  journal_hosts_fac = factor(`Size Guidelines if Journal Hosted Provided`,levels=c("a","b","c"),labels=journal_host_labels),
  copyright_licensing_fac = factor(`Copyright/Licensing Mark`,levels=c("a","b"),
                                   labels=c("Explicitly stated or mentioned","No Mention")),
  archival_retention_fac = factor(`Archival/Retention Mark`,levels=c("a","b"),
                                  labels=c("Explicitly stated","No Mention")),
  reproducibility_noted_fac = factor(`Reproducibility Noted Mark`,levels=c("a","b"),
                                  labels=c("Explicitly stated","No Mention"))
)

# create long data for plotting and dplyr
jdata_if_long  <- jdata%>%select(Journal,dsm:if_2015)%>%
  rename(`2013`=if_2013,`2014`=if_2014,`2015`=if_2015)%>%
  gather(year,impact_factor,`2013`:`2015`)
jdata_tc_long <- jdata%>%select(Journal,contains("total_cites"))%>%
  rename(`2013`=total_cites_2013,`2014`=total_cites_2014,`2015`=total_cites_2015)%>%
  gather(year,total_cites,`2013`:`2015`)
jdata_cit_long <- jdata%>%select(Journal,contains("citable_"))%>%
  rename(`2013`=citable_2013,`2014`=citable_2014,`2015`=citable_2015)%>%
  gather(year,citable_items,`2013`:`2015`)

jdata_long = left_join(jdata_if_long,jdata_tc_long)
jdata_long = left_join(jdata_long,jdata_cit_long)
```

# Data Summaries and Distributions

We have data from `r nrow(jdata)` journals, including impact factor (IF) and total number of citations in years 2013 and 2014.

Data Sharing Mark (DSM) has 6 categories, and we create a collapsed DSM variable that has two categories for "required" (DSM = 1 or 2) and "not required" (DSM = 3, 4, 5 or 6).

The labels for Data Sharing Mark (DSM) are as follows:

```{r dsmlabel}
kable(tibble("DSM Numeric Value"=1:6,"DSM Description Label"=dsm_labels))
```

## Summary of Continuous Data
```{r}
tmps = summary(jdata%>%select(if_2013:if_2015,
                              citable_2013:citable_2015,
                              total_cites_2013:total_cites_2015))
kable(tmps)
```


## TABLE METHODS X: Distribution of Impact Factors for Journal Titles

We include 2013 values in the manuscript and include 2014 values in the appendix.

This is a table categorizing journal impact factor (JIF) into ranges. The number of journals within each range and the percentage of the total number of journals is presented.

### 2013

```{r, results='asis'}
jdata$if_2013_cat = cut(jdata$if_2013,breaks = c(0,2,4,6,8,10,30,43),
                        include.lowest = TRUE,right = FALSE,
                        labels = c("<2","2-3.99","4-5.99","6-7.99","8-9.99","10-29.99","30-43")
)
jdata$if_2014_cat = cut(jdata$if_2014,breaks = c(0,2,4,6,8,10,30,43),
                        include.lowest = TRUE,right = FALSE,
                        c("<2","2-3.99","4-5.99","6-7.99","8-9.99","10-29.99","30-43")
)
jdata$if_2015_cat = cut(jdata$if_2015,breaks = c(0,2,4,6,8,10,30,43),
                        include.lowest = TRUE,right = FALSE,
                        c("<2","2-3.99","4-5.99","6-7.99","8-9.99","10-29.99","30-43")
)


jdata%>%group_by(if_2013_cat)%>%summarize("Number"=n())%>%mutate(
  "# (%)"=paste0(Number,"\t(",round(100*Number/nrow(jdata),1),"%)"))%>%select(-Number)%>%
  rename("IF 2013 Category"=if_2013_cat)%>%kable()
```

### 2014, 2015

```{r, results='asis'}
tmp = jdata%>%group_by(if_2014_cat)%>%summarize("Number"=n())%>%
  mutate(
  "2014 # (%)"=paste0(Number,"\t(",round(100*Number/nrow(jdata),1),"%)"))%>%select(-Number)%>%
  rename("IF Category"=if_2014_cat)
tmp2 = jdata%>%group_by(if_2015_cat)%>%summarize("Number"=n())%>%
  mutate(
  "2015 # (%)"=paste0(Number,"\t(",round(100*Number/nrow(jdata),1),"%)"))%>%select(-Number)%>%
  rename("IF Category"=if_2015_cat)
#replace NA
left_join(tmp,tmp2)%>%mutate(`IF Category`=ifelse(is.na(`IF Category`),
                                       "Not Reported",as.character(`IF Category`)))%>%kable()
```


## TABLE METHODS X: Distribution of Citable Items for Journal Titles by Category 


The total number of citable items per journal is categorized into ranges and the number of journals and percentage of total number of journals in each category is presented.

### 2013



```{r}

my_cut_fun <- function(x) {
  cut(x,
      breaks=c(0,100,500,1000,32000),
      labels=c("<100","100-500",
               "500-1000","1000-32000"),right=FALSE)
}

jdata = jdata%>%mutate(citable_2013_cat = my_cut_fun(jdata$citable_2013),
                       citable_2014_cat = my_cut_fun(jdata$citable_2014),
                       citable_2015_cat = my_cut_fun(jdata$citable_2015))


jdata%>%group_by(citable_2013_cat)%>%summarize("Number"=n())%>%
    mutate(
  "# (%)"=paste0(Number,"\t(",round(100*Number/nrow(jdata),1),"%)"))%>%select(-Number)%>%
  rename("Total Citable Items 2013 Category"=citable_2013_cat)%>%kable()

```

### 2014, 2015

```{r}
tmp1 = jdata%>%group_by(citable_2014_cat)%>%summarize("Number"=n())%>%
  rename(`Total Citable Items Category`=citable_2014_cat)%>%
  mutate(
    "2014 # (%)"=paste0(Number,"\t(",round(100*Number/nrow(jdata),1),"%)"))%>%select(-Number)
tmp2 = jdata%>%group_by(citable_2015_cat)%>%summarize("Number"=n())%>%
    rename(`Total Citable Items Category`=citable_2015_cat)%>%
  mutate(
    "2015 # (%)"=paste0(Number,"\t(",round(100*Number/nrow(jdata),1),"%)"))%>%select(-Number)

left_join(tmp1,tmp2)%>%mutate(`Total Citable Items Category`=ifelse(is.na(`Total Citable Items Category`),
                                         "Not Reported",
                                         as.character(`Total Citable Items Category`)))%>%kable()

```

## TABLE METHODS X: Distribution of Citable Items for Journal Titles

Min, 25th percentile, mean, median, 75th percentile, and max.

```{r}
tmp1 = jdata_long%>%group_by(year)%>%
  summarize_at(vars(citable_items),
               funs(Total="sum",Min="min",
                    Q25 = quantile(.,probs=.25),Mean="mean",Median="median",
                    Q75 = quantile(.,probs=.75),Max="max"),na.rm=T)%>%
  add_column("Journals"="With PLoS One",.after="year")
  
tmp2 = jdata_long%>%filter(!Journal=="PLoS One")%>%group_by(year)%>%
  summarize_at(vars(citable_items),
               funs(Total="sum",Min="min",
                    Q25 = quantile(.,probs=.25),Mean="mean",Median="median",
                    Q75 = quantile(.,probs=.75),Max="max"),na.rm=T) %>%
  add_column("Journals"="Remove PLoS One",.after="year") 
  
tmp = bind_rows(tmp1,tmp2)%>%arrange(Journals,year)
tmp = tmp%>%arrange(desc(Journals),year)
kable(tmp,caption="Citable Items by Year",digits=1)

```



## TABLE: Number of journals per data sharing mark (DSM)

```{r}
tmp = jdata%>%group_by(dsm,dsm_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata))

tmp1 = tmp%>%
  mutate("# Journals (%)" = paste0(N, "\t(",round(Percent,1),"%)"))%>%
  select(-N,-Percent)%>%
  rename("DSM"=dsm,"DSM Description"=dsm_fac)
kable(tmp1)
```

## TABLE: Publishing Volume by Data Sharing Mark


```{r}
tmp1 = jdata%>%group_by(dsm,dsm_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2013), 
            "MedianC"=median(citable_2013,na.rm=T),
            "PercentC" = 100*sum(citable_2013)/sum(jdata$citable_2013)
            )%>%
  mutate("# Journals (%)" = paste0(N, "\t(",round(Percent,1),"%)"),
         "Median # Citable Items per Journal 2013" = paste(MedianC," "), # need a space or kable malfunctions
         "# Citable Items 2013 (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC,-MedianC)%>%
  rename("DSM"=dsm,"DSM Description"=dsm_fac)

tmp2 = jdata%>%filter(!Journal=="PLoS One")%>%group_by(dsm,dsm_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2013), 
            "PercentC" = 100*sum(citable_2013)/sum(jdata$citable_2013[!jdata$Journal=="PLoS One"])
            )%>%
  mutate("# Citable Items 2013, Remove PLoS One (%)" = 
           paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC)%>%
  rename("DSM"=dsm,"DSM Description"=dsm_fac)

tmp = left_join(tmp1,tmp2)

tmp1 = jdata%>%group_by(dsm,dsm_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2014, na.rm=T), 
            "MedianC"=median(citable_2014,na.rm=T),
            "PercentC" = 100*sum(citable_2014,na.rm=T)/sum(jdata$citable_2014,na.rm=T)
            )%>%
  mutate("# Journals (%)" = paste0(N, "\t(",round(Percent,1),"%)"),
        "Median # Citable Items per Journal 2014" = paste(MedianC," "), # need a space or kable malfunctions
         "# Citable Items 2014 (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC,-MedianC)%>%
  rename("DSM"=dsm,"DSM Description"=dsm_fac)

tmp2 = jdata%>%filter(!Journal=="PLoS One")%>%group_by(dsm,dsm_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2014,na.rm=T), 
            "PercentC" = 
              100*sum(citable_2014,na.rm=T)/sum(jdata$citable_2014[!jdata$Journal=="PLoS One"],na.rm=T)
            )%>%
  mutate("# Citable Items 2014, Remove PLoS One (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC)%>%
  rename("DSM"=dsm,"DSM Description"=dsm_fac)

tmp3 = left_join(tmp1,tmp2)
tmp = left_join(tmp,tmp3)


tmp1 = jdata%>%group_by(dsm,dsm_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2015, na.rm=T), 
            "MedianC"=median(citable_2015,na.rm=T),
            "PercentC" = 100*sum(citable_2015,na.rm=T)/sum(jdata$citable_2015,na.rm=T)
            )%>%
  mutate("# Journals (%)" = paste0(N, "\t(",round(Percent,1),"%)"),
        "Median # Citable Items per Journal 2015" = paste(MedianC," "), # need a space or kable malfunctions
         "# Citable Items 2015 (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC,-MedianC)%>%
  rename("DSM"=dsm,"DSM Description"=dsm_fac)

tmp2 = jdata%>%filter(!Journal=="PLoS One")%>%group_by(dsm,dsm_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2015,na.rm=T), 
            "PercentC" = 
              100*sum(citable_2015,na.rm=T)/sum(jdata$citable_2015[!jdata$Journal=="PLoS One"],na.rm=T)
            )%>%
  mutate("# Citable Items 2015, Remove PLoS One (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC)%>%
  rename("DSM"=dsm,"DSM Description"=dsm_fac)

tmp3 = left_join(tmp1,tmp2)
tmp = left_join(tmp,tmp3)


#kable(tmp)

tmpt = data_frame("DSM Description"="Total",
           "# Journals (%)"=paste0(nrow(jdata)," (100%)"),
           "Median # Citable Items per Journal 2013"  =as.character(median(jdata$citable_2013)),
           "# Citable Items 2013 (%)"=paste0(sum(jdata$citable_2013)," (100%)"),
           "# Citable Items 2013, Remove PLoS One (%)"=
             paste0(sum(jdata$citable_2013[!jdata$Journal=="PLoS One"])," (100%)"),
           "Median # Citable Items per Journal 2014"  =as.character(median(jdata$citable_2014,na.rm=T)),
          "# Citable Items 2014 (%)"=paste0(sum(jdata$citable_2014,na.rm=T)," (100%)"),
           "# Citable Items 2014, Remove PLoS One (%)"=
             paste0(sum(jdata$citable_2014[!jdata$Journal=="PLoS One"],na.rm=T)," (100%)"),
            "Median # Citable Items per Journal 2015"  =as.character(median(jdata$citable_2015,na.rm=T)),
          "# Citable Items 2015 (%)"=paste0(sum(jdata$citable_2015,na.rm=T)," (100%)"),
           "# Citable Items 2015, Remove PLoS One (%)"=
             paste0(sum(jdata$citable_2015[!jdata$Journal=="PLoS One"],na.rm=T)," (100%)")
)
tmpt = bind_rows(tmp,tmpt)
kable(tmpt)

```


```{r barplot_pct,fig.width=9,fig.height=7,include=FALSE}
# 2013
tmp1 = jdata%>%group_by(dsm,dsm_fac)%>%
  summarize("N_journals"=n(),"Percent_journals"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2013), 
            "MedianC"=median(citable_2013,na.rm=T),
            "PercentC" = 100*sum(citable_2013)/sum(jdata$citable_2013)
            )
tmp1$year = "2013"
tmp1$plos = "with PLoS One"

tmp2 = jdata%>%filter(!Journal=="PLoS One")%>%group_by(dsm,dsm_fac)%>%
  summarize("N_journals"=n(),"Percent_journals"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2013), 
            "MedianC"=median(citable_2013,na.rm=T),
            "PercentC" = 100*sum(citable_2013)/sum(jdata$citable_2013[!jdata$Journal=="PLoS One"])
            )
tmp2$year = "2013"
tmp2$plos = "without PLoS One"

tmp = bind_rows(tmp1,tmp2)


# 2013
tmp1 = jdata%>%group_by(dsm,dsm_fac)%>%
  summarize("N_journals"=n(),"Percent_journals"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2014,na.rm=T), 
            "MedianC"=median(citable_2014,na.rm=T),
            "PercentC" = 100*sum(citable_2014,na.rm=T)/sum(jdata$citable_2014,na.rm=T)
            )
tmp1$year = "2014"
tmp1$plos = "with PLoS One"

tmp2 = jdata%>%filter(!Journal=="PLoS One")%>%group_by(dsm,dsm_fac)%>%
  summarize("N_journals"=n(),"Percent_journals"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2014,na.rm=T), 
            "MedianC"=median(citable_2014,na.rm=T),
            "PercentC" = 100*sum(citable_2014,na.rm=T)/sum(jdata$citable_2014[!jdata$Journal=="PLoS One"],na.rm=T))
tmp2$year = "2014"
tmp2$plos = "without PLoS One"

tmp3 = bind_rows(tmp1,tmp2)
tmp = bind_rows(tmp,tmp3)

tmp_n = tmp%>%select(dsm:Percent_journals,year,plos)%>%rename(N=N_journals,Percent=Percent_journals)
tmp_n$type = "Journal"
tmp_n = tmp_n%>%filter(plos=="with PLoS One") # basically the same
tmp_c = tmp%>%select(dsm,dsm_fac,Nc,PercentC,year,plos)%>%rename(N=Nc,Percent=PercentC)
tmp_c$type = "Citable Items"
tmplong = bind_rows(tmp_n,tmp_c)
tmplong = tmplong%>%unite(type2,type,plos,sep=" ",remove = FALSE)
tmplong$type2[tmplong$type=="Journal"] = "Journals"

tmplong = tmplong
tmplong = tmplong%>%filter(year=="2013") #similar by year
tmplong$type2 = factor(tmplong$type2,levels=c("Journals","Citable Items with PLoS One","Citable Items without PLoS One"))

ggplot(tmplong,aes(x=dsm,y=Percent,fill=dsm_fac,alpha=type2))+geom_bar(stat="identity",position = "dodge") +
  geom_text(aes(label=paste0(round(Percent,1),"%"),y=Percent),
            #alpha=1,
            position=position_dodge(width=0.9),vjust=1.4,size=2,color="black",show.legend = FALSE)+
  guides(alpha=guide_legend(title="Percent of",ncol=1),
         fill = guide_legend(title = "Data Sharing Mark", ncol = 1))+
#  facet_grid(year~.)+
  scale_alpha_manual(values = c(0.4, 0.6,0.9))+
  theme_minimal()+xlab("Data Sharing Mark (1-6)")+xlim(as.character(1:6))+ theme(legend.position = "bottom")

```


```{r barplot_pct_stacked,fig.width=9,fig.height=7,include=FALSE}

ggplot(tmplong,aes(x=type2,y=Percent,fill=dsm_fac))+geom_bar(stat="identity") +
  geom_text(aes(label=paste0(round(Percent,1),"%"),y=Percent),
            #alpha=1,
            position=position_stack(),vjust=1.4,size=3,color="black",show.legend = FALSE)+
  guides(alpha=guide_legend(title="Percent of",ncol=1),
         fill = guide_legend(title = "Data Sharing Mark", ncol = 1))+
  theme_minimal()+xlab("")+theme(axis.text.x = element_text(angle = 45, hjust = .9))
```

```{r barplot_pct_stacked_horiz,fig.width=9,fig.height=5}
tmplong$type1 = factor(tmplong$type2,levels=c("Citable Items without PLoS One","Citable Items with PLoS One","Journals"))

ggplot(tmplong,aes(x=type1,y=Percent,fill=dsm_fac))+geom_bar(stat="identity") +
  geom_text(aes(label=paste0(round(Percent,1),"%"),y=Percent),
            #alpha=1,
            position=position_stack(),hjust=1.2,size=3,color="black",show.legend = FALSE)+
  guides(alpha=guide_legend(title="Percent of",ncol=1),
         fill = guide_legend(title = "Data Sharing Mark", ncol = 1))+
  theme_minimal()+xlab("")+coord_flip()+ theme(legend.position = "bottom")
  

```

## FIGURE 1 - Percentage of journals per each data sharing mark (DSM).

The top bar shows the percentage of all journals for each data sharing mark. The middle bar
shows the percentage of citable items from each journal (including PLoS One) for each data
sharing mark. The lower bar shows the percentage of citable items for each journal
(excluding PLoS One) for each data sharing mark. Because of the journal PLoS One’s high
publishing activity, we analyzed the percentage of citable items for each data sharing mark
including and excluding PLoS One. The shades from dark to light represent DSM 1-6.

```{r fig1_barplot_pct_stacked_horiz,fig.width=9,fig.height=4,dpi=500}
tmplong$type1 = factor(tmplong$type2,levels=c("Citable Items without PLoS One","Citable Items with PLoS One","Journals"))
tmplong$dsm_fac_reverse = factor(tmplong$dsm_fac,levels=levels(tmplong$dsm_fac)[6:1])

ggplot(tmplong,aes(x=type1,y=Percent,fill=dsm_fac_reverse))+geom_bar(stat="identity") +
  geom_text(aes(label=paste0(round(Percent,1),"%"),y=Percent),
            #alpha=1,
            position=position_stack(),hjust=1.2,size=3,color="black",show.legend = FALSE)+
  #scale_fill_brewer(palette="Greys")+
  scale_fill_manual(values = RColorBrewer::brewer.pal(n=7, name="Blues")[1:6],
                    breaks=levels(tmplong$dsm_fac),
                    #labels=levels(tmplong$dsm_fac))+
                      labels=paste(1:6,levels(tmplong$dsm_fac)))+
  guides(fill = guide_legend(title = "Data Sharing Mark", ncol = 1))+
  theme_minimal()+xlab("")+coord_flip()+ theme(legend.position = "bottom")
  

```

## TABLE RESULTS X: Number of journals & Citable Items per data sharing 2 categories


```{r}

#summarize citable items
tmp1 = jdata%>%group_by(dsm2,dsm2_fac)%>%
  summarize("N"=n(),
            "Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2013), 
            "MedianC"=median(citable_2013,na.rm=T),
            "PercentC" = 100*sum(citable_2013)/sum(jdata$citable_2013)
            )%>%
  mutate("# Journals (%)" = paste0(N, "\t(",round(Percent,1),"%)"),
         "Median # Citable Items per Journal 2013" = paste(MedianC," "), # need a space or kable malfunctions
         "# Citable Items 2013 (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)"))%>%
  select(-N,-Nc,-Percent,-PercentC,-MedianC)%>%
  rename("DSM"=dsm2,"DSM Description"=dsm2_fac)

#summarize citable items, removing PLoS One
#do not need median since removing one journal will not change this except by .5
tmp2 = jdata%>%filter(!Journal=="PLoS One")%>%group_by(dsm2,dsm2_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2013), 
            "PercentC" = 100*sum(citable_2013)/sum(jdata$citable_2013[!jdata$Journal=="PLoS One"])
            )%>%
  mutate("# Citable Items 2013, Remove PLoS One (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC)%>%
  rename("DSM"=dsm2,"DSM Description"=dsm2_fac)

tmp = left_join(tmp1,tmp2)

tmp1 = jdata%>%group_by(dsm2,dsm2_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2014, na.rm=T), 
            "MedianC"=median(citable_2014,na.rm=T),
            "PercentC" = 100*sum(citable_2014,na.rm=T)/sum(jdata$citable_2014,na.rm=T)
            )%>%
  mutate("# Journals (%)" = paste0(N, "\t(",round(Percent,1),"%)"),
         "Median # Citable Items per Journal 2014" = paste(MedianC," "),
         "# Citable Items 2014 (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") 
         )%>%
  select(-N,-Nc,-Percent,-PercentC,-MedianC)%>%
  rename("DSM"=dsm2,"DSM Description"=dsm2_fac)

tmp2 = jdata%>%filter(!Journal=="PLoS One")%>%group_by(dsm2,dsm2_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2014,na.rm=T), 
            "PercentC" = 
              100*sum(citable_2014,na.rm=T)/sum(jdata$citable_2014[!jdata$Journal=="PLoS One"],na.rm=T)
            )%>%
  mutate("# Citable Items 2014, Remove PLoS One (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC)%>%
  rename("DSM"=dsm2,"DSM Description"=dsm2_fac)

tmp3 = left_join(tmp1,tmp2)
tmp = left_join(tmp,tmp3)


tmp1 = jdata%>%group_by(dsm2,dsm2_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2015, na.rm=T), 
            "MedianC"=median(citable_2015,na.rm=T),
            "PercentC" = 100*sum(citable_2015,na.rm=T)/sum(jdata$citable_2015,na.rm=T)
            )%>%
  mutate("# Journals (%)" = paste0(N, "\t(",round(Percent,1),"%)"),
         "Median # Citable Items per Journal 2015" = paste(MedianC," "),
         "# Citable Items 2015 (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") 
         )%>%
  select(-N,-Nc,-Percent,-PercentC,-MedianC)%>%
  rename("DSM"=dsm2,"DSM Description"=dsm2_fac)

tmp2 = jdata%>%filter(!Journal=="PLoS One")%>%group_by(dsm2,dsm2_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2015,na.rm=T), 
            "PercentC" = 
              100*sum(citable_2015,na.rm=T)/sum(jdata$citable_2015[!jdata$Journal=="PLoS One"],na.rm=T)
            )%>%
  mutate("# Citable Items 2015, Remove PLoS One (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC)%>%
  rename("DSM"=dsm2,"DSM Description"=dsm2_fac)

tmp3 = left_join(tmp1,tmp2)
tmp = left_join(tmp,tmp3)

tmp = tmp%>%arrange(desc(DSM))
tmp$DSM[1] = "DSM 1&2"
tmp$DSM[2] = "DSM 3-6"
kable(tmp)

```


## TABLE RESULTS X: Distribution of Open Access by Journal and Citable Item

```{r}
tmp1 = jdata%>%group_by(oam_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2013), 
            "MedianC"=median(citable_2013,na.rm=T),
            "PercentC" = 100*sum(citable_2013)/sum(jdata$citable_2013)
            )%>%
  mutate("# Journals (%)" = paste0(N, "\t(",round(Percent,1),"%)"),
         "Median # Citable Items per Journal 2013" = paste(MedianC," "),
         "# Citable Items 2013 (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC,-MedianC)%>%
  rename("Open Access"=oam_fac)

tmp2 = jdata%>%filter(!Journal=="PLoS One")%>%group_by(oam_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2013), 
            "PercentC" = 100*sum(citable_2013)/sum(jdata$citable_2013[!jdata$Journal=="PLoS One"])
            )%>%
  mutate("# Citable Items 2013, Remove PLoS One (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC)%>%
  rename("Open Access"=oam_fac)

tmp = left_join(tmp1,tmp2)



tmp1 = jdata%>%group_by(oam_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2014, na.rm =T), 
            "MedianC"=median(citable_2014,na.rm=T),
            "PercentC" = 100*sum(citable_2014, na.rm =T)/sum(jdata$citable_2014, na.rm=T)
            )%>%
  mutate("# Journals (%)" = paste0(N, "\t(",round(Percent,1),"%)"),
         "Median # Citable Items per Journal 2014" = paste(MedianC," "), # need  a space or kable malfunctions for some reason
         "# Citable Items 2014 (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC,-MedianC)%>%
  rename("Open Access"=oam_fac)

tmp2 = jdata%>%filter(!Journal=="PLoS One")%>%group_by(oam_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2014, na.rm =T), 
            "PercentC" = 100*sum(citable_2014, na.rm =T)/sum(jdata$citable_2014[!jdata$Journal=="PLoS One"], na.rm=T)
            )%>%
  mutate("# Citable Items 2014, Remove PLoS One (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC)%>%
  rename("Open Access"=oam_fac)

tmp3 = left_join(tmp1,tmp2)
tmp = left_join(tmp,tmp3)

tmp1 = jdata%>%group_by(oam_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2015, na.rm =T), 
            "MedianC"=median(citable_2015,na.rm=T),
            "PercentC" = 100*sum(citable_2015, na.rm =T)/sum(jdata$citable_2015, na.rm=T)
            )%>%
  mutate("# Journals (%)" = paste0(N, "\t(",round(Percent,1),"%)"),
         "Median # Citable Items per Journal 2015" = paste(MedianC," "), # need  a space or kable malfunctions for some reason
         "# Citable Items 2015 (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC,-MedianC)%>%
  rename("Open Access"=oam_fac)

tmp2 = jdata%>%filter(!Journal=="PLoS One")%>%group_by(oam_fac)%>%
  summarize("N"=n(),"Percent"=100*n()/nrow(jdata),
            "Nc"=sum(citable_2015, na.rm =T), 
            "PercentC" = 100*sum(citable_2015, na.rm =T)/sum(jdata$citable_2015[!jdata$Journal=="PLoS One"], na.rm=T)
            )%>%
  mutate("# Citable Items 2015, Remove PLoS One (%)" = paste0(Nc, "\t(",round(PercentC,1),"%)") )%>%
  select(-N,-Nc,-Percent,-PercentC)%>%
  rename("Open Access"=oam_fac)

tmp3 = left_join(tmp1,tmp2)
tmp = left_join(tmp,tmp3)


tmp = tmp%>%arrange(desc(`Open Access`))
kable(tmp)
```


## TABLE X: DSM by Open Access

Open Access Mark (OAM) has two categories, Subscription and Open Access.

Here we present a table of number of open access or subscription journals (number of citable items) within each DSM as well as the percent of journals (percent of citable items) that are open access within each DSM. For the DSM 1 we present number and percent citable items with and without PLoS One.

```{r, results='asis'}

# obtain open access counts at journal level
tmptable = with(jdata,table(dsm,oam_fac))
tmptable = cbind(tmptable,"% Open Access"=100*tmptable[,2]/rowSums(tmptable))
tmptable = as_data_frame(tmptable)
tmptable = tmptable%>%add_column()

# obtain open access counts at citable item level

jcit = left_join(
  left_join(jdata%>%group_by(dsm)%>%summarize("tot_cit"=sum(citable_2013)),
            jdata%>%filter(oam_fac=="Subscription")%>%group_by(dsm)%>%summarize("sub_cit"=sum(citable_2013))),
  jdata%>%filter(oam_fac=="Open Access")%>%group_by(dsm)%>%summarize(open_cit = sum(citable_2013))
)
jcit$open_cit[is.na(jcit$open_cit)] = 0
jcit$pct_open_cit = 100*jcit$open_cit/jcit$tot_cit

  
#without Plos One
open_cit = jdata%>%filter(Journal!="PLoS One",
               oam_fac=="Open Access",
               dsm==1)%>%summarize("open_cit"=sum(citable_2013))
tot_cit = jdata%>%filter(Journal!="PLoS One",
               dsm==1)%>%summarize("open_cit"=sum(citable_2013))
pct_cit = round(100*open_cit/tot_cit,2)


printtab = data_frame(
  "DSM"=dsm_labels,
  "Subscription"=paste0(tmptable$Subscription," (",jcit$sub_cit,")"),
  "Open Access"=paste0(tmptable$`Open Access`," (",jcit$open_cit,")"),
  "% Open Access"=paste0(round(tmptable$`% Open Access`,2),"% (",
                         round(jcit$pct_open_cit,2),"%)")
)
# add in without Plos One #s
printtab[1,3] = gsub(")",paste0("; ",open_cit,"*)"),printtab[1,3])
printtab[1,4] = gsub(")",paste0("; ",pct_cit,"%*)"),printtab[1,4])


kable(printtab)

```

*After removing PLoS One

## Table of 2 category DSM by Open Access

```{r}

# obtain open access counts at journal level
tmptable = with(jdata,table(dsm2,oam_fac))
tmptable = cbind(tmptable,"% Open Access"=100*tmptable[,2]/rowSums(tmptable))
tmptable = as_data_frame(tmptable)

# obtain open access counts at citable item level

jcit = left_join(
  left_join(jdata%>%group_by(dsm2)%>%summarize("tot_cit"=sum(citable_2013)),
            jdata%>%filter(oam_fac=="Subscription")%>%group_by(dsm2)%>%
              summarize("sub_cit"=sum(citable_2013))),
  jdata%>%filter(oam_fac=="Open Access")%>%group_by(dsm2)%>%summarize(open_cit = sum(citable_2013))
)
jcit$open_cit[is.na(jcit$open_cit)] = 0
jcit$pct_open_cit = 100*jcit$open_cit/jcit$tot_cit

  
#without Plos One
open_cit = jdata%>%filter(Journal!="PLoS One",
               oam_fac=="Open Access",
               dsm2==1)%>%summarize("open_cit"=sum(citable_2013))
tot_cit = jdata%>%filter(Journal!="PLoS One",
               dsm2==1)%>%summarize("open_cit"=sum(citable_2013))
pct_cit = round(100*open_cit/tot_cit,2)


printtab = data_frame(
  "DSM"=levels(jdata$dsm2_fac),
  "Subscription"=paste0(tmptable$Subscription," (",jcit$sub_cit,")"),
  "Open Access"=paste0(tmptable$`Open Access`," (",jcit$open_cit,")"),
  "% Open Access"=paste0(round(tmptable$`% Open Access`,2),"% (",
                         round(jcit$pct_open_cit,2),"%)")
)
printtab = printtab%>%arrange(desc(DSM))
# add in without Plos One #s
printtab[1,3] = gsub(")",paste0("; ",open_cit,"*)"),printtab[1,3])
printtab[1,4] = gsub(")",paste0("; ",pct_cit,"%*)"),printtab[1,4])


kable(printtab)

```

*After removing PLoS One



## TABLE X: Frequencies of Data Sharing Policies

```{r}

jdata$omic_data_sharing_fac = relevel(jdata$omic_data_sharing_fac,ref="No")
jdata$copyright_licensing_fac = relevel(jdata$copyright_licensing_fac,ref="No Mention")
jdata$archival_retention_fac = relevel(jdata$archival_retention_fac,ref="No Mention")
jdata$reproducibility_noted_fac = relevel(jdata$reproducibility_noted_fac,ref="No Mention")
varlist = c("omic_data_sharing_fac","rec_pref_sharing_fac","journal_hosts_fac","copyright_licensing_fac",
            "archival_retention_fac","reproducibility_noted_fac")
t1 = CreateTableOne(varlist,data = jdata)
x1 <- print(t1,printToggle = FALSE,noSpaces = FALSE)
tmp = rownames(x1)
tmptab = as_data_frame(x1)%>%add_column("Variable"=tmp,.before = "Overall")
tmptab$Variable = gsub("omic_data_sharing_fac","Omic Data Sharing Required",tmptab$Variable)
tmptab$Variable = gsub("rec_pref_sharing_fac","Recommended Sharing Method",tmptab$Variable)
tmptab$Variable = gsub("journal_hosts_fac","Journal Hosting Limit",tmptab$Variable)
tmptab$Variable = gsub("copyright_licensing_fac","Copyright Licensing of Data",tmptab$Variable)
tmptab$Variable = gsub("archival_retention_fac","Archival Retention Policy",tmptab$Variable)
tmptab$Variable = gsub("reproducibility_noted_fac","Reproducibility or Analogous Concepts Noted as Purpose of Data Policy",tmptab$Variable)


tmptab$Overall =  gsub(")","%)",tmptab$Overall)

kable(tmptab)
```

## TABLE: Recommended Data Sharing Method by DSM 1-5

```{r}
tmp = with(jdata,table(dsm_fac,rec_pref_sharing_fac))
tmp = cbind("dsm"=1:6,tmp)
tmp = tmp[-nrow(tmp),]
tmp = as_data_frame(tmp)
tmp[6,] = colSums(tmp)
tmp[6,1] = "Total"
tmppct = round(100*tmp[,-1]/rowSums(tmp[,-1]),2)

for(ii in 2:ncol(tmp)) {
  tmp[,ii] = paste0(unlist(tmp[,ii])," (",tmppct[,ii-1],"%)")
}
kable(tmp)

```

## FIGURE: Recommended Data Sharing Method by DSM 1-5

```{r fig2_tile_percent_recDSM_2, fig.width=9,fig.height=7,dpi=400}
tmp = with(jdata,table(dsm_fac,rec_pref_sharing_fac))
tmp = cbind("dsm"=1:6,tmp)
tmp = tmp[-nrow(tmp),]
tmp = as_data_frame(tmp)
tmp[6,] = colSums(tmp)
tmp[6,1] = "Total (DSM 1-5)"
tmppct = cbind(tmp[,1],100*tmp[,-1]/rowSums(tmp[,-1]))
tmpn = tmp

for(ii in 2:ncol(tmp)) {
  tmp[,ii] = paste0(unlist(tmp[,ii]),"\n(",round(tmppct[,ii],1),"%)")
}
tmplonglab = tmp%>%gather(key=type,value=label,-dsm)

tmplong = tmpn%>%gather(key=type,value=n,-dsm)
tmplong2 = tmppct%>%gather(key=type,value=Percent,-dsm)
tmplong = left_join(tmplong,tmplong2)
tmplong = left_join(tmplonglab,tmplong)
tmplong$type = factor(tmplong$type,
                      levels=c("Unspecified","Multiple methods equally recommended","By Reader Request to Authors",
                               "Journal Hosted","Public Online Repository"),
                      labels=c("Unspecified","Multiple Methods Equally Recommended","By Reader Request to Authors",
                               "Journal Hosted","Public Online Repository"))

ggplot(tmplong,aes(dsm,type,fill=Percent))+geom_tile(show.legend=TRUE)+geom_text(aes(label=label),color="white")+
  theme_minimal()+ylab("")+xlab("Data Sharing Mark")+scale_fill_continuous(limits=c(0,100))

```

## FIGURE 3: Recommended data sharing method by data sharing mark (DSM) 1-5.

The number (percent) of journals with each recommended data sharing method is represented by each tile, with darker shades denoting higher percentages of journals with the given data sharing method.

```{r fig3_tile_percent_recDSM, fig.width=9,fig.height=7,dpi=400}
tmplong$lowpct = factor(1*tmplong$Percent<25)
ggplot(tmplong,aes(dsm,type,fill=Percent))+
  geom_tile(show.legend=TRUE)+geom_text(aes(label=label,color=lowpct),show.legend = FALSE)+
  theme_minimal()+ylab("")+xlab("Data Sharing Mark")+
  scale_color_manual(values=c("White","Black"))+
  guides(color="none",label="none")+
  scale_fill_continuous(limits=c(0,100),low=RColorBrewer::brewer.pal(n=7, name="Blues")[1],high=RColorBrewer::brewer.pal(n=7, name="Blues")[7])

```






```{r}

tmptable = data.frame(table(jdata$omic_data_sharing_fac,useNA="ifany"))
colnames(tmptable) = c("Protein, Proteomic, Genomic Data Sharing Required with Deposit to Specific Data Banks","# (%)")
tmptable[,2] = paste0(tmptable[,2]," (",sprintf("%.1f",100*tmptable[,2]/sum(tmptable[,2]),1),"%)")
kable(tmptable)

tmptable = data.frame(table(jdata$rec_pref_sharing_fac,useNA="ifany"))
colnames(tmptable) = c("Recommended Sharing Method","# (%)")
tmptable[,2] = paste0(tmptable[,2]," (",sprintf("%.1f",100*tmptable[,2]/sum(tmptable[,2]),1),"%)")
kable(tmptable)

tmptable = data.frame(table(jdata$journal_hosts_fac,useNA="ifany"))
tmptable[,1] = as.character(tmptable[,1])
colnames(tmptable) = c("Journal Hosting Limit","# (%)")
tmptable[nrow(tmptable),1] ="N/A"
tmptable[,2] = paste0(tmptable[,2]," (",sprintf("%.1f",100*tmptable[,2]/sum(tmptable[,2]),1),"%)")
kable(tmptable)


tmptable = data.frame(table(jdata$copyright_licensing_fac,useNA="ifany"))
colnames(tmptable) = c("Copyright Licensing of Data","# (%)")
tmptable[,2] = paste0(tmptable[,2]," (",sprintf("%.1f",100*tmptable[,2]/sum(tmptable[,2]),1),"%)")
kable(tmptable)


tmptable = data.frame(table(jdata$archival_retention_fac,useNA="ifany"))
colnames(tmptable) = c("Archival Retention Policy","# (%)")
tmptable[,2] = paste0(tmptable[,2]," (",sprintf("%.1f",100*tmptable[,2]/sum(tmptable[,2]),1),"%)")
kable(tmptable)

tmptable = data.frame(table(jdata$reproducibility_noted_fac,useNA="ifany"))
colnames(tmptable) = c("Reproducibility or Analogous Concepts Noted as Purpose of Data Policy","# (%)")
tmptable[,2] = paste0(tmptable[,2]," (",sprintf("%.1f",100*tmptable[,2]/sum(tmptable[,2]),1),"%)")
kable(tmptable)

```

**Journal Hosting by Recommended Sharing Policy**

```{r}
tmptable = with(jdata,table(journal_hosts_fac,rec_pref_sharing_fac,useNA="ifany"))
rownames(tmptable)[nrow(tmptable)] = "N/A"
kable(tmptable)
```


## Distribution of Impact Factor - Normality Tests


First we must determine whether IF or some transformation of IF is normally distributed. We can visually assess this by plotting histograms and density plots.


```{r density_if, fig.show='hold', fig.width=5, fig.height=5}
# check normality
ggplot(jdata_long,aes(impact_factor,fill=year))+geom_density(alpha=0.4)+
  theme_minimal()+xlab("Impact Factor")+ggtitle("Density of IF by Year")
ggplot(jdata_long,aes(impact_factor,fill=year))+geom_histogram(alpha=0.4,bins=40)+
  theme_minimal()+xlab("Impact Factor")+ggtitle("Histogram of IF by Year")

# check normality of log
ggplot(jdata_long,aes(log(impact_factor),fill=year))+geom_density(alpha=0.4)+
  theme_minimal()+xlab("log(IF)")+ggtitle("Density of log(IF) by Year")
ggplot(jdata_long,aes(log(impact_factor),fill=year))+geom_histogram(alpha=0.4,bins=40)+
  theme_minimal()+xlab("Impact Factor")+ggtitle("Histogram of log(IF) by Year")

```

We can also test the Normality assumption with the Shapiro Wilk's Test. We can test IF by year, as well as log(IF). We can also test within each DSM group since the main assumption for ANOVA is that the outcome (IF) is normal within each group.

```{r shapirowilk}
norm_test1 = data.frame("Y"=c("IF 2013","IF 2014","log(IF) 2013","log(IF) 2014"),
           rbind(broom::tidy(shapiro.test(jdata$if_2013)),
                 broom::tidy(shapiro.test(jdata$if_2014)),
                 broom::tidy(shapiro.test(log(jdata$if_2013))),
                 broom::tidy(shapiro.test(log(jdata$if_2014)))
           ))
norm_test1$p.value = as.character(signif(norm_test1$p.value,2))
kable(norm_test1,digits=4,caption = "Shapiro Wilk Normality Test p-values of IF")

# check normality within dsm groups
norm_test = jdata%>%group_by(dsm_fac)%>%summarize(
  "if_2013.pvalue"=shapiro.test(if_2013)$p.value,
  "if_2014.pvalue"=shapiro.test(if_2014)$p.value,                                "log_if_2013.pvalue"=shapiro.test(log(if_2013))$p.value,
  "log_if_2014.pvalue"=shapiro.test(log(if_2014))$p.value
)
norm_test[,-1] = apply(norm_test[,-1] ,2,function(k) as.character(signif(k,2)))

kable(norm_test,digits=4,
      caption = "Shapiro Wilk Normality Test p-values of IF within DSM Groups")
```

The p-values from the Shapiro Wilk's test are mostly very significant (hence we reject the normality assumption) and in general Impact Factor appears to be quite skewed so **nonparametric tests** (Wilcoxon, Kruskal-Wallis) will be more appropriate. Taking the log helps somewhat but does not solve the problem. Furthermore log-IF is much less intepretable. **Hence, we use nonparametric tests to compare distributions.**

## Impact Factor vs. Data Sharing Mark

We wish to assess how Impact Factor differs between journals with different data sharing types.

### Summarize Impact Factor Within Groups

We can describe the minimum, maximum, mean, and median of Impact factor as well as total citations for each DSM type.



```{r impact_factor_dsm_table}
tmpsum = jdata_long%>%group_by(dsm,dsm_fac,year)%>%summarize(
  "Number of Journals"=n(),
  "min_IF"=min(impact_factor,na.rm=T),
  "mean_IF"=mean(impact_factor,na.rm=T),
  "median_IF"=median(impact_factor,na.rm=T),
  "max_IF"=max(impact_factor,na.rm=T),
  "min_TotalCites"=min(total_cites,na.rm=T),
  "mean_TotalCites"=mean(total_cites,na.rm=T),
  "median_TotalCites"=median(total_cites,na.rm=T),
  "max_TotalCites"=max(total_cites,na.rm=T)
)
tmpsum = tmpsum%>%rename(DSM=dsm_fac)
kable(tmpsum,digits=2)

tmpsum_IF_within_DSM = tmpsum
```



We can also visualize the distribution with a boxplot of IF by DSM type.

*Boxplot description*: The lower and upper "hinges" correspond to the first and third quartiles (the 25th and 75th percentiles). The upper whisker extends from the hinge to the highest value that is within 1.5 * IQR of the hinge, where IQR is the inter-quartile range, or distance between the first and third quartiles. The lower whisker extends from the hinge to the lowest value within 1.5 * IQR of the hinge. Data beyond the end of the whiskers are outliers and plotted as points (as specified by Tukey). (from `geom_boxplot` help page in the R `ggplot2` package.) The Triangle represents the mean.



```{r boxplot_if_dsm_1314,fig.width=9,fig.height=7}
ggplot(jdata_long%>%filter(year!=2015),aes(x=dsm,y=impact_factor,fill=dsm_fac,alpha=year))+
  geom_boxplot(position=position_dodge(.9))+theme_minimal()+
  scale_alpha_manual(values=c(0.3,0.7))+
  stat_summary(fun.y=mean,geom="point",position=position_dodge(.9),
               aes(x=dsm,y=impact_factor,shape="triangle"))+
  scale_shape_manual(values = 2)+
  guides(alpha="none",shape="none",fill=guide_legend(title="Data Sharing Mark",ncol=1))+
  xlab("Data Sharing Mark (1-6), Year (2013 light, 2014 dark)")+
  ylab("Impact Factor")+
  #ggtitle("Boxplot of Impact Factor by Data Sharing Mark")+
  scale_x_continuous(breaks=1:6)+
  theme(legend.position="bottom")
```

```{r boxplot_if_dsm,fig.width=9,fig.height=7}
ggplot(jdata_long,aes(x=dsm,y=impact_factor,fill=dsm_fac,alpha=year))+
  geom_boxplot(position=position_dodge(.9))+theme_minimal()+
  scale_alpha_manual(values=c(0.3,0.5,0.8))+
  stat_summary(fun.y=mean,geom="point",position=position_dodge(.9),
               aes(x=dsm,y=impact_factor,shape="triangle"))+
  scale_shape_manual(values = 2)+
  guides(alpha="none",shape="none",fill=guide_legend(title="Data Sharing Mark",ncol=1))+
  xlab("Data Sharing Mark (1-6), Year (2013-2015 light-dark)")+
  ylab("Impact Factor")+
  #ggtitle("Boxplot of Impact Factor by Data Sharing Mark")+
  scale_x_continuous(breaks=1:6)+
  theme(legend.position="bottom")
```

## FIGURE 2: Impact factors were higher for journals with the strongest data sharing policies (DSM 1) compared to journals with no mention of data sharing (DSM 6).

The median Impact Factor was calculated for the journals with each data sharing mark for
each report year (light color=2013, dark color=2014). The lower and upper hinges of the
boxplots represent the first and third quartiles of journal Impact Factor, the horizontal line
represents the median, the triangle represents the mean, and the upper and lower whiskers
extend from the hinge to the highest (lowest) value that is within 1.5 times the interquartile
range of the hinge, with journals outside this range represented as points.

```{r fig2_boxplot_if_dsm_1314,fig.width=9,fig.height=7,dpi=500}
tmplong = jdata_long%>%filter(year!=2015)
tmplong$dsm_fac_reverse = factor(tmplong$dsm_fac,levels=levels(tmplong$dsm_fac)[6:1])
ggplot(tmplong,aes(x=dsm,y=impact_factor,alpha=year,fill=dsm_fac_reverse))+
  geom_boxplot(position=position_dodge(.9))+theme_minimal()+
  scale_alpha_manual(values=c(0.7,0.9))+ 
  stat_summary(fun.y=mean,geom="point",position=position_dodge(.9),
               aes(x=dsm,y=impact_factor,shape="triangle"))+
  scale_shape_manual(values = 2)+
  scale_fill_manual(values = RColorBrewer::brewer.pal(n=7, name="Blues")[2:7],
                    #values=RColorBrewer::brewer.pal(n=6,name="Blues"),
                    breaks=levels(tmplong$dsm_fac),
                    labels=paste(1:6,levels(tmplong$dsm_fac)))+
  guides(alpha="none",shape="none",
         fill=guide_legend(title="Data Sharing Mark",ncol=1))+
  xlab("Data Sharing Mark (1-6), Year (2013 light, 2014 dark)")+
  ylab("Impact Factor")+
  #ggtitle("Boxplot of Impact Factor by Data Sharing Mark")+
  scale_x_continuous(breaks=1:6)+
  theme(legend.position="bottom")
```

```{r boxplot_citable_dsm_1314,fig.width=9,fig.height=7,include=FALSE}
ggplot(jdata_long%>%filter(!Journal=="PLoS One",year!=2015),aes(x=dsm,y=citable_items,fill=dsm_fac,alpha=year))+
  geom_boxplot(position=position_dodge(.9))+theme_minimal()+
  scale_alpha_manual(values=c(0.3,0.7))+
  stat_summary(fun.y=mean,geom="point",position=position_dodge(.9),
               aes(x=dsm,y=citable_items,shape="triangle"))+
  scale_shape_manual(values = 2)+
  guides(alpha="none",shape="none",fill="none")+
  xlab("Data Sharing Mark (1-6), Year (2013 light, 2014 dark)")+
  #ggtitle("Number of Citable Items per Journal by Data Sharing Mark")+
  ylab("Number of Citable Items per Journal")+scale_x_continuous(breaks=1:6)+
  theme(legend.position="bottom")
```

```{r boxplot_citable_dsm,fig.width=9,fig.height=7}
ggplot(jdata_long%>%filter(!Journal=="PLoS One"),aes(x=dsm,y=citable_items,fill=dsm_fac,alpha=year))+
  geom_boxplot(position=position_dodge(.9))+theme_minimal()+
  scale_alpha_manual(values=c(0.3,0.5,0.8))+
  stat_summary(fun.y=mean,geom="point",position=position_dodge(.9),
               aes(x=dsm,y=citable_items,shape="triangle"))+
  scale_shape_manual(values = 2)+
  guides(alpha="none",shape="none",fill="none")+
  xlab("Data Sharing Mark (1-6), Year (2013-2015 light-dark)")+
  #ggtitle("Number of Citable Items per Journal by Data Sharing Mark")+
  ylab("Number of Citable Items per Journal")+scale_x_continuous(breaks=1:6)+
  theme(legend.position="bottom")
```

Collapsing DSM into two categories:

```{r impact_factor_dsm_2cat_table}
tmpsum = jdata_long%>%group_by(dsm2_fac,year)%>%summarize(
  "Number of Journals"=n(),
  "min_IF"=min(impact_factor,na.rm=T),
  "mean_IF"=mean(impact_factor,na.rm=T),
  "median_IF"=median(impact_factor,na.rm=T),
  "max_IF"=max(impact_factor,na.rm=T),
  "min_TotalCites"=min(total_cites,na.rm=T),
  "mean_TotalCites"=mean(total_cites,na.rm=T),
  "median_TotalCites"=median(total_cites,na.rm=T),
  "max_TotalCites"=max(total_cites,na.rm=T)
)
tmpsum = tmpsum%>%rename(DSM=dsm2_fac)
kable(tmpsum,digits=2)

tmpsum_IF_within_DSM2 = tmpsum
```

```{r boxplot_if_dsm_2cat,fig.width=9,fig.height=5}
ggplot(jdata_long,aes(x=dsm2_fac_flip,y=impact_factor,fill=dsm2_fac_flip,alpha=year))+
  geom_boxplot(position=position_dodge(.9))+theme_minimal()+
  scale_alpha_manual(values=c(0.3,0.5,0.8))+
    stat_summary(fun.y=mean,geom="point",position=position_dodge(.9),
               aes(x=dsm2_fac,y=impact_factor,shape="triangle"))+
  scale_shape_manual(values = 2)+
  guides(alpha="none",shape="none",fill=guide_legend(title="Data Sharing"))+
  xlab("Data Sharing Required (No, Yes), Year (2013-2015 light-dark)")+
  #ggtitle("Boxplot of Impact Factor by Data Sharing Requirement")+
  ylab("Impact Factor")
```


```{r boxplot_citable_dsm_2cat,fig.width=9,fig.height=5}
ggplot(jdata_long%>%filter(!Journal=="PLoS One"),
       aes(x=dsm2_fac_flip,y=citable_items,fill=dsm2_fac_flip,alpha=year))+
  geom_boxplot(position=position_dodge(.9))+theme_minimal()+
    stat_summary(fun.y=mean,geom="point",position=position_dodge(.9),
               aes(x=dsm2_fac,y=citable_items,shape="triangle"))+
  scale_alpha_manual(values=c(0.3,0.5,0.8))+
  scale_shape_manual(values = 2)+
  guides(alpha="none",shape="none",fill=guide_legend(title="Data Sharing"))+
  xlab("Data Sharing Required (No, Yes), Year (2013-2015 light-dark)")+
  #ggtitle("Number of Citable Items per Journal by Data Sharing Requirement")+
  ylab("Number of Citable Items per Journal")
```

# Analysis: Journal Impact Factor, Open Access, Data Sharing Mark

## Methods

Continuous variables are summarized with medians and interquartile ranges (IQRs) denoting the 25th and 75th percentiles. Categorical variables are summarized with counts and percentages. The variables Impact Factor and total citable items are not normally distributed (Shapiro Wilk's Test p-values < 0.001) so medians are presented instead of means, and nonparametric methods are used for statistical tests.

The association of Impact Factor (IF) with 6-level data sharing mark (DSM) was tested with a nonparametric Kruskal-Wallis one-way analysis of variance (ANOVA) of IF in 2013 and 2014 with DSM as a grouping factor. Post-hoc pairwise two-sample Wilcoxon tests were used to determine whether the median IF for journals differ between the two level data sharing policy (required vs. not required) categories. P-values from the Wilcoxon tests were adjusted for multiple comparisons with the Holm procedure.

Pearson's chi-square test was used to test the association of data sharing policy (two levels: required vs not required) and open access status. Fisher's Exact Test was used to test the association of the 6-level DSM with open access status. Fisher’s Test was used as opposed to Chi-square test due to the low number of Open Access journals within some DSM categories. To examine the association of open access status and data sharing weighted by publishing volume we examined the number of citable items in each category and tested for the association of open access and data sharing with Pearson's chi-square test.

All statistical analyses were performed with `r R.Version()$version.string`. All code and data to reproduce these results can be found on github (https://github.com/OHSU-Ontology-Development-Group/DataSharingPolicies).


## Notes

- The Kruskal-Wallis test is a nonparametric version of ANOVA that tests whether the distribution of IF varies between DSM groups. 
- The Wilcoxon test is testing the difference in medians between two groups when the distributions of the outcome (IF in this case) are the same. Based on the boxplots and densities the distributions do look similar within DSM group so we are comfortable making inferences on the medians.

## IF ~ DSM

We perform a Kruskal Wallis test for the difference in median JIF between 6 category DSM.

We can also collapse DSM into two categories and perform a Wilcoxon test for difference in median JIF between required vs not required data sharing.

```{r if_kruskal}
k_if13 = with(jdata,kruskal.test(if_2013~dsm))
k_if14 = with(jdata,kruskal.test(if_2014~dsm))
k_if15 = with(jdata,kruskal.test(if_2015~dsm))

w_if13 = with(jdata,pairwise.wilcox.test(if_2013,dsm))
w_if14 = with(jdata,pairwise.wilcox.test(if_2014,dsm))
w_if15 = with(jdata,pairwise.wilcox.test(if_2015,dsm))

w_if13_2cat = with(jdata,wilcox.test(if_2013~dsm2))
w_if14_2cat = with(jdata,wilcox.test(if_2014~dsm2))
w_if15_2cat = with(jdata,wilcox.test(if_2015~dsm2))
```

```{r if_kruskal_res}
#kw results
kw_results = data.frame(Y=c("if_2013_6catDSM","if_2014_6catDSM","if_2015_6catDSM",
                            "if_2013_2catDSM","if_2014_2catDSM","if_2015_2catDSM"),
                        bind_rows(
  broom::tidy(k_if13),broom::tidy(k_if14),broom::tidy(k_if15),
  broom::tidy(w_if13_2cat),broom::tidy(w_if14_2cat),broom::tidy(w_if15_2cat)
  )
  )
kw_results$p.value = as.character(signif(kw_results$p.value,2))
kw_results = kw_results%>%rename("degrees of freedom"=parameter)
kable(kw_results,digits=2)
```

The pairwise Wilcoxon test p-values are below (by DSM group number) for 2013, 2014, and 2015.
The p-values are adjusted for multiple comparisons with the `r w_if13$p.adjust.method` method.

Significant p-values suggest that the median JIF is different for journals in the two DSM categories.

**Pairwise Wilcoxon p-values for IF 2013 between DSM:**
```{r wilcoxon_res_2013,  results='asis'}
signif_13 = signif(w_if13$p.value,2)
signifind = which(signif_13<.05,arr.ind=T)
emphasize.strong.cells(signifind)
signif_13[which(signif_13<.001,arr.ind=T)] = "< 0.001"
signif_13[which(is.na(signif_13),arr.ind=T)] = ""
pandoc.table(signif_13)
```

**Pairwise Wilcoxon p-values for IF 2014 between DSM:**
```{r wilcoxon_res_2014,  results='asis'}
signif_14 = signif(w_if14$p.value,2)
signifind = which(signif_14<.05,arr.ind=T)
emphasize.strong.cells(signifind)
signif_14[which(signif_14<.001,arr.ind=T)] = "< 0.001"
signif_14[which(is.na(signif_14),arr.ind=T)] = ""
pandoc.table(signif_14)
```

**Pairwise Wilcoxon p-values for IF 2015 between DSM:**
```{r wilcoxon_res_2015,  results='asis'}
signif_15 = signif(w_if15$p.value,2)
signifind = which(signif_15<.05,arr.ind=T)
emphasize.strong.cells(signifind)
signif_15[which(signif_15<.001,arr.ind=T)] = "< 0.001"
signif_15[which(is.na(signif_15),arr.ind=T)] = ""
pandoc.table(signif_15)
```




### Results

Impact factor is significantly associated with the six category data sharing mark (Kruskal-Wallis rank sum test, 5 df, p `r ifelse(k_if13$p.value<0.001,"< 0.001",round(k_if13$p.value,2))`, 2013 and 2014). Examining pairwise differences between DSMs we see that journals with DSM 1 have significantly higher JIF than journals with DSM 3, 4, 5, or 6 (Wilcoxon test, p `r signif_13["3","1"]`, `r signif_13["4","1"]`, `r signif_13["5","1"]`, `r signif_13["6","1"]`; 2013 data, 2014 similar). Journals with DSM 2 have significantly higher JIF than journals with DSM 3, 4, or 6 (Wilcoxon test, p = `r signif_13["3","2"]`,`r signif_13["4","2"]`, `r signif_13["6","2"]`; 2013 data, 2014 similar). Journals with DSM 5 have significantly higher JIF than journals with DSM 3, 4, and 6 (Wilcoxon test, p `r signif_13["5","3"]`, `r signif_13["5","4"]`, `r signif_13["6","5"]`; 2013 data, 2014 similar). In general, IF is not significantly different between DSM 1&2 and DSM 2&5, reflecting the similar JIF for journals with data sharing requirements, either full or partial sharing. After collapsing DSM into two categories, required (DSM 1-2) and not required (DSM 3-6) we still see a highly significant increase in JIF for journals with required data sharing (Wilcoxon Rank Sum Test, p `r ifelse(w_if13_2cat$p.value<0.001,"< 0.001",round(w_if13_2cat$p.value,2))`, 2013 and 2014 data). The median JIFs for DSM 1-6 are `r tmpsum_IF_within_DSM%>%filter(year==2013)%>%ungroup%>%select(median_IF)%>%unlist%>%round(2)%>%paste(collapse=", ")` in 2013 and `r tmpsum_IF_within_DSM%>%filter(year==2014)%>%ungroup%>%select(median_IF)%>%unlist%>%round(2)%>%paste(collapse=", ")` in 2014. The median JIFs for the collapsed two category DSM required and not required are `r tmpsum_IF_within_DSM2%>%filter(year==2013)%>%ungroup%>%arrange(desc(DSM))%>%select(median_IF)%>%unlist%>%round(2)%>%paste(collapse=", ")` in 2013 and
`r tmpsum_IF_within_DSM2%>%filter(year==2014)%>%ungroup%>%arrange(desc(DSM))%>%select(median_IF)%>%unlist%>%round(2)%>%paste(collapse=", ")` in 2014.

Results for 2015 are similar, except DSM 3 and 5 are no longer significantly different.

## Kruskal-Wallis test for year effect

A Kruskal-Wallis test was performed for year effect with impact factor and number of citable items. There were no significant differences in distribution impact factor nor citable items by year in the entire data set nor in subsets defined by DSM group. 

### Impact Factor

```{r}
tmp1 = broom::tidy(with(jdata_long,kruskal.test(impact_factor~as.factor(year))))

tmp2 = bind_rows(lapply(1:6, function(k){
  broom::tidy(with(jdata_long%>%filter(dsm==k),kruskal.test(impact_factor~as.factor(year))))
}))
tmp = bind_rows(tmp1,tmp2)
tmp = cbind("DSM"=c("All",paste0("DSM = ",1:6)),tmp)
tmp%>%kable(digits=3)
```

### Citable Items

```{r}
tmp1 = broom::tidy(with(jdata_long,kruskal.test(citable_items ~as.factor(year))))

tmp2 = bind_rows(lapply(1:6, function(k){
  broom::tidy(with(jdata_long%>%filter(dsm==k),kruskal.test(citable_items~as.factor(year))))
}))
tmp = bind_rows(tmp1,tmp2)
tmp = cbind("DSM"=c("All",paste0("DSM = ",1:6)),tmp)
tmp%>%kable(digits=3)
```



## Open Access vs. DSM


### Fisher's Exact Test for 6 category DSM vs. OAM

The Fisher's Exact Test for DSM vs. OAM tests for the independence of the categories of DSM and OAM (unordered).

Table of counts for 6 category DSM and two category OAM.
```{r oam_dsm_table}
tab1 = with(jdata,table(oam_fac,dsm))
kable(tab1)
```

Table of proportion with Open Access in each DSM category:
```{r oam_proportion_in_dsm}
tibble("DSM"=dsm_labels,"Proportion Open Access"=round(tab1[2,]/colSums(tab1),2))%>%kable()
```

We can test for the association of DSM and OAM with Fisher's Exact Test. The Test result is below:
```{r oam_dsm_fisher}
fishres = broom::tidy(fisher.test(tab1))
kable(fishres)
```

### Chi-Square test for 2 category DSM vs. OAM

Collapsing the categories into a 2x2 table makes the test hypothesis and result easier to interpret. When we collapse DSM into two categories (required vs. not required) there are more counts in each cell so we do not need to use the Fisher's Exact Test but instead can use a Chi-square test (commonly used for large samples).


The number of journals in each of the 2x2 table categories are below:
```{r oam_dsm2cat}
tab2 = with(jdata,table(oam_fac,dsm2_fac))
kable(tab2)
```

Table of proportion with Open Access in each DSM category:
```{r}
jdata%>%group_by(dsm2_fac)%>%summarize("Proportion  Open Access"=mean(oam))%>%
  rename("DSM"=dsm2_fac)%>%kable(digits=3)
```

Table of proportion data sharing required in each Open Access category:
```{r}
jdata%>%group_by(oam_fac)%>%summarize("Proportion  DSM Required"=mean(dsm2))%>%
  rename("OAM"=oam_fac)%>%kable(digits=3)
```


```{r oam_dsm2fac_chisq}
chires = broom::tidy(chisq.test(tab2))%>%rename(df=parameter)
kable(chires,digits=3)
```


### Results

The Fisher's Exact test is testing the hypothesis that open access status is associated with data sharing mark. The test is not significant (Fisher's Exact Test, p = `r signif(fishres$p.value,2)`) which suggests that the proportion of open access journals is not significantly differenct across data sharing mark categories.

The Chi-square test is testing the hypothesis that open access status is associated with data sharing requirement (two categories DSM 1-2 vs DSM 3-6). The test is not significant (Chi-square Test, df=1, p = `r round(chires$p.value,2)`) which suggests that journals with data sharing requirements are not any more likely to be open access than journals without data sharing requirement. Also, open access journals are not more likely to have data sharing requirements than subscription journals. This is further supported by the evidence that the proportion of open access journals that is similar for data sharing "required" vs. "non-required" journals.


# Publishing Volume

## All Journals

Here we determine how data sharing and open access are related once incorporating publishing volume. In this case, we are considering the "citable item" as the unit of measurement as opposed to journal. In other words, we ask, if we are given a citable item that is open access, is it more likely to have data sharing requirements than a citable item that is subscription based?

### Summary of Number of Citable Items by DSM and OAM

In 2013, the total number of citable items in the set of studied journals was `r jdata_long%>%filter(year==2013)%>%select(citable_items)%>%sum%>%as.character`, in 2014 it was
`r jdata_long%>%filter(year==2014)%>%select(citable_items)%>%na.omit%>%sum%>%as.character` and in 2015 it was
`r jdata_long%>%filter(year==2015)%>%select(citable_items)%>%na.omit%>%sum%>%as.character`.

Summary of citable items in 2013/2014/2015:
```{r citable_volume}



tmp = jdata_long%>%group_by(year,dsm)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))%>%
  rename(DSM=dsm)
kable(tmp,caption="Citable Items by DSM")

tmp = jdata_long%>%group_by(year,dsm2_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))%>%
  rename(DSM=dsm2_fac)
kable(tmp,caption="Citable Items by Required/Not Required DSM")

tmp = jdata_long%>%group_by(year,oam_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))%>%
  rename(OAM=oam_fac)
kable(tmp,caption="Citable Items by Open Access")
```

Summary with proportions:

```{r}
tmp = jdata_long%>%group_by(year,dsm)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T),
            "Proportion Citable Items Open Access"=
              sum(citable_items[oam==1],na.rm=T)/sum(citable_items,na.rm=T))%>%
  rename(DSM=dsm)
kable(tmp,digits=3)

tmp = jdata_long%>%group_by(year,dsm2_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T),
            "Proportion Citable Items Open Access"=
              sum(citable_items[oam==1],na.rm=T)/sum(citable_items,na.rm=T))%>%
  rename(DSM=dsm2_fac)
kable(tmp,digits=3)


tmp = jdata_long%>%group_by(year,oam_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T),
            "Proportion Citable Items Required DSM"=
              sum(citable_items[dsm2==1],na.rm=T)/sum(citable_items,na.rm=T))%>%
  rename(OAM=oam_fac)
kable(tmp,digits=3)
```


### Chi-square analysis - Citable Item

When we weight the number of open access journals and required data sharing journals by total citable items within each category there is a significant association between open access and data sharing requirement at the citable item level. That is, a citable item that is open access is much more likely to also have a data sharing requirement. This is mainly due to the fact that although the number of journals who have these open access or data sharing requirements is smaller than the number of journals that do not, the total citable articles within those journals is much larger. The p-values for the chi-square test at the citable item level are <2e-16, very significant.

```{r chisq_totalcitable,include=FALSE}
tmp = jdata_long%>%group_by(year,oam_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))
kable(tmp)

tmp = jdata_long%>%group_by(year,dsm2_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))
kable(tmp)

tmp = jdata_long%>%group_by(year,dsm2_fac,oam_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))
kable(tmp)
```

Chi-square test results:

```{r chisq_volume}
tmpm = data.frame(matrix(tmp$`Total Citable`[1:4],ncol=2))
colnames(tmpm) = levels(jdata$dsm2_fac)
rownames(tmpm) = levels(jdata$oam_fac)
kable(tmpm,caption="Citable Items by OAM and DSM, 2013")
tmpb = broom::tidy(chisq.test(tmpm))
if(tmpb$p.value < 2e-16) tmpb$p.value = "< 2e-16"
kable(tmpb)

tmpm = data.frame(matrix(tmp$`Total Citable`[1:4+4],ncol=2))
colnames(tmpm) = levels(jdata$dsm2_fac)
rownames(tmpm) = levels(jdata$oam_fac)
kable(tmpm,caption="Citable Items by OAM and DSM, 2014")
tmpb = broom::tidy(chisq.test(tmpm))
if(tmpb$p.value < 2e-16) tmpb$p.value = "< 2e-16"
kable(tmpb)

tmpm = data.frame(matrix(tmp$`Total Citable`[1:4+8],ncol=2))
colnames(tmpm) = levels(jdata$dsm2_fac)
rownames(tmpm) = levels(jdata$oam_fac)
kable(tmpm,caption="Citable Items by OAM and DSM, 2014")
tmpb = broom::tidy(chisq.test(tmpm))
if(tmpb$p.value < 2e-16) tmpb$p.value = "< 2e-16"
kable(tmpb)
```

Chi-square test results for DSM1-6 vs Open Access:

2013:
```{r}
tmp = jdata_long%>%group_by(year,dsm,oam_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))
tmp = add_row(as.data.frame(tmp),year=2013,dsm=4,oam_fac="Open Access",`Num Journals`=0,`Total Citable`=0)
tmp = add_row(as.data.frame(tmp),year=2014,dsm=4,oam_fac="Open Access",`Num Journals`=0,`Total Citable`=0)%>%arrange(year,dsm,oam_fac)

tmpm = data.frame(
  matrix(tmp%>%filter(year==2013)%>%ungroup%>%select(`Total Citable`)%>%unlist,nrow=2))
tmpb = broom::tidy(chisq.test(tmpm))
if(tmpb$p.value < 2e-16) tmpb$p.value = "< 2e-16"
kable(tmpb)
```

2014:
```{r}
tmpm = data.frame(
  matrix(tmp%>%filter(year==2014)%>%ungroup%>%select(`Total Citable`)%>%unlist,nrow=2))
tmpb = broom::tidy(chisq.test(tmpm))
if(tmpb$p.value < 2e-16) tmpb$p.value = "< 2e-16"
kable(tmpb)
```

2015:
```{r}
tmpm = data.frame(
  matrix(tmp%>%filter(year==2015)%>%ungroup%>%select(`Total Citable`)%>%unlist,nrow=2))
tmpb = broom::tidy(chisq.test(tmpm))
if(tmpb$p.value < 2e-16) tmpb$p.value = "< 2e-16"
kable(tmpb)
```

## Trend test for proportion of citable items with required data sharing within open access categories across years

```{r}
tmp = jdata_long%>%group_by(year,oam_fac,dsm2)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))

# within OA
tmpreq = tmp%>%filter(oam_fac=="Open Access",dsm2==1)%$%`Total Citable`
tmpnot = tmp%>%filter(oam_fac=="Open Access",dsm2==0)%$%`Total Citable`
tmptot = tmpreq + tmpnot

tmpres_oa = prop.test(tmpreq,tmptot)
tmpres_oa_trend = prop.trend.test(tmpreq,tmptot)


# within subscription
tmpreq = tmp%>%filter(oam_fac=="Subscription",dsm2==1)%$%`Total Citable`
tmpnot = tmp%>%filter(oam_fac=="Subscription",dsm2==0)%$%`Total Citable`
tmptot = tmpreq + tmpnot

tmpres_sub_trend = prop.trend.test(tmpreq,tmptot)
tmpres_sub = prop.test(tmpreq,tmptot)


tmpres = bind_rows(broom::tidy(tmpres_sub),broom::tidy(tmpres_oa))
tmpres_trend = bind_rows(broom::tidy(tmpres_sub_trend),broom::tidy(tmpres_oa_trend))

colnames(tmpres)[1:3] = c("2013","2014","2015")
tmpres_withPLOS = cbind(tmpres[,1:3],tmpres_trend)

```

```{r}
tmp = jdata_long%>%filter(!Journal=="PLoS One")%>%group_by(year,oam_fac,dsm2)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))

# within OA
tmpreq = tmp%>%filter(oam_fac=="Open Access",dsm2==1)%$%`Total Citable`
tmpnot = tmp%>%filter(oam_fac=="Open Access",dsm2==0)%$%`Total Citable`
tmptot = tmpreq + tmpnot

tmpres_oa = prop.test(tmpreq,tmptot)
tmpres_oa_trend = prop.trend.test(tmpreq,tmptot)


# within subscription
tmpreq = tmp%>%filter(oam_fac=="Subscription",dsm2==1)%$%`Total Citable`
tmpnot = tmp%>%filter(oam_fac=="Subscription",dsm2==0)%$%`Total Citable`
tmptot = tmpreq + tmpnot

tmpres_sub_trend = prop.trend.test(tmpreq,tmptot)
tmpres_sub = prop.test(tmpreq,tmptot)


tmpres = bind_rows(broom::tidy(tmpres_sub),broom::tidy(tmpres_oa))
tmpres_trend = bind_rows(broom::tidy(tmpres_sub_trend),broom::tidy(tmpres_oa_trend))

colnames(tmpres)[1:3] = c("2013","2014","2015")
tmpres = cbind(tmpres[,1:3],tmpres_trend)

```


When looking at the open access journals including PLoS One, we can see that the proportion of data sharing required citable items is *decreasing* from 2013 to 2015 within open access journals (`r paste(round(tmpres_withPLOS[2,1:3],2),collapse=", ")`,   p=`r signif(tmpres_withPLOS$p.value[2],2)`). 

After removing PLos One, we can see that the proportion of data sharing required citable items is *increasing* from 2013 to 2015 within open access journals (`r paste(round(tmpres[2,1:3],2),collapse=", ")`,   p=`r signif(tmpres$p.value[2],2)`).

However, this proportion does not increase within subscription journals (`r paste(round(tmpres[1,1:3],2),collapse=", ")`, p=`r signif(tmpres$p.value[1],2)`).

The results are below:

```{r}
kable(cbind("PLOS"=c("withPLOS","withPLOS","withoutPLOS","withoutPLOS"),bind_rows(tmpres_withPLOS,tmpres)),digits=3)
```


## Remove PLoS One

PLoS One may be skewing the results since it has such high volume, so we try the above analysis after removing this journal.

```{r remove_plosone}
jdata_long0 =  jdata_long # keep old jdata_long object
jdata_long = jdata_long%>%filter(!Journal=="PLoS One")
```

### Summary of Number of Citable Items by DSM and OAM

In 2013, the total number of citable items in the set of studied journals (removing Plos One) was `r jdata_long%>%filter(year==2013)%>%select(citable_items)%>%sum%>%as.character`, in 2014 it was
`r jdata_long%>%filter(year==2014)%>%select(citable_items)%>%na.omit%>%sum%>%as.character`, and in 2015 it was
`r jdata_long%>%filter(year==2015)%>%select(citable_items)%>%na.omit%>%sum%>%as.character`

Summary of citable items in 2013/2014/2015:
```{r citable_volume2}

tmp = jdata_long%>%group_by(year,dsm)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))%>%
  rename(DSM=dsm)
kable(tmp,caption="Citable Items by DSM")

tmp = jdata_long%>%group_by(year,dsm2_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))%>%
  rename(DSM=dsm2_fac)
kable(tmp,caption="Citable Items by Required/Not Required DSM")

tmp = jdata_long%>%group_by(year,oam_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))%>%
  rename(OAM=oam_fac)
kable(tmp,caption="Citable Items by Open Access")
```

Summary with proportions:
```{r}
tmp = jdata_long%>%group_by(year,dsm)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T),
            "Proportion Citable Items Open Access"=
              sum(citable_items[oam==1],na.rm=T)/sum(citable_items,na.rm=T))%>%
  rename(DSM=dsm)
kable(tmp,digits=3)

tmp = jdata_long%>%group_by(year,dsm2_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T),
            "Proportion Citable Items Open Access"=
              sum(citable_items[oam==1],na.rm=T)/sum(citable_items,na.rm=T))%>%
  rename(DSM=dsm2_fac)
kable(tmp,digits=3)


tmp = jdata_long%>%group_by(year,oam_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T),
            "Proportion Citable Items Required DSM"=
              sum(citable_items[dsm2==1],na.rm=T)/sum(citable_items,na.rm=T))%>%
  rename(OAM=oam_fac)
kable(tmp,digits=3)
```

### Chi-square analysis - Citable Item

The results are still highly significant, likely due to the large number of citable items, but an open access article is still more likely to have been published under data sharing requirements than a subscription article. In 2013, `r signif(tmp[2,5],3)*100`% of open access articles had data sharing requirements as opposed to `r signif(tmp[1,5],3)*100`% for subscription articles. In 2014 and 2015 it was `r signif(tmp[4,5],3)*100`% vs. `r signif(tmp[3,5],3)*100`% and `r signif(tmp[6,5],3)*100`% vs. `r signif(tmp[5,5],3)*100`%.

```{r chisq_totalcitable2}
tmp = jdata_long%>%group_by(year,oam_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))
kable(tmp)

tmp = jdata_long%>%group_by(year,dsm2_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))
kable(tmp)

tmp = jdata_long%>%group_by(year,dsm2_fac,oam_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))
kable(tmp)
```

Chi-squre test results: 
```{r chisq_volume2}
tmpm = data.frame(matrix(tmp$`Total Citable`[1:4],ncol=2))
colnames(tmpm) = levels(jdata$dsm2_fac)
rownames(tmpm) = levels(jdata$oam_fac)
kable(tmpm,caption="Citable Items by OAM and DSM, 2013")
tmpb = broom::tidy(chisq.test(tmpm))
if(tmpb$p.value < 2e-16) tmpb$p.value = "< 2e-16"
kable(tmpb)

tmpm = data.frame(matrix(tmp$`Total Citable`[1:4+4],ncol=2))
colnames(tmpm) = levels(jdata$dsm2_fac)
rownames(tmpm) = levels(jdata$oam_fac)
kable(tmpm,caption="Citable Items by OAM and DSM, 2014")
tmpb = broom::tidy(chisq.test(tmpm))
if(tmpb$p.value < 2e-16) tmpb$p.value = "< 2e-16"
kable(tmpb)

tmpm = data.frame(matrix(tmp$`Total Citable`[1:4+8],ncol=2))
colnames(tmpm) = levels(jdata$dsm2_fac)
rownames(tmpm) = levels(jdata$oam_fac)
kable(tmpm,caption="Citable Items by OAM and DSM, 2015")
tmpb = broom::tidy(chisq.test(tmpm))
if(tmpb$p.value < 2e-16) tmpb$p.value = "< 2e-16"
kable(tmpb)
```

Chi-square test results for DSM1-6 vs Open Access:

2013:
```{r}
tmp = jdata_long%>%group_by(year,dsm,oam_fac)%>%
  summarize("Num Journals"=n(),"Total Citable"=sum(citable_items,na.rm=T))
tmp = add_row(as.data.frame(tmp),year=2013,dsm=4,oam_fac="Open Access",`Num Journals`=0,`Total Citable`=0)
tmp = add_row(as.data.frame(tmp),year=2014,dsm=4,oam_fac="Open Access",`Num Journals`=0,`Total Citable`=0)%>%arrange(year,dsm,oam_fac)

tmpm = data.frame(
  matrix(tmp%>%filter(year==2013)%>%ungroup%>%select(`Total Citable`)%>%unlist,nrow=2))
tmpb = broom::tidy(chisq.test(tmpm))
if(tmpb$p.value < 2e-16) tmpb$p.value = "< 2e-16"
kable(tmpb)
```

2014:
```{r}
tmpm = data.frame(
  matrix(tmp%>%filter(year==2014)%>%ungroup%>%select(`Total Citable`)%>%unlist,nrow=2))
tmpb = broom::tidy(chisq.test(tmpm))
if(tmpb$p.value < 2e-16) tmpb$p.value = "< 2e-16"
kable(tmpb)
```

2015:
```{r}
tmpm = data.frame(
  matrix(tmp%>%filter(year==2015)%>%ungroup%>%select(`Total Citable`)%>%unlist,nrow=2))
tmpb = broom::tidy(chisq.test(tmpm))
if(tmpb$p.value < 2e-16) tmpb$p.value = "< 2e-16"
kable(tmpb)
```

## Results


Of the 318 journals examined, `r sum(jdata$dsm==1)` (`r round(mean(jdata$dsm==1)*100,1)`%) required data sharing as a condition of publication, barring exceptions and `r sum(jdata$dsm==2)` (`r round(mean(jdata$dsm==2)*100,1)`%) required data sharing but made no explicit statement regarding the effect on publication and editorial decisions. `r sum(jdata$dsm==3)` (`r round(mean(jdata$dsm==3)*100,1)`%) journals explicitly encouraged or addressed data sharing, but did not require it. And, `r sum(jdata$dsm==5)` (`r round(mean(jdata$dsm==5)*100,1)`%) journals only addressed data sharing for proteomic, genomic data, or other specific structural data.

In order to understand the potential influence of the policies on the published literature, we also evaluated the distribution of policies by publication volume. In 2013, the total number of citable items (papers) in the studied journals was `r jdata_long%>%filter(year==2013)%>%select(citable_items)%>%sum%>%as.character` (after removing PlOS One), in 2014 the total number of citable items was `r jdata_long%>%filter(year==2014)%>%select(citable_items)%>%sum(.,na.rm=T)%>%as.character`
and in 2015 the total number of citable items was `r jdata_long%>%filter(year==2015)%>%select(citable_items)%>%sum(.,na.rm=T)%>%as.character`,

with median numbers of citable items per journal of 
`r jdata_long%>%filter(year==2013)%>%select(citable_items)%>%unlist%>%median(na.rm=T)`,
`r jdata_long%>%filter(year==2014)%>%select(citable_items)%>%unlist%>%median(na.rm=T)`, and 
`r jdata_long%>%filter(year==2015)%>%select(citable_items)%>%unlist%>%median(na.rm=T)`, respectively. 



# Sample Representation of All Journals

We can compare sampled journals to non-sampled journals ("out of sample"") or to the entire set of journals ("All Journals").

```{r}
jdata_all <- read_excel("../Data-Sharing-Policies_2017-01-25.xlsx", sheet = "2014 ESI Data",
                    na = "N/A")

jdata_all = jdata_all%>%rename(Journal=`Full Journal Title`,
                               if_2014=`Journal Impact Factor`,citable_2014=`Citable Items`,
                               total_cites_2014=`Total Cites`)

jdata_sample = jdata%>%select(Journal,ISSN,if_2014,citable_2014,total_cites_2014)%>%add_column("in_sample"="Sample")
```

Note, these journals are in the sampled data but not in the "full data", but they have no data anyway:
```{r}
table(jdata_sample$Journal%in%jdata_all$Journal)
jdata_sample%>%filter(!Journal%in%jdata_all$Journal)%>%kable()
```

Number of journals in and out of sample (removing journals with missing data):

```{r}
# for comparing in sample to all journals
jdata_all_long = jdata_all
jdata_all_long$in_sample = "All Journals"
jdata_all_long = bind_rows(jdata_all_long,jdata_sample)
jdata_all_long = jdata_all_long%>%filter(!is.na(total_cites_2014))
table(jdata_all_long$in_sample)

# for comparing in sample to out of sample
jdata_all = left_join(jdata_all,jdata_sample%>%select(-ISSN))
jdata_all$in_sample[is.na(jdata_all$in_sample)] = "Out of sample"
jdata_all = jdata_all%>%filter(!is.na(Journal))
table(jdata_all$in_sample)

```

## Distribution of Impact Factor

```{r}
ggplot(jdata_all,aes(x=in_sample,y=if_2014,fill=in_sample))+geom_boxplot()+theme_minimal()+ggtitle("Distribution of Impact Factor")
ggplot(jdata_all_long,aes(x=in_sample,y=if_2014,fill=in_sample))+geom_boxplot()+theme_minimal()+ggtitle("Distribution of Impact Factor")

jdata_all_long%>%group_by(in_sample)%>%
  summarize_at(vars(if_2014),
               funs(Min="min",
                    Q25 = quantile(.,probs=.25),Mean="mean",Median="median",
                    Q75 = quantile(.,probs=.75),Max="max"),na.rm=T)%>%kable(digits=2)
```



## Distribution of Citable Items

```{r}

ggplot(jdata_all,aes(x=in_sample,y=citable_2014,fill=in_sample))+geom_boxplot()+theme_minimal()+ggtitle("Distribution of Citable Items")
ggplot(jdata_all%>%filter(citable_2014<20000),aes(x=in_sample,y=citable_2014,fill=in_sample))+geom_boxplot()+theme_minimal()+ggtitle("Distribution of Citable Items (remove Plos One)")

ggplot(jdata_all_long,aes(x=in_sample,y=citable_2014,fill=in_sample))+geom_boxplot()+theme_minimal()+ggtitle("Distribution of Citable Items")
ggplot(jdata_all_long%>%filter(citable_2014<20000),aes(x=in_sample,y=citable_2014,fill=in_sample))+geom_boxplot()+theme_minimal()+ggtitle("Distribution of Citable Items (remove Plos One)")


jdata_all_long%>%group_by(in_sample)%>%
  summarize_at(vars(citable_2014),
               funs(Total="sum",Min="min",
                    Q25 = quantile(.,probs=.25),Mean="mean",Median="median",
                    Q75 = quantile(.,probs=.75),Max="max"),na.rm=T)%>%kable(digits=2)
```



After removing PLoS One:

```{r}

jdata_all_long%>%filter(!Journal=="PLoS One")%>%group_by(in_sample)%>%
  summarize_at(vars(citable_2014),
               funs(Min="min",
                    Q25 = quantile(.,probs=.25),Mean="mean",Median="median",
                    Q75 = quantile(.,probs=.75),Max="max"),na.rm=T)%>%kable(digits=2)


```


## Distribution ofTotal Citations

```{r}
ggplot(jdata_all,aes(x=in_sample,y=total_cites_2014,fill=in_sample))+geom_boxplot()+theme_minimal()+ggtitle("Distribution of Total Citations")
ggplot(jdata_all_long,aes(x=in_sample,y=total_cites_2014,fill=in_sample))+geom_boxplot()+theme_minimal()+ggtitle("Distribution of Total Citations")


jdata_all_long%>%group_by(in_sample)%>%
  summarize_at(vars(total_cites_2014),
               funs(Total="sum",Min="min",
                    Q25 = quantile(.,probs=.25),Mean="mean",Median="median",
                    Q75 = quantile(.,probs=.75),Max="max"),na.rm=T)%>%kable(digits=2)
```



After removing PLoS One:

```{r}
jdata_all_long%>%filter(!Journal=="PLoS One")%>%group_by(in_sample)%>%
  summarize_at(vars(total_cites_2014),
               funs(Total="sum",Min="min",
                    Q25 = quantile(.,probs=.25),Mean="mean",Median="median",
                    Q75 = quantile(.,probs=.75),Max="max"),na.rm=T)%>%kable(digits=2)
```

# Publisher Comparison

Number of journals for each publisher, sorted:

```{r}
table(jdata$Publisher)%>%sort
```

Filter out publishers with 3 or fewer journals:

```{r}
tmp = table(jdata$Publisher)
tmp[tmp>3]%>%sort
tmpjournals = names(tmp[tmp>3])
```

Percent of sampled journals in each DSM category by publisher:

```{r}
library(janitor)
tmpdat = jdata%>%filter(Publisher%in%tmpjournals)
tmpout1 = tmpdat%>%tabyl(Publisher)
tmpout = tmpdat%>%crosstab(Publisher,dsm,"row")
tmpout2 = tmpout
tmpout2[,-1] = tmpout2[,-1]*100
tmpout2 = cbind(tmpout2,"# Journals"=tmpout1[,2])
kable(tmpout2,digits=1)
```

Summary of maximum percentages for each publisher:
```{r}
summary(apply(tmpout2[,-1],1,max))
```

Categorize into 2 level data sharing mark:

```{r}
tmpdat = jdata%>%filter(Publisher%in%tmpjournals)
tmpout = tmpdat%>%crosstab(Publisher,dsm2,"row")
tmpout1 = tmpdat%>%tabyl(Publisher)
tmpout2 = tmpout
tmpout2[,-1] = tmpout2[,-1]*100
tmpout2 = cbind(tmpout2,"# Journals"=tmpout1[,2])
kable(tmpout2,digits=1)
```

Summary of maximum percentages for each publisher:
```{r}
summary(apply(tmpout2[,-1],1,max))
```

# Curator Reliability

```{r}
cdata1 <- read_excel("../Data-Sharing-Policies_2017-01-25.xlsx",sheet = "Curator Data")
cdata0 <- read_excel("../Data-Sharing-Policies_2017-01-25.xlsx",sheet = "Data")
cdata1 <- cdata1[-1,]
cdata0 <- cdata0[-1,]



cdata1_all = cdata1%>%clean_names()%>%remove_empty_rows()%>%
  select(journal,open_access_mark,data_sharing_mark:reproducibility_noted_mark)
cdata0_all = cdata0%>%clean_names()%>%remove_empty_rows()%>%
  select(journal,open_access_mark,data_sharing_mark:reproducibility_noted_mark)

cdata0_all = data.frame(cdata0_all)

cdata0_all[which(cdata0_all=="a\n",arr.ind=T)] = "a"
cdata0_all[which(cdata0_all=="b\n",arr.ind=T)] = "b"


cdata1_all = data.frame(cdata1_all)

cdata1_all[which(cdata1_all=="a\n",arr.ind=T)] = "a"
cdata1_all[which(cdata1_all=="b\n",arr.ind=T)] = "b"

allvars = colnames(cdata1_all)[-1]
```

For each we have a table of agreement, the % agreement and Cohen's Kappa.


```{r}
restable = NULL
for(tmpvar in allvars) {
  tmp0 = cdata0_all%>%select(journal,contains(tmpvar))
  tmp1 = cdata1_all%>%select(journal,contains(tmpvar))
  colnames(tmp0)[2] = "var_orig"
  colnames(tmp1)[2] = "var_new"
  newdat = left_join(tmp1,tmp0)
  newdat = na.omit(newdat)
  mytab = table(newdat[,-1])
  print(tmpvar)
  print(mytab)
  
  res = data_frame("Score"=tmpvar,
                   "# Journals"=nrow(newdat),
                   "% Agreement"=100*sum(diag(mytab))/nrow(newdat),
                   "Cohen's Kappa"=kappa2(newdat[,-1])$value)
  restable = bind_rows(restable,res)
  
  
}
restable$Score = gsub("_"," ",restable$Score)

kable(restable,digits=3)
```


# Session Info

This analysis was performed in R/Rstudio using knitr with the following session:

```{r sessioninfo}
devtools::session_info()
```