---
title: 'Wood preservatives utilizing low-value olive oil production by-products: Analysis'
author: "| ***********\n| Michael Burnard, University of Primorska Muzejski trg 2,
  6000 Koper, Slovenia; InnoRenew CoE, Livade 6, 6310 Izola, Slovenia\n| Matthew Schwarzkopf,
  University of Primorska Muzejski trg 2, 6000 Koper, Slovenia; InnoRenew CoE, Livade
  6, 6310 Izola, Slovenia\n| Viacheslav Tverezovskiy, BioComposites Centre Bangor
  University, LL57 2UW, Bangor, Gwynedd, UK\n| Andreas Treu, Norwegian Institute of
  Bioeconomy Research Postboks 115, NO-1431 Ås, Norway \n| Miha Humar, Department
  of Wood Science and Technology, Biotechnical Faculty, University of Ljubljana, Rožina
  dolina, Cesta VIII/34, 1000 Ljubljana\n| Andreja Kutnar, University of Primorska
  Muzejski trg 2, 6000 Koper, Slovenia; InnoRenew CoE, Livade 6, 6310 Izola, Slovenia\n|
  ***********\n"
date: "16 February 2018"
output:
  pdf_document:
    df_print: kable
    toc_depth: 5
  html_notebook: default
  html_document:
    df_print: paged
    toc_depth: '5'
---
&nbsp;
\raggedright
\pagebreak
\tableofcontents
\newpage   

# Introduction   

Pine and beech wood specimens were treated with three types of preservative treatments and exposed to leaching and mini block fungal decay tests. The preservative treatments were developed from by-products of olive oil production and were intended to:

* Reduce leaching of the preservative out of the wood 
* Increase hydrophobicity
* Impart some level of antimicrobial performance

The four treatments were:

* Lampante Oil (LO), non-modified
* Direct Maleinization of LO (DM)
* Glycerolysis Maleinization of LO (GM)
* Water / untreated (Control)


## Set up the environment   

Load these packages to perform the analysis and create the visualisations.
```{r message=FALSE}
#basic packages
library(tidyverse)
library(scales)
library(zoo)
library(gridExtra)
library(ggalt)
```

# Leaching

Leaching tests consisted of 4 wet-dry cycles. The mass and dimensions of the specimens were measured and then submersed in distilled water for 48 hours. After the wet period their dimensions were measured. They were then dried in a climate chamber at 40 C for 24 hours. After the drying period, the mass and dimensions were measured again. This procedure was completed for 4 cycles.

Within this analysis we will refer to the treatments as LO, DM, GM, and Control (water).

## Load data set   

Load the most recent data, and remove the first column (just row numbers). Then take a look at its structure.
```{r}
d <- read.csv("data/OOTLeaching1MS-150417_v1_r0.csv", stringsAsFactors=FALSE)
d <- d[complete.cases(d),] # only weights for now.
glimpse(d)
```

## Clean, organise, and modify the data.

First, organise the data and add some extra data of interest.

Add a nicely named identifier column, *Cycle* that is equivelent to dates, and blocking variable.
```{r}
d$Cycle <- ifelse(d$Date == "11/10/16", "Cycle-0", 
                 ifelse(d$Date == "11/15/16", "Cycle-1",
                       ifelse(d$Date == "11/18/16", "Cycle-2",
                             ifelse(d$Date == "11/21/16", "Cycle-3",
                                   ifelse(d$Date == "11/28/16", "Cycle-4", NA)))))
d$Block <- ifelse(is.na(d$Cycle), "B", "A")
d$Cycle <- ifelse(d$Date == "1/5/17", "Cycle-0",
                  ifelse(d$Date == "1/8/17", "Cycle-1",
                         ifelse(d$Date == "1/11/17", "Cycle-2",
                                ifelse(d$Date == "1/14/17", "Cycle-3",
                                       ifelse(d$Date == "1/19/17", "Cycle-4", d$Cycle)))))

```

* *dffPrev* and *perDffPrev* are mass differences from the previous cycle.
```{r}
d <- d %>% arrange(Block, Specimen) %>% group_by(Specimen) %>% 
           mutate(dffPrev = Mass - lag(Mass, default=Mass[1]), 
          perDffPrev = (Mass - lag(Mass,default=Mass[1]))/Mass[1])
```

Now, to calculate the weight difference from the original value for each specimen after each cycle.

* *origMass* is a column that contains the original mass for each specimen.

The na.locf() function is from the zoo package, and fills in the NA values with the value from the previous row.

* *totalDff* is the difference between the weight at the current cycle compared to the original mass
* *perDffOrig* is same value as a percent of the original mass
```{r}
d$origMass <- ifelse(d$Cycle == "Cycle-0", d$Mass, NA)
d$origMass <- na.locf(d$origMass)

d <- d %>% group_by(Specimen) %>% mutate(totalDff = Mass - origMass,
           perDffOrig = (((Mass - origMass) / origMass)))
```

A little clean up... Split the Specimen name column into seperate values for modeling and summaries.
```{r}
d <- d %>% separate(Specimen, c("Species", "Treatment", "Test", "SpecNo"),
                sep="-", remove=FALSE, fill="right", extra="warn")

d$SpecNo <- ifelse(is.na(d$SpecNo), d$Test, d$SpecNo)
d <- d[,-4]
```

Now the data is ready for visualisation and further analysis.

\newpage
## Visualisation   

First, we create a second data set with only the rows we need. Then, we do some summary statistics to help us annotate the charts we'll make.
```{r}
d2 <- d %>% filter(!Cycle == "Cycle-0")
anno_text <- d2 %>% group_by(Cycle) %>% summarise(muPer=round(mean(perDffOrig),5),
                                                  muTot=round(mean(totalDff),4))

anno_text$muPer <- paste("Mean Diff: ", as.character(anno_text$muPer), "%")
anno_text$muTot <- paste("Mean Diff: ", as.character(anno_text$muTot), "g")
```

### Mass changes as a percent of original mass   
```{r fig.height=8, fig.width=7.5}
ggplot(data=d2, aes(x=Specimen, y=perDffOrig)) + theme_bw() +
  geom_hline(yintercept=0, colour="#CCCCCC") +
  geom_lollipop(point.colour="#4233E7", point.size=2) +
  facet_grid(Cycle~.) +
  scale_y_continuous(labels=percent) +
  geom_text(data=anno_text, aes(x=12,y=.015, label=muPer)) +
  labs(y="Weight difference (%)", x=NULL,
       title="Percent weight difference after drying phase",
       subtitle="Difference = (Cycle # - Original) / Original") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
```
\newpage

### Mass changes in grams   
```{r fig.height=8, fig.width=7.5}
ggplot(d2, aes(x=Specimen, y=totalDff)) + theme_bw() +
  geom_hline(yintercept=0, colour="#CCCCCC") +
  geom_lollipop(point.colour="#4233E7", point.size=2) +
  facet_grid(Cycle~.) +
  geom_text(data=anno_text, aes(x=12,y=.2, label=muTot)) +
  labs(y="Weight difference (g)", x=NULL,
       title="Weight difference after drying phase",
       subtitle="Difference = Cycle # - Original") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
```
\newpage
This is a convenient summary of weight change by species and treatment from cycle to cycle.
```{r}
leach.sum <- d %>% group_by(Treatment,Species,Cycle) %>% summarise(meanDff = mean(totalDff))
ggplot(leach.sum, aes(x=Cycle, y=meanDff, colour=Treatment, group=Treatment)) + theme_bw() +
  geom_line(size=0.5) +
  geom_point(size = 3) +
  scale_colour_manual(values=c("#D62F1F","#1E2DD6", "#1ED63D", "#AB1ED6")) +
  facet_grid(.~Species) +
  labs(y = "Absolute weight change [g]", x="")
leach.sum %>% spread(Cycle, meanDff)
```

## Statistics and models   

First, lets get just the data we need and take a look at it. The data is for the final Cycle only, so the differences are the total difference from the original state to the final state.
```{r}
ds <- d %>% filter(Cycle=="Cycle-4") %>% select(Specimen, Species, Treatment, Block, Mass,
                                                totalDff, SpecNo, perDffOrig, Cycle)
glimpse(ds)
```


### Some summary data   

By treatment and species.
```{r}
ds %>% group_by(Species, Treatment) %>% summarise(T.med=median(totalDff), T.sd=sd(totalDff),
                                       P.med=median(perDffOrig), P.sd=sd(perDffOrig), 
                                       N.obs=n())
```

\newpage
### Significance tests   

Assumptions were not met for parametric tests, so we use non-parametric alternatives.
Seperate by species (does not make sense to compare them this way) -- there is no expectation they would behave similarly.
```{r}
ds.b <- ds %>% filter(Species == "B")
ds.p <- ds %>% filter(Species == "P")
ds.b$Treatment <- factor(ds.b$Treatment)
ds.p$Treatment <- factor(ds.p$Treatment)
ds.b$Block <- factor(ds.b$Block)
ds.p$Block <- factor(ds.p$Block)
```
#### Beech by Block   
No difference detected.
```{r}
wilcox.test(perDffOrig ~ Block, data=ds.b)
```
#### Beech by Block    
No difference detected.
```{r}
wilcox.test(perDffOrig ~ Block, data=ds.p)
```
#### Beech by Treatment   
Treatment does seem to have an effect for beech speciemens.
```{r}
kruskal.test(perDffOrig ~ Treatment, data=ds.b)
```
#### Pine by treatment
Treatment does seem to have an effect for pine speciemens.
```{r}
kruskal.test(perDffOrig ~ Treatment, data=ds.p)
```

No apparent differences between blocks as a whole. There seem to be clear differences between treatments within species.

Now we test for differences between treatments within species.

#### Beech treatment comparisons   
To determine if there were differences between treatments, the effect size, and confidence intervals for the effect size, we use Wilcoxon rank sum tests for each pairing. We save the results to a dataframe for further inspection. We repeat the procedure the pine leaching data, in the following section.
```{r}
bl.df <- data_frame(Treat1=as.character("treat1"), Treat2=as.character("treat2"),
                    Loc.shift = 0, CI.lwr = 0, CI.upr = 0, p.value = 0)
b.cl <- ds.b %>% filter(Treatment %in% c("Control", "LO"))
b.cg <- ds.b %>% filter(Treatment %in% c("Control", "GM"))
b.cd <- ds.b %>% filter(Treatment %in% c("Control", "DM"))
b.lg <- ds.b %>% filter(Treatment %in% c("LO", "GM"))
b.ld <- ds.b %>% filter(Treatment %in% c("LO", "DM"))
b.dg <- ds.b %>% filter(Treatment %in% c("DM", "GM"))
x <- wilcox.test(perDffOrig ~ Treatment, data=b.cl, conf.int=TRUE, conf.level=0.95, 
                 exact=TRUE)
bl.df <- bind_rows(bl.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(perDffOrig ~ Treatment, data=b.cg, conf.int=TRUE, conf.level=0.95,
                 exact=TRUE)
bl.df <- bind_rows(bl.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(perDffOrig ~ Treatment, data=b.cd, conf.int=TRUE, conf.level=0.95, 
                 exact=TRUE)
bl.df <- bind_rows(bl.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("DM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(perDffOrig ~ Treatment, data=b.lg, conf.int=TRUE, conf.level=0.95, 
                 exact=TRUE)
bl.df <- bind_rows(bl.df, data_frame(Treat1=as.character("GM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(perDffOrig ~ Treatment, data=b.ld, conf.int=TRUE, conf.level=0.95, 
                 exact=TRUE)
bl.df <- bind_rows(bl.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(perDffOrig ~ Treatment, data=b.dg, conf.int=TRUE, conf.level=0.95, 
                 exact=TRUE) 
bl.df <- bind_rows(bl.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
bl.df <- bl.df[-1,]
```

Examine the results and plot. We make some new columns for absolute values to simplify reporting and plotting.
```{r}
bl.df$AbsValue <- abs(bl.df$Loc.shift)
bl.df$AbsCI.upr <- ifelse(bl.df$Loc.shift < 0, bl.df$CI.lwr * -1, bl.df$CI.upr)
bl.df$AbsCI.lwr <- ifelse(bl.df$Loc.shift < 0, bl.df$CI.upr * -1, bl.df$CI.lwr)
# create a nice label for the chart
bl.df$Label <- ifelse(bl.df$Loc.shift > 1, paste(bl.df$Treat1, " - ", bl.df$Treat2),
                      paste(bl.df$Treat2, " - ", bl.df$Treat1))
bl.df$Significant <- ifelse(bl.df$p.value < 0.05, "Yes", "No") #identify sig. differences
bl.df <- bl.df %>% arrange(AbsValue) # order by the magnitude of the effect size (step 1)
levels <- bl.df$Label # step 2
bl.df$Label <- factor(bl.df$Label, levels=levels) #step 3
bl.df %>% select(Label, AbsValue, AbsCI.lwr, AbsCI.upr, Significant, p.value) #inspect

#plot
ggplot(data=bl.df, aes(x=Label, y=AbsValue, ymin=AbsCI.lwr, ymax=AbsCI.upr, 
                       colour=Significant)) + theme_bw() +
  geom_pointrange() +
  geom_hline(yintercept=0, linetype=2) +
  scale_colour_manual(values=c("#D62F1F","#1E2DD6"), 
                      name=expression(paste("Significant (", alpha,  " < 0.05)"))) +
  scale_y_continuous(labels=percent, limits=(c(-0.01, 0.05)), 
                     breaks=seq(-0.01, 0.05, 0.01)) +
  coord_flip() +
  labs(title="Leaching: Treatment differences for Beech",
       subtitle="Median difference bewtween samples, with 95 % CIs",
       x="Comparison", y="Median weight loss difference ( % )") +
  theme(legend.position = "bottom")

```



#### Pine treatment comparisons   
 ** need to check abs ci vs normal ci and compre to beech **
```{r}
pl.df <- data_frame(Treat1=as.character("treat1"), Treat2=as.character("treat2"),
                    Loc.shift = 0, CI.lwr = 0, CI.upr = 0, p.value = 0)
p.cl <- ds.p %>% filter(Treatment %in% c("Control", "LO"))
p.cg <- ds.p %>% filter(Treatment %in% c("Control", "GM"))
p.cd <- ds.p %>% filter(Treatment %in% c("Control", "DM"))
p.lg <- ds.p %>% filter(Treatment %in% c("LO", "GM"))
p.ld <- ds.p %>% filter(Treatment %in% c("LO", "DM"))
p.dg <- ds.p %>% filter(Treatment %in% c("DM", "GM"))
x <- wilcox.test(perDffOrig ~ Treatment, data=p.cl, conf.int=TRUE, conf.level=0.95, 
                 exact=TRUE)
pl.df <- bind_rows(pl.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(perDffOrig ~ Treatment, data=p.cg, conf.int=TRUE, conf.level=0.95,
                 exact=TRUE)
pl.df <- bind_rows(pl.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(perDffOrig ~ Treatment, data=p.cd, conf.int=TRUE, conf.level=0.95, 
                 exact=TRUE)
pl.df <- bind_rows(pl.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("DM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(perDffOrig ~ Treatment, data=p.lg, conf.int=TRUE, conf.level=0.95, 
                 exact=TRUE)
pl.df <- bind_rows(pl.df, data_frame(Treat1=as.character("GM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(perDffOrig ~ Treatment, data=p.ld, conf.int=TRUE, conf.level=0.95, 
                 exact=TRUE)
pl.df <- bind_rows(pl.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(perDffOrig ~ Treatment, data=p.dg, conf.int=TRUE, conf.level=0.95, 
                 exact=TRUE)
pl.df <- bind_rows(pl.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
pl.df <- pl.df[-1,]
```
Examine the Results and plot.
```{r}
pl.df$AbsValue <- abs(pl.df$Loc.shift)
pl.df$AbsCI.upr <- ifelse(pl.df$Loc.shift < 0, pl.df$CI.lwr * -1, pl.df$CI.upr)
pl.df$AbsCI.lwr <- ifelse(pl.df$Loc.shift < 0, pl.df$CI.upr * -1, pl.df$CI.lwr)
# create a nice label for the chart
pl.df$Label <- ifelse(pl.df$Loc.shift > 1, paste(pl.df$Treat1, " - ", pl.df$Treat2),
                      paste(pl.df$Treat2, " - ", pl.df$Treat1))
pl.df$Significant <- ifelse(pl.df$p.value < 0.05, "Yes", "No") #identify sig. differences
pl.df <- pl.df %>% arrange(AbsValue) # order by the magnitude of the effect size (step 1)
levels <- pl.df$Label # step 2
pl.df$Label <- factor(pl.df$Label, levels=levels) #step 3
pl.df %>% select(Label, AbsValue, AbsCI.lwr, AbsCI.upr, Significant, p.value) #inspect

ggplot(data=pl.df, aes(x=Label, y=AbsValue, ymin=AbsCI.lwr, ymax=AbsCI.upr, 
                       colour=Significant)) + theme_bw() +
  geom_pointrange() +
  geom_hline(yintercept=0, linetype=2) +
  scale_colour_manual(values=c("#D62F1F","#1E2DD6"), 
                      name=expression(paste("Significant (", alpha,  " < 0.05)"))) +
  scale_y_continuous(labels=percent, limits=(c(-0.045, 0.045)),
                     breaks=seq(-0.045, 0.045, 0.015)) +
  coord_flip() +
  labs(title="Leaching: Treatment differences for Pine",
       subtitle="Median difference bewtween samples, with 95 % CIs",
       x="Comparison", y="Median weight change difference ( % ) ") +
  theme(legend.position = "bottom")

```

\newpage
# Fungi Results   

## Read and prepare the data   

```{r}
fg <- read.csv("data/OOTFungi1MS_11112017_v1_r0.csv", stringsAsFactors=FALSE)
fg <- fg %>% separate(SpecimenID, c("Species", "Treatment", "Test", "SpecNo"),
                sep="-", remove=FALSE, fill="right", extra="warn")

fg$Test <- "MB" #mini block

fg$SpecNo <- ifelse(fg$Treatment == "kontrola" & fg$Species == "P", seq(1,40,1),
                    ifelse(fg$Treatment == "kontrola" & fg$Species == "B", seq(1,40,1),
                           fg$Number))
fg$EWMW_per <- fg$EndWeightMinusWoodWeight / fg$WoodWeight

fg$Treatment <- ifelse(fg$Treatment == "kontrola", "Control", fg$Treatment)

fg$MassLoss <- (fg$InitialDryWeight - fg$PostInnocDryWeight) #Mass loss in grams
fg$InocLongName <- ifelse(fg$Inoculant == "gt", "Gloeopyllum trabeum", 
                   ifelse(fg$Inoculant == "pl", "Pleurotus ostreatus",
                   ifelse(fg$Inoculant == "pm", "Rhodonia placenta",
                   ifelse(fg$Inoculant == "tv", "Trametes versicolor", "Zed"))))

fg.p16 <- fg[ which(fg$Duration == "16" & fg$Species == "P"),]
fg.b16 <- fg[ which(fg$Duration == "16" & fg$Species == "B"),]

glimpse(fg)
```

```{r}
fg$Treatment = as.factor(fg$Treatment)
fg$Treatment = factor(fg$Treatment, levels = c("Control", "LO", "DM", "GM"))
fg.b16$Treatment <- factor(fg.b16$Treatment, levels = c("Control", "LO", "DM", "GM"))
fg.p16$Treatment <- factor(fg.p16$Treatment, levels = c("Control", "LO", "DM", "GM"))

ggplot(data=fg, aes(y=MassLoss, x=Treatment, colour=Treatment)) + theme_bw() +
  geom_boxplot() +
  scale_colour_manual(values=c("#D62F1F","#1E2DD6", "#1ED63D", "#AB1ED6")) +
  facet_grid(Duration~Species) +
  labs(y="Mass Loss (g)")
ggplot(data=fg, aes(y=RelativeMassLoss, x=Treatment, colour=Treatment)) + theme_bw() +
  geom_boxplot() +
  scale_colour_manual(values=c("#D62F1F","#1E2DD6", "#1ED63D", "#AB1ED6")) +
  scale_y_continuous(labels=percent) +
  facet_grid(Duration~Species) +
  labs(y="Relative Mass Loss (%)")

ggplot(data=fg.p16, aes(y=RelativeMassLoss, x=Treatment)) + theme_bw() +
  geom_boxplot() +
  facet_grid(.~InocLongName) +
  scale_colour_manual(values=c("#D62F1F","#1E2DD6", "#1ED63D", "#AB1ED6")) +
  scale_y_continuous(limits=c(0,.7)) +
  labs(y="Relative Mass Loss") +
  theme(axis.text.x = element_text(size=9, angle=0, colour="#000000", vjust=.5),
        axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12),
        axis.title = element_text(size=9), axis.text = element_text(size = 9, colour="#000000"),
        axis.text.y = element_text(size=9, angle=0, colour="#000000", vjust=.5), 
        strip.background = element_blank(), strip.text = element_text(size=12, face = "italic"))

ggplot(data=fg.b16, aes(y=RelativeMassLoss, x=Treatment)) + theme_bw() +
  geom_boxplot() +
  facet_grid(.~InocLongName) +
  scale_colour_manual(values=c("#D62F1F","#1E2DD6", "#1ED63D", "#AB1ED6")) +
  scale_y_continuous(limits=c(0,0.7)) +
  labs(y="Relative Mass Loss") +
  theme(axis.text.x = element_text(size=9, angle=0, colour="#000000", vjust=.5),
        axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12),
        axis.title = element_text(size=9), axis.text = element_text(size = 9, colour="#000000"),
        axis.text.y = element_text(size=9, angle=0, colour="#000000", vjust=.5), 
        strip.background = element_blank(), strip.text = element_text(size=12, face = "italic"))
```
## Data summary   

Our indicator of choice is Relative Mass Loss (InitialDryWeight - PostInoculantDryWeight) / InititalDryWeight,
```{r}
fg %>% group_by(Species, Treatment) %>% 
  summarise(RelativeMassLoss.mu = mean(RelativeMassLoss), 
            RelativeMassLoss.med = median(RelativeMassLoss),
            RelativeMassLoss.sd = sd(RelativeMassLoss),n = n())
```

## Comparisons   
Here we compare the difference in *RelativeMassLoss* between treatments for each species and Inoculant.   

### Beech :: Inoculant PL   

First, we compare responses to the PL Inoculant for beech at 16 weeks.
```{r}
cat("Beech :: PL")
bpl <- fg.b16 %>% filter(Inoculant == "pl")
kruskal.test(bpl$RelativeMassLoss, as.factor(bpl$Treatment))
```
The p-value indicates there is some difference between treatments, so we check them all with wilcoxian rank sum tests.
```{r}
bpl.df <- data_frame(Treat1=as.character("treat1"), Treat2=as.character("treat2"),
                    Loc.shift = 0, CI.lwr = 0, CI.upr = 0, p.value = 0)
bpl.cl <- bpl %>% filter(Treatment %in% c("Control", "LO"))
bpl.cd <- bpl %>% filter(Treatment %in% c("Control", "DM"))
bpl.cg <- bpl %>% filter(Treatment %in% c("Control", "GM"))
bpl.ld <- bpl %>% filter(Treatment %in% c("LO", "DM"))
bpl.lg <- bpl %>% filter(Treatment %in% c("LO", "GM"))
bpl.dg <- bpl %>% filter(Treatment %in% c("DM", "GM"))

x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=bpl.cl, conf.int=TRUE,
                 conf.level=0.95, exact=TRUE)
bpl.df <- bind_rows(bpl.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=bpl.cd, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
bpl.df <- bind_rows(bpl.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("DM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=bpl.cg, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
bpl.df <- bind_rows(bpl.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=bpl.ld, conf.int=TRUE,
                 conf.level=0.95, exact=TRUE)
bpl.df <- bind_rows(bpl.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=bpl.lg, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
bpl.df <- bind_rows(bpl.df, data_frame(Treat1=as.character("GM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=bpl.dg, conf.int=TRUE, 
                 conf.level=0.95,exact=TRUE)
bpl.df <- bind_rows(bpl.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
bpl.df <- bpl.df[-1,]
```
Examine the results and plot.
```{r}
bpl.df$AbsValue <- abs(bpl.df$Loc.shift)
bpl.df$AbsCI.upr <- ifelse(bpl.df$Loc.shift < 0, bpl.df$CI.lwr * -1, bpl.df$CI.upr)
bpl.df$AbsCI.lwr <- ifelse(bpl.df$Loc.shift < 0, bpl.df$CI.upr * -1, bpl.df$CI.lwr)
# create a nice label for the chart
bpl.df$Label <- ifelse(bpl.df$Loc.shift > 1, paste(bpl.df$Treat1, " - ", bpl.df$Treat2),
                      paste(bpl.df$Treat2, " - ", bpl.df$Treat1))
bpl.df$Significant <- ifelse(bpl.df$p.value < 0.05, "Yes", "No") #identify sig. differences
bpl.df <- bpl.df %>% arrange(AbsValue) # order by the magnitude of the effect size (step 1)
levels <- bpl.df$Label # step 2
bpl.df$Label <- factor(bpl.df$Label, levels=levels) #step 3
bpl.df %>% select(Label, AbsValue, AbsCI.lwr, AbsCI.upr, Significant, p.value) #inspect

ggplot(data=bpl.df, aes(x=Label, y=AbsValue, ymin=AbsCI.lwr, ymax=AbsCI.upr, 
                        colour=Significant)) + theme_bw() +
  geom_pointrange() +
  geom_hline(yintercept=0, linetype=2) +
  scale_colour_manual(values=c("#D62F1F","#1E2DD6"), 
                      name=expression(paste("Significant (", alpha,  " < 0.05)"))) +
  scale_y_continuous(labels=percent, limits=(c(-0.5, 0.5)), breaks=seq(-0.5, 0.5, 0.1)) +
  coord_flip() +
  labs(title="PL Inoculant, Beech",
       subtitle="Median difference bewtween samples, with 95 % CIs",
       x="Comparison", y="Median Relative Mass Loss difference ( % ) ") +
  theme(legend.position = "bottom")

```
   
### Beech :: Inoculant TV   

```{r}
btv <- fg.b16 %>% filter(Inoculant == "tv")
kruskal.test(btv$RelativeMassLoss, as.factor(btv$Treatment))
```
The p-value here suggests a difference between at least two of the treatments, so we check them all.
```{r}
btv.df <- data_frame(Treat1=as.character("treat1"), Treat2=as.character("treat2"),
                    Loc.shift = 0, CI.lwr = 0, CI.upr = 0, p.value = 0)
btv.cl <- btv %>% filter(Treatment %in% c("Control", "LO"))
btv.cd <- btv %>% filter(Treatment %in% c("Control", "DM"))
btv.cg <- btv %>% filter(Treatment %in% c("Control", "GM"))
btv.ld <- btv %>% filter(Treatment %in% c("LO", "DM"))
btv.lg <- btv %>% filter(Treatment %in% c("LO", "GM"))
btv.dg <- btv %>% filter(Treatment %in% c("DM", "GM"))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=btv.cl, conf.int=TRUE,
                 conf.level=0.95, exact=TRUE)
btv.df <- bind_rows(btv.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=btv.cd, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
btv.df <- bind_rows(btv.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("DM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=btv.cg, conf.int=TRUE, 
                 conf.level=0.95,  exact=TRUE)
btv.df <- bind_rows(btv.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=btv.ld, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
btv.df <- bind_rows(btv.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=btv.lg, conf.int=TRUE, 
                 conf.level=0.95,exact=TRUE)
btv.df <- bind_rows(btv.df, data_frame(Treat1=as.character("GM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=btv.dg, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
btv.df <- bind_rows(btv.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
btv.df <- btv.df[-1,]
```
Examine the results and plot.
```{r}
btv.df$AbsValue <- abs(btv.df$Loc.shift)
btv.df$AbsCI.upr <- ifelse(btv.df$Loc.shift < 0, btv.df$CI.lwr * -1, btv.df$CI.upr)
btv.df$AbsCI.lwr <- ifelse(btv.df$Loc.shift < 0, btv.df$CI.upr * -1, btv.df$CI.lwr)
# create a nice label for the chart
btv.df$Label <- ifelse(btv.df$Loc.shift > 1, paste(btv.df$Treat1, " - ", btv.df$Treat2),
                      paste(btv.df$Treat2, " - ", btv.df$Treat1))
btv.df$Significant <- ifelse(btv.df$p.value < 0.05, "Yes", "No") #identify sig. differences
btv.df <- btv.df %>% arrange(AbsValue) # order by the magnitude of the effect size (step 1)
levels <- btv.df$Label # step 2
btv.df$Label <- factor(btv.df$Label, levels=levels) #step 3
btv.df %>% select(Label, AbsValue, AbsCI.lwr, AbsCI.upr, Significant, p.value) #inspect

ggplot(data=btv.df, aes(x=Label, y=AbsValue, ymin=AbsCI.lwr, ymax=AbsCI.upr,
                        colour=Significant)) + theme_bw() +
  geom_pointrange() +
  geom_hline(yintercept=0, linetype=2) +
  scale_colour_manual(values=c("#D62F1F","#1E2DD6"), 
                      name=expression(paste("Significant (", alpha,  " < 0.05)"))) +
  scale_y_continuous(labels=percent, limits=(c(-0.5, 0.5)), breaks=seq(-0.5, 0.5, 0.1)) +
  coord_flip() +
  labs(title="TV Inoculant, Beech",
       subtitle="Median difference bewtween samples, with 95 % CIs",
       x="Comparison", y="Median Relative Mass Loss difference ( % ) ") +
  theme(legend.position = "bottom")

```
   
### Pine :: Inoculant GT    

```{r}
pgt <- fg.p16 %>% filter(Inoculant == "gt")
kruskal.test(pgt$RelativeMassLoss, as.factor(pgt$Treatment))
```

The p-value suggests we need to check for differences between pairs of treatments.

```{r}
pgt.df <- data_frame(Treat1=as.character("treat1"), Treat2=as.character("treat2"),
                    Loc.shift = 0, CI.lwr = 0, CI.upr = 0, p.value = 0)
pgt.cl <- pgt %>% filter(Treatment %in% c("Control", "LO"))
pgt.cd <- pgt %>% filter(Treatment %in% c("Control", "DM"))
pgt.cg <- pgt %>% filter(Treatment %in% c("Control", "GM"))
pgt.ld <- pgt %>% filter(Treatment %in% c("LO", "DM"))
pgt.lg <- pgt %>% filter(Treatment %in% c("LO", "GM"))
pgt.dg <- pgt %>% filter(Treatment %in% c("DM", "GM"))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=pgt.cl, conf.int=TRUE,
                 conf.level=0.95, exact=TRUE)
pgt.df <- bind_rows(pgt.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=pgt.cd, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
pgt.df <- bind_rows(pgt.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("DM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=pgt.cg, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
pgt.df <- bind_rows(pgt.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=pgt.ld, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
pgt.df <- bind_rows(pgt.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=pgt.lg, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
pgt.df <- bind_rows(pgt.df, data_frame(Treat1=as.character("GM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=pgt.dg, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
pgt.df <- bind_rows(pgt.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
pgt.df <- pgt.df[-1,]
```
Examine the results and plot.
```{r}
pgt.df$AbsValue <- abs(pgt.df$Loc.shift)
pgt.df$AbsCI.upr <- ifelse(pgt.df$Loc.shift < 0, pgt.df$CI.lwr * -1, pgt.df$CI.upr)
pgt.df$AbsCI.lwr <- ifelse(pgt.df$Loc.shift < 0, pgt.df$CI.upr * -1, pgt.df$CI.lwr)
# create a nice label for the chart
pgt.df$Label <- ifelse(pgt.df$Loc.shift > 1, paste(pgt.df$Treat1, " - ", pgt.df$Treat2),
                      paste(pgt.df$Treat2, " - ", pgt.df$Treat1))
pgt.df$Significant <- ifelse(pgt.df$p.value < 0.05, "Yes", "No") #identify sig. differences
pgt.df <- pgt.df %>% arrange(AbsValue) # order by the magnitude of the effect size (step 1)
levels <- pgt.df$Label # step 2
pgt.df$Label <- factor(pgt.df$Label, levels=levels) #step 3
pgt.df %>% select(Label, AbsValue, AbsCI.lwr, AbsCI.upr, Significant, p.value) #inspect

ggplot(data=pgt.df, aes(x=Label, y=AbsValue, ymin=AbsCI.lwr, ymax=AbsCI.upr, 
                        colour=Significant)) + theme_bw() +
  geom_pointrange() +
  geom_hline(yintercept=0, linetype=2) +
  scale_colour_manual(values=c("#D62F1F","#1E2DD6"), 
                      name=expression(paste("Significant (", alpha,  " < 0.05)"))) +
  scale_y_continuous(labels=percent, limits=(c(-0.5, 0.5)), breaks=seq(-0.5, 0.5, 0.1)) +
  coord_flip() +
  labs(title="GT Inoculant, Pine",
       subtitle="Median difference bewtween samples, with 95 % CIs",
       x="Comparison", y="Median Relative Mass Loss difference ( % ) ") +
  theme(legend.position = "bottom")

```
   
### Pine :: Inoculant PM   

```{r}
ppm <- fg.p16 %>% filter(Inoculant == "pm")
kruskal.test(pgt$RelativeMassLoss, as.factor(pgt$Treatment))
```

The p-value suggests we need to check for differences between pairs of treatments.
```{r}
ppm.df <- data_frame(Treat1=as.character("treat1"), Treat2=as.character("treat2"),
                    Loc.shift = 0, CI.lwr = 0, CI.upr = 0, p.value = 0)
ppm.cl <- ppm %>% filter(Treatment %in% c("Control", "LO"))
ppm.cd <- ppm %>% filter(Treatment %in% c("Control", "DM"))
ppm.cg <- ppm %>% filter(Treatment %in% c("Control", "GM"))
ppm.ld <- ppm %>% filter(Treatment %in% c("LO", "DM"))
ppm.lg <- ppm %>% filter(Treatment %in% c("LO", "GM"))
ppm.dg <- ppm %>% filter(Treatment %in% c("DM", "GM"))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=ppm.cl, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
ppm.df <- bind_rows(ppm.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=ppm.cd, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
ppm.df <- bind_rows(ppm.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("DM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=ppm.cg, conf.int=TRUE,
                 conf.level=0.95, exact=TRUE)
ppm.df <- bind_rows(ppm.df, data_frame(Treat1=as.character("Control"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=ppm.ld, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
ppm.df <- bind_rows(ppm.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=ppm.lg, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
ppm.df <- bind_rows(ppm.df, data_frame(Treat1=as.character("GM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(RelativeMassLoss~as.factor(Treatment), data=ppm.dg, conf.int=TRUE, 
                 conf.level=0.95, exact=TRUE)
ppm.df <- bind_rows(ppm.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
ppm.df <- ppm.df[-1,]
```
Examine the results and plot.
```{r}
ppm.df$AbsValue <- abs(ppm.df$Loc.shift)
ppm.df$AbsCI.upr <- ifelse(ppm.df$Loc.shift < 0, ppm.df$CI.lwr * -1, ppm.df$CI.upr)
ppm.df$AbsCI.lwr <- ifelse(ppm.df$Loc.shift < 0, ppm.df$CI.upr * -1, ppm.df$CI.lwr)
ppm.df$Label <- ifelse(ppm.df$Loc.shift > 1, paste(ppm.df$Treat1, " - ", ppm.df$Treat2),
                      paste(ppm.df$Treat2, " - ", ppm.df$Treat1))
ppm.df$Significant <- ifelse(ppm.df$p.value < 0.05, "Yes", "No") #identify sig. differences
ppm.df <- ppm.df %>% arrange(AbsValue) # order by the magnitude of the effect size (step 1)
levels <- ppm.df$Label # step 2
ppm.df$Label <- factor(ppm.df$Label, levels=levels) #step 3
ppm.df %>% select(Label, AbsValue, AbsCI.lwr, AbsCI.upr, Significant, p.value) #inspect

ggplot(data=ppm.df, aes(x=Label, y=AbsValue, ymin=AbsCI.lwr, ymax=AbsCI.upr,  
                        colour=Significant)) + theme_bw() +
  geom_pointrange() + 
  geom_hline(yintercept=0, linetype=2) +
  scale_colour_manual(values=c("#D62F1F","#1E2DD6"), 
                      name=expression(paste("Significant (", alpha,  " < 0.05)"))) +
  scale_y_continuous(labels=percent, limits=(c(-0.4, 0.4)), breaks=seq(-0.4, 0.4, 0.1)) +
  coord_flip() +
  labs(title="PM Inoculant, Pine",
       subtitle="Median difference bewtween samples, with 95 % CIs",
       x="Comparison", y="Median Relative Mass Loss difference ( % ) ") +
  theme(legend.position = "bottom")

```

## Oil Extraction   

Here we examine the amount of oil extracted by solvent extraction.
```{r}
oe <- fg %>% filter(!SpecimenID %in% c("B-kontrola", "P-kontrola"))

ggplot(data=oe, aes(x=Treatment, y=OilExtract_rel, colour=Treatment)) + theme_bw() +
  geom_boxplot() +
  scale_colour_manual(values=c("#D62F1F","#1E2DD6", "#1ED63D")) +
  facet_grid(~Species) +
  scale_y_continuous(labels=percent) +
  labs(y = "Relative Oil Extracted (%)")
```

## Fungal block test comparisons   

For each species we compare the results of the fungal block tests by treatment.   


### Beech

```{r}
fbb <- oe %>% filter(Species == "B")
kruskal.test(fbb$OilExtract_rel, as.factor(fbb$Treatment))
```
The Kruskal-Wallis test suggests a difference between some of the treatments. Now we need to compare each treatment individually. The following tests produce warnings because of ties, so exact p-values are not reported. (Warnings squelched).
```{r warning=FALSE}
fbb.df <- data_frame(Treat1=as.character("treat1"), Treat2=as.character("treat2"),
                    Loc.shift = 0, CI.lwr = 0, CI.upr = 0, p.value = 0)
fbb.ld <- fbb %>% filter(Treatment %in% c("LO", "DM"))
fbb.lg <- fbb %>% filter(Treatment %in% c("LO", "GM"))
fbb.dg <- fbb %>% filter(Treatment %in% c("DM", "GM"))

x <- wilcox.test(OilExtract_rel~as.factor(Treatment), data=fbb.ld, conf.int=TRUE, 
                 conf.level=0.95)
fbb.df <- bind_rows(fbb.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(OilExtract_rel~as.factor(Treatment), data=fbb.lg, conf.int=TRUE, 
                 conf.level=0.95)
fbb.df <- bind_rows(fbb.df, data_frame(Treat1=as.character("GM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(OilExtract_rel~as.factor(Treatment), data=fbb.dg, conf.int=TRUE, 
                 conf.level=0.95)
fbb.df <- bind_rows(fbb.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
fbb.df <- fbb.df[-1,]
fbb.df
```

```{r}
fbb.df$AbsValue <- abs(fbb.df$Loc.shift)
fbb.df$AbsCI.upr <- ifelse(fbb.df$Loc.shift < 0, fbb.df$CI.lwr * -1, fbb.df$CI.upr)
fbb.df$AbsCI.lwr <- ifelse(fbb.df$Loc.shift < 0, fbb.df$CI.upr * -1, fbb.df$CI.lwr)
fbb.df$Label <- ifelse(fbb.df$Loc.shift > 1, paste(fbb.df$Treat1, " - ", fbb.df$Treat2),
                      paste(fbb.df$Treat2, " - ", fbb.df$Treat1))
fbb.df$Significant <- ifelse(fbb.df$p.value < 0.05, "Yes", "No") #identify sig. differences
fbb.df <- fbb.df %>% arrange(AbsValue) # order by the magnitude of the effect size (step 1)
levels <- fbb.df$Label # step 2
fbb.df$Label <- factor(fbb.df$Label, levels=levels) #step 3
fbb.df %>% select(Label, AbsValue, AbsCI.lwr, AbsCI.upr, Significant, p.value) #inspect

ggplot(data=fbb.df, aes(x=Label, y=AbsValue, ymin=AbsCI.lwr, ymax=AbsCI.upr,  
                        colour=Significant)) + theme_bw() +
  geom_pointrange() +
  geom_hline(yintercept=0, linetype=2) +
  scale_colour_manual(values=c("#1E2DD6"), 
                      name=expression(paste("Significant (", alpha,  " < 0.05)"))) +
  scale_y_continuous(labels=percent, limits=(c(-0.4, 0.4)), breaks=seq(-0.4, 0.4, 0.1)) +
  geom_hline(yintercept=0, linetype=2) +
  coord_flip() +
  labs(title="Oil Extracted (Relative quantity), Beech",
       subtitle="Median difference bewtween samples, with 95 % CIs",
       x="Comparison", y="Median Oil Extracted Quantity ( % ) ") +
  theme(legend.position = "bottom")
```

### Pine     

```{r}
fbp <- oe %>% filter(Species == "P")
kruskal.test(fbp$OilExtract_rel, as.factor(fbp$Treatment))
```
The Kruskal-Wallis test suggests a difference between some of the treatments. Now we need to compare each treatment individually.
The following tests produce warnings because of ties, so exact p-values are not reported. (Warnings squelched).
```{r warning=FALSE}
fbp.df <- data_frame(Treat1=as.character("treat1"), Treat2=as.character("treat2"),
                    Loc.shift = 0, CI.lwr = 0, CI.upr = 0, p.value = 0)
fbp.ld <- fbp %>% filter(Treatment %in% c("LO", "DM"))
fbp.lg <- fbp %>% filter(Treatment %in% c("LO", "GM"))
fbp.dg <- fbp %>% filter(Treatment %in% c("DM", "GM"))

x <- wilcox.test(OilExtract_rel~as.factor(Treatment), data=fbp.ld, conf.int=TRUE, 
                 conf.level=0.95)
fbp.df <- bind_rows(fbp.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(OilExtract_rel~as.factor(Treatment), data=fbp.lg, conf.int=TRUE, 
                 conf.level=0.95)
fbp.df <- bind_rows(fbp.df, data_frame(Treat1=as.character("GM"), 
                                     Treat2=as.character("Lampante"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
x <- wilcox.test(OilExtract_rel~as.factor(Treatment), data=fbp.dg, conf.int=TRUE, 
                 conf.level=0.95)
fbp.df <- bind_rows(fbp.df, data_frame(Treat1=as.character("DM"), 
                                     Treat2=as.character("GM"), 
                                     Loc.shift=x$estimate, 
                                     CI.lwr=x$conf.int[1], CI.upr=x$conf.int[2], 
                                     p.value=x$p.value))
fbp.df <- fbp.df[-1,]
fbp.df
```

```{r}
fbp.df$AbsValue <- abs(fbp.df$Loc.shift)
fbp.df$AbsCI.upr <- ifelse(fbp.df$Loc.shift < 0, fbp.df$CI.lwr * -1, fbp.df$CI.upr)
fbp.df$AbsCI.lwr <- ifelse(fbp.df$Loc.shift < 0, fbp.df$CI.upr * -1, fbp.df$CI.lwr)
fbp.df$Label <- ifelse(fbp.df$Loc.shift > 1, paste(fbp.df$Treat1, " - ", fbp.df$Treat2),
                      paste(fbp.df$Treat2, " - ", fbp.df$Treat1))
fbp.df$Significant <- ifelse(fbp.df$p.value < 0.05, "Yes", "No") #identify sig. differences
fbp.df <- fbp.df %>% arrange(AbsValue) # order by the magnitude of the effect size (step 1)
levels <- fbp.df$Label # step 2
fbp.df$Label <- factor(fbp.df$Label, levels=levels) #step 3
fbp.df %>% select(Label, AbsValue, AbsCI.lwr, AbsCI.upr, Significant, p.value) #inspect

ggplot(data=fbp.df, aes(x=Label, y=AbsValue, ymin=AbsCI.lwr, ymax=AbsCI.upr,  
                        colour=Significant)) + theme_bw() +
  geom_pointrange() +
  geom_hline(yintercept=0, linetype=2) +
  scale_colour_manual(values=c("#1E2DD6"), 
                      name=expression(paste("Significant (", alpha,  " < 0.05)"))) +
  scale_y_continuous(labels=percent, limits=(c(-0.4, 0.4)), breaks=seq(-0.4, 0.4, 0.1)) +
  geom_hline(yintercept=0, linetype=2) +
  coord_flip() +
  labs(title="Oil Extracted (Relative quantity), Pine",
       subtitle="Median difference bewtween samples, with 95 % CIs",
       x="Comparison", y="Median Oil Extracted Quantity ( % ) ") +
  theme(legend.position = "bottom")
```

\newpage   


# Acknowledgements

The authors graciously thank COST Actions FP1303 and FP1407, the IO-0035 Slovenian Infrastructural program for their funding and support and the European Commission for funding the InnoRenew CoE project (Grant Agreement #739574) under the Horizon2020 Widespread-Teaming program. The authors also thank and acknowledge the technical support of Primož Habjan, Courtney Williamson, and Marko Posavčević.  

\newpage

# Environment
```{r}
sessionInfo()
```

