enter image description here I have different categories of environmental variables such as the biotic factor (temperature, windspeed), spatial factor (longitude and latitude), and soil properties (soil ph) and I would like to perform Canonical Correspondence analysis. I successfully performed the analysis but I don't know how to customize the plot. I want to customize the different categories of environmental variables by allotting them with different colours to differentiate spatial from biotic. In addition to this, how can I group the species based on a variable such as elevation etc? Thanks Below is my code
spe <- read.csv("C:/Microbiome/Result files/spe.csv", row.names=1)
env <- read.csv("C:/Microbiome/Result files/env.csv", row.names=1)
spatial <- read.csv("C://Microbiome/Result files/spatial.csv", row.names=1)
# Apply log+1 transformation to your species occurrences data (spe matrix) in order to correct for possible statistical errors associated to rare or very common species:
spelog <- decostand(spe, "log")
env <- na.omit(env) #This removes NA values in the dataframe, this is important to avoid errors
## Perform CCA. We need to specify that spe (species distribution matrix) is explained by env (environmental matrix).
ccamodel <- cca(spe~., env)
## To perform a Partial CCA, we need to use a third matrix (conditioning matrix). To do #that, we must first combine variables from the environmental ("env")
#and conditional ("spatial") matrices. We have to do that in order to later apply another function #("ordistep" function).
#After we have combined all variables, we will apply the Partial CCA using a formula where we specify the #response ("spe" matrix),
#the constraint variables (each of the variables from "env" matrix), and the conditioning variables ##c(variables from the conditioning matrix; in our case "spatial" matrix)
envspatial<-cbind(env,spatial)
nams <- names(envspatial)
partialccamodel <- formula(paste("spe ~", paste(nams[1: (length(envspatial)-(length(spatial)) )], collapse = " + "),
"+ Condition(", paste(nams[(length(envspatial)-(length(spatial)-1) ):length(envspatial)], collapse ="+"),")"))
partialccamodel<-cca(partialccamodel, envspatial)
# However, we also have to automatically select variables of "env" matrix that best explain "spe" matrix. We can do that by using a stepwise model from "ordistep" function. Let us do that with our "ccamodel" (not the partial cca).
finalmodel<- ordistep(ccamodel, scope=formula(ccamodel))
# Then, we can calculate Variance Inflation Factors (VIF) for each of the constraints (variables) from the "env" matrix (environmental matrix). If we find an environmental variable with VIF>10, we'll know that this variable presents colinearity with another or other variables. In that case, we would have to delete the variable from our initial dataset and redo all the analysis. In our example, no variable is redundant with each other (all of them have VIF<10).
vif.cca(finalmodel)
# Let us now fit the stepwise model (ordistep function) using our partial cca model ("partialccamodel"). Note that we will use
#X (longitude) and Y (latitude) variables from the previously created "envspatial" object as our #conditioning variables.
#After "vif.cca", note that "X" (spatial) variable has a value > 10, so one should consider to delete the #variable before the analysis.
partialccamodel <- formula(paste("spe ~", paste(nams[1: (length(envspatial)-(length(spatial)) )], collapse = " + "),"+ Condition(", paste(nams[(length(envspatial)-(length(spatial)-1) ):length(envspatial)], collapse ="+"),")"))
simplemodel<-cca(partialccamodel, envspatial)
finalmodelpartial<- ordistep(simplemodel, scope=formula(partialccamodel))
vif.cca(finalmodelpartial)
### Ok, let us now call "ccamodel" object (not the "partialccamodel") to see how to interpret results.
finalmodelpartial
#for color for plot
colvec <-c("magenta", "seagreen", "deepskyblue")
# Testing the significance of the CCA model:
anova.cca(finalmodel, perm.max=999)
# Testing the significance of terms (environmental variables):
anova.cca(finalmodel, by="terms", perm.max=999)
# Testing the significance of CCA axes (at least the first two or three should present a significant p value):
anova.cca(finalmodel, by="axis", perm.max=999)
### Finally, we may want to generate graphs in order to better understand results. To do that, we have to take a look
#at "xlim" (limits of x axis), "ylim" (limits of y axis), and "display" (if we want to observe species, #environmental gradients,
#and/or sites in the graph) arguments. For example, to show only species scores along CCA axes, use only #"sp" inside display" argument.
plot(finalmodel, xlim=c(-1.5,2), ylim=c(-1,1.5), display=c("species"), type = "text")
# If you want to show species and environmental gradients in our plot, use both "sp" and "cn" inside "display" argument,
plot(finalmodel, xlim=c(-3,3), ylim=c(-3,3), display=c("species","cn"), type = "text")
# To show species, environmental gradients, and sites, use "sp", "cn", and "wa" inside "display" argument,
plot(finalmodel, xlim=c(-4,4), ylim=c(-4,4), display=c("species","cn","wa"), type = "text", col= c("magenta", "seagreen", "deepskyblue"))
This code got me the diagram A but i want B
I have tried using different packages but all to no avail. CanonicalCorrespondenceanalysis