Reads MLSN data from GitHub and makes charts of Mehlich 3 P distribution

p_high_ph.R

# load packages, then read in data from Github

library(VGAM)
library(ggplot2)
library(cowplot)
library(dplyr)
library(RColorBrewer)

# generates an mlsn value for a vector using log logistic distribution
mlsn <- function(x) {
  fit.x <- vglm(x ~ 1, fisk)
  z <- Coef(fit.x)
  new.mlsn <- qfisk(0.1, z[1], z[2])
  return(round(new.mlsn, 0))
}

# read the full atc and pace data from GitHub
atc <- read.csv('https://raw.githubusercontent.com/micahwoods/2016_mlsn_paper/master/data/atc_samples.csv')

pace <- read.csv('https://raw.githubusercontent.com/micahwoods/2016_mlsn_paper/master/data/pace_samples.csv')

# only the columns I'm going to look at here
# makes it trivial to combine without complication 
atc2 <- select(atc, TECM3, pH, PM3)
pace2 <- select(pace, TECM3, pH, PM3)

# note the data set in case I want to check later
atc2$loc <- "atc"
pace2$loc <- "pace"

# combine the atc and pace data for CEC, pH, and Mehlich 3 P
both <- rbind.data.frame(atc2, pace2)

# do 4 curves
# mlsn normal
normal <- subset(both, pH >= 5.5 & pH <= 8.5 & TECM3 <= 6)
normal$type = "mlsn: pH >= 5.5 & pH <= 8.5 & TECM3 <= 6"

# mlsn data cut only by pH
uncut <- subset(both, pH >= 5.5 & pH <= 8.5)
uncut$type = "pH >= 5.5 & pH <= 8.5"

# pH 7.5 to 8 with the same CEC cut as mlsn
high_ph_cut <- subset(both, pH >= 7.5 & pH <= 8 & TECM3 <= 6)
high_ph_cut$type <- "pH >= 7.5 & pH <= 8 & TECM3 <= 6"

# ph 7.5 to 8 with no cut by CEC
high_ph_nocut <- subset(both, pH >= 7.5 & pH <= 8)
high_ph_nocut$type <- "pH >= 7.5 & pH <= 8"

# make a data frame to show the 0.1 value for each filtered data set
# so I can print that on a faceted ggplot
adf <- data.frame(xvalue = c(rep(250, 4)),
                   yvalue = c(rep(0.015, 4)),
                   type = c("mlsn: pH >= 5.5 & pH <= 8.5 & TECM3 <= 6",
                            "pH >= 5.5 & pH <= 8.5",
                            "pH >= 7.5 & pH <= 8 & TECM3 <= 6",
                            "pH >= 7.5 & pH <= 8" ),
                   value = c(mlsn(normal$PM3),
                             mlsn(uncut$PM3),
                             mlsn(high_ph_cut$PM3),
                             mlsn(high_ph_nocut$PM3)))

# combine the filtered data sets for plotting
d <- rbind.data.frame(normal, uncut, high_ph_cut, high_ph_nocut)

# make a chart that shows kernel density for each of the 4 ways I've cut the data set
p <- ggplot(data = d, aes(x = PM3))
yo1 <- p + theme_cowplot(font_family = "Roboto Condensed") +
  background_grid() +
  geom_vline(xintercept = 21, linetype = "dashed", alpha = 0.5) +
  geom_density(aes(fill = type, colour = type), alpha = 0.5) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  facet_wrap(~ type) +
  theme(legend.position = "none",
        plot.caption = element_text(size = 8)) +
  geom_label(data = adf, aes(x = xvalue, y = yvalue,
                            label = paste("0.1 level = ", value),
                            family = "Roboto Condensed",
             hjust = 0))

yo1

# same plot as above but cut the x axis short for detail
p <- ggplot(data = d, aes(x = PM3))
yo2 <- p + theme_cowplot(font_family = "Roboto Condensed") +
  background_grid() +
  geom_vline(xintercept = 21, linetype = "dashed", alpha = 0.5) +
  geom_density(aes(fill = type, colour = type), alpha = 0.5) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  facet_wrap(~ type) +
  scale_x_continuous(limits = c(0, 250)) +
  theme(legend.position = "none",
        plot.caption = element_text(size = 8)) +
  geom_label(data = adf, aes(x = xvalue, y = yvalue,
                            label = paste("0.1 level = ", value)),
             family = "Roboto Condensed",
            hjust = 1)

yo2

# save the plots
save_plot("~/Dropbox/charts/mlsnP_ph.png", yo1, base_asp = 1.78, scale = 1.25)
save_plot("~/Dropbox/charts/mlsnP_ph_250ppm_max.png", yo2, base_asp = 1.78, scale = 1.25)