Visualize effects for binomial models

Load packages, read data and source custom scripts

rm(list = ls())
library(mgcv)

#> Loading required package: nlme

#> This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.

library(dplyr)

#> 
#> Attaching package: 'dplyr'

#> The following object is masked from 'package:nlme':
#> 
#>     collapse

#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag

#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

library(purrr)
library(ggplot2)
library(colorspace)

path_proj <- day2day::git_path()
path_data <- file.path(path_proj, "data")
path_processed <- file.path(path_data, "processed")
path_modelled <- file.path(path_data, "modelled")

source(file.path(path_proj, "src", "68-mgcv-visualize.R"))
source(file.path(path_proj, "src", "63-tidy-names.R"))

land_cover <- fst::read.fst(file.path(path_processed, "uk_1km_dataframe_train_maskout.fst"))

models <- list(readRDS(file.path(path_modelled, "model_binomials_maskout_interact.rds")),
               readRDS(file.path(path_modelled, "model_binomials_maskout_simple.rds")))

models <- map(models, ~ mutate(., name = tidy_names_models(.$name)))

Custom function

plot_land_effects <- function(models, index, covars1, covars2, free_fill = FALSE,
                              rug = TRUE, rug_prop = 0.01, n = 80, n2 = 80) {
    cov_data <- map(
        models, ~ get_data_covs(.$model[[index]], covars1, model = .$name[[index]],
                                rug = rug, rug_prop = rug_prop, n = n, n2 = n2)
        ) %>% do.call(c, .)
    gg1 <- plot_gam_1d(cov_data) & guides(fill = "none")
    cov_data <- map(
        models, ~ get_data_covs(.$model[[index]], covars2, model = .$name[[index]],
                                rug = rug, rug_prop = rug_prop, n = n, n2 = n2)
        ) %>% do.call(c, .)

    gg2 <- plot_gam_2d(cov_data, free_fill = free_fill)
    return(list(gg1, gg2))
}

Arable

covars1 <- list("gdd", "elev", "slope_nb8")
covars2 <- list(c("min_tas", "max_tas"), c("smd", "sms"), c("elev", "slope_nb8"))
ggs <- plot_land_effects(models, 1, covars1, covars2, free_fill = TRUE)

#> Loading required package: patchwork

One dimensional effects.

ggs[[1]]

Two dimensional effects.

ggs[[2]] &
    theme(legend.position = "bottom", legend.key.width = unit(2, "lines"))

Wetland

covars1 <- list("gdd", "elev", "slope_nb8")
covars2 <- list(c("min_tas", "max_tas"), c("smd", "sms"), c("elev", "slope_nb8"))
ggs <- plot_land_effects(models, 2, covars1, covars2, free_fill = TRUE)

One dimensional effects.

ggs[[1]]

Two dimensional effects.

ggs[[2]] &
    theme(legend.position = "bottom", legend.key.width = unit(2, "lines"))

Improved grassland

covars1 <- list("gdd")
covars2 <- list(c("min_tas", "max_tas"), c("smd", "sms"))
ggs <- plot_land_effects(models, 3, covars1, covars2, free_fill = TRUE)

One dimensional effects.

ggs[[1]]

Two dimensional effects.

ggs[[2]] &
    theme(legend.position = "bottom", legend.key.width = unit(2, "lines"))

Forest

covars1 <- list("gdd")
covars2 <- list(c("min_tas", "max_tas"), c("smd", "sms"))
ggs <- plot_land_effects(models, 4, covars1, covars2, free_fill = TRUE)

One dimensional effects.

ggs[[1]]

Two dimensional effects.

ggs[[2]] &
    theme(legend.position = "bottom", legend.key.width = unit(2, "lines"))

Semi natural grassland

covars1 <- list("gdd")
covars2 <- list(c("min_tas", "max_tas"), c("smd", "sms"))
ggs <- plot_land_effects(models, 5, covars1, covars2, free_fill = TRUE)

One dimensional effects.

ggs[[1]]

Two dimensional effects.

ggs[[2]] &
    theme(legend.position = "bottom", legend.key.width = unit(2, "lines"))

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#> 105.657   1.450 107.288