knitr::opts_chunk$set(echo = TRUE)

Aim

Summarise the results of all GDM models listed in the table model_names and create a results table GDM_devexpl, as well as a table of all isplines in the table isplines_all_models. Serves as input for assessing the thresholds and comparing quality of the models, and the summary models across thresholds.

Requirements

run all the previous scripts to generate necessary output (see "the_analysis_hive.Rmd")

if(!exists("pathtodata")){
  sections_to_be_loaded <- c()
  source("vignettes/analysis_nonpublic.R")
}

Summarise deviance explained

Aim : Summarise deviance explained for all single predictors.

NOTE : here, potentially, we could include the maxsplines values as well. But they are already included in the table "isplines_all_models". No need to duplicate its calculation. If wanted, check version control of 16.12.22 for keywords MAJOR CHANGE and "maxsplines" in this script.

# Prepare table for deviance explained and variance of response by model
devexpl <- data.table(modelname = model_names[model_class %in% c("multifun", "singlefun") & lui == "LUI", modelname],
           GDM_dev_expl = as.numeric(NA))

for(mod in model_names[model_class %in% c("multifun", "singlefun") & lui == "LUI", modelname]){
  print(mod) # note : added this print statement because recently problems with readRDS() function (July 2023)
  model_names_selection <- model_names[which(model_names$modelname == mod), ]
  sections_to_be_loaded <- c("gdmoutput") # loads `gdmoutput` and `model_specs`
  source("vignettes/analysis_nonpublic.R")

  # add values to results table
  devexpl[modelname == mod, GDM_dev_expl := gdmoutput$explained]
}
rm(mod, model_names_selection, sections_to_be_loaded, gdmoutput, model_specs)

saveRDS(devexpl, file = paste(pathtoout, "/devexpl_all_models.RDS", sep = "")) # note this dataset has been called "model_results.CSV" in the past.
# is read in by analysis_nonpublic.R, "results" section into variable `model_results`)

Summarise all multi-threshold models

We aim to calculate an average spline over all thresholds. Therefore, we need the individual spline values for each threshold.

isplineExtract calculates predicted values for each predictor individually (the x and y values for isplines). The x values are just 200 samples from the predictor range (between predictor, e.g. LUI min and LUI max, take 200 equidistant values). The x values are the same across models, because all models use the same predictors (and thus the same range, i.e. min and max, of predicotrs). We can thus extract the x values of one predictor to form a results data table. We then go through each model and add the spline values for each of the model. For merging the models, we can use the x value of the predictor.

source("vignettes/analysis_nonpublic.R")

# # #
# PREPARE RESULTS TABLE
#
# take random model, only extract the X values and use as a basis for the results table.
mod <- "gdm_EFturnover_0.1_LUI"
model_names_selection <- model_names[which(model_names$modelname == mod), ]
sections_to_be_loaded <- c("gdmoutput") # loads `gdmoutput` and `model_specs`
source("vignettes/analysis_nonpublic.R")
# calculate maxsplines
exSplines <- gdm::isplineExtract(gdmoutput)
isplines_all_models <- melt(data.table(exSplines$x), measure.vars = colnames(exSplines$x),
     variable.name = "predictor", value.name = "Xpredictor")
rm(mod, model_names_selection, sections_to_be_loaded, exSplines)

# # #
# FILL RESULTS TABLE
#
# add all spline values from all modesl to this table --> form a huge table where
# columns are the models and rows are individual values of all predictors.
for(mod in model_names[model_class %in% c("multifun") & lui == "LUI", modelname]){
  model_names_selection <- model_names[which(model_names$modelname == mod), ]
  sections_to_be_loaded <- c("gdmoutput") # loads `gdmoutput` and `model_specs`
  source("vignettes/analysis_nonpublic.R")

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

  # get data for splines
  # x values : 200 values are created per predictor, with equal distances between min and max predictor
  # y values : the predicted values of the given predictor.
  # We extract both x and y values, merge them to a combined data table and use the x values for merging with the full results dataset.
  # get x and y values of current model in desired format
  Xpreds <- melt(data.table(exSplines$x), measure.vars = colnames(exSplines$x),
     variable.name = "predictor", value.name = "Xpredictor")
  Ypreds <- melt(data.table(exSplines$y), measure.vars = colnames(exSplines$x),
     variable.name = "predictor", value.name = mod)
  XYpreds <- cbind(Xpreds, Ypreds[, ..mod])
  #
  # Add the extracted x and y predictors to the results table
  isplines_all_models <- merge(isplines_all_models, XYpreds, by = c("predictor", "Xpredictor"))
}
rm(mod, model_names_selection, sections_to_be_loaded, exSplines, Xpreds, Ypreds, XYpreds)

# save the output
saveRDS(isplines_all_models, file = paste(pathtoout, "/isplines_all_threshold_models.RDS", sep = ""))


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