Combined.Rmd

---
title: "PairedData"
author: "Kendalynn A. Morris"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: yes
    toc_depth: '2'
    toc_float: yes
---

```{r setup, include=FALSE}
#load packages
library(dplyr)
library(lubridate)
library(ggplot2)
library(kableExtra)
library(ggbreak)
select <- dplyr::select

theme_set(theme_minimal() + theme(text = element_text(size = 11)))

# 'Pretty n' function to round a numeric value and print that # of digits 
pn <- function(x, n) {
  formatC(round(unlist(x), n),
          digits = n, format = "f")
}

# 'Clean p value' function to pretty-print p value(s), specifically
pclean <- function(x, digits = 3, printP = TRUE) {
  x <- as.vector(x)
  ltstring <- paste0("< 0.", paste(rep("0", digits - 1), 
                                   collapse = ""), "1")
  valstring <- ifelse(x < 10 ^ -digits, 
                      ltstring, pn(x, digits))
  if(printP) {
    paste("P", ifelse(x < 10 ^ -digits, 
                      valstring, paste("=", valstring)))
  } else {
    valstring
  }
}
```

Here the field and lab data are combined for comparison and graphing.

# Read in and combine data

```{r Field, echo = FALSE, warning = FALSE, message = TRUE}
f_dat <- read.csv("field-measurements/licorRTA.csv")
f_dat$date <- as_date(f_dat$date, "%m/%d/%Y", tz="EST")
f_dat$timestamp <- mdy_hm(f_dat$timestamp, tz="EST")
f_dat$Collar <- as.factor(f_dat$Collar)

f_dat$SWC <- replace(f_dat$SWC, f_dat$SWC < 0.10, NA)
f_dat %>% #select locations where gross rates were measured
  filter(Location %in% c("g_up", "g_low"),
         #select week when soils were collected
         week(timestamp) == 30) -> f_dat

f_dat$Origin <- recode_factor(f_dat$Origin,
                             "g_low" = "lowland",
                             "g_mid" = "midslope",
                             "g_up" = "upland",
                             "midstream" = "midstream",
                             "upsream" = "upstream")

#labelled column with current collar location in the same way
f_dat$Location <- recode_factor(f_dat$Location,
                               "g_low" = "lowland",
                               "g_mid" = "midslope",
                               "g_up" = "upland")

summary(f_dat)
```

```{r Lab, echo = FALSE, warning = FALSE, message = TRUE}
l_dat <- read.csv("labDat.csv")
l_dat %>%
  select("id","Origin","Location","round",
         "time_days","mt","mass","sm",
         "umolP","umolC",
         "P_rate","C_rate") %>%
  #ml in 130 ml jar divided by 1000 ml (1 L) for ppm
  mutate(ppm = mt / 0.13,
         minute = time_days * 24 * 60,
         Collar = as.factor(id)) -> l_dat

l_dat$Origin <- recode_factor(l_dat$Origin,
                             "g_low" = "lowland",
                             "g_mid" = "midslope",
                             "g_up" = "upland",
                             "midstream" = "midstream",
                             "upsream" = "upstream")

#labelled column with current collar location in the same way
l_dat$Location <- recode_factor(l_dat$Location,
                               "g_low" = "lowland",
                               "g_mid" = "midslope",
                               "g_up" = "upland")

summary(l_dat)
```

# Combine Data

```{r Merge baby Merge!}

l_dat %>%
  filter(round == "T4") %>%
  #assuming that flux from soil plug represents flux from 3.14 square centimeters
  #(surface area of the core collected)
  # current value * days/second * 1/3.14 * cm2/m2 * nmol/umol
  # (rate in umol per g dry soil per d) * dry mass *
  # 1/86400 * 1/3.14 * 1/0.0001 * 1000 = 36.9
  mutate(lab_net = (P_rate - C_rate) * mass * 36.9) %>%
  left_join(f_dat, by = c("Collar","Location","Origin")) -> pairedDat 
  

Olabs <- c("lowland", "midslope", "upland", "midstream", "upstream")
Llabs <- c("lowland", "upland")
```

```{r graphs, warning = FALSE, message = FALSE}

ggplot(pairedDat, aes(lab_net, FCH4, color = Origin)) +
    geom_point() +
    scale_fill_discrete(name = "Soil Origin",
                      labels = Olabs) +
    labs(x = "Lab net flux", y = "Field net flux",
         title = "Lab vs Field") +
  facet_wrap(Location~., scales = "free") +
    theme_bw() + theme(legend.position = "bottom")

pairedDat %>%
  group_by(Location, Origin) %>%
  mutate(n = n(),
         Field = median(FCH4),
         Field_se = sd(FCH4)/sqrt(n),
         Lab = median(lab_net),
         Lab_se = sd(lab_net)/sqrt(n),
         ) %>%
  mutate_if(is.numeric, round, digits=2) %>%
  select(Location, Origin,
         Field, Field_se,
         Lab, Lab_se, n) %>%
  unique(.) %>%
  arrange(Location, Origin) -> MethaneComparison

kable(MethaneComparison)

```

# Lab Net Rates

```{r Net Lab Rate Figure}

pairedDat$group <- "near zero"
pairedDat[pairedDat$lab_net < -0.5,]$group <- "consuming"
pairedDat[pairedDat$lab_net > 10,]$group <- "strongly producing"

pairedDat %>%
  select(group, lab_net, Collar) %>%
  right_join(l_dat, by = "Collar") %>%
  filter(round != "T5") %>%
  mutate(type = paste(Origin, "in", Location)) -> l_dat_fluxgroups

ggplot(l_dat_fluxgroups, aes(minute, ppm)) +
  geom_line(aes(group=Collar)) +
  geom_jitter(aes(shape = Location, color = lab_net), size = 3) #+
  #scale_y_break(c(20, 60))

ggplot(l_dat_fluxgroups, aes(lab_net, fill = type)) +
  xlab("Net Rate \n nmol per meter squared per second") +
  ggtitle("Lab Net Rates") +
  geom_histogram(color = "black") +
  scale_x_break(c(100, 900))

```

# Field Net Rates

```{r Net Field Rate Figure}

pairedDat %>%
  select(FCH4, Collar, Origin, Location) %>%
  mutate(type = paste(Origin, "in", Location)) -> f_dat_fluxgroups

ggplot(f_dat_fluxgroups, aes(FCH4, fill = type)) +
  xlab("Net Rate \n nmol per meter squared per second") +
  ggtitle("Field Net Rates") +
  geom_histogram(color = "black")

```