# Scripts for performing PRC analysis of soil chemistry through the profile (author Gregg Milligan).
rm(list=ls())

# Function for loading required packages.
pacFunc <- function(pkg){
" new.pkg <- pkg[!(pkg %in% installed.packages()[, ""Package""])]"
 if (length(new.pkg))
" install.packages(new.pkg, dependencies = TRUE)"
" sapply(pkg, require, character.only = TRUE)"
}

"packages <- c(""aqp"", ""scales"", ""lattice"", ""latticeExtra"", ""ggplot2"", ""vegan"", ""BiodiversityR"", ""plyr"", ""doMC"", ""doMC"", ""parallel"", ""foreach"", ""compiler"")"
pacFunc(packages)

# read data for both sites in from a tab delimited “.txt” file
" soil <- read.table(""soil.txt"", header = TRUE)"
# check heads
 head(soil)
#split dataframe into one for each site “york” and “Kirkby”
 attach(soil)
" york<-soil[site==""york"",]"
" kirkby<-soil[site==""kirkby"",]"
 dim(soil)
 dim(york)
 dim(Kirkby)
 detach(soil)


# Read in data; soil is York + Kirkby combined.
"york <- read.delim(""data/soil/soil/york.csv"", header = TRUE, sep = "","")"
"kirkby <- read.delim(""data/soil/soil/kirkby.csv"", header = TRUE, sep = "","")"
"soil <- read.delim(""data/soil/soil/soil.csv"", header = TRUE, sep = "","")"

# Code for running PRCs on soil chemistry data.

# Create Hellinger-transformed chemistry matrix from each of the data frames.
"yorkHell <- disttransform(york[,8:17], method = ""hellinger"")"
"kirkbyHell <- disttransform(kirkby[,8:17], method = ""hellinger"")"

# Create scaled and centered chemistry matrix from ech of the data frames.
"yorkCen <- scale(york[,8:17], center = TRUE, scale = TRUE)"
"kirkbyCen <- scale(kirkby[,8:17], center = TRUE, scale = TRUE)"

# Create predictor suite for use in the PRCs.
yorkDepth <- factor(york$depth)
"yorkTreat <- factor(york$trt, levels = c(""NT"", ""FIT""))"
kirkbyDepth <- factor(kirkby$depth)
"kirkbyTreat <- factor(kirkby$trt, levels = c(""NT"", ""FIT""))"

# PRC
# York
"yorkPRChell <- prc(yorkHell, yorkTreat, yorkDepth)"
"yorkPRCcen <- prc(yorkCen, yorkTreat, yorkDepth)"

# Kirkby
"kirkbyPRChell <- prc(kirkbyHell, kirkbyTreat, kirkbyDepth)"
"kirkbyPRCcen <- prc(kirkbyCen, kirkbyTreat, kirkbyDepth)"

# ANOVAs for PRC models.
#York
anova(yorkPRChell)
"anova(yorkPRChell, strata = yorkDepth, first = TRUE)"
anova(yorkPRCcen)
"anova(yorkPRCcen, strata = yorkDepth, first = TRUE)"

# Kirkby
anova(kirkbyPRChell)
"anova(kirkbyPRChell, strata = kirkbyDepth, first = TRUE)"
anova(kirkbyPRCcen)
"anova(kirkbyPRCcen, strata = kirkbyDepth, first = TRUE)"

# Plots for PRC.

# Backup the old par.
old.par <- par(no.readonly=T)

# York
"yorkPlothell <- plot(yorkPRChell, xlab = ""Depth (cm)"", ylab = ""Effect"", scaling = 3, species = TRUE, lwd = 2, cex = 1.2, legpos = NA)"
"yorkPlotcen <- plot(yorkPRCcen, xlab = ""Depth (cm)"", ylab = ""Effect"", scaling = 3, species = TRUE, lwd = 2, cex = 1.2, legpos = NA)"

# Kirkby
"kirkbyPRChell <- plot(kirkbyPRChell, xlab = ""Depth (cm)"", ylab = ""Effect"", scaling = 3, species = TRUE, lwd = 2, cex = 1.2, legpos = NA)"
"kirkbyPRCcen <- plot(kirkbyPRCcen, xlab = ""Depth (cm)"", ylab = ""Effect"", scaling = 3, species = TRUE, lwd = 2, cex = 1.2, legpos = NA)"

## Calculates explained variance per axis.
# York
yorkAxisVarHell <- yorkPRChell$CCA$eig/sum(yorkPRChell$CCA$eig)
yorkAxisVarCen <- yorkPRCcen$CCA$eig/sum(yorkPRCcen$CCA$eig)

# Kirkby
kirkbyAxisVarHell <- kirkbyPRChell$CCA$eig/sum(kirkbyPRChell$CCA$eig)
kirkbyAxisVarCen <- kirkbyPRCcen$CCA$eig/sum(kirkbyPRCcen$CCA$eig)

# Calculates P values for mulitvariate differences at each strata.
# York
# Hellinger-matrix
registerDoMC(cores = 4)
yorkOutHell <- NULL # Creates a NULL variable.
foreach(i = levels(yorkDepth)) %do% { # Depth-wise selection
" takeChemHell <- yorkHell[yorkDepth == i, ] # takeChemHell is an array of chemical-wise pRDA scores by depth."
 takeTreatHell <- yorkTreat[yorkDepth == i] # takeTreatHell is a factor of depth-wise treatments.
" yorkOutHell[[i]] <- anova(rda(takeChemHell ~ takeTreatHell), by = ""terms"", step = 1000) # Performs a term-wise ANOVA on take.spec as a function of take.treatment at each depth.out is an array of depth-wise ANOVA results."
}

"yorkPvalsHell <- sapply(yorkOutHell, function(x) x[1, 4]) # Grabs the P-values from out."

# Centered-matrix
registerDoMC(cores = 4)
yorkOutCen <- NULL
foreach(i = levels(yorkDepth)) %do% {
" takeChemCen <- yorkCen[yorkDepth == i, ]"
 takeTreatCen <- yorkTreat[yorkDepth == i]
" yorkOutCen[[i]] <- anova(rda(takeChemCen ~ takeTreatCen), by = ""terms"", step = 1000)"
}

"yorkPvalsCen <- sapply(yorkOutCen, function(x) x[1, 4])"

# Kirkby
# Hellinger-matrix
registerDoMC(cores = 4)
kirkbyOutHell <- NULL # Creates a NULL variable.
foreach(i = levels(kirkbyDepth)) %do% { # Depth-wise selection
" takeChemHell <- kirkbyHell[kirkbyDepth == i, ] # takeChemHell is an array of chemical-wise pRDA scores by depth."
 takeTreatHell <- kirkbyTreat[kirkbyDepth == i] # takeTreatHell is a factor of depth-wise treatments.
" kirkbyOutHell[[i]] <- anova(rda(takeChemHell ~ takeTreatHell), by = ""terms"", step = 1000) # Performs a term-wise ANOVA on take.spec as a function of take.treatment at each depth.out is an array of depth-wise ANOVA results."
}

"kirkbyPvalsHell <- sapply(kirkbyOutHell, function(x) x[1, 4]) # Grabs the P-values from out."

# Centered-matrix
registerDoMC(cores = 4)
kirkbyOutCen <- NULL
foreach(i = levels(kirkbyDepth)) %do% {
" takeChemCen <- kirkbyCen[kirkbyDepth == i, ]"
 takeTreatCen <- kirkbyTreat[kirkbyDepth == i]
" kirkbyOutCen[[i]] <- anova(rda(takeChemCen ~ takeTreatCen), by = ""terms"", step = 1000)"
}

"kirkbyPvalsCen <- sapply(kirkbyOutCen, function(x) x[1, 4])"

# t-test plots. Plot the relative change in each of the variables above in repsonse to depth.
## Calculates P-values for Hellinger-transformed community matrix.
# Grabs the P-values from out.
## Dunnett contrasts comparing FIT to NT.

"df <- data.frame(Treatment = Treatment, Depth = Depth)"
# Package for multiple comparisons
require(multcomp)

# Create empty object

out_willi <- NULL

"# Loop through time, compute PCA, extract scores and apply Williams test."
registerDoMC(cores = 4)
foreach(i = levels(Depth)) %do% {
" pca <- rda(cen[Depth == i, ]) # Compute PCA"
" pca_scores <- scores(pca, display = ""sites"", choices = 1) # Scores of first principle component"
" out_willi[[i]] <- summary(glht(aov(pca_scores ~ Treatment, data = df[Depth == i, ]), alternative = ""t"", linfct = mcp(Treatment = ""Dunnett"")))"
}

# extract p-values
"result <- lapply(out_willi, function(x) data.frame(comp = levels(df$Treatment)[-1],"
"pval = x$test$pvalues, sig = x$test$pvalues < 0.05))"
# shows the results of Williams-Test on PCA-scores for depth 1:
"result[[""1""]]"

## Relative change in C with depth.

"ratc.mean <- array(0,length(crat$ratc)/9) ## Create depth-length vector."
"ratc.mean[1] <- mean((crat$ratc)[crat$samp == ""1""]) ## Depth-wise mean change in soil C"
"ratc.mean[2] <- mean((crat$ratc)[crat$samp == ""2""])"
"ratc.mean[3] <- mean((crat$ratc)[crat$samp == ""3""])"
"ratc.mean[4] <- mean((crat$ratc)[crat$samp == ""4""])"
"ratc.mean[5] <- mean((crat$ratc)[crat$samp == ""5""])"
"ratc.mean[6] <- mean((crat$ratc)[crat$samp == ""6""])"
"ratc.mean[7] <- mean((crat$ratc)[crat$samp == ""7""])"
"ratc.mean[8] <- mean((crat$ratc)[crat$samp == ""8""])"

"d.1 <- crat$ratc[crat$samp == ""1""]"
"d.2 <- crat$ratc[crat$samp == ""2""]"
"d.3 <- crat$ratc[crat$samp == ""3""]"
"d.4 <- crat$ratc[crat$samp == ""4""]"
"d.5 <- crat$ratc[crat$samp == ""5""]"
"d.6 <- crat$ratc[crat$samp == ""6""]"
"d.7 <- crat$ratc[crat$samp == ""7""]"
"d.8 <- crat$ratc[crat$samp == ""8""]"

"ratc.asin <- array(0,length(crat$ratc)/9)"
ratc.asin[1] <- mean(sqrt(abs(d.1)) * asin(sign(d.1)))
ratc.asin[2] <- mean(sqrt(abs(d.2)) * asin(sign(d.2)))
ratc.asin[3] <- mean(sqrt(abs(d.3)) * asin(sign(d.3)))
ratc.asin[4] <- mean(sqrt(abs(d.4)) * asin(sign(d.4)))
ratc.asin[5] <- mean(sqrt(abs(d.5)) * asin(sign(d.5)))
ratc.asin[6] <- mean(sqrt(abs(d.6)) * asin(sign(d.6)))
ratc.asin[7] <- mean(sqrt(abs(d.7)) * asin(sign(d.7)))
ratc.asin[8] <- mean(sqrt(abs(d.8)) * asin(sign(d.8)))

t.test(sqrt(abs(d.1)) * asin(sign(d.1)))


"plot(crat$depth ~ crat$ratc, pch = 20, ylim = c(100,0), xlim = c(-0.5, 1.5), xlab = ""(FIT-NT)/NT"", ylab = ""Depth (cm)"", las = 1)"
"abline(v = c(0,1))"
par(new = TRUE)
"plot(unique(crat$depth) ~ ratc.mean, axes = F, xlab = """", ylab = """", ylim = c(100,0), xlim = c(-0.5, 1.5), col = ""red"", pch = 15, cex = 1.1)"


# Create dataframe containg relative changes in variables.

"relYork <- read.delim(""./data/soilRel/relYork.csv"", header = TRUE, sep = "","")"
"relKirkby <- read.delim(""./data/soilRel/relKirkby.csv"", header = TRUE, sep = "","")"

# York
"YmeanRelTOC <- tapply(relYork$toc, relYork$samp, mean)"
"YmeanRelN <- tapply(relYork$n, relYork$samp, mean)"
"YmeanRelCN <- tapply(relYork$cn, relYork$samp, mean)"
"YmeanRelpH <- tapply(relYork$ph, relYork$samp, mean)"
"YmeanRelNO3 <- tapply(relYork$no3, relYork$samp, mean)"
"YmeanRelNH4 <- tapply(relYork$nh4, relYork$samp, mean)"
"YmeanRelPO4 <- tapply(relYork$po4, relYork$samp, mean)"
"YmeanRelK <- tapply(relYork$k, relYork$samp, mean)"
"YmeanRelMg <- tapply(relYork$mg, relYork$samp, mean)"
"YmeanRelCa <- tapply(relYork$ca, relYork$samp, mean)"

"Yd1TOC <- relYork$toc[relYork$samp == ""1""]"
"Yd2TOC <- relYork$toc[relYork$samp == ""2""]"
"Yd3TOC <- relYork$toc[relYork$samp == ""3""]"
"Yd4TOC <- relYork$toc[relYork$samp == ""4""]"
"Yd5TOC <- relYork$toc[relYork$samp == ""5""]"
"Yd6TOC <- relYork$toc[relYork$samp == ""6""]"
"Yd7TOC <- relYork$toc[relYork$samp == ""7""]"
"Yd8TOC <- relYork$toc[relYork$samp == ""8""]"

"YtocAsin <- array(0, length(relYork$samp)/9)"
YtocAsin[1] <- mean(sqrt(abs(Yd1TOC)) * asin(sign(Yd1TOC)))
YtocAsin[2] <- mean(sqrt(abs(Yd2TOC)) * asin(sign(Yd2TOC)))
YtocAsin[3] <- mean(sqrt(abs(Yd3TOC)) * asin(sign(Yd3TOC)))
YtocAsin[4] <- mean(sqrt(abs(Yd4TOC)) * asin(sign(Yd4TOC)))
YtocAsin[5] <- mean(sqrt(abs(Yd5TOC)) * asin(sign(Yd5TOC)))
YtocAsin[6] <- mean(sqrt(abs(Yd6TOC)) * asin(sign(Yd6TOC)))
YtocAsin[7] <- mean(sqrt(abs(Yd7TOC)) * asin(sign(Yd7TOC)))
YtocAsin[8] <- mean(sqrt(abs(Yd8TOC)) * asin(sign(Yd8TOC)))

"plot(relYork$depth ~ relYork$toc, pch = 20, ylim = c(100,0), xlim = c(-0.5, 1.5), xlab = ""(FIT-NT)/NT"", ylab = ""Depth (cm)"", las = 1)"
"abline(v = c(0,1))"
par(new = TRUE)
"plot(unique(relYork$depth) ~ YmeanRelTOC, axes = F, xlab = """", ylab = """", ylim = c(100,0), xlim = c(-0.5, 1.5), col = ""red"", pch = 15, cex = 1.1)"

# N
"Yd1n <- relYork$n[relYork$samp == ""1""]"
"Yd2n <- relYork$n[relYork$samp == ""2""]"
"Yd3n <- relYork$n[relYork$samp == ""3""]"
"Yd4n <- relYork$n[relYork$samp == ""4""]"
"Yd5n <- relYork$n[relYork$samp == ""5""]"
"Yd6n <- relYork$n[relYork$samp == ""6""]"
"Yd7n <- relYork$n[relYork$samp == ""7""]"
"Yd8n <- relYork$n[relYork$samp == ""8""]"

"YnAsin <- array(0, length(relYork$samp)/9)"
YnAsin[1] <- mean(sqrt(abs(Yd1n)) * asin(sign(Yd1n)))
YnAsin[2] <- mean(sqrt(abs(Yd2n)) * asin(sign(Yd2n)))
YnAsin[3] <- mean(sqrt(abs(Yd3n)) * asin(sign(Yd3n)))
YnAsin[4] <- mean(sqrt(abs(Yd4n)) * asin(sign(Yd4n)))
YnAsin[5] <- mean(sqrt(abs(Yd5n)) * asin(sign(Yd5n)))
YnAsin[6] <- mean(sqrt(abs(Yd6n)) * asin(sign(Yd6n)))
YnAsin[7] <- mean(sqrt(abs(Yd7n)) * asin(sign(Yd7n)))
YnAsin[8] <- mean(sqrt(abs(Yd8n)) * asin(sign(Yd8n)))

"plot(relYork$depth ~ relYork$n, pch = 20, ylim = c(100,0), xlim = c(-1.5, 1.5), xlab = ""(FIT-NT)/NT"", ylab = ""Depth (cm)"", las = 1)"
"abline(v = c(0,1))"
par(new = TRUE)
"plot(unique(relYork$depth) ~ YmeanRelN, axes = F, xlab = """", ylab = """", ylim = c(100,0), xlim = c(-1.5, 1.5), col = ""red"", pch = 15, cex = 1.1)"

# CN
"Yd1cn <- relYork$n[relYork$samp == ""1""]"
"Yd2cn <- relYork$n[relYork$samp == ""2""]"
"Yd3cn <- relYork$n[relYork$samp == ""3""]"
"Yd4cn <- relYork$n[relYork$samp == ""4""]"
"Yd5cn <- relYork$n[relYork$samp == ""5""]"
"Yd6cn <- relYork$n[relYork$samp == ""6""]"
"Yd7cn <- relYork$n[relYork$samp == ""7""]"
"Yd8cn <- relYork$n[relYork$samp == ""8""]"

"YcnAsin <- array(0, length(relYork$samp)/9)"
YcnAsin[1] <- mean(sqrt(abs(Yd1cn)) * asin(sign(Yd1cn)))
YcnAsin[2] <- mean(sqrt(abs(Yd2cn)) * asin(sign(Yd2cn)))
YcnAsin[3] <- mean(sqrt(abs(Yd3cn)) * asin(sign(Yd3cn)))
YcnAsin[4] <- mean(sqrt(abs(Yd4cn)) * asin(sign(Yd4cn)))
YcnAsin[5] <- mean(sqrt(abs(Yd5cn)) * asin(sign(Yd5cn)))
YcnAsin[6] <- mean(sqrt(abs(Yd6cn)) * asin(sign(Yd6cn)))
YcnAsin[7] <- mean(sqrt(abs(Yd7cn)) * asin(sign(Yd7cn)))
YcnAsin[8] <- mean(sqrt(abs(Yd8cn)) * asin(sign(Yd8cn)))

"plot(relYork$depth ~ relYork$cn, pch = 20, ylim = c(100,0), xlim = c(-0.5, 1.5), xlab = ""(FIT-NT)/NT"", ylab = ""Depth (cm)"", las = 1)"
"abline(v = c(0,1))"
par(new = TRUE)
"plot(unique(relYork$depth) ~ YmeanRelCN, axes = F, xlab = """", ylab = """", ylim = c(100,0), xlim = c(-0.5, 1.5), col = ""red"", pch = 15, cex = 1.1)"

#pH
"Yd1pH <- relYork$n[relYork$samp == ""1""]"
"Yd2pH <- relYork$n[relYork$samp == ""2""]"
"Yd3pH <- relYork$n[relYork$samp == ""3""]"
"Yd4pH <- relYork$n[relYork$samp == ""4""]"
"Yd5pH <- relYork$n[relYork$samp == ""5""]"
"Yd6pH <- relYork$n[relYork$samp == ""6""]"
"Yd7pH <- relYork$n[relYork$samp == ""7""]"
"Yd8pH <- relYork$n[relYork$samp == ""8""]"

"YpHAsin <- array(0, length(relYork$samp)/9)"
YpHAsin[1] <- mean(sqrt(abs(Yd1pH)) * asin(sign(Yd1pH)))
YpHAsin[2] <- mean(sqrt(abs(Yd2pH)) * asin(sign(Yd2pH)))
YpHAsin[3] <- mean(sqrt(abs(Yd3pH)) * asin(sign(Yd3pH)))
YpHAsin[4] <- mean(sqrt(abs(Yd4pH)) * asin(sign(Yd4pH)))
YpHAsin[5] <- mean(sqrt(abs(Yd5pH)) * asin(sign(Yd5pH)))
YpHAsin[6] <- mean(sqrt(abs(Yd6pH)) * asin(sign(Yd6pH)))
YpHAsin[7] <- mean(sqrt(abs(Yd7pH)) * asin(sign(Yd7pH)))
YpHAsin[8] <- mean(sqrt(abs(Yd8pH)) * asin(sign(Yd8pH)))

"plot(relYork$depth ~ relYork$ph, pch = 20, ylim = c(100,0), xlim = c(-0.5, 1.5), xlab = ""(FIT-NT)/NT"", ylab = ""Depth (cm)"", las = 1)"
"abline(v = c(0,1))"
par(new = TRUE)
"plot(unique(relYork$depth) ~ YmeanRelpH, axes = F, xlab = """", ylab = """", ylim = c(100,0), xlim = c(-0.5, 1.5), col = ""red"", pch = 15, cex = 1.1)"

#PO4
"Yd1po4 <- relYork$n[relYork$samp == ""1""]"
"Yd2po4 <- relYork$n[relYork$samp == ""2""]"
"Yd3po4 <- relYork$n[relYork$samp == ""3""]"
"Yd4po4 <- relYork$n[relYork$samp == ""4""]"
"Yd5po4 <- relYork$n[relYork$samp == ""5""]"
"Yd6po4 <- relYork$n[relYork$samp == ""6""]"
"Yd7po4 <- relYork$n[relYork$samp == ""7""]"
"Yd8po4 <- relYork$n[relYork$samp == ""8""]"

"Ypo4Asin <- array(0, length(relYork$samp)/9)"
Ypo4Asin[1] <- mean(sqrt(abs(Yd1po4)) * asin(sign(Yd1po4)))
Ypo4Asin[2] <- mean(sqrt(abs(Yd2po4)) * asin(sign(Yd2po4)))
Ypo4Asin[3] <- mean(sqrt(abs(Yd3po4)) * asin(sign(Yd3po4)))
Ypo4Asin[4] <- mean(sqrt(abs(Yd4po4)) * asin(sign(Yd4po4)))
Ypo4Asin[5] <- mean(sqrt(abs(Yd5po4)) * asin(sign(Yd5po4)))
Ypo4Asin[6] <- mean(sqrt(abs(Yd6po4)) * asin(sign(Yd6po4)))
Ypo4Asin[7] <- mean(sqrt(abs(Yd7po4)) * asin(sign(Yd7po4)))
Ypo4Asin[8] <- mean(sqrt(abs(Yd8po4)) * asin(sign(Yd8po4)))

"plot(relYork$depth ~ relYork$po4, pch = 20, ylim = c(100,0), xlim = c(-5.0, 1.0), xlab = ""(FIT-NT)/NT"", ylab = ""Depth (cm)"", las = 1)"
"abline(v = c(0,1))"
par(new = TRUE)
"plot(unique(relYork$depth) ~ YmeanRelPO4, axes = F, xlab = """", ylab = """", ylim = c(100,0), xlim = c(-5.0, 1.0), col = ""red"", pch = 15, cex = 1.1)"

#NH4
"Yd1nh4 <- relYork$n[relYork$samp == ""1""]"
"Yd2nh4 <- relYork$n[relYork$samp == ""2""]"
"Yd3nh4 <- relYork$n[relYork$samp == ""3""]"
"Yd4nh4 <- relYork$n[relYork$samp == ""4""]"
"Yd5nh4 <- relYork$n[relYork$samp == ""5""]"
"Yd6nh4 <- relYork$n[relYork$samp == ""6""]"
"Yd7nh4 <- relYork$n[relYork$samp == ""7""]"
"Yd8nh4 <- relYork$n[relYork$samp == ""8""]"

"Ynh4Asin <- array(0, length(relYork$samp)/9)"
Ynh4Asin[1] <- mean(sqrt(abs(Yd1nh4)) * asin(sign(Yd1nh4)))
Ynh4Asin[2] <- mean(sqrt(abs(Yd2nh4)) * asin(sign(Yd2nh4)))
Ynh4Asin[3] <- mean(sqrt(abs(Yd3nh4)) * asin(sign(Yd3nh4)))
Ynh4Asin[4] <- mean(sqrt(abs(Yd4nh4)) * asin(sign(Yd4nh4)))
Ynh4Asin[5] <- mean(sqrt(abs(Yd5nh4)) * asin(sign(Yd5nh4)))
Ynh4Asin[6] <- mean(sqrt(abs(Yd6nh4)) * asin(sign(Yd6nh4)))
Ynh4Asin[7] <- mean(sqrt(abs(Yd7nh4)) * asin(sign(Yd7nh4)))
Ynh4Asin[8] <- mean(sqrt(abs(Yd8nh4)) * asin(sign(Yd8nh4)))

"plot(relYork$depth ~ relYork$nh4, pch = 20, ylim = c(100,0), xlim = c(-30.0, 30), xlab = ""(FIT-NT)/NT"", ylab = ""Depth (cm)"", las = 1)"
"abline(v = c(0,1))"
par(new = TRUE)
"plot(unique(relYork$depth) ~ YmeanRelNH4, axes = F, xlab = """", ylab = """", ylim = c(100,0), xlim = c(-30,30), col = ""red"", pch = 15, cex = 1.1)"

# Kirkby
"KmeanRelTOC <- tapply(relKirkby$toc, relKirkby$samp, mean)"
"KmeanRelN <- tapply(relKirkby$n, relKirkby$samp, mean)"
"KmeanRelCN <- tapply(relKirkby$cn, relKirkby$samp, mean)"
"KmeanRelpH <- tapply(relKirkby$ph, relKirkby$samp, mean)"
"KmeanRelNO3 <- tapply(relKirkby$no3, relKirkby$samp, mean)"
"KmeanRelNH4 <- tapply(relKirkby$nh4, relKirkby$samp, mean)"
"KmeanRelPO4 <- tapply(relKirkby$po4, relKirkby$samp, mean)"
"KmeanRelK <- tapply(relKirkby$k, relKirkby$samp, mean)"
"KmeanRelMg <- tapply(relKirkby$mg, relKirkby$samp, mean)"
"KmeanRelCa <- tapply(relKirkby$ca, relKirkby$samp, mean)"