knitr::opts_chunk$set(echo = TRUE)
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.
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") }
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`)
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 = ""))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.