rmarkdown-covid-19-transparent-analysis_short-version.Rmd

---
title: "Transparency of COVID-19-related research - Analysis"
author: Ahmad Sofi-Mahmudi, Eero Raittio, Sergio E. Uribe
date: "2022-02-11"
output: html_document
---

# Loading required packages

```{r}
pacman::p_load(dplyr,
               ggplot2,
               knitr,
               here,
               tidyr,
               europepmc,
               lubridate,
               ggpubr,
               epiR,
               gtsummary,
               forcats,
               nortest,
               caret,
               tibble,
               janitor)
```


# Loading the datasets

```{r}
opendata <- read.csv(here("data", "covid_transparency_opendata_final_new.csv"))
```


# Results
## Overall perspective
Number of all papers (open access and non open access), open-access papers and open-access percentage:
```{r}
kable(data.frame(hits_all = 258678, 
           hits_oa = 186279,
           oa_percentage = round((186279/258678)*100, 1)))
```

Undownloadable papers:
```{r}
c(number = 186279 - nrow(opendata), 
  percentage = round((186279 - nrow(opendata))/186279*100, 1))
```

Year of publication:

First, adding the real publication year and month to the datasets. The real publication year/month is the year/month the paper was first appear online stored in firstPublicationDate column.
```{r}
opendata <- opendata %>%
        mutate(year_firstpub = year(as.POSIXlt(firstPublicationDate,
                                               format = "%Y-%m-%d")),
               month_firstpub = month(as.POSIXlt(firstPublicationDate,
                                                 format = "%Y-%m-%d")))
```

Then:
```{r}
tabyl(opendata$year_firstpub) %>% adorn_pct_formatting()
```


Types of papers according to EPMC classification:
```{r}
kable(tabyl(opendata$type) %>%
              adorn_pct_formatting() %>%
              as.data.frame() %>%
              arrange(desc(n)))
```

Total number of journals:
```{r}
length(table(opendata$journalTitle))
```


Top journals:
```{r}

kable(head(tabyl(opendata$journalTitle) %>%
                 adorn_pct_formatting() %>%
                 as.data.frame() %>%
                 arrange(desc(n)), 10))
```

## Research papers
First, we subset the data to only include research papers:
```{r}
research <- opendata %>% filter(type == "research-article")
```

Number of open access papers per year:
```{r}
kable(tabyl(research$year_firstpub) %>% adorn_pct_formatting())
```

The number of journals in our dataset:
```{r}
length(table(research$journalTitle))
```

Top 10 journals with the highest number of research articles in our dataset, from high to low:
```{r}
kable(head(tabyl(research$journalTitle) %>%
                 adorn_pct_formatting() %>%
                 as.data.frame() %>%
                 arrange(desc(n)), 10))
```


The mean and the median of the number of citations to these references:
```{r}
kable(data.frame(Mean = round(mean(research$citedByCount),1),
           SD = round(sd(research$citedByCount), 2),
           Median = median(research$citedByCount),
           IQR = IQR(research$citedByCount)
           ))

```

Characteristics of the paper with the highest number of citations:
```{r}
kable(research[which.max(research$citedByCount),] %>% 
  select(citedByCount, pmid, pmcid, doi, title, authorString, journalTitle, year_firstpub))
```

The number and percentage of papers with a conflict of interest (CoI) disclosure:
```{r}
kable(data.frame(
  number = length(research$is_coi_pred[research$is_coi_pred == TRUE]),
           percentage = round(length(research$is_coi_pred[research$is_coi_pred == TRUE])/nrow(research)*100, 1)
           ))
```

Confidence interval for CoI:
```{r}

kable(round(epi.prev(pos = length(research$is_coi_pred[research$is_coi_pred == TRUE]),
         tested = nrow(research),
         se = 0.992,
         sp = 0.995)$ap, 
      1))
```

Comparing CoI disclosure rate in three years:
```{r}
kable(round(prop.table(table(research$is_coi_pred, research$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(research$is_coi_pred, research$year_firstpub)
```

The number and percentage of papers with a funding statement:
```{r}
kable(data.frame(number = length(research$is_fund_pred[research$is_fund_pred == TRUE]),
           percentage = round(length(research$is_fund_pred[research$is_fund_pred == TRUE])/nrow(research)*100, 1)
           ))
```

Confidence interval for funding statement:
```{r}

kable(round(epi.prev(pos = length(research$is_fund_pred[research$is_fund_pred == TRUE]),
         tested = nrow(research),
         se = 0.997,
         sp = 0.981)$ap, 
      1))
```

Comparing funding statement rate in three years:
```{r}
kable(round(prop.table(table(research$is_fund_pred, research$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(research$is_fund_pred, research$year_firstpub)
```

The number and percentage of papers that were registered beforehand:
```{r}
kable(data.frame(number = length(research$is_register_pred[research$is_register_pred == TRUE]),
           percentage = round(length(research$is_register_pred[research$is_register_pred == TRUE])/nrow(research)*100, 1)
           ))
```

Confidence interval for registration:
```{r}

kable(round(epi.prev(pos = length(research$is_register_pred[research$is_register_pred == TRUE]),
         tested = nrow(research),
         se = 0.955,
         sp = 0.997)$ap, 
      1))
```

Comparing registration rate in three years:
```{r}
kable(round(prop.table(table(research$is_register_pred, research$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(research$is_register_pred, research$year_firstpub)
```


The number and percentage of papers that shared data:
```{r}
kable(data.frame(number = length(research$is_open_data[research$is_open_data == TRUE]),
           percentage = round(length(research$is_open_data[research$is_open_data == TRUE])/nrow(research)*100, 1)
           ))
```

Confidence interval for data sharing:
```{r}

kable(round(epi.prev(pos = length(research$is_open_data[research$is_open_data == TRUE]),
         tested = nrow(research),
         se = 0.758,
         sp = 0.986)$ap, 
      1))
```

Comparing data sharing rate in three years:
```{r}
kable(round(prop.table(table(research$is_open_data, research$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(research$is_open_data, research$year_firstpub)
```

The number and percentage of papers that shared code:
```{r}
kable(data.frame(number = length(research$is_open_code[research$is_open_code == TRUE]),
           percentage = round(length(research$is_open_code[research$is_open_code == TRUE])/nrow(research)*100, 1)
           ))
```

Confidence interval for sharing code:
```{r}

kable(round(epi.prev(pos = length(research$is_open_code[research$is_open_code == TRUE]),
         tested = nrow(research),
         se = 0.587,
         sp = 0.997)$ap, 
      1))
```

Comparing sharing code rate in three years:
```{r}
kable(round(prop.table(table(research$is_open_code, research$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(research$is_open_code, research$year_firstpub)
```

Sum of five indicators:
```{r}
research %>% mutate(sum = rowSums(research %>% select(is_coi_pred, is_register_pred, is_fund_pred, is_open_data, is_open_code))) %>% select(sum) %>% tbl_summary()
```



## Randomized controlled trials
First, we subset the data to only include RCTs using is_rct variable:
```{r}
rct <- opendata %>% filter(is_rct == TRUE)
```

Total number of RCTs:
```{r}
nrow(rct)
```


Number of open access papers per year:
```{r}
kable(tabyl(rct$year_firstpub) %>% adorn_pct_formatting())
```

The number of journals in our dataset:
```{r}
length(table(rct$journalTitle))
```
Top 10 journals with the highest number of articles in our dataset, from high to low:
```{r}
kable(head(tabyl(rct$journalTitle) %>%
                 adorn_pct_formatting() %>%
                 as.data.frame() %>%
                 arrange(desc(n)), 10))
```


The mean and the median of the number of citations to these references:
```{r}
kable(data.frame(Mean = round(mean(rct$citedByCount),1),
           SD = round(sd(rct$citedByCount), 2),
           Median = median(rct$citedByCount),
           IQR = IQR(rct$citedByCount)
           ))

```

Characteristics of the paper with the highest number of citations:
```{r}
kable(rct[which.max(rct$citedByCount),] %>% 
  select(citedByCount, pmid, pmcid, doi, title, authorString, journalTitle, year_firstpub))
```

The number and percentage of papers with a conflict of interest (CoI) disclosure:
```{r}
kable(data.frame(
  number = length(rct$is_coi_pred[rct$is_coi_pred == TRUE]),
           percentage = round(length(rct$is_coi_pred[rct$is_coi_pred == TRUE])/nrow(rct)*100, 1)
           ))
```

Confidence interval for CoI:
```{r}

kable(round(epi.prev(pos = length(rct$is_coi_pred[rct$is_coi_pred == TRUE]),
         tested = nrow(rct),
         se = 0.992,
         sp = 0.995)$ap, 
      1))
```


The number and percentage of papers with a funding statement:
```{r}
kable(data.frame(number = length(rct$is_fund_pred[rct$is_fund_pred == TRUE]),
           percentage = round(length(rct$is_fund_pred[rct$is_fund_pred == TRUE])/nrow(rct)*100, 1)
           ))
```

Confidence interval for funding statement:
```{r}

kable(round(epi.prev(pos = length(rct$is_fund_pred[rct$is_fund_pred == TRUE]),
         tested = nrow(rct),
         se = 0.997,
         sp = 0.981)$ap, 
      1))
```


The number and percentage of papers that were registered beforehand:
```{r}
kable(data.frame(number = length(rct$is_register_pred[rct$is_register_pred == TRUE]),
           percentage = round(length(rct$is_register_pred[rct$is_register_pred == TRUE])/nrow(rct)*100, 1)
           ))
```

Confidence interval for registration:
```{r}

kable(round(epi.prev(pos = length(rct$is_register_pred[rct$is_register_pred == TRUE]),
         tested = nrow(rct),
         se = 0.955,
         sp = 0.997)$ap, 
      1))
```



The number and percentage of papers that shared data:
```{r}
kable(data.frame(number = length(rct$is_open_data[rct$is_open_data == TRUE]),
           percentage = round(length(rct$is_open_data[rct$is_open_data == TRUE])/nrow(rct)*100, 1)
           ))
```

Confidence interval for data sharing:
```{r}

kable(round(epi.prev(pos = length(rct$is_open_data[rct$is_open_data == TRUE]),
         tested = nrow(rct),
         se = 0.758,
         sp = 0.986)$ap, 
      1))
```


The number and percentage of papers that shared code:
```{r}
kable(data.frame(number = length(rct$is_open_code[rct$is_open_code == TRUE]),
           percentage = round(length(rct$is_open_code[rct$is_open_code == TRUE])/nrow(rct)*100, 1)
           ))
```

Confidence interval for sharing code:
```{r}

kable(round(epi.prev(pos = length(rct$is_open_code[rct$is_open_code == TRUE]),
         tested = nrow(rct),
         se = 0.587,
         sp = 0.997)$ap, 
      1))
```


### Comparing years
Comparing CoI disclosure rate in three years:
```{r}
kable(round(prop.table(table(rct$is_coi_pred, rct$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(rct$is_coi_pred, rct$year_firstpub)
```

Comparing funding statement rate in three years:
```{r}
kable(round(prop.table(table(rct$is_fund_pred, rct$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(rct$is_fund_pred, rct$year_firstpub)
```

Comparing registration rate in three years:
```{r}
kable(round(prop.table(table(rct$is_register_pred, rct$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(rct$is_register_pred, rct$year_firstpub)
```


Comparing data sharing rate in three years:
```{r}
kable(round(prop.table(table(rct$is_open_data, rct$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(rct$is_open_data, rct$year_firstpub)
```

Comparing sharing code rate in three years:
```{r}
kable(round(prop.table(table(rct$is_open_code, rct$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(rct$is_open_code, rct$year_firstpub)
```


Sum of five indicators:
```{r}
rct %>% mutate(sum = rowSums(rct %>% select(is_coi_pred, is_register_pred, is_fund_pred, is_open_data, is_open_code))) %>% select(sum) %>% tbl_summary()
```


## Review papers
First, we subset the data to only include review papers using is_review variable:
```{r}
review <- opendata %>% filter(is_review == TRUE)
```

Total number of reviews:
```{r}
nrow(review)
```

Number of open access papers per year:
```{r}
kable(tabyl(review$year_firstpub) %>% adorn_pct_formatting())
```

The number of journals in our dataset:
```{r}
length(table(review$journalTitle))
```
Top 10 journals with the highest number of articles in our dataset, from high to low:
```{r}
kable(head(tabyl(review$journalTitle) %>%
                 adorn_pct_formatting() %>%
                 as.data.frame() %>%
                 arrange(desc(n)), 10))
```


The mean and the median of the number of citations to these references:
```{r}
kable(data.frame(Mean = round(mean(review$citedByCount),1),
           SD = round(sd(review$citedByCount), 2),
           Median = median(review$citedByCount),
           IQR = IQR(review$citedByCount)
           ))

```

Characteristics of the paper with the highest number of citations:
```{r}
kable(review[which.max(review$citedByCount),] %>% 
  select(citedByCount, pmid, pmcid, doi, title, authorString, journalTitle, year_firstpub))
```

The number and percentage of papers with a conflict of interest (CoI) disclosure:
```{r}
kable(data.frame(
  number = length(review$is_coi_pred[review$is_coi_pred == TRUE]),
           percentage = round(length(review$is_coi_pred[review$is_coi_pred == TRUE])/nrow(review)*100, 1)
           ))
```

Confidence interval for CoI:
```{r}

kable(round(epi.prev(pos = length(review$is_coi_pred[review$is_coi_pred == TRUE]),
         tested = nrow(review),
         se = 0.992,
         sp = 0.995)$ap, 
      1))
```


The number and percentage of papers with a funding statement:
```{r}
kable(data.frame(number = length(review$is_fund_pred[review$is_fund_pred == TRUE]),
           percentage = round(length(review$is_fund_pred[review$is_fund_pred == TRUE])/nrow(review)*100, 1)
           ))
```

Confidence interval for funding statement:
```{r}

kable(round(epi.prev(pos = length(review$is_fund_pred[review$is_fund_pred == TRUE]),
         tested = nrow(review),
         se = 0.997,
         sp = 0.981)$ap, 
      1))
```


The number and percentage of papers that were registered beforehand:
```{r}
kable(data.frame(number = length(review$is_register_pred[review$is_register_pred == TRUE]),
           percentage = round(length(review$is_register_pred[review$is_register_pred == TRUE])/nrow(review)*100, 1)
           ))
```

Confidence interval for registration:
```{r}

kable(round(epi.prev(pos = length(review$is_register_pred[review$is_register_pred == TRUE]),
         tested = nrow(review),
         se = 0.955,
         sp = 0.997)$ap, 
      1))
```


The number and percentage of papers that shared data:
```{r}
kable(data.frame(number = length(review$is_open_data[review$is_open_data == TRUE]),
           percentage = round(length(review$is_open_data[review$is_open_data == TRUE])/nrow(review)*100, 1)
           ))
```

Confidence interval for data sharing:
```{r}

kable(round(epi.prev(pos = length(review$is_open_data[review$is_open_data == TRUE]),
         tested = nrow(review),
         se = 0.758,
         sp = 0.986)$ap, 
      1))
```


The number and percentage of papers that shared code:
```{r}
kable(data.frame(number = length(review$is_open_code[review$is_open_code == TRUE]),
           percentage = round(length(review$is_open_code[review$is_open_code == TRUE])/nrow(review)*100, 1)
           ))
```

Confidence interval for sharing code:
```{r}

kable(round(epi.prev(pos = length(review$is_open_code[review$is_open_code == TRUE]),
         tested = nrow(review),
         se = 0.587,
         sp = 0.997)$ap, 
      1))
```

### Comparing years
Comparing CoI disclosure rate in three years:
```{r}
kable(round(prop.table(table(review$is_coi_pred, review$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(review$is_coi_pred, review$year_firstpub)
```

Comparing funding statement rate in three years:
```{r}
kable(round(prop.table(table(review$is_fund_pred, review$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(review$is_fund_pred, review$year_firstpub)
```

Comparing registration rate in three years:
```{r}
kable(round(prop.table(table(review$is_register_pred, review$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(review$is_register_pred, review$year_firstpub)
```


Comparing data sharing rate in three years:
```{r}
kable(round(prop.table(table(review$is_open_data, review$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(review$is_open_data, review$year_firstpub)
```

Comparing sharing code rate in three years:
```{r}
kable(round(prop.table(table(review$is_open_code, review$year_firstpub), 2)*100, 1))
```

Chi-square test:
```{r}
chisq.test(review$is_open_code, review$year_firstpub)
```


Sum of five indicators:
```{r}
review %>% mutate(sum = rowSums(review %>% select(is_coi_pred, is_register_pred, is_fund_pred, is_open_data, is_open_code))) %>% select(sum) %>% tbl_summary() 
```



# Tables
## Journal-related differences in transparency (Table 1)

Table 1A
```{r}
set.seed(100)
research %>% 
  select(journalTitle,is_coi_pred,is_fund_pred,is_register_pred,is_open_data,is_open_code) %>%
        mutate(journalTitle = fct_lump(journalTitle,n=5)) %>%
  tbl_summary(by=journalTitle, percent="column",
              label = c(is_coi_pred ~ "COI disclosure", is_fund_pred ~ "Funding disclosure", is_register_pred ~ "Protocol registration", is_open_data ~ "Data sharing", is_open_code ~ "Code sharing")) %>%
  add_p(test.args=all_tests("fisher.test") ~ list(simulate.p.value=TRUE)) %>%
        as_flex_table()
```

Table 1B
```{r}
set.seed(100)
rct %>% 
  select(journalTitle,is_coi_pred,is_fund_pred,is_register_pred,is_open_data,is_open_code) %>%
        mutate(journalTitle = fct_lump(journalTitle,n=5)) %>%
  tbl_summary(by=journalTitle, percent="column",
              label = c(is_coi_pred ~ "COI disclosure", is_fund_pred ~ "Funding disclosure", is_register_pred ~ "Protocol registration", is_open_data ~ "Data sharing", is_open_code ~ "Code sharing")) %>%
  add_p(test.args=all_tests("fisher.test") ~ list(simulate.p.value=TRUE)) %>%
        as_flex_table()
```

Table 1C
```{r}
set.seed(100)
review %>% 
  select(journalTitle,is_coi_pred,is_fund_pred,is_register_pred,is_open_data,is_open_code) %>%
        mutate(journalTitle = fct_lump(journalTitle,n=5)) %>%
  tbl_summary(by=journalTitle, percent="column",
              label = c(is_coi_pred ~ "COI disclosure", is_fund_pred ~ "Funding disclosure", is_register_pred ~ "Protocol registration", is_open_data ~ "Data sharing", is_open_code ~ "Code sharing")) %>%
  add_p(test.args=all_tests("fisher.test") ~ list(simulate.p.value=TRUE)) %>%
        as_flex_table()
```


## Publisher-related differences in transparency (Table 2)

Table 2A
```{r}
set.seed(100)
research %>% 
  select(scimago_publisher,is_coi_pred,is_fund_pred,is_register_pred,is_open_data,is_open_code) %>%
        mutate(scimago_publisher = fct_lump(scimago_publisher,n=5)) %>%
  tbl_summary(by=scimago_publisher, percent="column",
              label = c(is_coi_pred ~ "COI disclosure", is_fund_pred ~ "Funding disclosure", is_register_pred ~ "Protocol registration", is_open_data ~ "Data sharing", is_open_code ~ "Code sharing")) %>%
  add_p(test.args=all_tests("fisher.test") ~ list(simulate.p.value=TRUE)) %>%
        as_flex_table()
```


Table 2B
```{r}
set.seed(100)
rct %>% 
  select(scimago_publisher,is_coi_pred,is_fund_pred,is_register_pred,is_open_data,is_open_code) %>%
        mutate(scimago_publisher = fct_lump(scimago_publisher,n=5)) %>%
  tbl_summary(by=scimago_publisher, percent="column",
              label = c(is_coi_pred ~ "COI disclosure", is_fund_pred ~ "Funding disclosure", is_register_pred ~ "Protocol registration", is_open_data ~ "Data sharing", is_open_code ~ "Code sharing")) %>%
  add_p(test.args=all_tests("fisher.test") ~ list(simulate.p.value=TRUE)) %>%
        as_flex_table()
```

Table 2C
```{r}
set.seed(100)
review %>% 
  select(scimago_publisher,is_coi_pred,is_fund_pred,is_register_pred,is_open_data,is_open_code) %>%
        mutate(scimago_publisher = fct_lump(scimago_publisher,n=5)) %>%
  tbl_summary(by=scimago_publisher, percent="column",
              label = c(is_coi_pred ~ "COI disclosure", is_fund_pred ~ "Funding disclosure", is_register_pred ~ "Protocol registration", is_open_data ~ "Data sharing", is_open_code ~ "Code sharing")) %>%
  add_p(test.args=all_tests("fisher.test") ~ list(simulate.p.value=TRUE)) %>%
        as_flex_table()
```


## Journal Impact Factor (JIF) and citation count analyses (Table 3)

Table 3A
```{r}
set.seed(100)

map(
  c("is_coi_pred", "is_fund_pred", "is_register_pred", "is_open_data", "is_open_code"),
  ~ tbl_summary(
      research,
      by = all_of(.x),
      include = citedByCount,
      label = list(citedByCount = .x)
    ) %>% add_p()
) %>% tbl_stack()

```


Table 3B
```{r}
set.seed(100)

map(
  c("is_coi_pred", "is_fund_pred", "is_register_pred", "is_open_data", "is_open_code"),
  ~ tbl_summary(
      research,
      by = all_of(.x),
      include = jif2020,
      label = list(jif2020 = .x)
    ) %>% add_p()
) %>% tbl_stack()

```


# Figure 1

## Figure 1A
### Figure 1A - left: proportion of research papers that adhered to each transparency indicator
```{r}

prop_research <- research %>%
        summarise("COI Disclosure" = sum(is_coi_pred == TRUE),
                  "Funding disclosure" = sum(is_fund_pred == TRUE),
                  "Protocol registration" = sum(is_register_pred == TRUE),
                  "Data sharing" = sum(is_open_data == TRUE),
                  "Code sharing" = sum(is_open_code == TRUE)) %>%
        t() %>%
        as.data.frame() %>%
        rownames_to_column(var = "indicator") %>%
        mutate(percentage = round(V1/nrow(research)*100, 1))

fig1A_left <- prop_research %>% ggplot() +
        aes(x = reorder(indicator, V1),
            y = V1,
            fill = indicator,
            col = indicator) +
        geom_col(alpha = 0.75) +
        geom_text(aes(label = percentage), hjust = -0.1, size = 3) +
        scale_y_continuous(
                name = "\nNumber of research articles",
                breaks = seq(0, 10E4, 2.5E4),
                labels = seq(0, 10E4, 2.5E4),
                expand = expansion(mult = c(0, 0.2))
                ) +
        coord_flip() +
        xlab(NULL) +
        theme(axis.title.x=element_blank(),
                legend.position = "none", 
                panel.grid.major.y = element_blank()
        )
```

### Figure 1A - right: monthly trend for each indicator
For the line graph, first, we add the real publication year and month to the datasets. The real publication year/month is the year/month the paper was first appear online stored in firstPublicationDate column.
```{r}
research$YearMonth <- format(as.Date(research$firstPublicationDate), "%Y-%m")
```


```{r}

monthly_research <- 
        research %>% 
        select(YearMonth,
               is_coi_pred,
               is_fund_pred,
               is_register_pred,
               is_open_data,
               is_open_code) %>%
        gather("indicator", "value", -YearMonth) %>%
        count(YearMonth, indicator, value) %>%
        mutate(indicator = recode(indicator,
                                  is_coi_pred = "COI Disclosure",
                                  is_fund_pred = "Funding disclosure",
                                  is_register_pred = "Protocol registration",
                                  is_open_data = "Data sharing",
                                  is_open_code = "Code sharing")) %>%
        complete(indicator, value, YearMonth, fill = list(n = 0)) %>%
        group_by(YearMonth, indicator) %>% 
        mutate(p = n / sum(n)) %>%
        filter(value) %>%
        ungroup()

fig1A_right <-
        monthly_research %>% 
        ggplot() +
        aes(x = YearMonth, 
            y = p,
            group = indicator,
            color = indicator) +
        geom_line(size = 0.75) +
        scale_y_continuous(limits = c(0, 1), 
                           labels = scales::percent) +
        scale_color_discrete(name = NULL) +
        scale_fill_discrete(breaks = c("COI Disclosure",
                                       "Funding disclosure",
                                       "Protocol registration",
                                       "Data sharing",
                                       "Code sharing")) +
        labs(y = NULL, 
             x = NULL) +
        theme(panel.grid.minor = element_blank(),
              legend.position = "none",
              axis.text.x = element_blank(),
              legend.key.size = unit(0.2, 'cm')
              )

fig1A <- ggarrange(fig1A_left, fig1A_right,
                    ncol = 2, nrow = 1, common.legend = F)

# tiff("Figure.tiff", width = 14, height = 7, units = "in", res = 300)
fig1A
#dev.off()

```

## Figure 1B
### Figure 1B - left: proportion of RCT papers that adhered to each transparency indicator
```{r}

prop_rct <- rct %>%
        summarise("COI Disclosure" = sum(is_coi_pred == TRUE),
                  "Funding disclosure" = sum(is_fund_pred == TRUE),
                  "Protocol registration" = sum(is_register_pred == TRUE),
                  "Data sharing" = sum(is_open_data == TRUE),
                  "Code sharing" = sum(is_open_code == TRUE)) %>%
        t() %>%
        as.data.frame() %>%
        rownames_to_column(var = "indicator") %>%
        mutate(percentage = round(V1/nrow(rct)*100, 1))

fig1B_left <- prop_rct %>% ggplot() +
        aes(x = reorder(indicator, V1),
            y = V1,
            fill = indicator,
            col = indicator) +
        geom_col(alpha = 0.75) +
        geom_text(aes(label = percentage), hjust = -0.1, size = 3) +
        scale_y_continuous(
                name = "\nNumber of RCT articles (in hundreds)",
                breaks = seq(0, 10E2, 2.5E2),
                labels = seq(0, 10E2, 2.5E2),
                expand = expansion(mult = c(0, 0.2))
                ) +
        coord_flip() +
        xlab(NULL) +
        theme(axis.title.x=element_blank(),
                legend.position = "none", 
                panel.grid.major.y = element_blank()
        )
```

### Figure 1B - right: monthly trend for each indicator
For the line graph, first, we add the real publication year and month to the datasets. The real publication year/month is the year/month the paper was first appear online stored in firstPublicationDate column.
```{r}
rct$YearMonth <- format(as.Date(rct$firstPublicationDate), "%Y-%m")
```


```{r}

monthly_rct <- 
        rct %>% 
        select(YearMonth,
               is_coi_pred,
               is_fund_pred,
               is_register_pred,
               is_open_data,
               is_open_code) %>%
        gather("indicator", "value", -YearMonth) %>%
        count(YearMonth, indicator, value) %>%
        mutate(indicator = recode(indicator,
                                  is_coi_pred = "COI Disclosure",
                                  is_fund_pred = "Funding disclosure",
                                  is_register_pred = "Protocol registration",
                                  is_open_data = "Data sharing",
                                  is_open_code = "Code sharing")) %>%
        complete(indicator, value, YearMonth, fill = list(n = 0)) %>%
        group_by(YearMonth, indicator) %>% 
        mutate(p = n / sum(n)) %>%
        filter(value) %>%
        ungroup()

fig1B_right <-
        monthly_rct %>% 
        ggplot() +
        aes(x = YearMonth, 
            y = p,
            group = indicator,
            color = indicator) +
        geom_line(size = 0.75) +
        scale_y_continuous(limits = c(0, 1), 
                           labels = scales::percent) +
        scale_color_discrete(name = NULL) +
        scale_fill_discrete(breaks = c("COI Disclosure",
                                       "Funding disclosure",
                                       "Protocol registration",
                                       "Data sharing",
                                       "Code sharing")) +
        labs(y = NULL, 
             x = NULL) +
        theme(panel.grid.minor = element_blank(),
              legend.position = "none",
              axis.text.x = element_blank(),
              legend.key.size = unit(0.2, 'cm'),
              )

fig1B <- ggarrange(fig1B_left, fig1B_right,
                    ncol = 2, nrow = 1, 
                    align = "hv", common.legend = F)

# tiff("Figure.tiff", width = 14, height = 7, units = "in", res = 300)
fig1B
#dev.off()

```


## Figure 1C
### Figure 1C - left: proportion of review papers that adhered to each transparency indicator
```{r}

prop_review <- review %>%
        summarise("COI Disclosure" = sum(is_coi_pred == TRUE),
                  "Funding disclosure" = sum(is_fund_pred == TRUE),
                  "Protocol registration" = sum(is_register_pred == TRUE),
                  "Data sharing" = sum(is_open_data == TRUE),
                  "Code sharing" = sum(is_open_code == TRUE)) %>%
        t() %>%
        as.data.frame() %>%
        rownames_to_column(var = "indicator") %>%
        mutate(percentage = round(V1/nrow(review)*100, 1))

fig1C_left <- prop_review %>% ggplot() +
        aes(x = reorder(indicator, V1),
            y = V1,
            fill = indicator,
            col = indicator) +
        geom_col(alpha = 0.75) +
        geom_text(aes(label = percentage), hjust = -0.1, size = 3) +
        scale_y_continuous(
                name = "\nNumber of papers",
                breaks = seq(0, 20E3, 5E3),
                labels = seq(0, 20E3, 5E3),
                expand = expansion(mult = c(0, 0.2))
                ) +
        coord_flip() +
        xlab(NULL) +
        theme(
                legend.position = "none", 
                panel.grid.major.y = element_blank()
        )
```

### Figure 1C - right: monthly trend for each indicator
For the line graph, first, we add the real publication year and month to the datasets. The real publication year/month is the year/month the paper was first appear online stored in firstPublicationDate column.
```{r}
review$YearMonth <- format(as.Date(review$firstPublicationDate), "%Y-%m")
```


```{r}

monthly_review <- 
        review %>% 
        select(YearMonth,
               is_coi_pred,
               is_fund_pred,
               is_register_pred,
               is_open_data,
               is_open_code) %>%
        gather("indicator", "value", -YearMonth) %>%
        count(YearMonth, indicator, value) %>%
        mutate(indicator = recode(indicator,
                                  is_coi_pred = "COI Disclosure",
                                  is_fund_pred = "Funding disclosure",
                                  is_register_pred = "Protocol registration",
                                  is_open_data = "Data sharing",
                                  is_open_code = "Code sharing")) %>%
        complete(indicator, value, YearMonth, fill = list(n = 0)) %>%
        group_by(YearMonth, indicator) %>% 
        mutate(p = n / sum(n)) %>%
        filter(value) %>%
        ungroup()

fig1C_right <-
        monthly_review %>% 
        ggplot() +
        aes(x = YearMonth, 
            y = p,
            group = indicator,
            color = indicator) +
        geom_line(size = 0.75) +
        scale_y_continuous(limits = c(0, 1), 
                           labels = scales::percent) +
        scale_color_discrete(name = NULL) +
        scale_fill_discrete(breaks = c("COI Disclosure",
                                       "Funding disclosure",
                                       "Protocol registration",
                                       "Data sharing",
                                       "Code sharing")) +
        labs(y = NULL, 
             x = "\nMonth") +
        theme(panel.grid.minor = element_blank(),
              legend.position = c(0.5, 0.4),
              axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1),
              legend.key.size = unit(0.05, 'cm')
              )

fig1C <- ggarrange(fig1C_left, fig1C_right,
                    ncol = 2, nrow = 1, 
                    align = "hv", common.legend = F)

# tiff("Figure.tiff", width = 14, height = 7, units = "in", res = 300)
fig1C
#dev.off()

```


## Final Figure 1

```{r}
fig1 <- ggarrange(
        ggarrange(fig1A_left, fig1A_right, ncol = 2),
        ggarrange(fig1B_left, fig1B_right, ncol = 2),
        ggarrange(fig1C_left, fig1C_right, ncol = 2),
        nrow = 3,
        labels = c("A - Research articles", "B - RCTs", "C - Reviews"),
        align = "hv", 
        common.legend = F,
        hjust = -0.05,
        vjust = 1.1
)

# tiff("Figure 1.tiff", width = 15, height = 10, units = "in", res = 300)
fig1
#dev.off()
```


# Appendix 1

Citations to article and journal impact factor by transparency practices for RCTs
Table S1A
```{r}
set.seed(100)

map(
  c("is_coi_pred", "is_fund_pred", "is_register_pred", "is_open_data", "is_open_code"),
  ~ tbl_summary(
      rct,
      by = all_of(.x),
      include = citedByCount,
      label = list(citedByCount = .x)
    ) %>% add_p()
) %>% tbl_stack()

```


Table S1B
```{r}
set.seed(100)

map(
  c("is_coi_pred", "is_fund_pred", "is_register_pred", "is_open_data", "is_open_code"),
  ~ tbl_summary(
      rct,
      by = all_of(.x),
      include = jif2020,
      label = list(jif2020 = .x)
    ) %>% add_p()
) %>% tbl_stack()

```

Citations to article and journal impact factor by transparency practices for reviews
Table S2A
```{r}
set.seed(100)

map(
  c("is_coi_pred", "is_fund_pred", "is_register_pred", "is_open_data", "is_open_code"),
  ~ tbl_summary(
      review,
      by = all_of(.x),
      include = citedByCount,
      label = list(citedByCount = .x)
    ) %>% add_p()
) %>% tbl_stack()

```


Table S2B
```{r}
set.seed(100)

map(
  c("is_coi_pred", "is_fund_pred", "is_register_pred", "is_open_data", "is_open_code"),
  ~ tbl_summary(
      review,
      by = all_of(.x),
      include = jif2020,
      label = list(jif2020 = .x)
    ) %>% add_p()
) %>% tbl_stack()

```

## Journals

Total number of papers published in journals with an SJR:
```{r}
table(is.na(opendata$sjr2020))
```

Top five journals in SJR:
```{r}
tapply(opendata$sjr2020, opendata$journalTitle, max) %>% as.data.frame() %>% arrange(desc(.)) %>% head(5)

opendata$sjr2020c = ifelse(opendata$sjr2020 < 26.304,"Other", opendata$journalTitle)
```


Table S1
```{r}
set.seed(100)
opendata %>% 
  select(sjr2020c,is_coi_pred,is_fund_pred,is_register_pred,is_open_data,is_open_code) %>%
  tbl_summary(by=sjr2020c, percent="column",
              label = c(is_coi_pred ~ "COI disclosure", is_fund_pred ~ "Funding disclosure", is_register_pred ~ "Protocol registration", is_open_data ~ "Data sharing", is_open_code ~ "Code sharing")) %>%
  add_p(test.args=all_tests("fisher.test") ~ list(simulate.p.value=TRUE)) %>%
        as_flex_table()
```


Top five journals in H-index:
```{r}
tapply(opendata$scimago_hindex, opendata$journalTitle, max) %>% as.data.frame() %>% arrange(desc(.)) %>% head(5)

opendata$hindexc = ifelse(opendata$scimago_hindex < 771,"Other", opendata$journalTitle)
```

Table S2
```{r}
set.seed(100)
opendata %>% 
  select(hindexc,is_coi_pred,is_fund_pred,is_register_pred,is_open_data,is_open_code) %>%
  tbl_summary(by=hindexc, percent="column",
              label = c(is_coi_pred ~ "COI disclosure", is_fund_pred ~ "Funding disclosure", is_register_pred ~ "Protocol registration", is_open_data ~ "Data sharing", is_open_code ~ "Code sharing")) %>%
  add_p(test.args=all_tests("fisher.test") ~ list(simulate.p.value=TRUE)) %>%
        as_flex_table()
```

## Appendix 3:
### Research articles:
#### Correct p-values (considering all levels)
```{r}
research %>% 
  select(journalTitle,is_coi_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
research %>% 
  select(journalTitle,is_fund_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
research %>% 
  select(journalTitle,is_register_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
research %>% 
  select(journalTitle,is_open_data) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
research %>% 
  select(journalTitle,is_open_code) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
```

### Publisher-related differences in transparency
Correct p-values (considering all levels)
```{r}
research %>% 
  select(scimago_publisher,is_coi_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
research %>% 
  select(scimago_publisher,is_fund_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
research %>% 
  select(scimago_publisher,is_register_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
research %>% 
  select(scimago_publisher,is_open_data) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
research %>% 
  select(scimago_publisher,is_open_code) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
```

## Appendix 4:
### Review articles:
#### Correct p-values (considering all levels)
```{r}
review %>% 
  select(journalTitle,is_coi_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
review %>% 
  select(journalTitle,is_fund_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
review %>% 
  select(journalTitle,is_register_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
review %>% 
  select(journalTitle,is_open_data) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
review %>% 
  select(journalTitle,is_open_code) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
```

### Publisher-related differences in transparency
Correct p-values (considering all levels)
```{r}
review %>% 
  select(scimago_publisher,is_coi_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
review %>% 
  select(scimago_publisher,is_fund_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
review %>% 
  select(scimago_publisher,is_register_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
review %>% 
  select(scimago_publisher,is_open_data) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
review %>% 
  select(scimago_publisher,is_open_code) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
```


## Appendix 5:
### Other articles:
#### Correct p-values (considering all levels)
```{r}
other %>% 
  select(journalTitle,is_coi_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
other %>% 
  select(journalTitle,is_fund_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
other %>% 
  select(journalTitle,is_open_data) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
other %>% 
  select(journalTitle,is_open_code) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
```

### Publisher-related differences in transparency
Correct p-values (considering all levels)
```{r}
other %>% 
  select(scimago_publisher,is_coi_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
other %>% 
  select(scimago_publisher,is_fund_pred) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
other %>% 
  select(scimago_publisher,is_open_data) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
other %>% 
  select(scimago_publisher,is_open_code) %>%
        table() %>%
        fisher.test(simulate.p.value=TRUE)
```