gcbb_manuscript.Rmd

---
title: "Analyses, figures, and tables for BETYdb manuscript"
output: github_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(cache = TRUE, autodep=TRUE, echo = FALSE, message = FALSE, results='hide')
```

```sh
git checkout tags/1.4.7 -b 1.4.7
./scripts/build.sh
```

## Trait Meta Analysis

1. Estimate parameters for multiple crops
2. Compare sources of variability Site, Treatment, Greenhouse
3. For one site / two species, look at fetilizer effects

## Assess potential productivity: 

* use meta-analysis estimates of traits to predict yield

### model validation

* compare BioCro predictions with observations

# Results

## Figure: Web Interface

* Data Entry Workflow
* Accessing Data
* Visualizing Model output

# Code

```{r load-packages}
library(betydb)
library(pander)
require(knitr)
library(data.table)
library(scales)
library(dplyr)


theme_set(theme_bw(base_size = 10, base_family = 'Times'))

```

## Summary of content

```{r record-counts, eval = FALSE}
library(dplyr)
library(PEcAn.utils)
d <- list(host = 'ebi-forecast.igb.illinois.edu',
          dbname = 'ebi_production',
          user = 'bety',
          password = 'bety')

bety <- src_postgres(host = d$host, user = d$user, password = d$password, dbname = d$dbname)

traits_yields_view <- tbl(
  bety, sql("select * from traits_and_yields_view_private where checked > -1")) %>% collect

```

```{r glopnet, eval = FALSE}


## bring in Glopnet
library(gdata)
library(dplyr)
#download.file("http://www.nature.com/nature/journal/v428/n6985/extref/nature02403-s2.xls", "wright-2004.xls")
glopnet.data <- read.xls("wright-2004.xls", skip = 7)

g <- glopnet.data %>% mutate(leaf_longevity = 12 * 10^(log.LL),  SLA =  1000/ (10^log.LMA), GS = 10^log.Gs,
                             leafN = 10^log.Nmass, leafP = 10^log.Pmass, leaf_respiration_rate_m2 = 10^log.Rdarea) %>% 
  dplyr::select(scientificname = Species, leaf_longevity, SLA, GS, leafN, leafP, leaf_respiration_rate_m2)

library(tidyr)
gg <- g %>% gather(trait, mean, -scientificname, na.rm = TRUE) %>% mutate(citation_id = 160, result_type = 'traits', author = 'Wright',access_level = 2, checked = 1)

```

```{r summary, eval=FALSE}
bety_glopnet <- tbl(
  bety, sql("select * from traits_and_yields_view_private where author = 'Wright'")) %>% collect 


traits_yields_view <- rbind(traits_yields_view, bety_glopnet) %>% data.table


data_summary <- traits_yields_view
save(data_summary, file ="data/data_summary.RData")
save(traits_yields_view, file ="data/traits_yields_view.RData")
write.csv(data_summary, file = 'data/data_summary.csv', row.names = FALSE)
write.csv(traits_yields_view, file = 'data/traits_yields_view.csv', row.names = FALSE)
```

```{r}
load('data/data_summary.RData')

.d <- data_summary[,list(ntraits = sum(grepl('traits', result_type)), nyields = sum(grepl('yield', result_type)), total = sum(!is.na(result_type))), by = genus][order(total, decreasing = TRUE)]

rbind(.d[!genus == "" ][1:55], .d[, list(genus = "Total", ntraits = sum(ntraits), nyields = sum(nyields), total = sum(total))]) %>% kable

rbind(.d[!genus == "" ][1:15], .d[, list(genus = "Total", ntraits = sum(ntraits), nyields = sum(nyields), total = sum(total))]) %>% kable



```

### Status of records

```{r record-status}
load('data/traits_yields_view.RData')
## Number of records checked

traits_yields_view[,list(n = length(n)),by='checked'] %>% kable
## Number of public records
traits_yields_view[,list(n = length(n)),by='access_level'] %>% kable

traits_yields_view[checked >=0 & (!author == 'Wright') & access_level < 4,list(n = length(mean)),by='name,author'] %>% kable
## Number of public bioenergy records

# Bioenergy Genera
genera <- c("Miscanthus", "Panicum", "Populus", "Saccharum", "Pinus", "Salix", "Agave", "Robinia", "Pinus")

### evaluate proportion public v private
.t <- traits_yields_view[genus %in% genera,list(n = length(n)),by='access_level,result_type'] 
.t %>% kable 
.t %>% group_by(access_level,result_type) %>% summarise(r = result_type,a = access_level,x = 100*n / sum(n))

traits_yields_view[genus %in% genera ,list(n = length(n)),by='access_level'] %>% kable
traits_yields_view[genus %in% genera, list(n = length(n)),by='access_level,result_type'] %>% kable

```

```{r biofuel_traits}

traits <-  c("Ayield", "leafN", "LAI", "SLA", "LMA", "root_respiration_rate", "stemN", "rootN", "root2shoot", "Vcmax", "leaf_turnover_rate", "Jmax", "c2n_leaf")
genera_x_traits <- expand.grid(genus = genera, trait = c(traits))

T <-  traits_yields_view %>% collect %>% setDT
trait_summary <- T[genus %in% genera & trait %in% traits, list(n = length(mean)), by = 'trait,genus']  


#data(trait.dictionary, package = "PEcAn.utils")
#write.csv(trait.dictionary, 'data/trait.dictionary.csv', row.names = FALSE)
trait.dictionary <- read.csv('data/trait.dictionary.csv')
trait.dictionary <- trait.dictionary %>% dplyr::select(trait = id, label = figid) %>% rbind(data.frame(trait = "Ayield", label = "Yield"))

trait_summary_full <- genera_x_traits  %>% 
  full_join(trait_summary, by = c('genus', 'trait')) %>%
  left_join(trait.dictionary, by = 'trait')%>%  
  filter(genus !='Agave')

```


```{r trait-plot}
#devtools::install_github("infotroph/ggplotTicks")
library(ggplotTicks)
gcbb_plot_theme <- function(p){ 
  require(ggplotTicks)
  newp <- p +
#  scale_y_continuous(limits = c(0, 60), breaks = c(0, 20, 40, 60)) +
  theme(text = element_text(family = 'Times', size = 10, color = 'black'),
        line = element_line(color = 'black'),
        axis.title = element_text(vjust = 0.1),
        panel.grid = element_blank(),
        axis.ticks.length = unit(-0.1, 'cm'),
        axis.text.x = element_text(size = 8, margin=margin(2,2,2,2,"mm")),
        axis.text.y = element_text(size = 8, margin=margin(2,2,2,2,"mm")),
        #panel.grid = element_blank(),
        #panel.border = element_blank(),
        strip.background = element_blank(), #http://docs.ggplot2.org/current/facet_wrap.html
        strip.text = element_text(size = 8))
  newp <- mirror_ticks(newp)
  return(newp)
}

trait_plot <- ggplot(data = trait_summary_full) +
  geom_bar(aes(genus, n, fill = label), stat = "identity", position = "dodge", width = 0.75) +
  scale_fill_manual(values=c(brewer_pal(palette="Paired")(11), "black")) +
  guides(fill = guide_legend(reverse=TRUE, title="Trait" )) +
  ylab("Number of Records") + xlab("") +
  scale_y_log10(breaks = c(1, 10, 100, 1000, 10000), limits = c(1,7500)) +
  scale_x_discrete(limits = c("Pinus", "Populus", "Salix", "Robinia", "Saccharum", "Panicum", "Miscanthus")) +
  coord_flip() +
  theme_bw()

trait_plot2 <- trait_plot +
  theme(text = element_text(family = 'Times', size = 10, color = 'black'),
        axis.title = element_text(vjust = 0.1),
        panel.grid = element_blank(),
        axis.ticks.length = unit(-0.1, 'cm'),
        axis.text.x = element_text(size = 8, margin=margin(5,5,5,5,"mm")),
        axis.text.y = element_text(size = 8, margin=margin(5,5,5,5,"mm")),
        axis.ticks.y = element_blank(),
        axis.line = element_line(colour = 'black')) +
  geom_vline(xintercept = 0) 

trait_plot2.1 <- trait_plot2 + theme(legend.key.size = unit(2, 'mm'),
                    legend.key = element_rect(size = 0),
                    legend.title = element_blank(), 
                    legend.position = c(0.85,0.4))

trait_plot3 <- mirror_ticks(trait_plot2.1)
trait_plot3
#print(trait_plot)
ggsave(plot = trait_plot3,
       filename = 'figures/trait_records_by_genus.svg', 
       width = 169, height = 100, dpi = 300, units = 'mm')


```

```{r trait-data-histogram}

T <- traits_yields_view %>% 
  filter(!is.na(mean) & checked >= 0) %>% setDT

leafN <- T %>% filter(trait == 'c2n_leaf') %>% mutate(trait = 'leafN', mean = 48/mean)
yield <- T %>% filter(trait == 'Ayield') %>% mutate(trait = 'Yield')
sla1 <- T %>% filter(trait == 'SLA')
sla2 <- T %>% filter(trait == 'LMA') %>% mutate(mean = 1/mean)
vmax <- T %>% filter(trait == 'Vcmax') %>% mutate(trait = 'Vmax')

T <- rbind(T, leafN, yield, sla1, sla2, vmax)

histogram_traits <- c("Yield", "leafN", "LAI", "SLA","Vmax", "Jmax")
T <- T %>% filter(trait %in% histogram_traits)

trait_labels <- list(Jmax = expression(atop(Jmax~~at~~25*degree*C,
                                             group("(",mu*mol~~CO[2]~~m^-2~~s^-1,")"))),
                     LAI = "Leaf Area Index\n(ratio)",
                     leafN = "Leaf Nitrogen\n(%)",
                     SLA = expression(atop(Specific~~Leaf~~Area,
                                           group("(",m^2~~kg^-1,")"))),
                     Vmax = expression(atop(Vcmax~~or~~Vpmax~~at~~25*degree*C,
                                            group("(",mu*mol~~CO[2]~~m^-2~~s^-1,")"))),
                     Yield =  expression(Yield~~"(Mg "*ha^"-1"*yr^"-1"*")"))

label_plots <- list()
for(l in names(trait_labels)){
  label_plots[[l]] <- qplot(1, xlab = trait_labels[[l]])
  
}
 

T <- T %>% mutate(genus = ifelse(genus %in% genera, genus, 'Other'))
T$genus <- factor(T$genus, 
                       levels = c("Agave", "Miscanthus", "Panicum", "Pinus", "Populus", "Robinia", "Saccharum", "Salix", "Other"),
                       labels = c("Agave", "Miscanthus", "Panicum (Switchgrass)", "Pinus (Pine)", "Populus (Poplar)", "Robinia (Locust)", "Saccharum (Sugarcane)", "Salix (Willow)", "Other"))


trait_hist <- ggplot(data = T) +
  geom_histogram(bins = 25, aes(mean, fill = genus, color = genus)) +
  scale_color_manual(values = c(brewer_pal(type = "div", palette = "Set1")(8), 'grey80')) +
  scale_fill_manual(values = c(brewer_pal(type = "div", palette = "Set1")(8), 'grey80')) +
  ylab("Count") +
  xlab("Observed mean") +
  facet_wrap(~trait, scales = 'free', switch = "x")


trait_hist2 <- trait_hist +  
  theme(text = element_text(family = 'Times', size = 10, color = 'black'),
        line = element_line(color = 'black'),
        axis.title = element_text(vjust = 0.1),
        panel.grid = element_blank(),
        axis.ticks.length = unit(-0.1, 'cm'),
        axis.text.x = element_text(margin=margin(2,2,2,2,"mm")),
        axis.text.y = element_text(margin=margin(2,2,2,2,"mm")),
        legend.key.size = unit(3, 'mm'),
        legend.key = element_rect(size = 0),
        legend.title = element_blank(),
        #panel.grid = element_blank(),
        #panel.border = element_blank(),
        strip.background = element_blank())

trait_hist3 <- mirror_ticks(trait_hist2, allPanels = TRUE)

ggplot() + 
  geom_text(data = trait_labels, aes(x=20, y=wt, 
      label = lab2), size = 3, parse = T)


bigplot <- cowplot::plot_grid(plotlist = c(list(trait_hist3), label_plots), ncol = 1)
ggsave(plot = bigplot, 
      filename = 'figures/trait_histograms_by_genus.svg', 
      width = 169, dpi = 300, units = 'mm')


```
## Trait Meta-analysis

```{r ma-data-import, echo=FALSE,message=FALSE,results='hide'}
## Import traits
library(PEcAn.all)
logger.setQuitOnSevere(severeQuits = FALSE)
settings <- read.settings("inputs/biofuelspecies.xml")

settings$pfts <- get.trait.data(pfts = settings$pfts, 
                                modeltype = "BIOCRO", 
                                dbfiles = settings$database$dbfiles, 
                                database = settings$database$bety,   
                                forceupdate = FALSE, 
                                trait.names = TRUE)

```


```{r trait-meta-analysis, eval = FALSE}
require(plyr)
settings$database$bety$write = TRUE
run.meta.analysis(pfts = settings$pfts, iterations = 10000, 
                  dbfiles = settings$database$dbfiles, 
                  database = settings$database$bety)  

get.parameter.samples(pfts=settings$pfts)

```

### Summarize Trait Data

```{r data-contents, eval = FALSE}
library(data.table)
tmpdata <- list()
settings$outdir <- file.path("~/dev/gcbb_manuscript/outputs/")
for (i in seq(settings$pfts)){
  settings$pfts[[i]]$outdir <- file.path("~/dev/gcbb_manuscript/outputs/pft", settings$pfts[[i]]$name)
  pft <- settings$pfts[[i]]
  tmpdata[[pft$name]] <- cbind(pft = pft$name, fread(file.path(pft$outdir, "trait.data.csv")))
}
biofuel.traitdata <- data.table(do.call("rbind", tmpdata))

setnames(biofuel.traitdata, "vname", "variable")

## convert leafN to c2n 
biofuel.traitdata[variable == "leafN", `:=` (mean = 48/mean, vname = 'c2n_leaf')]


pander::pander(biofuel.traitdata[,list(n.obs = sum(n, is.na(n), na.rm = TRUE), n.records = length(n),  
                                       mean = mean(mean), SE = sqrt(sum(stat^2, na.rm = TRUE))), 
                                 by = c("variable", "pft")][order(variable)],
               split.tables = 200, round = 2)

save(settings, biofuel.traitdata, file = 'data/data_contents.RData')
```

### Trait Meta-analysis summary

#### Define Priors

```{r trait-summary}
load('data/data_contents.RData')
priors <- NULL
for(pft in settings$pfts){
  pftname <- pft$name
  load(file.path(settings$outdir, "pft", pftname, "prior.distns.Rdata"))
  for(trait in rownames(prior.distns)){
    tmp <- do.call(paste0("q", prior.distns[trait,"distn"]), list(c(0.025, 0.5, 0.975), prior.distns[trait,]$parama, prior.distns[trait,]$paramb))
    priors <- rbind(priors, data.table(pft=gsub("biocro.", "", pftname), variable=trait, median = tmp[2], lcl=tmp[1], ucl=tmp[3]))
  
  }
}

```

### Traits

```{r meta-analysis-plot-dataframe}
library(cowplot)
library(jagstools)
ma.output <- list()
priors.list <- list()
for(pft in settings$pfts){
  load(file.path(pft$outdir, "trait.mcmc.Rdata"))
  ma.output[[pft$name]] <- trait.mcmc
  load(file.path(pft$outdir, "prior.distns.Rdata"))
  priors.list[[pft$name]] <- prior.distns
  
}

trait.summary <- NULL
prior.summary <- NULL
for(pft in settings$pfts){
  for(trait in c("leaf_respiration_rate_m2", "stomatal_slope.BB", "SLA", "Vcmax")){
    if(trait %in% names(ma.output[[pft$name]])){
      x <- ma.output[[pft$name]][[trait]]
      z <- signif(jagsresults(x = x, params = c('beta.o', 'beta.site', 'beta.trt', 'beta.ghs', 'sd.y', 'sd.site', 'sd.trt')), 4)
      d <- data.table(type = 'posterior', pft = pft$name, variable = trait, param = rownames(z), z) 
      trait.summary <- rbind(trait.summary, d)
      rm(x); rm(z)
    } 
    prior.distns <- priors.list[[pft$name]]
    x1 <-  do.call(paste0("r", prior.distns[trait,"distn"]), list(10000, prior.distns[trait,]$parama, prior.distns[trait,]$paramb))
    z1 <- cbind(type = 'prior', signif(data.table(mean = mean(x1), sd = sd(x1), t(quantile(x1, c(0.025, 0.25, 0.5, 0.75, 0.975)))), 4))
    rownames(z1) <- 'beta.o'
    prior.summary <- rbind(prior.summary, data.table(pft = pft$name, variable = trait, param = rownames(z1), z1))
  }
  rm(x1); rm(z1)
}

t <- rbind(trait.summary, prior.summary)
t[, `:=`(prior = type == 'prior')]
t[, `:=`(posterior = type == 'posterior')]

#vars <- c("SLA", "Vcmax",  "leaf_respiration_rate_m2", "stomatal_slope.BB")#, "c2n_leaf")#, "cuticular_cond", "extinction_coefficient_diffuse", "growth_respiration_coefficient", "Jmax")

vars <- c("SLA", "Vcmax",  "leaf_respiration_rate_m2")#, "stomatal_slope.BB", "c2n_leaf")#, "cuticular_cond", "extinction_coefficient_diffuse", "growth_respiration_coefficient", "Jmax")

```

```{r ma_plots, fig.width=8, fig.height=1.5}
ma.plots <- list()

ts1 <- t[!(pft %in% c('saof', "Miscanthus_x_giganteus")) & param == 'beta.o' & !(variable == 'stomatal_slope.BB') & !prior]

n_traits <- biofuel.traitdata[variable %in% vars & pft %in% c("misp", "pavi", "salix", 'populus'), list(n = length(mean)), by = 'variable,pft']

ts2 <- merge(ts1, n_traits, by = c('variable', 'pft'))

for(var in vars){
  ma.plots[[var]] <- ggplot(data = ts2[variable == var], aes(pft, y = mean)) + 
    coord_flip() +
    geom_pointrange(aes(ymin = mean-sd, ymax =  mean + sd)) +
    geom_text(aes(pft, y = mean, label = paste0("(", n, ")")), vjust = 2, family = 'Times') +
    expand_limits(x = 0, y = 0) + 
    coord_flip() + xlab("") +
    scale_x_discrete(breaks=c("salix", "populus", "pavi", "misp"),
                     labels=c("Willow", "Poplar", "Switchgrass", "Miscanthus")) +
    theme(legend.position="none")
  
}

ma.plots$SLA <- ma.plots$SLA + ylab(expression(atop("SLA", group("(",m^2~kg^-1,")")))) 
#ma.plots$Vcmax <- 

ma.plots$Vcmax <- ma.plots$Vcmax + ylab(expression(atop(Vcmax~~or~~Vpmax~~at~~25*degree*C,group("(",mu*mol~~CO[2]~~m^-2~~s^-1,")"))))

ma.plots$leaf_respiration_rate_m2 <- ma.plots$leaf_respiration_rate_m2 + ylab(expression(atop(Leaf~~Respiration~~Rate~~at~~25*degree*C,group("(",mu*mol~~CO[2]~~m^-2~~s^-1,")"))))
#ma.plots$stomatal_slope.BB <- ma.plots$stomatal_slope.BB + ylab("Stomatal Slope \n Ball Berry Model (ratio)")

ma.plots2 <- list()
for(i in 1:length(ma.plots)){
  ma.plots2[[i]] <- ma.plots[[i]] +
    theme(axis.text.x = element_text(size = 10, margin = unit(c(0.7,0,0,0), "cm")),
          axis.text.y = element_blank(), 
          axis.title.x = element_text(size = 10),
          axis.ticks.y = element_blank(), 
          axis.ticks.length = unit(-0.2, 'cm'),
          panel.grid = element_blank())
}

labels <- ggplot() + 
  geom_text(aes(x = factor(1:4), y = 0, label = as.factor(c("Willow", "Poplar", "Switchgrass", "Miscanthus"))), hjust = 1, vjust = -0.5, position = position_dodge(width = 1.5)) + 
  coord_flip() +
  scale_y_continuous(limits = c(0, 0))+ ylab('genus\n newline') +
  xlab(label = "  \n  ") +
  theme(axis.text.x = element_text(color = 'white', size = 16),
        axis.text.y = element_blank(),
        axis.title = element_text(color = 'white', size = 18),
        axis.ticks = element_line(color = 'white'),
        axis.line.x = element_line(color = 'white'),
        axis.line.y = element_line(color = 'white'),
        panel.grid = element_blank())

ma_plotlist <- c(list(labels = labels), ma.plots2)
ma_plotlist2 <- lapply(ma_plotlist, mirror_ticks)
ma_plotgrid <- plot_grid(plotlist = ma_plotlist2, nrow = 1)#, rel_widths = c(0.4, 0.2, 0.2, 0.2))

ggsave(plot = ma_plotgrid,
       filename = 'figures/ma_plotgrid.svg', 
       width = 225, height = 75, dpi = 300, units = 'mm')


```

```{r}

sla_nsites <- ts1[,list(nsites = sum(grepl("beta.site", param))), by = pft]
sla_sd <- ts1[variable == "SLA" & grepl('sd.site', param)]

sla_sd <- merge(sla_sd, sla_nsites, by = 'pft')

ggplot(data = sla_sd, aes(x = pft, y = mean)) + 
  coord_flip() +
  geom_pointrange(aes(ymin = `25%`, ymax = `75%`), size = 1.2) +
  geom_text(aes(label = paste0("(", nsites, ")")), vjust = -3, size = 5) +
  expand_limits(y = 0) + 
  coord_flip() + xlab("") + ylab("SD among sites")
  scale_x_discrete(breaks=c("salix", "populus", "pavi", "misp"),
                     labels=c("Willow", "Poplar", "Switchgrass", "Miscanthus")) +
  theme(axis.text = element_text(size = 16),
        axis.title = element_text(size = 18))


```
### Meta-analysis of Miscanthus and Switchgrass yields

#### Extract, subset Miscanthus and Switchgrass yield data

* removed data with no fert or planting date
* only Miscanthus and Switchgrass
* only >= 3 years after planting
* following Heaton et al 2008 http://onlinelibrary.wiley.com/doi/10.1111/j.1365-2486.2008.01662.x/pdf

**heaton model**

> yield = harvestdate + precip + gdd + species*(harvestdate + precip + gdd)


```{r}

detach("package:raster", unload=TRUE, force = TRUE)
detach("package:plyr", unload=TRUE, force = TRUE)
library(dplyr)
library(data.table)
library(ggplot2)
## connection to database
d <- list(host = 'ebi-forecast.igb.illinois.edu',
          dbname = 'ebi_production',
          user = 'bety',
          password = 'bety')
bety <- src_postgres(host = d$host, user = d$user, password = d$password, dbname = d$dbname)

## query and join tables
species <- tbl(bety, 'species') %>% 
  dplyr::select(id, scientificname, genus) %>% 
  dplyr::rename(specie_id = id)

sites <- tbl(bety, sql(
  paste("select id as site_id, st_y(st_centroid(sites.geometry)) AS lat,",
        "st_x(st_centroid(sites.geometry)) AS lon,",
        " sitename, city, country from sites"))
  )

citations <- tbl(bety, 'citations') %>%
  select(citation_id = id, author, year, title)

yields <- tbl(bety, 'yields') %>%
  select(id, date, mean, n, statname, stat, site_id, specie_id, treatment_id, citation_id, cultivar_id) %>% 
  left_join(species, by = 'specie_id') %>%
  left_join(sites, by = 'site_id') %>% 
  left_join(citations, by = 'citation_id')

managements_treatments <- tbl(bety, 'managements_treatments') %>%
  select(treatment_id, management_id)

treatments <- tbl(bety, 'treatments') %>% 
  dplyr::mutate(treatment_id = id) %>% 
  dplyr::select(treatment_id, name, definition, control)

managements <- tbl(bety, 'managements') %>%
  filter(mgmttype %in% c('fertilizer_N', 'fertilizer_N_rate', 'planting', 'irrigation')) %>%
  dplyr::mutate(management_id = id) %>%
  dplyr::select(management_id, date, mgmttype, level, units) %>%
  left_join(managements_treatments, by = 'management_id') %>%
  left_join(treatments, by = 'treatment_id') 


nitrogen <- managements %>% 
  filter(mgmttype == "fertilizer_N_rate") %>%
  select(treatment_id, nrate = level)

planting <- managements %>% filter(mgmttype == "planting") %>%
  select(treatment_id, planting_date = date)

planting_rate <- managements %>% filter(mgmttype == "planting") %>%
  select(treatment_id, planting_date = date, planting_density = level) 

irrigation <- managements %>% 
  filter(mgmttype == 'irrigation') 

irrigation_rate <- irrigation %>% 
  filter(units == 'mm', !is.na(treatment_id)) %>% 
  group_by(treatment_id, year = sql("extract(year from date)"), units) %>% 
  summarise(irrig.mm = sum(level)) %>% 
  group_by(treatment_id) %>% 
  summarise(irrig.mm.y = mean(irrig.mm))

irrigation_boolean <- irrigation %>%
  collect %>%   
  group_by(treatment_id) %>% 
  mutate(irrig = as.logical(mean(level))) %>% 
  select(treatment_id, irrig = irrig)

irrigation_all <- irrigation_boolean %>%
  full_join(irrigation_rate, copy = TRUE, by = 'treatment_id')

grass_yields <- yields %>% 
  filter(genus %in% c('Miscanthus', 'Panicum')) %>%
  left_join(nitrogen, by = 'treatment_id') %>% 
  #left_join(planting, by = 'treatment_id') %>% 
  left_join(planting_rate, by = 'treatment_id') %>% 
  left_join(irrigation_all, by = 'treatment_id', copy = TRUE) %>% 
  collect %>% 
  mutate(age = year(date)- year(planting_date),
         nrate = ifelse(is.na(nrate), 0, nrate),
         SE = ifelse(statname == "SE", stat, ifelse(statname == 'SD', stat / sqrt(n), NA)),
         continent = ifelse(lon < -30, 'united_states', ifelse(lon < 75, 'europe', 'asia'))) %>% 
  filter(!duplicated(.))

save(grass_yields, file = "~/dev/bety_manuscript/data/grass_yields.RData")
# load("~/dev/bety_manuscript/data/grass_yields.RData")
g <- grass_yields %>% setDT 
if(nrow(g[is.na(age) & age < 0])>0) stop()


ggplot(data = grass_yields, aes(lon,lat)) + 
  geom_point(aes(color = genus, size = mean), 
             alpha = 0.1)
```

### Differences in Data since Heaton et al 2004

```{r heaton-2004-data}

heaton_authors <- c("Beale", "Beuch", "Ercoli", "Himken", "Jorgensen", 
                    "Lewandowski", "Schwarz", "Muir", "Ma", "Reynolds", "Sanderson",
                    "Madakadze", "Clifton-Brown")
heaton_authors[!heaton_authors %in% unique(grass_yields$author)]
### hmmm....

grass_yields %>% filter(genus == 'Miscanthus' & lon < -20) %>% summarise(n = length(genus))

```


### Climate parms

#### Daymet

```{r daymet-metdata, eval=FALSE}

#install_bitbucket("khufkens/daymetr")
library(DaymetR)
library(lubridate)
library(data.table)
y <- grass_yields %>% filter(!is.na(site_id), !is.na(lat)) %>% 
  select(site_id, lat, lon, date) %>% 
  distinct %>% data.frame %>% data.table 

y[,`:=`(year = year(ymd(date)))]


daymetsites <- y[year %in% 1980:2014 & lon < -40 & lat > 0 ,
                 list(start = min(year), end = max(year)), by = 'lat,lon,site_id']
save(daymetsites, file = 'daymetsites.RData')
#system("rsync daymetsites.RData dlebauer:")

daymet <- list()
idx <- 0
for(siteid in unique(daymetsites$site_id)){
  idx <- idx+1 
  s <- daymetsites[site_id == siteid]
  file <- paste(s$site_id, s$start, paste0(s$end, ".csv"), sep = '_')
  
  if(!file.exists(file)){
    download.daymet(site = s$site_id, lat = s$lat, lon = s$lon, start_yr = s$start, end_yr = s$end)
  }
  
  if(file.info(file)$size > 0){
    met <- fread(file)
  } else {
    met <- data.table(year = NA, yday = NA, `tmax (deg c)` = NA, 
    `tmin (deg c)` = NA, `dayl (s)` = NA, `prcp (mm/day)` = NA, 
    `srad (W/m^2)` = NA, `swe (kg/m^2)` = NA, `vp (Pa)` = NA)
  }

  daymet[[idx]] <- cbind(site_id = siteid, met) 

  save(daymet, file = 'daymet.RData')
}

load("~/dev/bety_manuscript/data/daymet.RData")

daymet <- rbindlist(daymet)
setnames(daymet, 
         c("yday", "tmax (deg c)", "tmin (deg c)", "dayl (s)", "prcp (mm/day)", "srad (W/m^2)", "swe (kg/m^2)",  "vp (Pa)"),
         c('doy', 'tmax', 'tmin', 'day1', 'precip', 'srad', 'swe', 'vp'))
daymet <- daymet[!is.na(tmax)]         

gd <- daymet[,list(gd = ifelse( sum(tmin,tmax)/2 > 10, (tmax + tmin) / 2 - 10, 0)), by = 'site_id,year,doy']
gdd <- gd[,list(daymet_gdd = sum(gd)), by = 'site_id,year']
p <- daymet[,list(daymet_precip = sum(precip)), by = 'site_id,year']
summer_p <- daymet[doy %in% yday(ymd(paste0(year, "-04-01"))):yday(ymd(paste0(year,"-09-30"))), list(summer_precip = sum(precip)),by = 'site_id,year']

daymet_vars <- merge(gdd, p, by = c('site_id', 'year'))
daymet_vars <- merge(daymet_vars, summer_p, by = c('site_id', 'year'))

save(daymet_vars, file = "~/dev/bety_manuscript/data/daymet_vars.RData")
```

### CruNCEP

```{r cruncep-metdata, eval=FALSE}

y2 <- y[year %in% 1950:2015 ,list(start = min(year), end = max(year)), by = 'lat,lon,site_id']

### Cruncep goes only to 2010. needto get data from pre-2010
y2[end == 2011]$end <- 2010 
y2[start == 2011]$start <- 2010 

write.csv(y2, '~/dev/bety_manuscript/data/sites2.csv')
system("rsync ~/dev/bety_manuscript/data/sites2.csv biocluster:")

## on biocluster
library(data.table)
library(PEcAn.data.atmosphere)
library(PEcAn.all)
library(ncdf4)
library(udunits2)

sites <- fread("~/sites2.csv")
met.nc <- nc_open("/home/groups/ebimodeling/met/cruncep/all.nc")
sitemet <- list()
for(siteid in unique(sites$site_id)){
    logger.info("extracting data for site_id", siteid)
    sitemet[[as.character(siteid)]] <- 
        sites[site_id == siteid,
              load.cfmet(met.nc = met.nc, lat = lat, lon = lon,
                         start.date=paste0(start, "-01-01"),
                         end.date = paste0(end, "-12-31"))]
    save(sitemet, file = '~/sitemet.RData')
}

system("rsync -routi biocluster:sitemet.RData ~/dev/bety_manuscript/data/")


## on ebi-forecast
load("~/dev/bety_manuscript/data/sitemet.RData")
sitemet[["281"]] <- NULL

library(udunits2)
library(lubridate)
library(PEcAn.data.atmosphere)
library(data.table)
z <- lapply(sitemet, 
  function(x) {
    PEcAn.data.atmosphere::cfmet.downscale.subdaily(x, output.dt = 1)
  }
)

zz <- lapply(z,
             function(x){
               as.data.table(PEcAn.BIOCRO::cf2biocro(x))
               })

for(siteid in names(zz)){
  zz[[siteid]][, `:=`(site_id = as.numeric(siteid))]
}

cruncepmet <- rbindlist(zz)


gd <- cruncepmet[,list(gd = ifelse(mean(Temp) > 10, (max(Temp) + min(Temp)) / 2  - 10, 0)), by = 'site_id,year,doy']
gdd <- gd[,list(cruncep_gdd = sum(gd)), by = 'site_id,year']

p <- cruncepmet[,list(cruncep_precip = sum(precip)), by = 'site_id,year']
summer_p <- cruncepmet[doy %in% yday(ymd(paste0(year, "-04-01"))):yday(ymd(paste0(year,"-09-31"))), list(summer_precip = sum(precip)),by = 'site_id,year']

cruncep_vars <- merge(gdd, p, by = c("site_id", "year"))
cruncep_vars <- merge(cruncep_vars, summer_p, by = c("site_id", "year"))
save(cruncep_vars, file = "~/dev/bety_manuscript/data/cruncep_vars.RData")
```

#### Compare Cruncep and Daymet



```{r}
load("~/dev/bety_manuscript/data/cruncep_vars.RData")
load("~/dev/bety_manuscript/data/daymet_vars.RData")
```


```{r}
met <- cruncep_vars %>% inner_join(daymet_vars, by = c('site_id','year'))

ggplot(data = met, aes(cruncep_gdd, daymet_gdd), alpha = 0.5) + 
  geom_point() +   
  geom_smooth(method = 'lm', formula = y ~ x -1, size = 0.9) +
  scale_x_continuous(limits = c(0, 4500)) +
  scale_y_continuous(limits = c(0, 4500))

ggplot(data = met, aes(cruncep_precip, daymet_precip), alpha = 0.3) + 
  geom_point() + 
  geom_smooth(method = 'lm', formula = y ~ x -1, size = 0.9) + 
  scale_x_continuous(limits = c(0, 2100)) +
  scale_y_continuous(limits = c(0, 2100))
ggplot(data = met) + 
  geom_point(aes(summer_precip.x, summer_precip.y, group = site_id))

met$site <- as.factor(met$site_id)
met$year <- as.factor(met$year)

precip_model <- lm(cruncep_precip ~ daymet_precip -1, data = met)
gdd_model <- lm(cruncep_gdd ~ daymet_gdd -1, data = met)
summer_model <- lm(summer_precip.x ~ summer_precip.y -1, data = met)
summary(precip_model)
summary(gdd_model)
summary(summer_model)

```

#### Fill in missing Cruncep data w/ Daymet

Use linear model to describe relationship between met variables derived from predict Cruncep from Daymet. Then use these models to estimate what Cruncep would have generated based on values derived from DayMet (this extends dataset to end of 2011) 

```{r fill-missing-daymet}

setnames(daymet_vars, 'summer_precip', 'summer_precip.y')
daymet_vars[,`:=`(cruncep_precip = predict(precip_model, newdata = daymet_vars), 
                  cruncep_gdd = predict(gdd_model, newdata = daymet_vars),
                  summer_precip = predict(summer_model, newdata = daymet_vars))]
daymet_est <- daymet_vars[,list(site_id, year, cruncep_precip, cruncep_gdd, summer_precip.y)] %>% 
  anti_join(cruncep_vars, by = c('site_id','year'))

setnames(daymet_est, 'summer_precip.y', 'summer_precip')
met <- rbind(cruncep_vars, daymet_est)
save(met, file = "~/dev/bety_manuscript/data/met.RData")
```

## Yield vs Climate, Fert Analysis

#### Merge Met vars w/ yield, fertilization, and irrigation 

```{r merge-met-yield}
load("~/dev/bety_manuscript/data/met.RData")
load("~/dev/bety_manuscript/data/grass_yields.RData")
library(data.table)
library(lubridate)
library(dplyr)
gg <- grass_yields %>% filter(!is.na(site_id), !is.na(lat)) %>% 
  dplyr::select(site_id, lat, lon, date, genus, scientificname, lat, lon, sitename, mean, n, SE,
         nrate, planting_date, irrig, irrig.mm.y, age, continent) %>% 
  mutate(year = year(ymd(date))) %>% 
  left_join(met, by = c('site_id', 'year'))


library(lme4)

gg <- gg %>% filter(!duplicated(gg))
ggg <- gg %>% 
  mutate(gdd = cruncep_gdd, precip = cruncep_precip, summer_precip = summer_precip, nrate = nrate, year = year, abslat = abs(lat)) %>% 
  mutate(sitename = ifelse(grepl("EBI Energy farm", sitename), "EBI Energy farm", sitename),
         continent = ifelse(lon < -30, "North America", ifelse(lon > 100, "China", "Europe"))) 

gggg <- ggg %>% filter(gdd < 2900 & gdd > 970 & nrate <= 205 & age < 11 & gdd > 600 & age > 2)

## 'no site' = no site random effects, as with original Heaton paper
misp_lm <- lm(mean ~ nrate + summer_precip + gdd + age, data = gggg, 
              subset = genus == "Miscanthus")


continent_misp_lm <- lm(mean ~ continent* nrate + continent* summer_precip + continent* gdd + continent* age, data = gggg, 
              subset = genus == "Miscanthus" & lat < 100)

pavi_lm <- lm(mean ~ nrate + summer_precip + gdd + age, data = gggg, 
              subset = genus == "Panicum")
all_lm <- lm(mean ~ genus * nrate + genus * summer_precip + genus * gdd + genus * age, data = gggg)
misp_re <- lme4::lmer(mean ~ nrate + summer_precip + gdd  + age + (1 | sitename / year), 
                  data = gggg, subset = genus == "Miscanthus")
pavi_re <- lme4::lmer(mean ~ nrate + summer_precip + gdd  + age + (1 | sitename / year), 
                  data = gggg, subset = genus == "Panicum")

all_re <- lme4::lmer(mean ~ nrate + summer_precip + gdd  + age + genus:nrate + genus:summer_precip + genus:gdd  + genus:age + (1 | sitename / year), 
                  data = gggg)

misp_continent_re <- lme4::lmer(mean ~ nrate + summer_precip + gdd  + age + continent:nrate + continent:summer_precip + continent:gdd  + continent:age + (1 | sitename / year), 
                  data = gggg, subset = genus == 'Miscanthus')
pavi_continent_re <- lme4::lmer(mean ~ nrate + summer_precip + gdd  + age + continent:nrate + continent:summer_precip + continent:gdd  + continent:age + (1 | sitename / year), 
                  data = gggg, subset = genus == 'Panicum' & continent %in% c("Europe", "North America") )

library(texreg)
library(gridExtra)

## Heaton et al 2004 looked at summer (Apr-Sept) precip
## Wang et al 2010 looked at annual precip

```

#### Miscanthus Switchgrass Yield Meta-analysis

```{r misp-pavi-yield-ma}

#library(lmerTest)
summary(pavi_re)
summary(misp_re)

save(pavi_re, misp_re, file = '~/dev/gcbb_manuscript/data/misp_pavi_yield_ma.RData')
texreg(l = list(misp_re, pavi_re),single.row = TRUE,
       digits = 4, dcolumn = TRUE,
       custom.model.names = c('Miscanthus', 'Switchgrass'), custom.coef.names = c('Intercept', 'Nitrogen Fertilization', 'Summer Precipitation', 'Growing Degree Days', "Stand Age"))

my_coefs <- function(mod){
  a <- data.frame(coef(summary(mod)), stringsAsFactors = FALSE)
  b <- cbind(parameter = rownames(a), a, 
             p.z = round(2 * (1 - pnorm(abs(a[,'t.value']))), 4))
  return(b)
}

#pavi_lm_coef <- my_coefs(pavi_lm) # lm
pavi_re_coef <- my_coefs(pavi_re) # lme4 random effect of site
#misp_lm_coef <- my_coefs(misp_lm)
misp_re_coef <- my_coefs(misp_re)

require(knitr)
rbind(misp_re_coef,
  pavi_re_coef) %>% 
  kable('markdown', row.names = FALSE, digits = 4)
#  kable(format = 'latex', row.names = FALSE, digits = 4, booktabs = TRUE)

```

```{r using-effects-package}
library(effects)
pavi_effects <- lapply(allEffects(pavi_re), 
                       function(x) data.frame(genus = 'Panicum', as.data.frame(x)))
misp_effects <- lapply(allEffects(misp_re), 
                       function(x) data.frame(genus = 'Miscanthus', as.data.frame(x)))
           
eff <- list(nrate = rbind(pavi_effects$nrate, misp_effects$nrate),
            summer_precip = rbind(pavi_effects$summer_precip, misp_effects$summer_precip),
            gdd = rbind(pavi_effects$gdd, misp_effects$gdd),
            age = rbind(pavi_effects$age, misp_effects$age))
```

```{r yield_ma_plots, fig.width=8, fig.height=2}
p1 <-  ggplot(data = ggg, aes(x = nrate, color = genus)) +
  geom_point(aes(x = jitter(nrate, 20), y = mean), alpha = 0.25, size = 0.25) +
  geom_smooth(data = eff$nrate, aes(y = fit), method = 'lm', fullrange = FALSE, se = FALSE) +
  #geom_smooth(data = eff$nrate, aes(y = fit + se), method = 'lm', size= 0.5, linetype = 2, fullrange = FALSE) +
  #geom_smooth(data = eff$nrate, aes(y = fit - se), method = 'lm', size= 0.5, linetype = 2, fullrange = FALSE) + 
  ylab(expression(Yield~~"(Mg "*ha^"-1"*yr^"-1"*")")) +
  xlab(expression("Nitrogen Fertilization Rate"~~"(kg "*ha^"-1"*yr^"-1"*")")) + 
  xlim(0,250)

p2 <-  ggplot(data = ggg, aes(x = summer_precip, color = genus)) +
  geom_point(aes(y = mean), alpha = 0.25, size = 0.25) +
  geom_smooth(data = eff$summer_precip, aes(y = fit), method = 'lm', fullrange = FALSE, se = FALSE) +
  ylab(" ") +
  xlab(expression("Growing Season Precipitation"~~"(mm "*yr^"-1"*")"))

p3 <-  ggplot(data = ggg, aes(x = gdd, color = genus)) +
  geom_point(aes(y = mean), alpha = 0.25, size = 0.25) +
  geom_smooth(data = eff$gdd, aes(y = fit), method = 'lm', fullrange = FALSE, se = FALSE, size = 0.5) +
  ylab(expression(Yield~~"(Mg "*ha^"-1"*yr^"-1"*")")) + 
  scale_x_continuous(name = expression(paste("Growing Degree Days", ~~(degree*C))))

p4 <- ggplot(data = ggg, aes(x = age, color = genus)) +
  geom_point(aes(x = jitter(age, 0.5), y = mean), alpha = 0.25, size = 0.25) +
  geom_smooth(data = eff$age, aes(y = fit), method = 'lm', fullrange = FALSE, se = FALSE, size = 0.5) +
  ylab(" ") +
  scale_x_continuous(name = 'Stand Age (yr)', limits = c(0, 10), breaks = c(0,3,6,9)) 

p <- list(p1, p2, p3)

p <- lapply(p, function(x) { 
  x + theme(legend.position = 'none') +
    scale_color_discrete()})
p[[4]] <- p4 + 
  theme(legend.key.size = unit(3, 'mm'),
        legend.background = element_rect(color = NULL),
        legend.key = element_rect(size = 0),
        legend.title = element_blank(),
        legend.position = c(0.7,0.9))

pp <- lapply(p, gcbb_plot_theme)
ppp <- lapply(pp, function(x) x + scale_colour_brewer(palette = "Set1", labels = c('Miscanthus', 'Panicum (Switchgrass)')))

pppp <- cowplot::plot_grid(plotlist = ppp)
#ppp
cowplot::ggsave(plot = pppp,
      filename = 'figures/yield_ma_plots.svg', 
      width = 169, height = 100, dpi = 300, units = 'mm')


```