# Loading necessary libraries
library(data.table)
library(magrittr)
library(spdep)
library(adespatial)
library(ade4)
library(adegraphics)
library(maptools)


setwd("") # set to directory with project files

coordinates <- read.delim2("coordinates.csv", header = T) %>% as.data.table()
coords <- coordinates[,.(x, y)] %>% as.matrix

# select a distance based spatial matrix based on maximising adjusted RČ
candidates <- listw.candidates(coords, nb = "dnear", d2=seq(from = 20, to = 100, by = 10))
listwSelected <- listw.select(x = coords,candidates = candidates) #a maximum distance of 100m generates MEMs with maximised adjRČ

# plot the best model
plot(listwSelected$best$MEM.select, SpORcoords = coords, nb = candidates$Dnear100_Linear)


# load the measured data
measurements <- read.delim2("data.csv", header = T, row.names = 1)

# add the spatial predictors to the data table
measurements$MEM1 <- listwSelected$best$MEM.select$MEM1
measurements$MEM2 <- listwSelected$best$MEM.select$MEM2

# detrend the dependent variables based on the generated spatial predictors
for(var in c("Ss", "Sf", "G", "R", "Na", "Ne", "Ho", "He", "F")) {
  measurements[paste0(var, "_detrended")] <- lm(as.formula(paste0(var, " ~ MEM1 + MEM2")), data = measurements) %>% resid()
}

# FOR THE FOLLOWING BAYESIAN MULTIPLE REGRESSION, UTILITY FUNCTIONS 
# PROVIDED BY KRUSCHKE (2015) WERE USED, WHICH WERE NOT REPRODUCED HERE
# AS SOURCE.
# THEY CAN BE OBTAINED FROM: https://sites.google.com/site/doingbayesiandataanalysis/software-installation
# THE BAYESIAN MODEL WAS ADJUSTED TO ESTIMATE THE REGRESSION PARAMETERS FROM
# A WEAKLY INFORMED T-DISTRIBUTION WITH PRECISION OF 4 AND NORMALITY PARAMETER
# OF 5. 4 INDEPENDENT CHAINS WERE SAMPLED.

# load utility functions and model
source("DBDA2E-utilities.R")
source("Jags-Ymet-XmetMulti-Mrobust.R")


# settings for MCMC chains; notation as used in Kruschke (2015)
myData = measurements
xName = c("El", "P", "K")
numSavedSteps=20000 
thinSteps=5
graphFileType = "png" 


# do MCMC sampling, create diagnostic plots and summarise results for all dependent variables
for(var in c("G", "Na", "Ne", "He")) {
  
  yName = paste0(var, "_detrended")
  fileNameRoot = paste0(var, "-Jags-") 
  
  # generate the MCMC chain
  mcmcCoda = genMCMC( data=myData , xName=xName , yName=yName , 
                      numSavedSteps=numSavedSteps , thinSteps=thinSteps , 
                      saveName=fileNameRoot, nChains = 4)
  # display diagnostics of chain, for specified parameters
  parameterNames = varnames(mcmcCoda)
  for ( parName in parameterNames ) {
    diagMCMC( codaObject=mcmcCoda , parName=parName , 
              saveName=fileNameRoot , saveType=graphFileType )
    graphics.off()
  }

  # get summary statistics of chain
  summaryInfo = smryMCMC( mcmcCoda , 
                          saveName=fileNameRoot  )
  show(summaryInfo)
  # display posterior information
  plotMCMC( mcmcCoda , data=myData , xName=xName , yName=yName , 
            pairsPlot=TRUE , showCurve=FALSE ,
            saveName=fileNameRoot , saveType=graphFileType )

  graphics.off()
}