# Scripts for performing generic PRC analysis of soil chemistry through the profile through time using the vegan package (author Gregg Milligan).
"#Requires a “samp” dataframe containing:  1. A factor column named “trt”  containing treatment in the form of “NT” for no tillage (control); “FIT” for full inversion tillage (time zero); “FIT+1” … “FIT+n” for each time period following.  2. A numeric/integer column named “Depth” containing the depth of the soil strata, e.g. the mid-point of a 10cm auger sample."
# Requires a “chem” dataframe containing soil chemistry data with each column containg one variable.

rm(list=ls())
 #Clear workspace.

# 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(""scales"", ""lattice"", ""latticeExtra"", ""ggplot2"", ""vegan"", ""BiodiversityR"", ""plyr"", ""doMC"", ""doMC"", ""parallel"", ""foreach"", ""compiler"")"
pacFunc(packages)

## Load data.
"samp <- read.table(""<replace with sampling data file>"", header = TRUE) # Load dataframe containing sampling information."
"chem <- read.table(""<replace with soil chemistry data file>"", header = TRUE) # Load soil chemistry data frame."

str(samp) # ensure “trt” is a factor; “Depth” is integer.
str(chem) # ensure all are numeric.

"chemHellinger <- disttransform(chem, method = ""hellinger"") # Hellinger transform chem dataframe."

Depth <- factor(samp$depth) # Create depth factor for us ein PRC.
"Treatment ? factor(samp$trt, levels = c(""NT"", ""FIT"", ""FIT+1"", ""FIT+2"", ""FIT+3"", ""FIT+4"", ""FIT+5"", ""FIT+6"", ""FIT+7"", ""FIT+8"", ""FIT+9"", ""FIT+10"") #  Create treatment factor for us ein PRC. ***delete/add post-Fit period as required***"

"soilPRC <- prc(chemHellinger, Treatment, Depth) # Preform PRC."

anv0 <- anova(soilPRC) # ANOVA for model fit.
"anv1 <- anova(soilPRC, strata = Depth, first = TRUE) # ANOVA for first axis."

"prcPlot <- plot(soilPRC, xlab = ""Depth (cm)"", ylab = ""Effect"", scaling = 3, species = TRUE, lwd = 2, cex = 1.2, legpos = NA) # Plot of PRC diagram."

# Calculates explained variance per axis.
AxisVar <- soilPRC$CCA$eig/sum(soilPRC$CCA$eig)

# Calculates P values for mulitvariate differences at each strata.
registerDoMC(cores = 4)
Out <- NULL # Creates a NULL variable.
foreach(i = levels(Depth)) %do% { # Depth-wise selection
"    takeChem <- chem[Depth == i, ] # takeChem is an array of chemical-wise pRDA scores by depth."
    takeTreat <- Treatment[Depth == i] # takeTreat is a factor of depth-wise treatments.
"    Out[[i]] <- anova(rda(takeChem ~ takeTreat), 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."
}

"Pvals <- sapply(Out, function(x) x[1, 4]) # Grabs the P-values from out. Replace ""i"" in ""x[1,i]"" with depth strata required."