plotting.Rmd

## Load libraries
```{r}
library(tidyverse)
library(ggpubr)
library(gridExtra)
library(superb)

```

## Read in data
```{r}

#tidyData <- read.csv("chemistry-data-tidy.csv")

# Read in the same data but with nutrients for the dry treatment reported as mg/kg wet weight
tidyDataDry <- read.csv("chemistry-data-tidy-dryCorrected.csv")

```

```{r}
str(tidyDataDry)

```


## Convert weeks to factor so it can be used to group the plots
```{r}
tidyDataDry$Weeks <- as.factor(tidyDataDry$Weeks)
str(tidyDataDry)
```

```{r}

tidyDataDry<- tidyDataDry %>% 
  rename(
   "NH4 (mg/kg)" = NH4..mg.kg.,
   "K (mg/kg)" = K..mg.kg.,
   "NO2 (mg/kg)" = NO2..mg.kg.,
   "NO3 (mg/kg)" = NO3..mg.kg.,
   "Moisture (%)" = Moisture....
  )

head(tidyDataDry)

```
## Replacing negative values in the moisture and nitrate columns with 0. 
```{r}

tidyDataDry <- tidyDataDry %>% 
  mutate(across(c(`Moisture (%)`, `NO3 (mg/kg)`:`NH4 (mg/kg)`), ~ ifelse(.x<0, 0, .x)))


tail(tidyDataDry, 20)
```


## Make frequency historgrams
```{r}

for (i in 3:7) {
  plot1 <- ggplot(tidyDataDry, aes(tidyDataDry[,i])) +
    geom_histogram() +
    labs(x = colnames(tidyDataDry[i]))
  print(plot1)
  
}
```

Make boxplots
```{r}
plot.list <- list()
for(i in 4:7){
  plot <- ggplot(tidyDataDry, aes(x = Weeks, y = tidyDataDry[,i])) +
          stat_boxplot(aes(Weeks, tidyDataDry[,i]), geom = "errorbar") +
          geom_boxplot(aes(Weeks, tidyDataDry[,i]), outlier.shape = NA, coef = 0) +
          ylab(colnames(tidyDataDry[i])) +
          facet_wrap(.~Treatment) +
          theme(text = element_text(family = "Arial", size = 16), panel.background = element_rect(fill = NA, colour = 'black')) +
          scale_x_discrete(labels = c("Fresh control", "7", "24")) +
          scale_y_continuous(limits = c(min(tidyDataDry[,i]), max(tidyDataDry[,i]))) 
  plot.list[[i]] <- plot
 #ggsave(file = paste0("boxplot-by-treatment-", colnames(tidyData[i]), ".jpeg"), plot = plot, height = 7, width = 9, units = "in")
}

no3.plot <- plot.list[[4]]
no2.plot <- plot.list[[5]]
nh4.plot <- plot.list[[6]]
k.plot <- plot.list[[7]]


```

Plotting pH separately so that y axis can be adjusted in order to remove white space on plot
```{r}

pH_plot <- ggplot(tidyDataDry, aes(x = Weeks, y = pH)) +
          stat_boxplot(aes(Weeks, pH), geom = "errorbar") +
          geom_boxplot(aes(Weeks, pH), outlier.shape = NA, coef = 0) +
          ylab("pH") +
          facet_wrap(.~Treatment) +
          theme(text = element_text(family = "Arial", size = 16), 
                panel.background = element_rect(fill = NA, colour ='black')) +
          scale_x_discrete(labels = c("Fresh control", "7", "24")) +
          scale_y_continuous(limits = c(5.3, 6.7))

ph.plot <- pH_plot + 
  showSignificance(c(2,3), 6.6, -0.07, "***", panel = list(Treatment = "Dry")) +
  showSignificance(c(1,2.9), 5.6, +0.09, "**", panel = list(Treatment = "Dry")) +
  showSignificance(c(1,2), 6.67, -0.07, "***", panel = list(Treatment = "Fridge")) +
  showSignificance(c(2.2,3), 6.64, -0.07, "***", panel = list(Treatment = "Fridge")) +
  showSignificance(c(1,3), 5.6, +0.07, "***", panel = list(Treatment = "Fridge")) +
  showSignificance(c(2,3), 6.6, -0.07, "***", panel = list(Treatment = "Frozen"))
  
ph.plot

ggsave(file = "pH-boxplot.jpeg", pH.plot, height = 7, width = 9, units = "in")



```

```{r}




```



## Plotting Nitrogen line graph (fig. 7)
Select only nitrogen data and assign to new dataframe
```{r}
N_data <- dryData %>%
  select(Sample.ID, `NO3 (mg/kg)`, `NO2 (mg/kg)`, `NH4 (mg/kg)`, Treatment, Weeks)

head(N_data)
```

```{r}

N_data[!complete.cases(N_data), ]

```

```{r}

N_data <- N_data %>%
  drop_na()

N_data[!complete.cases(N_data), ]
```

## Make data into long format for plotting on line graph
```{r}

N_data_long <- N_data %>%
  pivot_longer(cols = `NO3 (mg/kg)`:`NH4 (mg/kg)`,
               names_to = "Nitrogen compound",
               values_to = "Concentration (mg/kg)")

head(N_data_long, 20)
```
## Calculate summary statistics
```{r}

summaryN <- N_data_long %>%
  group_by(Treatment, Weeks, `Nitrogen compound`) %>%
  summarise(
    meanConc = mean(`Concentration (mg/kg)`),
    sdConc = sd(`Concentration (mg/kg)`)
  )

summaryN
```

```{r}

summaryN[!complete.cases(summaryN), ]

```

## Make the plot (fig. 7)
```{r}


N_plot <- ggplot(summaryN, aes(x = factor(Weeks), y = meanConc, color = factor(`Nitrogen compound`), group = `Nitrogen compound`)) +
  geom_line() +
  geom_point() +
  facet_wrap(~Treatment) +
  theme(text = element_text(family = "Arial", size = 16), panel.grid.major = element_blank(), 
           panel.grid.minor = element_blank(), legend.title = element_text(size = 12), 
           legend.position = "bottom") +
  labs(x = "Time (weeks)", y = "Mean concentration (mg/kg)", shape = "Nitrogen form") +
  scale_color_discrete(labels = c("NH4", "NO2", "NO3"))

  

N_plot + scale_color_manual(values = c("#E69F00", "#009E73",  "#56B4E9"))
```

Save the plot as jpeg
```{r}

ggsave(filename = "Nitrogen-linegraph.jpeg", plot = N_plot,  width = 9, unit = "in")

```