3_index.figures.Rmd

---
title: "index.figures"
author: "Simon Topp"
date: "7/6/2018"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(RColorBrewer)
library(scales)
library(tidyverse)
library(gridExtra)
library(grid)
library(lattice)
library(viridis)
library(knitr)
library(kableExtra)
library(feather)
```

# Bring in the munged data set and do a little more cleaning

```{r}
# Read in munged index.  Filter out incomplete entries.
munge <- read_feather('out/IndexMungedFull.feather') %>%  
  filter(!is.na(Model))

#Do a little more munging
#Reorder Models
munge$Model <- factor(munge$Model, levels = c('Empirical', 'Semi.Empirical', 'Semi.Analytical', 'Machine.Learning', 'Mixed', 'Product'))

#Reorder Parameters
munge$Parameter <- factor(munge$Parameter, levels = c('Chl', 'TSS', 'Clarity', 'CDOM', 'Nutrients', 'Turbidity', 'Carbon.Other', 'Cyanobacteria', 'Sediments.Other', 'Metals', 'Chromaticity', 'Trophic', 'Other'), labels = c('Chl', 'TSS', 'Clarity', 'CDOM', 'Nutrients', 'Turbidity', 'Other Carbon', 'Cyanobacteria', 'Other Sediments', 'Metals', 'Chromaticity', 'Trophic State', 'Other'))

#Reorder Paramter Group
munge$Parameter.Grp <- factor(munge$Parameter.Grp, 
                              levels = c("Algae","Carbon", "Clarity", "Suspended.Sediment","Turbidity", "Nutrients", "Other"))

# Factor coverage
munge$Coverage.KM <- factor(munge$Coverage.KM, levels = c('0','1','2','3','4','5','6'), labels = c('unknown','10^1','10^2', '10^3', '10^4', '10^5', '10^6'))

# Reorder category
munge$Cat <- factor(munge$Cat, levels = c('Methods', 'Methods w/ Pattern App', 'Trends/Patterns', 'Applied'))

# Reorder waterbodies
munge$Waterbody <- factor(munge$Waterbody, levels = c('Lakes', 'Rivers', 'Deltas', 'Estuaries', 'Multiple'))

#Create Filtered Set With only 1 entry per paper.
filt <- munge %>%
  select(-c(Parameter, Landscape.var, Num.Models, Best.AC,Error.fit1, Error.fit2,Error.fit3,Error.val1,Error.val2,Error.val3,Fit.metric1, Fit.metric2,Fit.metric3,Val.metric1,Val.metric2,Val.metric3, Atm.Comp)) %>%
  distinct(ID, .keep_all = T)

```

# Figures from the actual index.

```{r}
### Four panel figure 

##### For version with mixed years as labels, create mixed year variable where
# early years are binned together.   Otherwise leave this commented out.
######

filt <- filt %>% mutate(Year.Bin = paste0(trunc(Year/5) * 5 - 1900,'-', trunc(Year/5) * 5 + 4 - 1900),
                        Year.Mixed = ifelse(Year < 2000, Year.Bin, Year),
                        Year.Mixed = factor(Year.Mixed, levels = unique(Year.Mixed[order(Year)], ordered = T)),
                        Year.Bin = paste0(trunc(Year/5) * 5,'-', trunc(Year/5) * 5 + 4))


munge <- munge %>% mutate(Year.Bin = paste0(trunc(Year/5) * 5 - 1900,'-', trunc(Year/5) * 5 + 4 - 1900),
                        Year.Mixed = ifelse(Year < 2000, Year.Bin, Year),
                        Year.Mixed = factor(Year.Mixed, levels = unique(Year.Mixed[order(Year)], ordered = T)),
                        Year.Bin = paste0(trunc(Year/5) * 5,'-', trunc(Year/5) * 5 + 4))


# Make function for individual plots
viz <- function(data, x, fill, stack, label){
        if(stack == 'stack'){
        plot <-ggplot(data, aes(x = x,fill = fill)) +
                geom_bar(position = 'fill', alpha = 0.8) +
                scale_y_continuous(labels = percent) +
                scale_fill_brewer(palette = 'Paired') +
                scale_x_discrete(breaks = c('75-79','85-89','95-99', paste0(seq(2001,2018,2)))) +
                theme_bw() + 
                labs(x = 'Years', y = 'Count', fill = label) + 
                theme(axis.text.x = element_text(angle = 90, 
                                                 hjust = 1,
                                                 vjust = 0.5),
                legend.text=element_text(size=9),
                axis.title = element_blank(),
                legend.position = 'top',
                legend.title=element_text(size=10),
                legend.key.size = unit(0.5,'cm')) +
                guides(fill = guide_legend(nrow = 2, title.position = 'top'))
        } 
        if(stack == 'fill'){
        plot <- ggplot(data, aes(x = x, fill = fill)) + 
        geom_density(position = 'fill', alpha = 0.6, adjust = 1/2) + 
        scale_y_continuous(labels = percent) +
        scale_x_continuous(breaks = pretty_breaks(n=10)) +
        scale_fill_brewer(palette = 'Paired') +
        theme_bw() + 
        labs(x = 'Year', y = 'Distribution', fill = label) + 
        theme(axis.text.x = element_text(angle = 90, 
                                         hjust = 1, 
                                         vjust = 0.5),
              legend.text=element_text(size=9),
              axis.title = element_blank(),
              legend.position = 'top',
              legend.title=element_text(size=10),
              legend.key.size = unit(0.5,'cm')) +
        guides(fill = guide_legend(nrow = 2, title.position = 'top'))
        }
      if(stack == 'bar'){
        plot <-ggplot(data, aes(x = x,fill = fill)) +
                geom_bar(position = 'stack', alpha = 0.8) +
                scale_fill_brewer(palette = 'Paired') +
                scale_x_discrete(breaks = pretty_breaks(n=8)) +
                theme_bw() + 
                labs(x = 'Years', y = 'Count', fill = label) + 
                theme(axis.text.x = element_text(angle = 90,
                                                 hjust = 1,
                                                 vjust = 0.5),
                legend.text=element_text(size=9),
                axis.title = element_blank(),
                legend.position = 'top',
                legend.title=element_text(size=10),
                legend.key.size = unit(0.5,'cm')) +
                guides(fill = guide_legend(nrow = 2, title.position = 'top'))
        } 
  return(plot)
  }


####Yearly Distributions of Publication Characteristics (Figure 5)

### Frequency Distribution Version
## Parameter Groups
p1 <- viz(munge, munge$Year, munge$Parameter.Grp, 'fill', 'Modelled Parameters') +
  scale_fill_brewer(breaks = levels(munge$Parameter.Grp), labels = c("Algae", "Carbon", "Clarity", "Suspended\nSediment", "Turbidity", "Nutrients", "Other"),
                     palette = 'Paired')

## Modelling Approach
p2 <- viz(filt, filt$Year, filt$Model, 'fill', 'Modelling Approach')

## Study Coverage
cov <- filt %>% filter(Coverage.KM != 'unknown')
p3 <- viz(cov, cov$Year, cov$Coverage.KM, 'fill', expression(paste('Coverage Scale (km'^'2',')')))

## Study Time Spand
p4 <- viz(filt, filt$Year, filt$Time.Bin, 'fill', 'Study Duration')

g <- grid.arrange(
  grobs = list(p1,p2,p3,p4),
  layout_matrix = rbind(c(1, 2),
                        c(3, 4)),
  left = 'Distribution',
  bottom = 'Year'
)
 
ggplot2::ggsave(filename='figures/YearDist.png',plot = g, width=6.5,height=6, units='in', dpi=250)


#############Stacked Bars Version Figure
## Parameter Groups
p1 <- viz(munge, as.factor(munge$Year.Mixed), munge$Cat, 'stack', 'Publication Focus') +
  scale_fill_brewer(palette = 'Paired') 

p1 <- viz(munge, as.factor(munge$Year.Mixed), munge$Parameter.Grp, 'stack', 'Modelled Parameters') +
  scale_fill_brewer(breaks = levels(munge$Parameter.Grp), labels = c("Algae", "Carbon", "Clarity", "Suspended\nSediment", "Turbidity", "Nutrients", "Other"),
                     palette = 'Paired') 

## Modelling Approach
filt <- filt %>% filter(Year > 1984)
p2 <- viz(filt, filt$Year.Mixed, filt$Model, 'stack', 'Modelling Approach')  +
  scale_fill_brewer(breaks = levels(munge$Model),
                     palette = 'Paired')

## Experiment with continueous Gradient
cc <- scales::seq_gradient_pal("light blue", "red", "Lab")(seq(0,1,length.out=5))
## Study Coverage
cov <- filt %>% filter(Coverage.KM != 'unknown')
levels(cov$Coverage.KM) = c("unknown", "< 10^3","< 10^3", "10^3 < 10^4", "10^3 < 10^5", "10^5 < 10^6", "> 10^6")
p3 <- viz(cov, as.factor(cov$Year.Mixed), cov$Coverage.KM, 'stack', expression(paste('Coverage Scale (km'^'2',')'))) + scale_fill_manual(values = cc)

## Experiment with continueous Gradient
cc2 <- scales::seq_gradient_pal("light blue", "red", "Lab")(seq(0,1,length.out=6))
## Study Duration
p4 <- viz(filt, as.factor(filt$Year.Mixed), filt$Time.Bin, 'stack', 'Study Duration') + scale_fill_manual(values = cc2)

g <- grid.arrange(
  grobs = list(p2,p3,p4),
  nrow = 1,
  #layout_matrix = rbind(c(2, NA),
  #                      c(3, 4)),
  left = 'Distribution',
  bottom = 'Year'
)

ggsave(filename='figures/YearDistStackYearMixedv3.png',plot = g, width=12,height=3, units='in', dpi=250)

viz(filt, as.factor(filt$Year.Mixed), filt$Cat, 'stack', 'Analysis Level')

#############Bar Version
## Parameter Group
p1 <- viz(munge, as.factor(munge$Year), munge$Parameter.Grp, 'bar', 'Modelled Parameters') +
  scale_fill_brewer(breaks = levels(munge$Parameter.Grp), labels = c("Algae", "Carbon", "Clarity", "Suspended\nSediment", "Turbidity", "Nutrients", "Other"),
                     palette = 'Paired')

## Mondelling Appraoch
p2 <- viz(filt, as.factor(filt$Year), filt$Model, 'bar', 'Modelling Approach')  +
  scale_fill_brewer(breaks = levels(munge$Model),
                     palette = 'Paired')

## Study Coverage
cov <- filt %>% filter(Coverage.KM != 'unknown')
p3 <- viz(cov, as.factor(cov$Year), cov$Coverage.KM, 'bar', expression(paste('Coverage Scale (km'^'2',')')))

## Study Duration
p4 <- viz(filt, as.factor(filt$Year), filt$Time.Bin, 'bar', 'Study Duration')

g <- grid.arrange(
  grobs = list(p1,p3,p2,p4),
  layout_matrix = rbind(c(1, 2),
                        c(3, 4)),
  left = 'Distribution',
  bottom = 'Year'
)
ggsave(filename='figures/YearDistBar.png',plot = g, width=6.5,height=6, units='in', dpi=250)


levels(filt$Cat)[levels(filt$Cat) == 'Applied'] = 'Water Quality Science'
## Causal fig (Figure 7)
p1 <- ggplot(filt, aes(x = Year, fill = Cat)) + 
        geom_density(position = 'fill', alpha = 0.6, adjust = 1/2) + 
        scale_y_continuous(labels = percent) +
        scale_x_continuous(breaks = pretty_breaks(n=10)) +
        scale_fill_brewer(palette = 'Paired') +
        #scale_fill_viridis_d() +
        theme_bw() + 
        labs(x = 'Year', y = 'Distribution', fill = 'Analysis Level') + 
        theme(axis.text.x = element_text(angle = 90, 
                                         hjust = 1, 
                                         vjust = 0.5),
              legend.text=element_text(size=7),
              axis.text = element_text(size = 7),
              legend.position = 'top',
              legend.title=element_text(size=9),
              legend.key.size = unit(0.5,'cm')) +
        guides(fill = guide_legend(nrow = 2, title.position = 'top'))

ggsave(filename='figures/CausalFigViridis.png', width=4,height=2, units='in', dpi=250)

check <- filt %>%
  mutate(check = 1) %>%
  group_by(Cat) %>%
  arrange(Year) %>%
  mutate(cumsum = cumsum(check)) %>%
  ungroup()

p2<-ggplot(check, aes(x = Year, y = cumsum,color = Cat)) + 
        geom_line() + 
        scale_x_continuous(breaks = pretty_breaks(n=10)) +
        scale_color_brewer(palette = 'Paired') +
        #scale_color_viridis_d() +
        theme_bw() + 
        labs(x = 'Year', y = 'Cumulative Sum', color = 'Analysis Level') + 
        theme(axis.text.x = element_text(angle = 90, 
                                         hjust = 1, 
                                         vjust = 0.5),
              legend.position = 'none',
              axis.text = element_text(size = 7))
              #legend.position = 'right',
              #legend.title=element_text(size=9),
              #legend.key.size = unit(0.5,'cm')) #+
        #guides(fill = guide_legend(nrow = 2, title.position = 'top'))

g <- grid.arrange(p1,p2, ncol = 1, heights = c(2,1.5))
ggsave(filename='figures/CausalFig2Panel.png',plot = g, width=3.5,height=5, units='in', dpi=250)


### Yearly Publications and Parameter Publications by Waterbody (Figure 3)

p1 <- ggplot(munge, aes(x= Parameter, y= ..count.., fill = Waterbody)) + 
  geom_bar()+
  scale_fill_manual(values = c('#a6cee3', '#fdbf6f', '#1f78b4', '#b2df8a', '#33a02c')) +
  labs(x = 'Parameter', fill = 'Waterbody\nType', y = 'Count') +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
              legend.text=element_text(size=7),
              legend.title=element_text(size=8),
              legend.key.size = unit(0.5,'cm'),
              #axis.title.y = element_blank(),
              axis.text=element_text(size=7),
               axis.title=element_text(size=8)
        )
        
p2 <- ggplot(filt, aes(x = Year.Bin, y = ..count.., fill = Waterbody)) + 
  geom_bar() +
  scale_fill_manual(values = c('#a6cee3', '#fdbf6f', '#1f78b4', '#b2df8a', '#33a02c')) +
  labs(y = 'Publication Count', x = 'Year') +
  labs(x = 'Time Period', y = 'Publication Count') +
  theme_bw()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        legend.position = 'none',
        axis.text=element_text(size=6),
        axis.title=element_text(size=8),
        axis.title.x = element_text(margin = margin(t = 18)))
              

g <- grid.arrange(p2, p1, ncol = 2, widths = c(1.5,2))

ggsave(filename='figures/YearParam2Panel.png',plot = g, width=6,height=2.8, units='in', dpi=250)
```

## Correlation table

```{r}
library(Hmisc)
library(GGally)
library(corrplot)

df <- filt %>% 
  mutate(Coverage.KM = as.numeric(Coverage.KM),
         Cat = as.numeric(Cat)) %>%
  dplyr::select(Pub.Year = Year, Study.Duration = Years, Study.Scale = Coverage.KM, Study.Category = Cat) %>% as.matrix()

corcos <- rcorr(df, type = 'pearson')
# Correlation Coefficients
write.csv(corcos$r, file = 'out/cormatrix.csv')

# P values
write.csv(corcos$P, file = 'out/corpval.csv')

## Correlation plot
corrplot.mixed(corcos$r, order="hclust", tl.pos = 'd',  
         p.mat = corcos$P, sig.level = c(0.1,0.05, 0.01), insig = "label_sig", pch.cex = 0.7)
```

# Spatial Figure (Figure 6)

```{r}
library(ggmap)
library(maptools)
library(maps)
library(sf)
library(rgdal)
library(magrittr)
library(viridis)
library(rnaturalearth)
# Convert filtered index into spatial object and remove studies without explicit lat/long
filt.sf <- filt %>% 
  filter(!is.na(Lat)) %>% 
  st_as_sf(.,coords=c('Long','Lat'),crs=4326) %>%
  st_transform(.,54030)

#Check world v2
world <- rnaturalearth::ne_download(scale = 110,
                                       type = "land",
                                       category = "physical",
                                       destdir = tempdir(),
                                       load = TRUE,
                                       returnclass = "sf") %>%
    st_transform(.,54030) %>%
    st_buffer(.,0)
  

# Pull in world map and reproject
data("wrld_simpl")
world = st_as_sf(wrld_simpl) %>%
    st_transform(.,54030) %>%
    st_buffer(.,0)

# Set figure options
fill <- 'gray60'
color <- NA

# Check bounding box
st_bbox(world)

## Category Spatial (not used)
ggplot() + geom_sf(data = world, aes(), size = 0.2, fill = fill, color = color) + 
  geom_sf(data = filt.sf, aes(colour = Cat), size = .2,alpha =0.8) + 
  coord_sf(ylim = c(-6000000, 8338578), xlim = c(-16500000, 16874713)) +
  theme_bw() +
  scale_colour_viridis_d(name = 'Analysis \n Level') +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.grid.major = element_line(color = 'white'),
        legend.position = 'none')

ggsave(filename='figures/CatSpatial.png', width=3,height=1.5, units='in', dpi=350)

# Publication Year Spatial
p1 <- ggplot() + geom_sf(data = world, aes(), fill= fill, color = color) + 
  geom_sf(data = filt.sf, aes(colour = Year),size = 0.7, alpha =0.9) + 
  coord_sf(ylim = c(-5500000, 8338578), xlim = c(-13000000, 16874713)) +
  theme_bw() +
  scale_colour_viridis_c(name = 'Publication\nYear') +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.grid.major = element_line(color = 'white'),
        plot.margin=margin(0,0,0,0, "mm"),
        legend.margin = margin(0,0,0,0, "mm"),
        legend.text=element_text(size=7),
              axis.title = element_blank(),
              legend.position = 'bottom',
              legend.title=element_text(size=9),
              legend.key.size = unit(0.5,'cm'))

#Study Duration
p2<- ggplot() + geom_sf(data = world, aes(), fill = fill, color = color) + 
  geom_sf(data = filt.sf, aes(colour = Time.Bin),size = 0.7, alpha =0.9) + 
  coord_sf(ylim = c(-5500000, 8338578), xlim = c(-13000000, 16874713)) +
  theme_bw() +
  scale_colour_viridis_d(name = 'Time Span') +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.grid.major = element_line(color = 'white'),
        plot.margin=margin(0,0,0,0, "mm"),
        legend.margin = margin(0,0,0,0, "mm"),
        legend.text=element_text(size=7),
              axis.title = element_blank(),
              legend.position = 'bottom',
              legend.title=element_text(size=9),
              legend.key.size = unit(0.5,'cm'))


#Study Coverage
p3 <-ggplot() + geom_sf(data = world, aes(), fill = fill, color = color) + 
  geom_sf(data = filt.sf %>% filter(Coverage.KM != 'unknown'), aes(colour = Coverage.KM),size = .7, alpha =0.9) +
  coord_sf(ylim = c(-5500000, 8338578), xlim = c(-13000000, 16874713)) +
  theme_bw() +
  scale_colour_viridis_d(name = expression(paste('Coverage\nScale (km'^'2',')'))) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        panel.grid.major = element_line(color = 'white'),
        legend.margin = margin(0,0,0,0, "mm"),
        legend.text=element_text(size=7),
              axis.title = element_blank(),
              legend.position = 'bottom',
              legend.title=element_text(size=9),
              legend.key.size = unit(0.5,'cm'))

g <- grid.arrange(p1,p2,p3, nrow = 3)

ggsave(filename='figures/SpatialSimpleWorld.png', plot = g, width=3.5,height=6, units='in', dpi=350)
```

# Figure category distribution plot

```{r}

## Extraneous Figs Distribution Plots

df <- filt %>%
  mutate(fig.total = Figs.m + Figs.v + Figs.t) %>%
  gather(c(Figs.m:Figs.t), key = 'fig.type', value = 'fig.count') %>%
  na.omit() %>%
  mutate(fig.type = factor(fig.type, levels = c('Figs.m', 'Figs.v', 'Figs.t'), labels = c('Background', 'Validation', 'Applied')))

ggplot(df, aes(x = Cat, y = fig.count/fig.total)) +
  geom_boxplot(aes(colour = fig.type)) + 
  scale_x_discrete(labels = c('Methods', 'Methods w/\nPattern App', 'Trends/\nPatterns', 'Applied')) +
  labs(x = 'Paper Category', y = 'Figure Proportion', 
       colour = 'Figure Type') #+ 
  theme(axis.text.x = element_text(angle = 90, 
                                   vjust = 0.5,
                                   hjust = 1))

ggsave('figures/FigCountCat.png', width = 5, height = 3, units = 'in', dpi = 300)
```


# Some general summary tables
```{r}
## Country summary
country.sum <- filt %>% 
  group_by(Location) %>% 
  summarize(scale = mean(as.numeric(Coverage.KM)),
            duration = mean(Years),
            year = mean(Year),
            count = n(),
            prop = count/length(filt$ID))

# Applied Paper Summary
causal.sum <- filt %>%
  filter(Cat == 'Applied') %>%
  group_by(Location, Cat) %>% 
  dplyr::summarize(scale = mean(as.numeric(Coverage.KM)),
            duration = mean(Years),
            year = mean(Year),
            count = n(),
            prop = count/length(filt$ID))


# Summary for period before and after 2008
papes.sum <- filt %>% 
  mutate(period = ifelse(Year < 2008, 1, 2),
         Papes = 1) %>%
  spread(Cat, Papes) %>%
  group_by(period) %>%
  dplyr::summarize(scale = mean(as.numeric(Coverage.KM)),
            duration = mean(Years),
            year = mean(Year),
            count = n(),
            prop = count/length(filt$ID),
            methods = sum(Methods, na.rm = T),
            methodsPatterns = sum(`Methods w/ Pattern App`, na.rm = T),
            patterns = sum(`Trends/Patterns`, na.rm = T),
            applied = sum(`Applied`, na.rm = T),
            m.p = methods/count,
            mp.p = methodsPatterns/count,
            p.p = patterns/count,
            a.p = applied/count,
            ls = length(Sensor[grepl(x = Sensor,pattern = 'Landsat') == T])/count)


```

# Analyze the Error Metrics

```{r}

## Super annoying way to make long error data frame

# Pull and organize fit errors
ef <- munge %>%
  mutate(model.id = row_number()) %>%
  select(ID, model.id,Year,Year.Bin, Parameter, Parameter.Grp, Model, Error.fit1:Error.fit4) %>%
  gather(Error.fit1:Error.fit4, key = metric.num, value = Error) %>%
  mutate(type = 'Fit')

ef$fit[ef$metric.num == 'Error.fit1'] <- 1
ef$fit[ef$metric.num == 'Error.fit2'] <- 2
ef$fit[ef$metric.num == 'Error.fit3'] <- 3
ef$fit[ef$metric.num == 'Error.fit4'] <- 4
ef <- ef %>% mutate(UniqueID = paste0(model.id,'.', fit))
 
efm <- munge %>%
  mutate(model.id = row_number(),
         type = 'Fit') %>%
  select(ID, model.id, type, Year, Parameter, Parameter.Grp, Fit.metric1:Fit.metric4) %>%
  gather(Fit.metric1:Fit.metric4, key = metric.num, value = Metric)
efm$fit[efm$metric.num == 'Fit.metric1'] <- 1
efm$fit[efm$metric.num == 'Fit.metric2'] <- 2
efm$fit[efm$metric.num == 'Fit.metric3'] <- 3
efm$fit[efm$metric.num == 'Fit.metric4'] <- 4
efm <- efm %>% mutate(UniqueID = paste0(model.id, '.',fit)) 

error.fit <- left_join(ef, efm %>% select(Metric, UniqueID), by = 'UniqueID')

## Pull and organize validation errors
ev <- munge %>%
  mutate(model.id = row_number(),
         type = 'Val') %>%
  select(ID, model.id,type, Year, Year.Bin, Parameter, Parameter.Grp, Error.val1:Error.val4) %>%
  gather(Error.val1:Error.val4, key = metric.num, value = Error)
ev$fit[ev$metric.num == 'Error.val1'] <- 1
ev$fit[ev$metric.num == 'Error.val2'] <- 2
ev$fit[ev$metric.num == 'Error.val3'] <- 3
ev$fit[ev$metric.num == 'Error.val4'] <- 4
ev <- ev %>% mutate(UniqueID = paste0(model.id,'.', fit))
 
evm <- munge %>%
  mutate(model.id = row_number(),
         type = 'Val') %>%
  select(ID, model.id,Year,type, Parameter, Parameter.Grp, Val.metric1:Val.metric4) %>%
  gather(Val.metric1:Val.metric4, key = val.num, value = Metric)
evm$fit[evm$val.num == 'Val.metric1'] <- 1
evm$fit[evm$val.num == 'Val.metric2'] <- 2
evm$fit[evm$val.num == 'Val.metric3'] <- 3
evm$fit[evm$val.num == 'Val.metric4'] <- 4
evm <- evm %>% mutate(UniqueID = paste0(model.id, '.',fit)) 

error.val <- left_join(ev, evm %>% select(Metric, UniqueID), by = 'UniqueID')

# Combine into master error dataframe
error.long <- bind_rows(error.fit,error.val) %>%
  mutate(Metric = trimws(Metric)) %>%
  select(-metric.num)

error.long$Metric[error.long$Metric == ''] = 'none'
error.long$Error[error.long$Error == ''] = 'none'


## Check out the results
unique(error.long$Metric)[order(unique(error.long$Metric))]

### Pull the R2 Values to make error time series figure
r <- error.long %>% 
  filter(Metric == 'R2')%>%
  mutate(Error = as.numeric(Error)) %>%
  #na.omit() %>%
  filter(Error < 1) #Removes 1 data entry error, remove later

# Error Time Series (Figure 4)
ggplot(r,aes(x = Year, y = Error, color = type)) + geom_point(size = .7) + 
  geom_smooth(se = F, size = .7, method = 'lm') + 
  scale_color_discrete(name = 'Error\nType', breaks = c('Fit', 'Val'), labels = c('Fit', 'Validation')) +
  scale_x_continuous(breaks = pretty_breaks(n = 12)) +
  labs(y = expression(paste('Coefficient of\nDetermination (R'^2,')'))) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, 
                                   hjust = 1),
        #legend.position = 'bottom',
        plot.margin = margin(.7,.7,.7,.7, unit = 'cm'))

ggsave('figures/R2linearv2.png', width = 4, height = 3, units = 'in')

## Error trends
lmFit <- lm(Error ~ Year, r %>% filter(type == 'Fit'))
lmVal <- lm(Error ~ Year, r %>% filter(type == 'Val'))
summary(lmFit)
summary(lmVal)
confint(lmFit, 'Year', level=0.95)
confint(lmVal, 'Year', level=0.95)
##Metric Counts
#see the various reported metrics
unique(error.long$Metric)[order(unique(error.long$Metric))]

### Go through error metrics and combine to get total counts for different metrics.

# Fit counts
r.f = length(grep(x = error.fit$Metric, pattern = 'r2', ignore.case = T))
rmse.f = length(grep(x = error.fit$Metric, pattern = c('rmse|nrms|rmsd'), ignore.case = T))
se.f = length(grep(x = error.fit$Metric, pattern =
                     c('^se$|^mse$|^nse$|^sdv$|^see$'), ignore.case = T))
re.f = length(grep(x = error.fit$Metric, pattern =
                     c('^\\%|^sre$|^rpd$|^sdv$|^see$|^re|Relative|$mre|mrad|^mpe|lre|^average'), ignore.case = T))

mae.f = length(grep(x = error.fit$Metric, pattern =
                     c('mare|mape|^mae|^madp'), ignore.case = T))
bias.f = length(grep(x = error.fit$Metric, pattern =
                     c('bias|mnb'), ignore.case = T))

#Val Counts
r.v = length(grep(x = error.val$Metric, pattern = 'r2', ignore.case = T))
rmse.v = length(grep(x = error.val$Metric, pattern = c('rmse|nrms|rmsd'), ignore.case = T))
se.v = length(grep(x = error.val$Metric, pattern =
                     c('^se$|^mse$|^nse$|^sdv$|^see$'), ignore.case = T))
re.v = length(grep(x = error.val$Metric, pattern =
                     c('^\\%|^sre$|^rpd$|^sdv$|^see$|^re|Relative|$mre|mrad|^mpe|lre|^average'), ignore.case = T))

mae.v = length(grep(x = error.val$Metric, pattern =
                     c('mare|mape|^mae|^madp'), ignore.case = T))
bias.v = length(grep(x = error.val$Metric, pattern =
                     c('bias|mnb'), ignore.case = T))

counts <- data.frame('Metric' = c(rep(c('R2', 'RMSE', 'Standard Error', 'Relative Error', 'Mean Absolute Error', 'Bias'),2)), 
                     'Count' = 
                       c(r.f,rmse.f, se.f, re.f, mae.f, bias.f, r.v, rmse.v, se.v, re.v, mae.v, bias.v),
                     'type' = c(rep('Fit',6),rep('Val',6)))


ggplot(counts, aes(x = Metric, y = Count, fill = type)) + geom_col(position = 'dodge') +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

## Collapse terms list to get general idea of the various ways each metric is reported

rmse.t = error.long$Metric[grep(x = error.long$Metric, pattern = c('rmse|nrms|rmsd'), ignore.case = T)] %>%
  unique()
se.t = error.long$Metric[grep(x = error.long$Metric,  pattern =
c('^se$|^mse$|^nse$|^sdv$|^see$'), ignore.case = T)] %>%
  unique()
re.t = error.long$Metric[grep(x = error.long$Metric, pattern = c('^\\%|^sre$|^rpd$|^sdv$|^see$|^re|Relative|$mre|mrad|^mpe|lre|^average'), ignore.case = T)] %>%
  unique()
mae.t = error.long$Metric[grep(x = error.long$Metric, pattern = c('mare|mape|^mae|^madp'), ignore.case = T)] %>%
  unique()
bias.t = error.long$Metric[grep(x = error.long$Metric, pattern = c('bias|mnb'), ignore.case = T)] %>%
  unique()

error.metrics <- data.frame('Metric' = c('RMSE', 'Standard Error', 'Relative error', 'MAE', 'BIAS'), 'Terms' = c(paste0(rmse.t, collapse = ', '), paste0(se.t, collapse = ', '), paste0(re.t, collapse = ', '), paste0(mae.t, collapse = ', '), paste0(bias.t, collapse = ', ')))

kable(error.metrics) %>%
  kable_styling()
```

# Some extraneous single panel figures 


```{r}
##### Extraneous single panel Figures

##### Just made out of curiousity, not used in the manuscript

viz(filt, filt$Years, filt$Coverage.KM, 'fill', 'Study\nDuration')

viz(filt, filt$Year, filt$Cat, 'fill', 'Study Category')

viz(filt, as.numeric(filt$Coverage.KM), filt$Parameter.Grp, 'fill', 'Study Category')
ggsave('figures/yearsCat.png', height = 3, width = 3, units = 'in', dpi = 300)

viz(filt, filt$Year, filt$Cat, 'fill', 'Study Category')

viz(filt, filt$Year, filt)

ggplot(filt, aes(x = Cat, y = as.numeric(Coverage.KM), colour = Parameter.Grp)) + geom_boxplot()

```


```{r}
x <- filt %>%
  group_by(Year) %>%
  summarise(count = n())
##Fit a trend line to the data
# Linear
lm <- lm(x$count~x$Year)
# Log
lm2 <- lm(log(x$count)~x$Year)
summary(lm2)
## Calculate doubling time
log(2)/lm2$coefficients[2]

ggplot(x) + geom_point(aes(x = Year, y = exp(lm2$fitted.values))) + geom_point(aes(x = Year, y = count), colour = 'red')

# Average Annual Increase
(x$count[x$Year == max(x$Year)]/x$count[x$Year == min(x$Year)])^(1/(max(x$Year) - min(x$Year))) - 1


#Check 1980-2012 to compare to Bornman and Mutz 2015
check <- x %>% filter(Year >1979, Year <2013)
(x$count[x$Year == 2012]/x$count[x$Year == 1980])^(1/33) - 1
lmcheck<- lm(log(check$count)~check$Year)
#Doubling time
log(2)/lmcheck$coefficients[2]


```