---
title: "isotope posterior analyses_90CI"
author: "Kim O'Keefe"
date: "2/28/2019"
output: html_document
---

## This file contains isotope analyses that accompany: 
## O’Keefe K, Nippert JB, McCulloh KA. 2019. Plant water uptake along a diversity gradient provides evidence for complementarity in hydrological niches. Oikos 10.1111/oik.06529
## Included are posterior comparisons for SIAR isotope mixing model output.

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
Set working directory and read input files (produced earlier by siar):

```{r}
setwd("")
# Choose files for "data_1_2" and "data_4_16" objects for combined month/individual month SIAR output comparisons
# Comparing SIAR output with all months combined
# data_1_2 <-read.table('Isotope_Posterior_All_1_2', header=TRUE)
# data_4_16 <-read.table('Isotope_Posterior_All_4_16', header=TRUE)

# Comparing SIAR output with June only
data_1_2 <-read.table('Isotope_Posterior_June_A', header=TRUE)
data_4_16 <-read.table('Isotope_Posterior_June_B', header=TRUE)

# Comparing SIAR output with July only
# data_1_2 <-read.table('Isotope_Posterior_July_A', header=TRUE)
# data_4_16 <-read.table('Isotope_Posterior_July_B', header=TRUE)

# Comparing SIAR output with August only
# data_1_2 <-read.table('Isotope_Posterior_Aug_A', header=TRUE)
# data_4_16 <-read.table('Isotope_Posterior_Aug_B', header=TRUE)

```

# helper functions to get posterior information

```{r}
posteriorsummary=function(x){
  return(list(posteriorpositive = sum(x>0)/length(x),
              interval90credibility = quantile(x, c(0.05,0.95)),
              median = median(x)))
}
# posteriorsummary((-100):101) # just to doulbe-check
```
```{r}
comparecolumns=function(x,y){
  return(list(diff = posteriorsummary(y-x),
              relativediff = posteriorsummary((y-x)/x),
              twofolddiff = posteriorsummary(log2(y/x))))
}
# comparecolumns(1:100,2:101) #Input vectors always need to be positive
```
```{r}
getcolumn = function(diversity, depth, species){
  depthcontribution = 3 # default: if "snow" or anything other than shallow or deep
  if (depth == "shallow"){ depthcontribution = 1
  } else {if (depth == "deep") { depthcontribution = 2}}
  speciesskip = 5 * (species -1)
  # skip 5 columns for each species before the one of interest
  firstdatafile = diversity %in% c("monoculture","2sp") # true if we should use
  # data_1_2, false if we should use data_4_16
  diversityskip = 0 # default: works for monoculture and 4sp treatment
  if (diversity %in% c("2sp", "16sp")){
    diversityskip = 25
  }
  index = depthcontribution + speciesskip + diversityskip
  if (firstdatafile){
    values = data_1_2[,index]
  } else { values = data_4_16[,index] }
  return(list(
    index = index,
    values = values,
    name = ifelse(firstdatafile, names(data_1_2)[index], names(data_4_16)[index])
    ))
}
# getcolumn("16sp", depth="deep", species=4) # to check
```

# loop through species comparisons

set up empty data frame that will contain the results:

```{r}
ncomparisons = 10*3*4 # 10 pairs of species, 3 depths, 4 diversity trts
speciescomp = data.frame(
  diversity=rep(NA,ncomparisons),
  depth = rep(NA, ncomparisons),
  species1 = rep(NA, ncomparisons),
  species2 = rep(NA, ncomparisons)
)
speciescomp$columnname1 = NA
speciescomp$columnname2 = NA
speciescomp$diff.posteriorpositive = NA
speciescomp$diff.interval90low = NA
speciescomp$diff.interval90high = NA
speciescomp$diff.median = NA
speciescomp$twofolddiff.posteriorpositive = NA
speciescomp$twofolddiff.interval90low = NA
speciescomp$twofolddiff.interval90high = NA
speciescomp$twofolddiff.median = NA
# add other column for relative difference
```

then go and run all these comparisons:

```{r}
irow = 1 # index of row for result data frame
for(diversity in c("monoculture","2sp","4sp","16sp")){ # "monoculture"){ #
  for(depth in c("shallow","deep","snow")){
    for(sp1 in 1:4){ # 4){ # 
      sp1data = getcolumn(diversity, depth, sp1)
      for(sp2 in (sp1+1):5){
        speciescomp$diversity[irow] = diversity
        speciescomp$depth[irow] = depth
        speciescomp$species1[irow] = sp1
        speciescomp$species2[irow] = sp2
        sp2data = getcolumn(diversity, depth, sp2)
        speciescomp$columnname1[irow] = sp1data$name
        speciescomp$columnname2[irow] = sp2data$name
        result = comparecolumns(sp1data$values,sp2data$values)
        speciescomp$diff.posteriorpositive[irow] = result$diff$posteriorpositive
        speciescomp$diff.interval90low[irow] = result$diff$interval90credibility[1]
        speciescomp$diff.interval90high[irow] = result$diff$interval90credibility[2]
        speciescomp$diff.median[irow] = result$diff$median
        speciescomp$twofolddiff.posteriorpositive[irow] = result$twofolddiff$posteriorpositive
        speciescomp$twofolddiff.interval90low[irow] = result$twofolddiff$interval90credibility[1]
        speciescomp$twofolddiff.interval90high[irow] = result$twofolddiff$interval90credibility[2]
        speciescomp$twofolddiff.median[irow] = result$twofolddiff$median
        irow = irow+1
      }
    }
  }
}
speciescomp
write.csv(speciescomp, "June_speciesComparisons_90CI.csv", quote=F, row.names=F)
```

# compare diversity with monoculture

```{r}
# create empty data frame to contain results
ncomparisons = 3*5*3 # 3 depths, 5 species, 3 non-monoculture diversity trts
divcomp = data.frame(
  depth = rep(NA, ncomparisons),
  species = rep(NA, ncomparisons),
  diversity=rep(NA,ncomparisons)
)
divcomp$mononame = NA
divcomp$diversename = NA
divcomp$diff.posteriorpositive = NA
divcomp$diff.interval90low = NA
divcomp$diff.interval90high = NA
divcomp$diff.median = NA
divcomp$twofolddiff.posteriorpositive = NA
divcomp$twofolddiff.interval90low = NA
divcomp$twofolddiff.interval90high = NA
divcomp$twofolddiff.median = NA
# add other column for relative difference

# now go for the loop
irow=1
for(depth in c("shallow","deep","snow")){
  for (species in 1:5){
    monodata = getcolumn("monoculture", depth, species)
    for(diversity in c("2sp","4sp","16sp")){
      divdata = getcolumn(diversity, depth, species)
      divcomp$depth[irow] = depth
      divcomp$species[irow] = species
      divcomp$diversity[irow] = diversity
      divcomp$mononame[irow] = monodata$name
      divcomp$diversename[irow] = divdata$name
      result = comparecolumns(monodata$values,divdata$values)
      divcomp$diff.posteriorpositive[irow] = result$diff$posteriorpositive
      divcomp$diff.interval90low[irow] = result$diff$interval90credibility[1]
      divcomp$diff.interval90high[irow] = result$diff$interval90credibility[2]
      divcomp$diff.median[irow] = result$diff$median
      divcomp$twofolddiff.posteriorpositive[irow] = result$twofolddiff$posteriorpositive
      divcomp$twofolddiff.interval90low[irow] = result$twofolddiff$interval90credibility[1]
      divcomp$twofolddiff.interval90high[irow] = result$twofolddiff$interval90credibility[2]
      divcomp$twofolddiff.median[irow] = result$twofolddiff$median
      irow = irow+1
    }}}
divcomp
write.csv(divcomp, "June_diversityComparisons_90CI.csv", quote=F, row.names=F)
```