Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Reads MLSN data from GitHub and makes charts of Mehlich 3 P distribution
# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment