############################################################
############################################################
### ABLE 2025 - Generating a Simulated Dataset for Community
### Ecology Analyses 
############################################################

### SCRIPT FOR WORKSHOP EXERCISES

###########################################################
### SETTING UP
############################################################

## First, log into posit.cloud and initialize the RStudio
## workspace.

## Follow along with the provided guide to load both the
## RScript (Script.R) and demo dataset (Demo.txt) into
## your workspace.

## View your dataset to ensure that it loaded correctly:

View(Demo)

str(Demo)

## The first line of the str() output should read:

# "'data.frame':	27 obs. of  38 variables:"

## If it does not, please raise your hand and the instructor
## will assist you shortly.

## Once the data  is loaded, use this script to visualize the
## demo community data and run a multivariate analysis. To 
## complete these tasks, you will need to install and load
## the 'vegan' package:

install.packages("vegan")
library(vegan)

############################################################
### ORDINATION
############################################################

## In order to run a Non-metric Multidimensional Scaling
## (NMDS) ordination, you will need to subset, or section
## your data so that you have an object that only contains
## the species abundance information. This is called the
## sample x species matrix.

Samp <- Demo[ , 9:38]
View(Samp)

## Carefully review the above lines.

## When dealing with data frames or matrices in R, the data
## is held/indexed in DataFrame[rows, columns].

## You are using the index to pull every single row (the space
## following '[' ) in the matrix and only columns 9 to 38 (the
## range following the comma).

## Now you can run a NMDS ordination using the Bray-Curtis
## dissimilarity values.

Ord1 <- metaMDS(Samp, distance = "bray", k = 2, trymax = 20, 
                autotransform = TRUE, noshare = 0.1)

## After running the NMDS ordination, you can visualize the
## results using standard plotting practices.

## You might be interested in seeing if the community structure
## within each Terrain is different.

## To do this, you'll need to make an object representing the
## Terrain each sample is in.

## Since the Terrain variable (Demo$Terrain) is character data
## (strings of text), you need to convert it to a factor in
## order to use it correctly in the subsequent code:

Terra <- as.factor(Demo$Terrain)

## Double check the levels within the Terrain variable: 

length(unique(Demo$Terrain))

## Review the output in the console, which should report the
## value 5 - there are 5 levels or groups within this variable.
## You'll need the same number of colours as there are levels.

## Create a colour palette with the appropriate number of
## colours (one for each level):

ColPal <- c("seagreen", "dodgerblue", "firebrick2", 
            "lightgoldenrod1", "gray40")

## Apply the palette to your factor object and create a new,
## combined object:

ColPalTer <- ColPal[Terra]

## Now plot the NMDS ordination with the points coloured by
## Terrain (note that you should run each function in sequence):

plot(Ord1, display = c("sites", "species"), 
     choices = c(1, 2), type = "n", shrink = FALSE)
points(Ord1, display = "sites", choices = c(1, 2), 
       shrink = FALSE, pch = 16, col = ColPalTer)
legend("bottomright",
       legend = c("Forest", "Island", "Mountain", "Plains",
                  "Swamp"), 
       col = c("seagreen", "dodgerblue", "firebrick2", 
               "lightgoldenrod1", "gray40"),
       pch = 16, cex = 0.5)

############################################################
### ANALYSIS
############################################################

## The above ordination helps you visualize how similar the
## sites/samples are in terms of their species composition:
## i.e. WHICH species they have, and in what relative abundances.

## Samples that are close together in the ordination have low
## dissimilarity meaning they have similar animal communities.

## It SEEMS like at least 2 of the groups have quite different
## community composition based on their separation from each 
## other in the ordination diagram.

## BUT, you will need a way to test this statistically ...

## To do this, you can use something called PERMANOVA - which is
## short for:

## Permutational Multivariate Analysis of Variance 

## Use the code below to run a PERMANOVA to test for differences
## in community composition between groups:

PermAOV1 <- adonis2(Samp ~ Demo$Terrain, permutations = 999,
                    method = "bray", strata = NULL)

## Notice the equation notation used for the first argument -
## you're looking at sample composition as a function of Terrain.

## Take a look at the summary output, which will print in the
## console:

PermAOV1

## This tells you whether or not there is a difference between
## some of the groups/samples/sites, but not which groups!

## How can you determine exactly which groups are different?

## To do this, you need to create a function that conducts pairwise
## comparisons between groups.

## Highlight and run all the code below, starting on ~ line 164
## and ending with the } on ~ line 197:

pairwise.adonis <- function(x, factors, sim.function = 'vegdist', sim.method = 'bray',
                            p.adjust.m = 'bonferroni')
{
  library(vegan)
  
  co = combn(unique(as.character(factors)), 2)
  pairs = c()
  F.Model =c()
  R2 = c()
  p.value = c()
  
  for(elem in 1:ncol(co)){
    if(sim.function == 'daisy'){
      library(cluster); x1 = daisy(x[factors %in% c(co[1,elem], co[2,elem]), ], metric = sim.method)
    } else{x1 = vegdist(x[factors %in% c(co[1, elem], co[2, elem]), ], method = sim.method)}
    
    ad = adonis(x1 ~ factors[factors %in% c(co[1, elem], co[2, elem])] );
    pairs = c(pairs, paste(co[1, elem],'vs',co[2, elem]));
    F.Model = c (F.Model, ad$aov.tab[1, 4]);
    R2 = c(R2, ad$aov.tab[1, 5]);
    p.value = c(p.value, ad$aov.tab[1, 6])
  }
  p.adjusted = p.adjust(p.value, method = p.adjust.m)
  sig = c(rep('',length(p.adjusted)))
  sig[p.adjusted <= 0.05] <-'.'
  sig[p.adjusted <= 0.01] <-'*'
  sig[p.adjusted <= 0.001] <-'**'
  sig[p.adjusted <= 0.0001] <-'***'
  
  pairw.res = data.frame(pairs, F.Model, R2, p.value, p.adjusted, sig)
  print("Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1")
  return(pairw.res)
  
} 

## Here, you are creating a new function called pairwise.adonis()
## which will conduct pairwise comparisons of within vs. between
## group variation in composition for each PAIR of groups, where
## your sample x species matrix is selected using the x = argument 
## and the groups you want to compare are selected using the 
## factors = argument.

## Now use this function to run pairwise PERMANOVA tests on your 
## simulated community data:

Pairwise1 <- pairwise.adonis(x = Samp, factors = Demo$Terrain,
                             sim.function = 'vegdist',
                             sim.method = 'bray', p.adjust.m = 'fdr')

## Note that this may produce a list of errors describing the
## depreciation of 'adonis' and recommends using adonis2 - you
## can safely ignore these errors.

## Call the results to see which groups are significantly different
## from each other:

Pairwise1