In allanecology/BetaDivMultifun: Testing effects of multitrophic beta-diversity and land-use intensity (LUI) on ecosystem multifunctionality

GDM Predictions

Run predictions from the GDMs in order to visualise the results (analogously to the effects package for mixed effects models).

Requirements

From the file analysis_nonpublic.R.

sections_to_be_loaded <- c("gdminput", "gdmoutput")
funs <- "EFdistance"
compon_lui <- "LUI"
source("vignettes/analysis_nonpublic.R")
library(gdm)

EFdist

exSplines <- gdm::isplineExtract(gdmoutput)
maxsplines <- apply(exSplines$y, 2, max, na.rm = TRUE)

gdminput$predicted.distance<- predict(gdmoutput, gdminput) # is actually the same as gdmoutput$predicted

dtgdminput <- data.table::data.table(gdminput)

# check : plot pred vs. obs
plot(gdminput$predicted.distance, gdminput$distance)

# visualise results
ggplot2::ggplot(gdminput, aes( x = s1.LUI, y = s2.LUI, col = predicted.distance)) +
  geom_point()
# still quite messy. Try to clean out before plotting

# selecting ecologically interesting pairs of LUI dissim :
# - low-low
# - medium- medium
# - high-high
lowlui <- dtgdminput[s1.LUI < 1 & s2.LUI < 1 | s1.LUI > 1 & s2.LUI < 1, ]
ggplot2::ggplot(dtgdminput, aes( x = s1.LUI, y = s2.LUI, col = predicted.distance)) +
  geom_point()

plot(abs(gdminput$s1.LUI - gdminput$s2.LUI), gdminput$predicted.distance)

# just checking how it's in linear model.
# we need to plot the residuals of LUI in order to get the exact model values...
mod <- lm(distance ~ s1.LUI + s2.LUI, data = gdminput)
plot(predict(mod), gdminput$distance)
plot(predict(mod), gdminput$s2.LUI)

# maybe need to work with the partial effects? Yes, probably!
plot(exSplines[[1]][,"LUI"], exSplines[[2]][,"LUI"]) # second item : partial ecological distance
length(exSplines[[1]][,"LUI"]) # why just 200 values? either too much or too few... ??
# don't forget the intercept!

# what is a rug plot in gdm.plot() ?

From Vignette

# As in Section 1, we first fit a gdm using raster layers as predictors
# Load data from the gdm package
gdmExpData <- southwest
rastFile <- system.file("./extdata/swBioclims.grd", package="gdm")
# create a raster stack
envRast <- stack(rastFile)
# Create a 'species list' data input using the gdm package data
sppTab <- gdmExpData[, c("species", "site", "Lat", "Long")]
# prepare the gdm input table, using rasters as predictors
sitePairRast <- formatsitepair(bioData=sppTab, bioFormat=2, XColumn="Long", YColumn="Lat", sppColumn="species", siteColumn="site", predData=envRast)

# Remove any site-pairs containing NAs for the extracted raster-based predictors
sitePairRast <- na.omit(sitePairRast)
# fit the GDM
gdmRastMod <- gdm(data=sitePairRast,geo=TRUE)
# Generate transformed predictor rasters, based on the raw raster predictor layers
# and the fitted gdm
transRasts <- gdm.transform(model=gdmRastMod, data=envRast)

# Get the data from the gdm transformed rasters as a table
rastDat <- na.omit(getValues(transRasts))
# The PCA can be fit on a sample of grid cells if the rasters are large
rastDat <- sampleRandom(transRasts, 50000)
# perform the principle components analysis
pcaSamp <- prcomp(rastDat)
# Predict the first three principle components for every cell in the rasters
# note the use of the 'index' argument
pcaRast <- predict(transRasts, pcaSamp, index=1:3)
# scale the PCA rasters to make full use of the colour spectrum
pcaRast[[1]] <- (pcaRast[[1]]-pcaRast[[1]]@data@min) /
(pcaRast[[1]]@data@max-pcaRast[[1]]@data@min)*255
pcaRast[[2]] <- (pcaRast[[2]]-pcaRast[[2]]@data@min) /
(pcaRast[[2]]@data@max-pcaRast[[2]]@data@min)*255
pcaRast[[3]] <- (pcaRast[[3]]-pcaRast[[3]]@data@min) /
(pcaRast[[3]]@data@max-pcaRast[[3]]@data@min)*255
# Plot the three PCA rasters simultaneously, each representing a different colour
# (red, green, blue)
plotRGB(pcaRast, r=1, g=2, b=3)

What does the transform function do?

# use predictors
test <- merge(LUI[, 1:2], predictors[, 1:2], by = "Plotn")

gdm.transform(gdmoutput, test)

allanecology/BetaDivMultifun documentation built on Nov. 9, 2023, 8:47 p.m.