From 99f93c9bd9402e23d621ff2d92d0e3d8d66a5f75 Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Thu, 21 Sep 2023 13:32:54 +0100 Subject: [PATCH 01/21] WIP: Looking at adding prediction mapping --- NAMESPACE | 3 --- R/model_parse.R | 6 ++++++ R/util_coords.R | 2 +- R/util_paths.R | 2 +- R/util_shiny.R | 2 +- man/get_tmpdir.Rd | 1 + vignettes/covid.Rmd | 17 +++++++++++++++++ 7 files changed, 27 insertions(+), 6 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 49236911..39f576d4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,10 +1,7 @@ # Generated by roxygen2: do not edit by hand export(clear_caches) -export(get_busy_spinner) -export(get_tmpdir) export(get_tutorial_datapath) -export(has_coords) export(interactive_priors) export(latlong_to_utm) export(load_tutorial_data) diff --git a/R/model_parse.R b/R/model_parse.R index e02cfe81..78511975 100644 --- a/R/model_parse.R +++ b/R/model_parse.R @@ -6,9 +6,15 @@ #' @return list #' @keywords internal parse_model_output_bru <- function(model_output, measurement_data) { + cat(nrow(measurement_data), nrow(model_output$summary.fitted.values)) + return(NULL) + fitted_mean_post <- model_output$summary.fitted.values$mean[seq_len(nrow(measurement_data))] fitted_sd_post <- model_output$summary.fitted.values$mean[seq_len(nrow(measurement_data))] + pred.25 <- inlabru_model$summary.fitted.values$`0.025quant`[1:nrow(covid19_data)] + pred.975 <- inlabru_model$summary.fitted.values$`0.975quant`[1:nrow(covid19_data)] + mean_post <- model_output$summary.random$f$mean sd_post <- model_output$summary.random$f$sd fixed_mean <- model_output$summary.fixed$mean diff --git a/R/util_coords.R b/R/util_coords.R index 8024a0dc..07c7a649 100644 --- a/R/util_coords.R +++ b/R/util_coords.R @@ -28,7 +28,7 @@ latlong_to_utm <- function(lat, lon) { #' @param spatial_data #' #' @return bool -#' @export +#' @keywords internal has_coords <- function(spatial_data) { tryCatch( { diff --git a/R/util_paths.R b/R/util_paths.R index fea0de9a..f737f958 100644 --- a/R/util_paths.R +++ b/R/util_paths.R @@ -95,7 +95,7 @@ clear_caches <- function() { #' Get the system temporary directory #' #' @return NULL -#' @export +#' @keywords internal get_tmpdir <- function() { fs::path_dir(fs::path_temp()) } diff --git a/R/util_shiny.R b/R/util_shiny.R index f3d734f8..0d733d7e 100644 --- a/R/util_shiny.R +++ b/R/util_shiny.R @@ -3,7 +3,7 @@ #' TODO - add this gif to the package #' #' @return busyspinner -#' @export +#' @keywords internal get_busy_spinner <- function() { if (curl::has_internet()) { busy_spinner <- shinybusy::add_busy_gif("https://raw.githubusercontent.com/4DModeller/logo/main/4DMlogo_loading.gif", height = 100, width = 100, position = "top-right") diff --git a/man/get_tmpdir.Rd b/man/get_tmpdir.Rd index 15acca21..7ccccaa6 100644 --- a/man/get_tmpdir.Rd +++ b/man/get_tmpdir.Rd @@ -9,3 +9,4 @@ get_tmpdir() \description{ Get the system temporary directory } +\keyword{internal} diff --git a/vignettes/covid.Rmd b/vignettes/covid.Rmd index d483c1d7..a8eeb721 100644 --- a/vignettes/covid.Rmd +++ b/vignettes/covid.Rmd @@ -28,7 +28,9 @@ We first load INLA and then retrieve the data from the fdmr example data store. ```{r} library(INLA) +library(magrittr) ``` + ```{r} fdmr::retrieve_tutorial_data(dataset = "covid") ``` @@ -412,6 +414,21 @@ sp_data@data$ave.risk <- average_risk_by_nb$ave.risk domain <- sp_data@data$ave.risk +fdmr::plot_map( + polygon_data = sp_data, + domain = domain, + palette = "Reds", + legend_title = "Risk", + add_scale_bar = TRUE, + polygon_fill_opacity = 1, + polygon_line_colour = "transparent" +) <- + dplyr::group_by(predictions, MSOA11CD) %>% dplyr::summarize(ave.risk = mean(pred.mean)) + +sp_data@data$ave.risk <- average_risk_by_nb$ave.risk + +domain <- sp_data@data$ave.risk + fdmr::plot_map( polygon_data = sp_data, domain = domain, From 3a4ceb366e539bc88fc8f853356b017b5ca2a9a3 Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Wed, 27 Sep 2023 13:36:43 +0100 Subject: [PATCH 02/21] Remove duplicated code --- vignettes/covid.Rmd | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/vignettes/covid.Rmd b/vignettes/covid.Rmd index a8eeb721..aa26d13a 100644 --- a/vignettes/covid.Rmd +++ b/vignettes/covid.Rmd @@ -414,21 +414,6 @@ sp_data@data$ave.risk <- average_risk_by_nb$ave.risk domain <- sp_data@data$ave.risk -fdmr::plot_map( - polygon_data = sp_data, - domain = domain, - palette = "Reds", - legend_title = "Risk", - add_scale_bar = TRUE, - polygon_fill_opacity = 1, - polygon_line_colour = "transparent" -) <- - dplyr::group_by(predictions, MSOA11CD) %>% dplyr::summarize(ave.risk = mean(pred.mean)) - -sp_data@data$ave.risk <- average_risk_by_nb$ave.risk - -domain <- sp_data@data$ave.risk - fdmr::plot_map( polygon_data = sp_data, domain = domain, From 26705b8957230e0f4570a10935941852988b66ab Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Fri, 29 Sep 2023 10:42:07 +0100 Subject: [PATCH 03/21] WIP: First pass at Shiny model viewer --- NAMESPACE | 2 +- R/shiny_modelviewer.R | 90 +++++++++++++++++++++++++++++++++++++++++++ vignettes/hydro.Rmd | 6 ++- 3 files changed, 96 insertions(+), 2 deletions(-) create mode 100644 R/shiny_modelviewer.R diff --git a/NAMESPACE b/NAMESPACE index 7972db2c..f733d69d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,7 +1,7 @@ # Generated by roxygen2: do not edit by hand export(clear_caches) -export(get_tmpdir) +export(create_prediction_field) export(get_tutorial_datapath) export(interactive_priors) export(latlong_to_utm) diff --git a/R/shiny_modelviewer.R b/R/shiny_modelviewer.R new file mode 100644 index 00000000..49c0b499 --- /dev/null +++ b/R/shiny_modelviewer.R @@ -0,0 +1,90 @@ +#' Parse inlabru model output +#' +#' @param +#' @importFrom magrittr %>% +#' +#' @return shiny::app +#' @keywords internal +model_viewer_shiny <- function(model_output, measurement_data, mesh) { + busy_spinner <- get_busy_spinner() + + parsed_model_output <- parse_model_output(model_output = model_output, measurement_data = measurement_data) + + plot_choices <- c("Range", "Stdev", "AR(1)", "Boxplot", "Density", "DIC") + + parsed_names <- names(parsed_model_output) + map_vars <- parsed_names[!parsed_names %in% c("dic", "pars")] + + ui <- shiny::fluidPage( + busy_spinner, + shiny::headerPanel(title = "Model viewer"), + shiny::sidebarLayout( + shiny::sidebarPanel( + shiny::selectInput( + inputId = "plot_type", + label = "Model type:", + choices = c("inlabru"), + selected = "inlabru" + ), + ), + shiny::mainPanel( + shiny::tabsetPanel( + type = "tabs", + shiny::tabPanel( + "Plot", + shiny::h2("Plot output"), + shiny::selectInput(inputId = "plot_type", label = "Plot type:", choices = plot_choices, selected = plot_choices[1]), + shiny::plotOutput(outputId = "plot_model_out") + ), + shiny::tabPanel( + "Map", + shiny::selectInput(inputId = "select_run_map", label = "Select run:", choices = c()), + shiny::selectInput(inputId = "map_var_a", label = "Variable a:", choices = map_vars), + shiny::selectInput(inputId = "map_var_b", label = "Variable b:", choices = map_vars), + leaflet::leafletOutput(outputId = "map_out") + ), + shiny::tabPanel( + "Help", + shiny::h3("Help"), + ) + ) + ) + ) + ) + + # Define server logic required to draw a histogram + server <- function(input, output, session) { + map_raster <- shiny::reactive({ + pred_field <- create_prediction_field( + var_a = parsed_model_output[[input$map_var_a]], + var_b = parsed_model_output[[input$map_var_b]], + mesh = mesh + ) + + crs <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" + create_raster(dataframe = pred_field, crs = sp::CRS(crs)) + }) + + output$map_out <- leaflet::renderLeaflet({ + leaflet::leaflet() %>% + leaflet::addTiles(group = "OSM") %>% + leaflet::addRasterImage(map_raster(), opacity = 0.9, group = "Raster") + }) + # Run the application + shiny::shinyApp(ui = ui, server = server) + } +} + +#' Mesh building shiny app. Creates and visualises a mesh from some spatial data. +#' +#' @param spatial_data Spatial datas +#' +#' @return shiny::app +#' @export +model_viewer <- function(spatial_data, measurement_data, mesh) { + shiny::runApp(model_viewer_shiny( + spatial_data = spatial_data, + measurement_data = measurement_data, + mesh = mesh + )) +} diff --git a/vignettes/hydro.Rmd b/vignettes/hydro.Rmd index 10f5acd4..19eaad20 100644 --- a/vignettes/hydro.Rmd +++ b/vignettes/hydro.Rmd @@ -315,7 +315,11 @@ inlabru_model2 <- inlabru::bru(formula2, ) ``` -Now it's time to see its output +Now it's time to see its output, we can use the `fdmr::model_viewer` and pass in our model output, measurement data and mesh. This Shiny app helps analyse model output by automatically creating plots and a map of predictions for your data. + +```{r eval=FALSE} +fdmr::model_viewer() +``` ```{r error=TRUE} # TODO : this could be where the model_eval function I recommend could be used From 7e9444f5c05629c0a7885e6f0d6c975700c8b01b Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Fri, 29 Sep 2023 12:16:44 +0100 Subject: [PATCH 04/21] WIP: Adding better parsing --- R/model_parse.R | 13 +++++-------- R/shiny_modelviewer.R | 25 +++++++++++++------------ vignettes/hydro.Rmd | 20 ++++++++++---------- 3 files changed, 28 insertions(+), 30 deletions(-) diff --git a/R/model_parse.R b/R/model_parse.R index 2b3618d2..7f90c947 100644 --- a/R/model_parse.R +++ b/R/model_parse.R @@ -6,9 +6,6 @@ #' @return list #' @keywords internal parse_model_output_bru <- function(model_output, measurement_data) { - cat(nrow(measurement_data), nrow(model_output$summary.fitted.values)) - return(NULL) - fitted_mean_post <- model_output$summary.fitted.values$mean[seq_len(nrow(measurement_data))] fitted_sd_post <- model_output$summary.fitted.values$sd[seq_len(nrow(measurement_data))] @@ -53,9 +50,9 @@ parse_model_output <- function(model_output, measurement_data, model_type = "inl #' Create a prediction field from the parsed model output and the mesh #' -#' @param var_a -#' @param var_b -#' @param mesh +#' @param var_a +#' @param var_b +#' @param mesh #' #' @return data.frame #' @export @@ -64,11 +61,11 @@ create_prediction_field <- function(var_a, var_b, mesh) { xy_grid <- base::expand.grid(mod_proj$x, mod_proj$y) A_proj <- INLA::inla.spde.make.A(mesh = mesh, loc = as.matrix(xy_grid)) - z <- base::exp(base::as.numeric(A_proj %*% var_a[1:mesh$n]) + base::sum(var_b)) + z <- base::exp(base::as.numeric(A_proj %*% var_a[1:mesh$n]) + base::sum(var_b)) base::data.frame(x = xy_grid[, 1], y = xy_grid[, 2], z = z) } create_raster <- function(dataframe, crs) { raster::rasterFromXYZ(dataframe, crs = crs) -} \ No newline at end of file +} diff --git a/R/shiny_modelviewer.R b/R/shiny_modelviewer.R index 49c0b499..d028c28d 100644 --- a/R/shiny_modelviewer.R +++ b/R/shiny_modelviewer.R @@ -5,16 +5,19 @@ #' #' @return shiny::app #' @keywords internal -model_viewer_shiny <- function(model_output, measurement_data, mesh) { +model_viewer_shiny <- function(model_output, mesh, measurement_data) { busy_spinner <- get_busy_spinner() - parsed_model_output <- parse_model_output(model_output = model_output, measurement_data = measurement_data) + parsed_model_output <- parse_model_output( + model_output = model_output, + measurement_data = measurement_data + ) plot_choices <- c("Range", "Stdev", "AR(1)", "Boxplot", "Density", "DIC") parsed_names <- names(parsed_model_output) map_vars <- parsed_names[!parsed_names %in% c("dic", "pars")] - + browser() ui <- shiny::fluidPage( busy_spinner, shiny::headerPanel(title = "Model viewer"), @@ -38,7 +41,6 @@ model_viewer_shiny <- function(model_output, measurement_data, mesh) { ), shiny::tabPanel( "Map", - shiny::selectInput(inputId = "select_run_map", label = "Select run:", choices = c()), shiny::selectInput(inputId = "map_var_a", label = "Variable a:", choices = map_vars), shiny::selectInput(inputId = "map_var_b", label = "Variable b:", choices = map_vars), leaflet::leafletOutput(outputId = "map_out") @@ -71,20 +73,19 @@ model_viewer_shiny <- function(model_output, measurement_data, mesh) { leaflet::addRasterImage(map_raster(), opacity = 0.9, group = "Raster") }) # Run the application - shiny::shinyApp(ui = ui, server = server) } + + shiny::shinyApp(ui = ui, server = server) } #' Mesh building shiny app. Creates and visualises a mesh from some spatial data. #' -#' @param spatial_data Spatial datas +#' @param model_output INLA model output +#' @param mesh INLA mesh +#' @param measurement_data Measurement data #' #' @return shiny::app #' @export -model_viewer <- function(spatial_data, measurement_data, mesh) { - shiny::runApp(model_viewer_shiny( - spatial_data = spatial_data, - measurement_data = measurement_data, - mesh = mesh - )) +model_viewer <- function(model_output, mesh, measurement_data) { + shiny::runApp(model_viewer_shiny(model_output = model_output, mesh = mesh, measurement_data = measurement_data)) } diff --git a/vignettes/hydro.Rmd b/vignettes/hydro.Rmd index 19eaad20..f616d0f7 100644 --- a/vignettes/hydro.Rmd +++ b/vignettes/hydro.Rmd @@ -298,27 +298,27 @@ formula2 <- streamflow ~ 0 + Intercept + precip + ```{r error=TRUE, modelrun} -inlabru_model1 <- inlabru::bru(formula1, +model_out_a <- inlabru::bru(formula1, data = inla_data, family = "gaussian", options = list( - verbose = TRUE + verbose = FALSE ) ) -inlabru_model2 <- inlabru::bru(formula2, - data = inla_data, - family = "gaussian", - options = list( - verbose = TRUE - ) -) +# inlabru_model2 <- inlabru::bru(formula2, +# data = inla_data, +# family = "gaussian", +# options = list( +# verbose = TRUE +# ) +# ) ``` Now it's time to see its output, we can use the `fdmr::model_viewer` and pass in our model output, measurement data and mesh. This Shiny app helps analyse model output by automatically creating plots and a map of predictions for your data. ```{r eval=FALSE} -fdmr::model_viewer() +fdmr::model_viewer(model_output = model_out_a, mesh = mesh, measurement_data = inla_data) ``` ```{r error=TRUE} From 75ad316811f036371a342c7db7b1546a15de4e52 Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Fri, 29 Sep 2023 12:25:08 +0100 Subject: [PATCH 05/21] Remove unused variable read --- R/model_parse.R | 3 --- 1 file changed, 3 deletions(-) diff --git a/R/model_parse.R b/R/model_parse.R index 7f90c947..3e947e97 100644 --- a/R/model_parse.R +++ b/R/model_parse.R @@ -9,9 +9,6 @@ parse_model_output_bru <- function(model_output, measurement_data) { fitted_mean_post <- model_output$summary.fitted.values$mean[seq_len(nrow(measurement_data))] fitted_sd_post <- model_output$summary.fitted.values$sd[seq_len(nrow(measurement_data))] - pred.25 <- inlabru_model$summary.fitted.values$`0.025quant`[1:nrow(covid19_data)] - pred.975 <- inlabru_model$summary.fitted.values$`0.975quant`[1:nrow(covid19_data)] - mean_post <- model_output$summary.random$f$mean sd_post <- model_output$summary.random$f$sd fixed_mean <- model_output$summary.fixed$mean From 8c8066d0d787de069d70dfca57200040bcd649e1 Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Fri, 29 Sep 2023 14:39:32 +0100 Subject: [PATCH 06/21] WIP: Looking at change of behaviour since moving from rgdal --- R/model_parse.R | 6 +++--- R/shiny_modelviewer.R | 10 +++++++--- vignettes/hydro.Rmd | 17 +++++++++-------- 3 files changed, 19 insertions(+), 14 deletions(-) diff --git a/R/model_parse.R b/R/model_parse.R index 3e947e97..5120c187 100644 --- a/R/model_parse.R +++ b/R/model_parse.R @@ -47,9 +47,9 @@ parse_model_output <- function(model_output, measurement_data, model_type = "inl #' Create a prediction field from the parsed model output and the mesh #' -#' @param var_a -#' @param var_b -#' @param mesh +#' @param var_a Variable a data +#' @param var_b Variable b data +#' @param mesh INLA mesh #' #' @return data.frame #' @export diff --git a/R/shiny_modelviewer.R b/R/shiny_modelviewer.R index d028c28d..5e77bcaa 100644 --- a/R/shiny_modelviewer.R +++ b/R/shiny_modelviewer.R @@ -1,6 +1,9 @@ #' Parse inlabru model output #' -#' @param +#' @param model_output INLA model output +#' @param mesh INLA mesh +#' @param measurement_data Measurement data +#' #' @importFrom magrittr %>% #' #' @return shiny::app @@ -17,7 +20,7 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { parsed_names <- names(parsed_model_output) map_vars <- parsed_names[!parsed_names %in% c("dic", "pars")] - browser() + ui <- shiny::fluidPage( busy_spinner, shiny::headerPanel(title = "Model viewer"), @@ -64,7 +67,8 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { ) crs <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" - create_raster(dataframe = pred_field, crs = sp::CRS(crs)) + raster = create_raster(dataframe = pred_field, crs = sp::CRS(crs)) + browser() }) output$map_out <- leaflet::renderLeaflet({ diff --git a/vignettes/hydro.Rmd b/vignettes/hydro.Rmd index 5d9226a8..0b334ccd 100644 --- a/vignettes/hydro.Rmd +++ b/vignettes/hydro.Rmd @@ -28,7 +28,8 @@ First we will look at the location of the dam by making a map and plotting the s ```{r error=TRUE,warning=FALSE} # TODO : here is a good example for the map function norway_polygon_path <- fdmr::get_tutorial_datapath(dataset = "hydro", filename = "Kvilldal_Catch_Boundary.geojson") -norway_polygon <- sf::read_sf(norway_polygon_path) %>% sf::st_zm() +#norway_polygon <- sf::read_sf(norway_polygon_path) %>% sf::st_zm() +norway_polygon <- rgdal::readOGR(norway_polygon_path) sfc <- sf::st_transform(norway_polygon, crs = "+proj=longlat +datum=WGS84") ``` We want to mark the position of the dam and the stream gauges on the map so we'll create a `data.frame` to pass into the `plot_map` @@ -301,13 +302,13 @@ model_out_a <- inlabru::bru(formula1, ) ) -# inlabru_model2 <- inlabru::bru(formula2, -# data = inla_data, -# family = "gaussian", -# options = list( -# verbose = TRUE -# ) -# ) +inlabru_model2 <- inlabru::bru(formula2, + data = inla_data, + family = "gaussian", + options = list( + verbose = TRUE + ) +) ``` Now it's time to see its output, we can use the `fdmr::model_viewer` and pass in our model output, measurement data and mesh. This Shiny app helps analyse model output by automatically creating plots and a map of predictions for your data. From 53bf111d976836c793a878be17a4e5d77113452d Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Fri, 6 Oct 2023 09:31:32 +0100 Subject: [PATCH 07/21] WIP: Adding model viewer work --- R/model_parse.R | 15 +++++++-- R/plot_priors.R | 5 ++- R/shiny_modelviewer.R | 55 +++++++++++++++++++++++++++----- vignettes/data_preprocessing.Rmd | 2 +- vignettes/hydro.Rmd | 36 +++++++++++++++------ 5 files changed, 91 insertions(+), 22 deletions(-) diff --git a/R/model_parse.R b/R/model_parse.R index 5120c187..6845a8fd 100644 --- a/R/model_parse.R +++ b/R/model_parse.R @@ -50,11 +50,22 @@ parse_model_output <- function(model_output, measurement_data, model_type = "inl #' @param var_a Variable a data #' @param var_b Variable b data #' @param mesh INLA mesh +#' @param crs CRS as a proj4string #' #' @return data.frame #' @export -create_prediction_field <- function(var_a, var_b, mesh) { - mod_proj <- fmesher::fm_evaluator(mesh) +create_prediction_field <- function(var_a, var_b, mesh, crs = NULL) { + if (is.null(crs)) { + read_crs <- mesh$crs$input + if (is.null(read_crs)) { + warning("Cannot read CRS from mesh, pass proj4string to crs otherwise we'll assume crs = +proj=longlat +datum=WGS84") + crs <- "+proj=longlat +datum=WGS84" + } else { + crs <- read_crs + } + } + + mod_proj <- fmesher::fm_evaluator(mesh, crs = crs) xy_grid <- base::expand.grid(mod_proj$x, mod_proj$y) A_proj <- INLA::inla.spde.make.A(mesh = mesh, loc = as.matrix(xy_grid)) diff --git a/R/plot_priors.R b/R/plot_priors.R index 40a34772..40384e87 100644 --- a/R/plot_priors.R +++ b/R/plot_priors.R @@ -12,10 +12,13 @@ plot_line_comparison <- function(data, to_plot, title) { return("No pars data.") } - ggplot2::ggplot(single_df, ggplot2::aes(x = x, y = y, color = Run)) + + plot <- ggplot2::ggplot(single_df, ggplot2::aes(x = x, y = y, color = Run)) + ggplot2::geom_line() + ggplot2::ggtitle(title) + ggplot2::theme(text = ggplot2::element_text(size = 16)) + + browser() + plot } diff --git a/R/shiny_modelviewer.R b/R/shiny_modelviewer.R index 5e77bcaa..da1f2f49 100644 --- a/R/shiny_modelviewer.R +++ b/R/shiny_modelviewer.R @@ -3,7 +3,7 @@ #' @param model_output INLA model output #' @param mesh INLA mesh #' @param measurement_data Measurement data -#' +#' #' @importFrom magrittr %>% #' #' @return shiny::app @@ -16,11 +16,17 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { measurement_data = measurement_data ) + crs <- mesh$crs$input + if (is.null(crs)) { + warning("Unable to read CRS from mesh, using default WGS84.") + crs <- "+proj=longlat +datum=WGS84" + } + plot_choices <- c("Range", "Stdev", "AR(1)", "Boxplot", "Density", "DIC") parsed_names <- names(parsed_model_output) map_vars <- parsed_names[!parsed_names %in% c("dic", "pars")] - + ui <- shiny::fluidPage( busy_spinner, shiny::headerPanel(title = "Model viewer"), @@ -44,8 +50,8 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { ), shiny::tabPanel( "Map", - shiny::selectInput(inputId = "map_var_a", label = "Variable a:", choices = map_vars), - shiny::selectInput(inputId = "map_var_b", label = "Variable b:", choices = map_vars), + shiny::selectInput(inputId = "map_var_a", label = "Variable a:", choices = map_vars, selected = "mean_post"), + shiny::selectInput(inputId = "map_var_b", label = "Variable b:", choices = map_vars, selected = "fixed_mean"), leaflet::leafletOutput(outputId = "map_out") ), shiny::tabPanel( @@ -66,9 +72,8 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { mesh = mesh ) - crs <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" - raster = create_raster(dataframe = pred_field, crs = sp::CRS(crs)) - browser() + raster <- create_raster(dataframe = pred_field, crs = crs) + raster }) output$map_out <- leaflet::renderLeaflet({ @@ -76,7 +81,41 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { leaflet::addTiles(group = "OSM") %>% leaflet::addRasterImage(map_raster(), opacity = 0.9, group = "Raster") }) - # Run the application + + model_plot <- shiny::eventReactive(input$plot_type, ignoreNULL = FALSE, { + if (input$plot_type == "Range") { + return(plot_line_comparison( + data = parsed_model_output, + to_plot = "Range for f", + title = "Range" + )) + } else if (input$plot_type == "Stdev") { + return(plot_line_comparison( + data = parsed_model_output, + to_plot = "Stdev for f", + title = "Marginal standard deviation" + )) + } else if (input$plot_type == "AR(1)") { + return(plot_line_comparison( + data = parsed_model_output, + to_plot = "GroupRho for f", + title = "AR(1)" + )) + } else if (input$plot_type == "Boxplot") { + return(plot_priors_boxplot(data = parsed_model_output)) + } else if (input$plot_type == "Density") { + return(plot_priors_density( + data = parsed_model_output, + measurement_data = measurement_data + )) + } else if (input$plot_type == "DIC") { + return(plot_dic(data = parsed_model_output)) + } + }) + + output$plot_model_out <- shiny::renderPlot({ + model_plot() + }) } shiny::shinyApp(ui = ui, server = server) diff --git a/vignettes/data_preprocessing.Rmd b/vignettes/data_preprocessing.Rmd index 33bfd275..70d9f8c4 100644 --- a/vignettes/data_preprocessing.Rmd +++ b/vignettes/data_preprocessing.Rmd @@ -41,7 +41,7 @@ class(sp_data) Then we retrieve the projection attributes of the shapefile, i.e., `sp_data`, and transform it from its original coordinate reference system (CRS) to a new CRS, which is the World Geodetic System 1984 (WGS84). Finally, we convert it to a `SpatialPolygonsDataFrame` using [`sf::as_Spatial`](https://r-spatial.github.io/sf/reference/coerce-methods.html). ```{r transCRS} -sp_data <- sf::st_transform(sp_data, sf::st_crs("+proj=longlat +datum=WGS84 +no_defs")) +sp_data <- sf::st_transform(sp_data, sf::st_crs("+proj=longlat +datum=WGS84")) ``` In the COVID-19 tutorial, the raw COVID-19 infections data and the related covariate data were obtained from the official UK Government COVID-19 dashboard and the Office for National Statistics (ONS). The data were initially downloaded, organised and saved in a CSV file format. The CSV file can be imported into R using the `utils::read.csv()` function. diff --git a/vignettes/hydro.Rmd b/vignettes/hydro.Rmd index 0b334ccd..5bdae5e1 100644 --- a/vignettes/hydro.Rmd +++ b/vignettes/hydro.Rmd @@ -30,6 +30,17 @@ First we will look at the location of the dam by making a map and plotting the s norway_polygon_path <- fdmr::get_tutorial_datapath(dataset = "hydro", filename = "Kvilldal_Catch_Boundary.geojson") #norway_polygon <- sf::read_sf(norway_polygon_path) %>% sf::st_zm() norway_polygon <- rgdal::readOGR(norway_polygon_path) +sfc <- sf::st_transform(norway_polygon, crs = "+proj=longlat +datum=WGS84") +``` + +```{r} +norway_polygon_path <- fdmr::get_tutorial_datapath(dataset = "hydro", filename = "Kvilldal_Catch_Boundary.geojson") +norway_polygon <- rgdal::readOGR(norway_polygon_path) +norway_polygon <- sf::st_as_sf(norway_polygon, + coords = c("longitude", "latitude"), + crs = "+proj=utm +zone=32" +) + sfc <- sf::st_transform(norway_polygon, crs = "+proj=longlat +datum=WGS84") ``` We want to mark the position of the dam and the stream gauges on the map so we'll create a `data.frame` to pass into the `plot_map` @@ -193,7 +204,7 @@ get_stream_data <- function(row) { pixel_dist_gauge$streamflow <- apply(pixel_dist_gauge, 1, get_stream_data) ``` -Now we create a new `SpatialPointDataFrame` called `inla_data` for model fitting. The data frame contains the response and predictor data at each spatial location and time point. +Now we create a new `SpatialPointsDataFrame` called `inla_data` for model fitting. The `data.frame` contains the response and predictor data at each spatial location and time point. ```{r error=TRUE,inladat} inla_data <- merge(pixel_dist_gauge, raster_df, by = c("time", "x", "y")) @@ -218,7 +229,6 @@ In order to model this with the 4D-Modeller we need to: Create the triangulated mesh of the study region and have a quick look at it. ```{r error=TRUE,warning=FALSE, fig.width=7, fig.height=7} -# TODO: it would be nice to map the mesh onto the dataset using the mapmaker function i am suggesting e <- era5_precip_cropped@extent resolution <- raster::res(era5_precip_cropped) sr <- "+proj=utm +zone=32" @@ -227,12 +237,9 @@ xres <- resolution[1] * 2 yres <- resolution[2] * 2 xy <- sp::coordinates(era5_precip_cropped) - -# Maybe this makes the mesh smaller xy <- xy[seq(1, nrow(xy), by = 2), ] + colnames(xy) <- c("LONG", "LAT") -# mesh <- INLA::inla.mesh.2d(loc=xy, max.edge=c(xres*1000, xres*10000), cutoff=2000, crs=sr) -# mesh <- INLA::inla.mesh.2d(loc=xy, max.edge=c(xres*1, xres*100), cutoff=2000, crs=sr) mesh <- INLA::inla.mesh.2d(loc = xy, max.edge = c(xres * 1, xres * 100), cutoff = 75, crs = sr) @@ -294,7 +301,7 @@ formula2 <- streamflow ~ 0 + Intercept + precip + ```{r error=TRUE, modelrun} -model_out_a <- inlabru::bru(formula1, +inlabru_model1 <- inlabru::bru(formula1, data = inla_data, family = "gaussian", options = list( @@ -306,19 +313,28 @@ inlabru_model2 <- inlabru::bru(formula2, data = inla_data, family = "gaussian", options = list( - verbose = TRUE + verbose = FALSE ) ) ``` +As `create_prediction_field` expects us to have coordinates in degrees of latitude and longitude we need to transform `inla_data` the `SpatialPointsDataFrame` we pass to INLA so that we can easily plot it. We first take the extent of the map + +```{r} +parsed_output = parse_model_output(inlabru_model2, inla_data) +``` + +```{r} +pred_field <- create_prediction_field(parsed_model_output$mean_post, parsed_model_output$fixed_mean, mesh) +``` + Now it's time to see its output, we can use the `fdmr::model_viewer` and pass in our model output, measurement data and mesh. This Shiny app helps analyse model output by automatically creating plots and a map of predictions for your data. ```{r eval=FALSE} -fdmr::model_viewer(model_output = model_out_a, mesh = mesh, measurement_data = inla_data) +fdmr::model_viewer(model_output = inlabru_model2, mesh = mesh, measurement_data = inla_data) ``` ```{r error=TRUE} -# TODO : this could be where the model_eval function I recommend could be used summary(inlabru_model1) ``` ```{r error=TRUE} From e210314bb115ae60e2e5210045a24d8842908b3a Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Fri, 6 Oct 2023 14:04:55 +0100 Subject: [PATCH 08/21] Plots working --- R/model_parse.R | 1 - R/plot_priors.R | 5 +---- R/shiny_modelviewer.R | 17 ++++++++++------- vignettes/hydro.Rmd | 1 + 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/R/model_parse.R b/R/model_parse.R index 6845a8fd..64336691 100644 --- a/R/model_parse.R +++ b/R/model_parse.R @@ -68,7 +68,6 @@ create_prediction_field <- function(var_a, var_b, mesh, crs = NULL) { mod_proj <- fmesher::fm_evaluator(mesh, crs = crs) xy_grid <- base::expand.grid(mod_proj$x, mod_proj$y) A_proj <- INLA::inla.spde.make.A(mesh = mesh, loc = as.matrix(xy_grid)) - z <- base::exp(base::as.numeric(A_proj %*% var_a[1:mesh$n]) + base::sum(var_b)) base::data.frame(x = xy_grid[, 1], y = xy_grid[, 2], z = z) } diff --git a/R/plot_priors.R b/R/plot_priors.R index 40384e87..40a34772 100644 --- a/R/plot_priors.R +++ b/R/plot_priors.R @@ -12,13 +12,10 @@ plot_line_comparison <- function(data, to_plot, title) { return("No pars data.") } - plot <- ggplot2::ggplot(single_df, ggplot2::aes(x = x, y = y, color = Run)) + + ggplot2::ggplot(single_df, ggplot2::aes(x = x, y = y, color = Run)) + ggplot2::geom_line() + ggplot2::ggtitle(title) + ggplot2::theme(text = ggplot2::element_text(size = 16)) - - browser() - plot } diff --git a/R/shiny_modelviewer.R b/R/shiny_modelviewer.R index da1f2f49..fb2d4696 100644 --- a/R/shiny_modelviewer.R +++ b/R/shiny_modelviewer.R @@ -16,6 +16,9 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { measurement_data = measurement_data ) + # The comparison functions expect a list of lists + parsed_modeloutput_plots <- list(parsed_model_output) + crs <- mesh$crs$input if (is.null(crs)) { warning("Unable to read CRS from mesh, using default WGS84.") @@ -33,7 +36,7 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { shiny::sidebarLayout( shiny::sidebarPanel( shiny::selectInput( - inputId = "plot_type", + inputId = "model_type", label = "Model type:", choices = c("inlabru"), selected = "inlabru" @@ -85,31 +88,31 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { model_plot <- shiny::eventReactive(input$plot_type, ignoreNULL = FALSE, { if (input$plot_type == "Range") { return(plot_line_comparison( - data = parsed_model_output, + data = parsed_modeloutput_plots, to_plot = "Range for f", title = "Range" )) } else if (input$plot_type == "Stdev") { return(plot_line_comparison( - data = parsed_model_output, + data = parsed_modeloutput_plots, to_plot = "Stdev for f", title = "Marginal standard deviation" )) } else if (input$plot_type == "AR(1)") { return(plot_line_comparison( - data = parsed_model_output, + data = parsed_modeloutput_plots, to_plot = "GroupRho for f", title = "AR(1)" )) } else if (input$plot_type == "Boxplot") { - return(plot_priors_boxplot(data = parsed_model_output)) + return(plot_priors_boxplot(data = parsed_modeloutput_plots)) } else if (input$plot_type == "Density") { return(plot_priors_density( - data = parsed_model_output, + data = parsed_modeloutput_plots, measurement_data = measurement_data )) } else if (input$plot_type == "DIC") { - return(plot_dic(data = parsed_model_output)) + return(plot_dic(data = parsed_modeloutput_plots)) } }) diff --git a/vignettes/hydro.Rmd b/vignettes/hydro.Rmd index 5bdae5e1..f3e83900 100644 --- a/vignettes/hydro.Rmd +++ b/vignettes/hydro.Rmd @@ -325,6 +325,7 @@ parsed_output = parse_model_output(inlabru_model2, inla_data) ```{r} pred_field <- create_prediction_field(parsed_model_output$mean_post, parsed_model_output$fixed_mean, mesh) +pred_field ``` From 3020700848b49457265ee97e96c170652dff3fb4 Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Fri, 6 Oct 2023 15:45:20 +0100 Subject: [PATCH 09/21] Adding man files --- man/create_prediction_field.Rd | 23 +++++++++++++++++++++++ man/get_busy_spinner.Rd | 19 +++++++++++++++++++ man/has_coords.Rd | 20 ++++++++++++++++++++ man/model_viewer.Rd | 21 +++++++++++++++++++++ man/model_viewer_shiny.Rd | 22 ++++++++++++++++++++++ 5 files changed, 105 insertions(+) create mode 100644 man/create_prediction_field.Rd create mode 100644 man/get_busy_spinner.Rd create mode 100644 man/has_coords.Rd create mode 100644 man/model_viewer.Rd create mode 100644 man/model_viewer_shiny.Rd diff --git a/man/create_prediction_field.Rd b/man/create_prediction_field.Rd new file mode 100644 index 00000000..709de1b8 --- /dev/null +++ b/man/create_prediction_field.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/model_parse.R +\name{create_prediction_field} +\alias{create_prediction_field} +\title{Create a prediction field from the parsed model output and the mesh} +\usage{ +create_prediction_field(var_a, var_b, mesh, crs = NULL) +} +\arguments{ +\item{var_a}{Variable a data} + +\item{var_b}{Variable b data} + +\item{mesh}{INLA mesh} + +\item{crs}{CRS as a proj4string} +} +\value{ +data.frame +} +\description{ +Create a prediction field from the parsed model output and the mesh +} diff --git a/man/get_busy_spinner.Rd b/man/get_busy_spinner.Rd new file mode 100644 index 00000000..919a9d9f --- /dev/null +++ b/man/get_busy_spinner.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/util_shiny.R +\name{get_busy_spinner} +\alias{get_busy_spinner} +\title{Get a busy spinner for our Shiny apps. +If we have an internet connection we download the spinner +TODO - add this gif to the package} +\usage{ +get_busy_spinner() +} +\value{ +busyspinner +} +\description{ +Get a busy spinner for our Shiny apps. +If we have an internet connection we download the spinner +TODO - add this gif to the package +} +\keyword{internal} diff --git a/man/has_coords.Rd b/man/has_coords.Rd new file mode 100644 index 00000000..cd5dc91e --- /dev/null +++ b/man/has_coords.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/util_coords.R +\name{has_coords} +\alias{has_coords} +\title{Check if spatial data has coordinates set using sp::coordinates() +Returns TRUE if coordinates have been set} +\usage{ +has_coords(spatial_data) +} +\arguments{ +\item{spatial_data}{} +} +\value{ +bool +} +\description{ +Check if spatial data has coordinates set using sp::coordinates() +Returns TRUE if coordinates have been set +} +\keyword{internal} diff --git a/man/model_viewer.Rd b/man/model_viewer.Rd new file mode 100644 index 00000000..5947add3 --- /dev/null +++ b/man/model_viewer.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shiny_modelviewer.R +\name{model_viewer} +\alias{model_viewer} +\title{Mesh building shiny app. Creates and visualises a mesh from some spatial data.} +\usage{ +model_viewer(model_output, mesh, measurement_data) +} +\arguments{ +\item{model_output}{INLA model output} + +\item{mesh}{INLA mesh} + +\item{measurement_data}{Measurement data} +} +\value{ +shiny::app +} +\description{ +Mesh building shiny app. Creates and visualises a mesh from some spatial data. +} diff --git a/man/model_viewer_shiny.Rd b/man/model_viewer_shiny.Rd new file mode 100644 index 00000000..7bd6353b --- /dev/null +++ b/man/model_viewer_shiny.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shiny_modelviewer.R +\name{model_viewer_shiny} +\alias{model_viewer_shiny} +\title{Parse inlabru model output} +\usage{ +model_viewer_shiny(model_output, mesh, measurement_data) +} +\arguments{ +\item{model_output}{INLA model output} + +\item{mesh}{INLA mesh} + +\item{measurement_data}{Measurement data} +} +\value{ +shiny::app +} +\description{ +Parse inlabru model output +} +\keyword{internal} From 3a6fcc5221b40a583ff2cdd0ee3fad5fecd97bbe Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Fri, 6 Oct 2023 15:52:18 +0100 Subject: [PATCH 10/21] Update namespace --- NAMESPACE | 2 +- man/create_prediction_field.Rd | 18 +++++++++++++----- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index c38bd3f8..f733d69d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,7 +2,6 @@ export(clear_caches) export(create_prediction_field) -export(get_tmpdir) export(get_tutorial_datapath) export(interactive_priors) export(latlong_to_utm) @@ -10,6 +9,7 @@ export(load_tutorial_data) export(mesh_builder) export(mesh_checker) export(mesh_to_spatial) +export(model_viewer) export(numbers_only) export(parse_model_output) export(plot_barchart) diff --git a/man/create_prediction_field.Rd b/man/create_prediction_field.Rd index 709de1b8..5ad702f1 100644 --- a/man/create_prediction_field.Rd +++ b/man/create_prediction_field.Rd @@ -4,16 +4,24 @@ \alias{create_prediction_field} \title{Create a prediction field from the parsed model output and the mesh} \usage{ -create_prediction_field(var_a, var_b, mesh, crs = NULL) +create_prediction_field( + mesh, + plot_type = "predicted_mean_fields", + data_type = "poisson", + var_a = NULL, + var_b = NULL +) } \arguments{ -\item{var_a}{Variable a data} +\item{mesh}{INLA mesh} -\item{var_b}{Variable b data} +\item{plot_type}{Type of plot to create, "predicted_mean_fields" etc} -\item{mesh}{INLA mesh} +\item{data_type}{Type of data, "poisson" etc} + +\item{var_a}{Data for variable a, required for "predicted_mean_fields" and "random_effect_fields"} -\item{crs}{CRS as a proj4string} +\item{var_b}{Data for variable b, required for "predicted_mean_fields"} } \value{ data.frame From 8dfc79c44489d06a783edaecd16663dbc07a0662 Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Sat, 7 Oct 2023 10:22:05 +0100 Subject: [PATCH 11/21] Adding legend and correct colours to raster --- R/model_parse.R | 7 +---- R/shiny_modelviewer.R | 65 ++++++++++++++++++++++++------------------- R/shiny_priors.R | 35 +++++++++++++++++------ vignettes/hydro.Rmd | 22 +++++++-------- 4 files changed, 76 insertions(+), 53 deletions(-) diff --git a/R/model_parse.R b/R/model_parse.R index f1104c69..9bc99b42 100644 --- a/R/model_parse.R +++ b/R/model_parse.R @@ -92,9 +92,4 @@ create_prediction_field <- function(mesh, } base::data.frame(x = xy_grid[, 1], y = xy_grid[, 2], z = z) -} - - -create_raster <- function(dataframe, crs) { - raster::rasterFromXYZ(dataframe, crs = crs) -} +} \ No newline at end of file diff --git a/R/shiny_modelviewer.R b/R/shiny_modelviewer.R index 1c5f0748..cc3e68cb 100644 --- a/R/shiny_modelviewer.R +++ b/R/shiny_modelviewer.R @@ -20,8 +20,8 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { parsed_modeloutput_plots <- list(parsed_model_output) crs <- mesh$crs$input - if (is.null(crs)) { - warning("Unable to read CRS from mesh, using default WGS84.") + if (is.null(crs) || is.na(crs)) { + warning("Cannot read CRS from mesh, using default CRS = +proj=longlat +datum=WGS84") crs <- "+proj=longlat +datum=WGS84" } @@ -30,25 +30,23 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { ui <- shiny::fluidPage( busy_spinner, shiny::headerPanel(title = "Model viewer"), - shiny::sidebarLayout( - shiny::tabsetPanel( - type = "tabs", - shiny::tabPanel( - "Plots", - shiny::h2("Plot output"), - shiny::selectInput(inputId = "plot_type", label = "Plot type:", choices = plot_choices, selected = plot_choices[1]), - shiny::plotOutput(outputId = "plot_model_out") - ), - shiny::tabPanel( - "Map", - shiny::selectInput(inputId = "map_var_a", label = "Variable a:", choices = map_vars, selected = "mean_post"), - shiny::selectInput(inputId = "map_var_b", label = "Variable b:", choices = map_vars, selected = "fixed_mean"), - leaflet::leafletOutput(outputId = "map_out") - ), - shiny::tabPanel( - "Help", - shiny::h3("Help"), - ) + shiny::tabsetPanel( + type = "tabs", + shiny::tabPanel( + "Plots", + shiny::h2("Plot output"), + shiny::selectInput(inputId = "plot_type", label = "Plot type:", choices = plot_choices, selected = plot_choices[1]), + shiny::plotOutput(outputId = "plot_model_out") + ), + shiny::tabPanel( + "Map", + shiny::selectInput(inputId = "map_plot_type", label = "Plot type", choices = c("Predicted mean fields", "Random effect fields"), selected = "Predicted mean fields"), + shiny::selectInput(inputId = "map_data_type", label = "Data type", choices = c("Poisson", "Gaussian"), selected = "Poisson"), + leaflet::leafletOutput(outputId = "map_out") + ), + shiny::tabPanel( + "Help", + shiny::h3("Help"), ) ) ) @@ -56,14 +54,25 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { # Define server logic required to draw a histogram server <- function(input, output, session) { map_raster <- shiny::reactive({ - pred_field <- create_prediction_field( - var_a = parsed_model_output[[input$map_var_a]], - var_b = parsed_model_output[[input$map_var_b]], - mesh = mesh - ) + data_type <- tolower(input$map_data_type) + if (input$map_plot_type == "Predicted mean fields") { + pred_field <- create_prediction_field( + mesh = mesh, + plot_type = "predicted_mean_fields", + data_type = data_type, + var_a = parsed_model_output[["mean_post"]], + var_b = parsed_model_output[["fixed_mean"]] + ) + } else { + pred_field <- create_prediction_field( + mesh = mesh, + plot_type = "random_effect_fields", + data_type = data_type, + var_a = parsed_model_output[["mean_post"]] + ) + } - raster <- create_raster(dataframe = pred_field, crs = crs) - raster + raster::rasterFromXYZ(pred_field, crs = crs) }) output$map_out <- leaflet::renderLeaflet({ diff --git a/R/shiny_priors.R b/R/shiny_priors.R index 1e7e0e51..00cfc419 100644 --- a/R/shiny_priors.R +++ b/R/shiny_priors.R @@ -22,6 +22,18 @@ priors_shiny <- function(spatial_data, stop("Please make sure you have set coordinates on spatial_data using sp::coordinates.") } + spatial_crs <- sp::proj4string(spatial_data) + mesh_crs <- mesh$crs$input + + if (is.na(mesh_crs) && is.na(spatial_crs)) { + warning("Cannot read CRS from mesh or spatial_data, using default CRS = +proj=longlat +datum=WGS84") + crs <- "+proj=longlat +datum=WGS84" + } else if (is.na(mesh_crs)) { + crs <- spatial_crs + } else { + crs <- mesh_crs + } + # Text for priors help prior_range_text <- "A length 2 vector, with (range0, Prange) specifying that P(ρ < ρ_0)=p_ρ, where ρ is the spatial range of the random field." @@ -368,7 +380,7 @@ priors_shiny <- function(spatial_data, rownames = TRUE ) - map_raster <- shiny::reactive({ + prediction_field <- shiny::reactive({ if (length(model_vals$parsed_outputs) == 0) { return() } @@ -377,7 +389,7 @@ priors_shiny <- function(spatial_data, # c("Predicted mean fields", "Random effect fields") data_type <- tolower(input$map_data_type) if (input$map_plot_type == "Predicted mean fields") { - pred_field <- create_prediction_field( + create_prediction_field( mesh = mesh, plot_type = "predicted_mean_fields", data_type = data_type, @@ -385,15 +397,22 @@ priors_shiny <- function(spatial_data, var_b = data[["fixed_mean"]] ) } else { - pred_field <- create_prediction_field( + create_prediction_field( mesh = mesh, plot_type = "random_effect_fields", data_type = data_type, var_a = data[["mean_post"]] ) } + }) + + map_raster <- shiny::reactive({ + raster::rasterFromXYZ(prediction_field(), crs = crs) + }) - create_raster(dataframe = pred_field, crs = sp::proj4string(spatial_data)) + map_colours <- shiny::reactive({ + range <- prediction_field()[["z"]] + leaflet::colorNumeric(palette = "viridis", domain = range, reverse = FALSE) }) output$map_out <- leaflet::renderLeaflet({ @@ -401,10 +420,10 @@ priors_shiny <- function(spatial_data, return() } - m <- leaflet::leaflet() - m <- leaflet::addTiles(m, group = "OSM") - m <- leaflet::addRasterImage(m, map_raster(), opacity = 0.9, group = "Raster") - m + leaflet::leaflet() %>% + leaflet::addTiles(group = "OSM") %>% + leaflet::addRasterImage(map_raster(), colors=map_colours(), opacity = 0.9, group = "Raster") %>% + leaflet::addLegend(position = "topright", pal = map_colours(), values = prediction_field()[["z"]]) }) model_plot <- shiny::eventReactive(input$plot_type, ignoreNULL = FALSE, { diff --git a/vignettes/hydro.Rmd b/vignettes/hydro.Rmd index f3e83900..159e83f3 100644 --- a/vignettes/hydro.Rmd +++ b/vignettes/hydro.Rmd @@ -33,16 +33,16 @@ norway_polygon <- rgdal::readOGR(norway_polygon_path) sfc <- sf::st_transform(norway_polygon, crs = "+proj=longlat +datum=WGS84") ``` -```{r} -norway_polygon_path <- fdmr::get_tutorial_datapath(dataset = "hydro", filename = "Kvilldal_Catch_Boundary.geojson") -norway_polygon <- rgdal::readOGR(norway_polygon_path) -norway_polygon <- sf::st_as_sf(norway_polygon, - coords = c("longitude", "latitude"), - crs = "+proj=utm +zone=32" -) - -sfc <- sf::st_transform(norway_polygon, crs = "+proj=longlat +datum=WGS84") -``` +# ```{r} +# norway_polygon_path <- fdmr::get_tutorial_datapath(dataset = "hydro", filename = "Kvilldal_Catch_Boundary.geojson") +# norway_polygon <- rgdal::readOGR(norway_polygon_path) +# norway_polygon <- sf::st_as_sf(norway_polygon, +# coords = c("longitude", "latitude"), +# crs = "+proj=utm +zone=32" +# ) + +# sfc <- sf::st_transform(norway_polygon, crs = "+proj=longlat +datum=WGS84") +# ``` We want to mark the position of the dam and the stream gauges on the map so we'll create a `data.frame` to pass into the `plot_map` function with the names and positions. We'll create a list for each of the points and then use the (`dplyr::bind_rows`)[https://dplyr.tidyverse.org/reference/bind.html] function to create a `data.frame` from these. @@ -324,7 +324,7 @@ parsed_output = parse_model_output(inlabru_model2, inla_data) ``` ```{r} -pred_field <- create_prediction_field(parsed_model_output$mean_post, parsed_model_output$fixed_mean, mesh) +pred_field <- create_prediction_field(mesh, plot_type = "predicted_mean_fields", data_type = "poisson", parsed_model_output$mean_post, parsed_model_output$fixed_mean) pred_field ``` From dab636eee2ba6114c64b491efb1c91816c3f2522 Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Sat, 7 Oct 2023 18:24:45 +0100 Subject: [PATCH 12/21] Basic model viz working --- R/model_parse.R | 10 ++--- R/shiny_mapper.R | 11 ++--- R/shiny_modelviewer.R | 94 ++++++++++++++++++++++++++++++++++++------- R/shiny_priors.R | 39 +++++++++++++----- vignettes/hydro.Rmd | 15 +++++-- 5 files changed, 128 insertions(+), 41 deletions(-) diff --git a/R/model_parse.R b/R/model_parse.R index 9bc99b42..0af9d2e7 100644 --- a/R/model_parse.R +++ b/R/model_parse.R @@ -66,6 +66,7 @@ create_prediction_field <- function(mesh, stop("Invalid plot type, select from ", valid_plots) } + data_type <- tolower(data_type) valid_data_types <- c("poisson", "gaussian") if (!(data_type %in% valid_data_types)) { stop("Invalid data type, select from ", valid_data_types) @@ -80,11 +81,8 @@ create_prediction_field <- function(mesh, A_proj <- INLA::inla.spde.make.A(mesh = mesh, loc = as.matrix(xy_grid)) if (plot_type == "predicted_mean_fields") { - if (data_type == "poisson") { - z <- base::exp(base::as.numeric(A_proj %*% var_a[1:mesh$n]) + base::sum(var_b)) - } else { - z <- base::as.numeric(A_proj %*% var_a[1:mesh$n]) + base::sum(var_b) - } + z <- base::as.numeric(A_proj %*% var_a[1:mesh$n]) + base::sum(var_b) + if (data_type == "poisson") z <- base::exp(z) } else { # We get an error here as we only have 265 items # z <- var_a[1:mesh$n] @@ -92,4 +90,4 @@ create_prediction_field <- function(mesh, } base::data.frame(x = xy_grid[, 1], y = xy_grid[, 2], z = z) -} \ No newline at end of file +} diff --git a/R/shiny_mapper.R b/R/shiny_mapper.R index 6689fa06..72141c1d 100644 --- a/R/shiny_mapper.R +++ b/R/shiny_mapper.R @@ -66,13 +66,13 @@ raster_mapping_app <- function(raster_data = NULL, polygon_data = NULL, date_for inputId = "colour_category", label = "Palette type", choices = c("Sequential", "Diverging", "Qualitative", "Viridis"), - selected = "seq" + selected = "Viridis" ), shiny::selectInput( inputId = "colour_scheme", label = "Color Scheme", choices = default_colours, - selected = "Blues" + selected = "viridis" ), shiny::sliderInput( inputId = "raster_opacity", @@ -138,6 +138,10 @@ raster_mapping_app <- function(raster_data = NULL, polygon_data = NULL, date_for colours }) + colour_scheme <- shiny::reactive({ + input$colour_scheme + }) + raster_opacity <- shiny::reactive({ input$raster_opacity }) @@ -146,9 +150,6 @@ raster_mapping_app <- function(raster_data = NULL, polygon_data = NULL, date_for input$polygon_opacity }) - colour_scheme <- shiny::reactive({ - input$colour_scheme - }) colour_palette <- shiny::reactive({ if (is.null(palette)) { diff --git a/R/shiny_modelviewer.R b/R/shiny_modelviewer.R index cc3e68cb..2ccef2fd 100644 --- a/R/shiny_modelviewer.R +++ b/R/shiny_modelviewer.R @@ -3,14 +3,26 @@ #' @param model_output INLA model output #' @param mesh INLA mesh #' @param measurement_data Measurement data +#' @param data_type Type of data, Poisson or Gaussian #' #' @importFrom magrittr %>% #' #' @return shiny::app #' @keywords internal -model_viewer_shiny <- function(model_output, mesh, measurement_data) { +model_viewer_shiny <- function(model_output, mesh, measurement_data, data_type) { busy_spinner <- get_busy_spinner() + crs <- mesh$crs$input + if (is.null(crs) || is.na(crs)) { + warning("Cannot read CRS from mesh, using default CRS = +proj=longlat +datum=WGS84") + crs <- "+proj=longlat +datum=WGS84" + } + + data_type <- stringr::str_to_title(data_type) + if (!(data_type %in% c("Poisson", "Gaussian"))) { + stop("data_type must be one of Poisson or Gaussian") + } + parsed_model_output <- parse_model_output( model_output = model_output, measurement_data = measurement_data @@ -19,11 +31,8 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { # The comparison plotting functions expect a list of lists parsed_modeloutput_plots <- list(parsed_model_output) - crs <- mesh$crs$input - if (is.null(crs) || is.na(crs)) { - warning("Cannot read CRS from mesh, using default CRS = +proj=longlat +datum=WGS84") - crs <- "+proj=longlat +datum=WGS84" - } + brewer_palettes <- RColorBrewer::brewer.pal.info + default_colours <- rownames(brewer_palettes[brewer_palettes$cat == "seq", ]) plot_choices <- c("Range", "Stdev", "AR(1)", "Boxplot", "Density", "DIC") @@ -40,8 +49,27 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { ), shiny::tabPanel( "Map", - shiny::selectInput(inputId = "map_plot_type", label = "Plot type", choices = c("Predicted mean fields", "Random effect fields"), selected = "Predicted mean fields"), - shiny::selectInput(inputId = "map_data_type", label = "Data type", choices = c("Poisson", "Gaussian"), selected = "Poisson"), + shiny::fluidRow( + shiny::column( + 6, + shiny::selectInput(inputId = "map_plot_type", label = "Plot type", choices = c("Predicted mean fields", "Random effect fields"), selected = "Predicted mean fields"), + shiny::selectInput(inputId = "map_data_type", label = "Data type", choices = c("Poisson", "Gaussian"), selected = data_type), + ), + shiny::column( + 6, + shiny::selectInput( + inputId = "colour_category", + label = "Palette type", + choices = c("Sequential", "Diverging", "Qualitative", "Viridis"), + selected = "Viridis" + ), + shiny::selectInput( + inputId = "colour_scheme", + label = "Color Scheme", + choices = default_colours, + ), + ) + ), leaflet::leafletOutput(outputId = "map_out") ), shiny::tabPanel( @@ -53,10 +81,30 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { # Define server logic required to draw a histogram server <- function(input, output, session) { - map_raster <- shiny::reactive({ + category_colours <- shiny::reactive({ + if (input$colour_category == "Viridis") { + colours <- c("viridis", "magma", "inferno", "plasma") + } else { + palettes_mapping <- list("Sequential" = "seq", "Diverging" = "div", "Qualitative" = "qual") + chosen_cat <- palettes_mapping[input$colour_category] + colours <- rownames(subset(RColorBrewer::brewer.pal.info, category %in% chosen_cat)) + } + colours + }) + + colour_scheme <- shiny::reactive({ + input$colour_scheme + }) + + shiny::observe({ + shiny::updateSelectInput(session, inputId = "colour_scheme", label = "Colours", choices = category_colours()) + }) + + + prediction_field <- shiny::reactive({ data_type <- tolower(input$map_data_type) if (input$map_plot_type == "Predicted mean fields") { - pred_field <- create_prediction_field( + create_prediction_field( mesh = mesh, plot_type = "predicted_mean_fields", data_type = data_type, @@ -64,21 +112,36 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { var_b = parsed_model_output[["fixed_mean"]] ) } else { - pred_field <- create_prediction_field( + create_prediction_field( mesh = mesh, plot_type = "random_effect_fields", data_type = data_type, var_a = parsed_model_output[["mean_post"]] ) } + }) - raster::rasterFromXYZ(pred_field, crs = crs) + z_values <- shiny::reactive({ + prediction_field()[["z"]] + }) + + map_raster <- shiny::reactive({ + raster::rasterFromXYZ(prediction_field(), crs = crs) + }) + + map_colours <- shiny::reactive({ + leaflet::colorNumeric(palette = colour_scheme(), domain = z_values(), reverse = FALSE) }) output$map_out <- leaflet::renderLeaflet({ + if (is.null(map_raster())) { + return() + } + leaflet::leaflet() %>% leaflet::addTiles(group = "OSM") %>% - leaflet::addRasterImage(map_raster(), opacity = 0.9, group = "Raster") + leaflet::addRasterImage(map_raster(), colors = map_colours(), opacity = 0.9, group = "Raster") %>% + leaflet::addLegend(position = "topright", pal = map_colours(), values = z_values()) }) model_plot <- shiny::eventReactive(input$plot_type, ignoreNULL = FALSE, { @@ -125,9 +188,10 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data) { #' @param model_output INLA model output #' @param mesh INLA mesh #' @param measurement_data Measurement data +#' @param data_type Type of data, Poisson or Gaussian #' #' @return shiny::app #' @export -model_viewer <- function(model_output, mesh, measurement_data) { - shiny::runApp(model_viewer_shiny(model_output = model_output, mesh = mesh, measurement_data = measurement_data)) +model_viewer <- function(model_output, mesh, measurement_data, data_type = "Poisson") { + shiny::runApp(model_viewer_shiny(model_output = model_output, mesh = mesh, measurement_data = measurement_data, data_type = data_type)) } diff --git a/R/shiny_priors.R b/R/shiny_priors.R index 00cfc419..494d4792 100644 --- a/R/shiny_priors.R +++ b/R/shiny_priors.R @@ -380,19 +380,33 @@ priors_shiny <- function(spatial_data, rownames = TRUE ) - prediction_field <- shiny::reactive({ - if (length(model_vals$parsed_outputs) == 0) { - return() + category_colours <- shiny::reactive({ + if (input$colour_category == "Viridis") { + colours <- c("viridis", "magma", "inferno", "plasma") + } else { + palettes_mapping <- list("Sequential" = "seq", "Diverging" = "div", "Qualitative" = "qual") + chosen_cat <- palettes_mapping[input$colour_category] + colours <- rownames(subset(RColorBrewer::brewer.pal.info, category %in% chosen_cat)) } + colours + }) + + colour_scheme <- shiny::reactive({ + input$colour_scheme + }) + shiny::observe({ + shiny::updateSelectInput(session, inputId = "colour_scheme", label = "Colours", choices = category_colours()) + }) + + + prediction_field <- shiny::reactive({ data <- model_vals$parsed_outputs[[input$select_run_map]] - # c("Predicted mean fields", "Random effect fields") - data_type <- tolower(input$map_data_type) if (input$map_plot_type == "Predicted mean fields") { create_prediction_field( mesh = mesh, plot_type = "predicted_mean_fields", - data_type = data_type, + data_type = input$data_type, var_a = data[["mean_post"]], var_b = data[["fixed_mean"]] ) @@ -400,19 +414,22 @@ priors_shiny <- function(spatial_data, create_prediction_field( mesh = mesh, plot_type = "random_effect_fields", - data_type = data_type, + data_type = input$data_type, var_a = data[["mean_post"]] ) } }) + z_values <- shiny::reactive({ + prediction_field()[["z"]] + }) + map_raster <- shiny::reactive({ raster::rasterFromXYZ(prediction_field(), crs = crs) }) map_colours <- shiny::reactive({ - range <- prediction_field()[["z"]] - leaflet::colorNumeric(palette = "viridis", domain = range, reverse = FALSE) + leaflet::colorNumeric(palette = colour_scheme(), domain = z_values(), reverse = FALSE) }) output$map_out <- leaflet::renderLeaflet({ @@ -422,8 +439,8 @@ priors_shiny <- function(spatial_data, leaflet::leaflet() %>% leaflet::addTiles(group = "OSM") %>% - leaflet::addRasterImage(map_raster(), colors=map_colours(), opacity = 0.9, group = "Raster") %>% - leaflet::addLegend(position = "topright", pal = map_colours(), values = prediction_field()[["z"]]) + leaflet::addRasterImage(map_raster(), colors = map_colours(), opacity = 0.9, group = "Raster") %>% + leaflet::addLegend(position = "topright", pal = map_colours(), values = z_values()) }) model_plot <- shiny::eventReactive(input$plot_type, ignoreNULL = FALSE, { diff --git a/vignettes/hydro.Rmd b/vignettes/hydro.Rmd index 159e83f3..73619f09 100644 --- a/vignettes/hydro.Rmd +++ b/vignettes/hydro.Rmd @@ -27,9 +27,16 @@ First we will look at the location of the dam by making a map and plotting the s ```{r error=TRUE,warning=FALSE} # TODO : here is a good example for the map function -norway_polygon_path <- fdmr::get_tutorial_datapath(dataset = "hydro", filename = "Kvilldal_Catch_Boundary.geojson") +#norway_polygon_path <- fdmr::get_tutorial_datapath(dataset = "hydro", filename = "Kvilldal_Catch_Boundary.geojson") #norway_polygon <- sf::read_sf(norway_polygon_path) %>% sf::st_zm() +norway_polygon_path <- fdmr::get_tutorial_datapath(dataset = "hydro", filename = "Kvilldal_Catch_Boundary.geojson") norway_polygon <- rgdal::readOGR(norway_polygon_path) +norway_polygon <- sf::st_as_sf(norway_polygon, + coords = c("longitude", "latitude"), + crs = "+proj=utm +zone=32" +) + + sfc <- sf::st_transform(norway_polygon, crs = "+proj=longlat +datum=WGS84") ``` @@ -324,15 +331,15 @@ parsed_output = parse_model_output(inlabru_model2, inla_data) ``` ```{r} -pred_field <- create_prediction_field(mesh, plot_type = "predicted_mean_fields", data_type = "poisson", parsed_model_output$mean_post, parsed_model_output$fixed_mean) +pred_field <- create_prediction_field(mesh, plot_type = "predicted_mean_fields", data_type = "gaussian", parsed_output$mean_post, parsed_output$fixed_mean) pred_field ``` -Now it's time to see its output, we can use the `fdmr::model_viewer` and pass in our model output, measurement data and mesh. This Shiny app helps analyse model output by automatically creating plots and a map of predictions for your data. +Now it's time to see its output, we can use the `fdmr::model_viewer` and pass in our model output, mesh and measurement data. This Shiny app helps analyse model output by automatically creating plots and a map of predictions for your data. Make sure you have a CRS set on your mesh. If not CRS can be read a default of `"+proj=longlat +datum=WGS84"` will be used which may lead to errors with data using UTM as we've done here. ```{r eval=FALSE} -fdmr::model_viewer(model_output = inlabru_model2, mesh = mesh, measurement_data = inla_data) +fdmr::model_viewer(model_output = inlabru_model2, mesh = mesh, measurement_data = inla_data, data_type = "Gaussian") ``` ```{r error=TRUE} From 44efcd3825a8b07d4e1b85820f09ecc17fcd5568 Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Sat, 7 Oct 2023 18:41:27 +0100 Subject: [PATCH 13/21] Adding ability to select Poisson or Gaussian data to priors app --- R/shiny_priors.R | 53 +++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 46 insertions(+), 7 deletions(-) diff --git a/R/shiny_priors.R b/R/shiny_priors.R index 494d4792..8f7e74e0 100644 --- a/R/shiny_priors.R +++ b/R/shiny_priors.R @@ -4,6 +4,7 @@ #' @param measurement_data Measurement data #' @param time_variable Time variable in measurement_data #' @param mesh INLA mesh +#' @param data_distribution Data distribution, Poisson or Gaussian #' @param log_folder Folder to write out logs #' #' @importFrom INLA f @@ -14,9 +15,16 @@ priors_shiny <- function(spatial_data, measurement_data, time_variable, mesh, + data_distribution = "Poisson", log_folder = NULL) { future::plan(future::multisession()) + if (!(data_distribution %in% c("Poisson", "Gaussian"))) { + stop("We only support Poisson and Gaussian data") + } + + data_distribution_lower <- tolower(data_distribution) + got_coords <- has_coords(spatial_data = spatial_data) if (!got_coords) { stop("Please make sure you have set coordinates on spatial_data using sp::coordinates.") @@ -34,6 +42,9 @@ priors_shiny <- function(spatial_data, crs <- mesh_crs } + brewer_palettes <- RColorBrewer::brewer.pal.info + default_colours <- rownames(brewer_palettes[brewer_palettes$cat == "seq", ]) + # Text for priors help prior_range_text <- "A length 2 vector, with (range0, Prange) specifying that P(ρ < ρ_0)=p_ρ, where ρ is the spatial range of the random field." @@ -171,9 +182,27 @@ priors_shiny <- function(spatial_data, ), shiny::tabPanel( "Map", - shiny::selectInput(inputId = "select_run_map", label = "Select run:", choices = c()), - shiny::selectInput(inputId = "map_plot_type", label = "Plot type", choices = c("Predicted mean fields", "Random effect fields"), selected = "Predicted mean fields"), - shiny::selectInput(inputId = "map_data_type", label = "Data type", choices = c("Poisson", "Gaussian"), selected = "Poisson"), + shiny::fluidRow( + shiny::column( + 6, + shiny::selectInput(inputId = "map_plot_type", label = "Plot type", choices = c("Predicted mean fields", "Random effect fields"), selected = "Predicted mean fields"), + shiny::selectInput(inputId = "map_data_type", label = "Data type", choices = c("Poisson", "Gaussian"), selected = "Poisson"), + ), + shiny::column( + 6, + shiny::selectInput( + inputId = "colour_category", + label = "Palette type", + choices = c("Sequential", "Diverging", "Qualitative", "Viridis"), + selected = "Viridis" + ), + shiny::selectInput( + inputId = "colour_scheme", + label = "Color Scheme", + choices = default_colours, + ), + ) + ), leaflet::leafletOutput(outputId = "map_out") ), shiny::tabPanel( @@ -302,15 +331,20 @@ priors_shiny <- function(spatial_data, formula_local <- inla_formula() measurement_data_local <- measurement_data + family_control <- NULL + if (data_distribution_lower == "poisson") { + family_control <- list(link = "log") + } + promise <- promises::future_promise( { # Without loading INLA here we get errors require("INLA") inlabru::bru(formula_local, data = measurement_data_local, - family = "poisson", + family = data_distribution_lower, E = measurement_data_local[[exposure_param_local]], - control.family = list(link = "log"), + control.family = family_control, options = list( verbose = FALSE ) @@ -491,6 +525,11 @@ priors_shiny <- function(spatial_data, params <- model_vals$run_params[[input$select_run_code]] + family_control_str <- "NULL" + if (data_distribution_lower == "poisson") { + family_control_str <- "control.family = list(link = 'log')," + } + paste0( "spde <- INLA::inla.spde2.pcmatern( mesh = mesh, @@ -504,9 +543,9 @@ priors_shiny <- function(spatial_data, )", "\n\n", paste0("model_output <- inlabru::bru(formula, data = measurement_data, - family = 'poisson', + family = ", data_distribution_lower, ", E = measurement_data[[", input$exposure_param, "]], - control.family = list(link = 'log'), + control.family = ", family_control_str, " options = list( verbose = FALSE ) From 52c644514f58a33f4374549a277b01fa9073411f Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Sat, 7 Oct 2023 19:02:31 +0100 Subject: [PATCH 14/21] Combine observes --- R/shiny_priors.R | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/R/shiny_priors.R b/R/shiny_priors.R index 8f7e74e0..a8a49c49 100644 --- a/R/shiny_priors.R +++ b/R/shiny_priors.R @@ -253,6 +253,7 @@ priors_shiny <- function(spatial_data, }) shiny::observe({ + shiny::updateSelectInput(session, inputId = "colour_scheme", label = "Colours", choices = category_colours()) shiny::updateSelectInput(session = session, inputId = "select_run_map", choices = run_names()) shiny::updateSelectInput(session = session, inputId = "select_run_code", choices = run_names()) }) @@ -429,11 +430,6 @@ priors_shiny <- function(spatial_data, input$colour_scheme }) - shiny::observe({ - shiny::updateSelectInput(session, inputId = "colour_scheme", label = "Colours", choices = category_colours()) - }) - - prediction_field <- shiny::reactive({ data <- model_vals$parsed_outputs[[input$select_run_map]] if (input$map_plot_type == "Predicted mean fields") { From 159ca2c7ac4beac6341f29964ff1087b67bfe5cd Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Sat, 7 Oct 2023 19:18:32 +0100 Subject: [PATCH 15/21] Removing test section of tut --- vignettes/hydro.Rmd | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/vignettes/hydro.Rmd b/vignettes/hydro.Rmd index 73619f09..d590d421 100644 --- a/vignettes/hydro.Rmd +++ b/vignettes/hydro.Rmd @@ -324,18 +324,6 @@ inlabru_model2 <- inlabru::bru(formula2, ) ) ``` -As `create_prediction_field` expects us to have coordinates in degrees of latitude and longitude we need to transform `inla_data` the `SpatialPointsDataFrame` we pass to INLA so that we can easily plot it. We first take the extent of the map - -```{r} -parsed_output = parse_model_output(inlabru_model2, inla_data) -``` - -```{r} -pred_field <- create_prediction_field(mesh, plot_type = "predicted_mean_fields", data_type = "gaussian", parsed_output$mean_post, parsed_output$fixed_mean) -pred_field -``` - - Now it's time to see its output, we can use the `fdmr::model_viewer` and pass in our model output, mesh and measurement data. This Shiny app helps analyse model output by automatically creating plots and a map of predictions for your data. Make sure you have a CRS set on your mesh. If not CRS can be read a default of `"+proj=longlat +datum=WGS84"` will be used which may lead to errors with data using UTM as we've done here. ```{r eval=FALSE} From ee1adc033d49df2d7c74fc14e090bd1d63b9ea88 Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Sun, 8 Oct 2023 10:23:40 +0100 Subject: [PATCH 16/21] Add changelog entry [skip ci] --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 63a84096..7bc47511 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Checks on the data types being passed into the `mesh_builder` tool - [PR #101](https://github.com/4DModeller/fdmr/pull/101) - Ability to plot either polygon or point data on Leaflet map of `mesh_builder` tool - [PR #101](https://github.com/4DModeller/fdmr/pull/101) - The ability to plot model predictions on a `leaflet` map in the our [Interactive priors Shiny app](https://4dmodeller.github.io/fdmr/articles/priors_app.html) - [PR #147](https://github.com/4DModeller/fdmr/pull/147) +- A new Shiny app to parse and plot INLA model output, letting users easily view model parameters and predictions on a map - [PR #158](https://github.com/4DModeller/fdmr/pull/158) ### Fixed From a8d0796a732dfeb74cafb3974281114a0db40c2e Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Sun, 8 Oct 2023 10:29:26 +0100 Subject: [PATCH 17/21] Correct text on crs [skip ci] --- vignettes/hydro.Rmd | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/vignettes/hydro.Rmd b/vignettes/hydro.Rmd index d590d421..6dcdd889 100644 --- a/vignettes/hydro.Rmd +++ b/vignettes/hydro.Rmd @@ -238,7 +238,7 @@ Create the triangulated mesh of the study region and have a quick look at it. ```{r error=TRUE,warning=FALSE, fig.width=7, fig.height=7} e <- era5_precip_cropped@extent resolution <- raster::res(era5_precip_cropped) -sr <- "+proj=utm +zone=32" +crs <- "+proj=utm +zone=32" xres <- resolution[1] * 2 yres <- resolution[2] * 2 @@ -248,7 +248,7 @@ xy <- xy[seq(1, nrow(xy), by = 2), ] colnames(xy) <- c("LONG", "LAT") -mesh <- INLA::inla.mesh.2d(loc = xy, max.edge = c(xres * 1, xres * 100), cutoff = 75, crs = sr) +mesh <- INLA::inla.mesh.2d(loc = xy, max.edge = c(xres * 1, xres * 100), cutoff = 75, crs = crs) plot(mesh) ``` @@ -324,7 +324,7 @@ inlabru_model2 <- inlabru::bru(formula2, ) ) ``` -Now it's time to see its output, we can use the `fdmr::model_viewer` and pass in our model output, mesh and measurement data. This Shiny app helps analyse model output by automatically creating plots and a map of predictions for your data. Make sure you have a CRS set on your mesh. If not CRS can be read a default of `"+proj=longlat +datum=WGS84"` will be used which may lead to errors with data using UTM as we've done here. +Now it's time to see its output, we can use the `fdmr::model_viewer` and pass in our model output, mesh and measurement data. This Shiny app helps analyse model output by automatically creating plots and a map of predictions for your data. Make sure you have a CRS set on your mesh. If no CRS can be read a default of `"+proj=longlat +datum=WGS84"` will be used which may lead to errors with data using UTM as we've done here. ```{r eval=FALSE} fdmr::model_viewer(model_output = inlabru_model2, mesh = mesh, measurement_data = inla_data, data_type = "Gaussian") From 203f47bf6ecf3a67dd1d40bb8114151708048bb5 Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Mon, 9 Oct 2023 16:11:28 +0100 Subject: [PATCH 18/21] Renam variable for data distribution type --- R/shiny_modelviewer.R | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/R/shiny_modelviewer.R b/R/shiny_modelviewer.R index 2ccef2fd..637c902d 100644 --- a/R/shiny_modelviewer.R +++ b/R/shiny_modelviewer.R @@ -3,13 +3,13 @@ #' @param model_output INLA model output #' @param mesh INLA mesh #' @param measurement_data Measurement data -#' @param data_type Type of data, Poisson or Gaussian +#' @param data_dist Type of data, Poisson or Gaussian #' #' @importFrom magrittr %>% #' #' @return shiny::app #' @keywords internal -model_viewer_shiny <- function(model_output, mesh, measurement_data, data_type) { +model_viewer_shiny <- function(model_output, mesh, measurement_data, data_dist) { busy_spinner <- get_busy_spinner() crs <- mesh$crs$input @@ -18,9 +18,9 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data, data_type) crs <- "+proj=longlat +datum=WGS84" } - data_type <- stringr::str_to_title(data_type) - if (!(data_type %in% c("Poisson", "Gaussian"))) { - stop("data_type must be one of Poisson or Gaussian") + data_dist <- stringr::str_to_title(data_dist) + if (!(data_dist %in% c("Poisson", "Gaussian"))) { + stop("data_dist must be one of Poisson or Gaussian") } parsed_model_output <- parse_model_output( @@ -53,7 +53,7 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data, data_type) shiny::column( 6, shiny::selectInput(inputId = "map_plot_type", label = "Plot type", choices = c("Predicted mean fields", "Random effect fields"), selected = "Predicted mean fields"), - shiny::selectInput(inputId = "map_data_type", label = "Data type", choices = c("Poisson", "Gaussian"), selected = data_type), + shiny::selectInput(inputId = "map_data_type", label = "Data type", choices = c("Poisson", "Gaussian"), selected = data_dist), ), shiny::column( 6, @@ -102,12 +102,12 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data, data_type) prediction_field <- shiny::reactive({ - data_type <- tolower(input$map_data_type) + data_dist <- tolower(input$map_data_type) if (input$map_plot_type == "Predicted mean fields") { create_prediction_field( mesh = mesh, plot_type = "predicted_mean_fields", - data_type = data_type, + data_dist = data_dist, var_a = parsed_model_output[["mean_post"]], var_b = parsed_model_output[["fixed_mean"]] ) @@ -115,7 +115,7 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data, data_type) create_prediction_field( mesh = mesh, plot_type = "random_effect_fields", - data_type = data_type, + data_dist = data_dist, var_a = parsed_model_output[["mean_post"]] ) } @@ -188,10 +188,10 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data, data_type) #' @param model_output INLA model output #' @param mesh INLA mesh #' @param measurement_data Measurement data -#' @param data_type Type of data, Poisson or Gaussian +#' @param data_dist Type of data, Poisson or Gaussian #' #' @return shiny::app #' @export -model_viewer <- function(model_output, mesh, measurement_data, data_type = "Poisson") { - shiny::runApp(model_viewer_shiny(model_output = model_output, mesh = mesh, measurement_data = measurement_data, data_type = data_type)) +model_viewer <- function(model_output, mesh, measurement_data, data_dist = "Poisson") { + shiny::runApp(model_viewer_shiny(model_output = model_output, mesh = mesh, measurement_data = measurement_data, data_dist = data_dist)) } From 0b6fcc1d1c9a3126ed61f4114fdad5ee59ce86d6 Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Mon, 9 Oct 2023 16:26:57 +0100 Subject: [PATCH 19/21] Fix CRS checks, rename arg for data distribution and pass through --- R/model_parse.R | 3 ++- R/shiny_modelviewer.R | 18 +++++++++--------- R/shiny_priors.R | 24 +++++++++++------------- vignettes/hydro.Rmd | 4 +++- 4 files changed, 25 insertions(+), 24 deletions(-) diff --git a/R/model_parse.R b/R/model_parse.R index 657e27ea..01d04284 100644 --- a/R/model_parse.R +++ b/R/model_parse.R @@ -66,6 +66,7 @@ create_prediction_field <- function(mesh, stop("Invalid plot type, select from ", valid_plots) } + data_dist <- tolower(data_dist) valid_data_dists <- c("poisson", "gaussian") if (!(data_dist %in% valid_data_dists)) { stop("Invalid data type, select from ", valid_data_dists) @@ -81,7 +82,7 @@ create_prediction_field <- function(mesh, if (plot_type == "predicted_mean_fields") { z <- base::as.numeric(A_proj %*% var_a[1:mesh$n]) + base::sum(var_b) - if (data_type == "poisson") z <- base::exp(z) + if (data_dist == "poisson") z <- base::exp(z) } else { # We get an error here as we only have 265 items # z <- var_a[1:mesh$n] diff --git a/R/shiny_modelviewer.R b/R/shiny_modelviewer.R index 637c902d..285e5604 100644 --- a/R/shiny_modelviewer.R +++ b/R/shiny_modelviewer.R @@ -3,13 +3,13 @@ #' @param model_output INLA model output #' @param mesh INLA mesh #' @param measurement_data Measurement data -#' @param data_dist Type of data, Poisson or Gaussian +#' @param data_distribution Type of data, Poisson or Gaussian #' #' @importFrom magrittr %>% #' #' @return shiny::app #' @keywords internal -model_viewer_shiny <- function(model_output, mesh, measurement_data, data_dist) { +model_viewer_shiny <- function(model_output, mesh, measurement_data, data_distribution) { busy_spinner <- get_busy_spinner() crs <- mesh$crs$input @@ -18,9 +18,9 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data, data_dist) crs <- "+proj=longlat +datum=WGS84" } - data_dist <- stringr::str_to_title(data_dist) - if (!(data_dist %in% c("Poisson", "Gaussian"))) { - stop("data_dist must be one of Poisson or Gaussian") + data_distribution <- stringr::str_to_title(data_distribution) + if (!(data_distribution %in% c("Poisson", "Gaussian"))) { + stop("data_distribution must be one of Poisson or Gaussian") } parsed_model_output <- parse_model_output( @@ -53,7 +53,7 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data, data_dist) shiny::column( 6, shiny::selectInput(inputId = "map_plot_type", label = "Plot type", choices = c("Predicted mean fields", "Random effect fields"), selected = "Predicted mean fields"), - shiny::selectInput(inputId = "map_data_type", label = "Data type", choices = c("Poisson", "Gaussian"), selected = data_dist), + shiny::selectInput(inputId = "map_data_type", label = "Data type", choices = c("Poisson", "Gaussian"), selected = data_distribution), ), shiny::column( 6, @@ -188,10 +188,10 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data, data_dist) #' @param model_output INLA model output #' @param mesh INLA mesh #' @param measurement_data Measurement data -#' @param data_dist Type of data, Poisson or Gaussian +#' @param data_distribution Type of data, Poisson or Gaussian #' #' @return shiny::app #' @export -model_viewer <- function(model_output, mesh, measurement_data, data_dist = "Poisson") { - shiny::runApp(model_viewer_shiny(model_output = model_output, mesh = mesh, measurement_data = measurement_data, data_dist = data_dist)) +model_viewer <- function(model_output, mesh, measurement_data, data_distribution = "Poisson") { + shiny::runApp(model_viewer_shiny(model_output = model_output, mesh = mesh, measurement_data = measurement_data, data_distribution = data_distribution)) } diff --git a/R/shiny_priors.R b/R/shiny_priors.R index 010e777e..3afd0a78 100644 --- a/R/shiny_priors.R +++ b/R/shiny_priors.R @@ -23,8 +23,6 @@ priors_shiny <- function(spatial_data, stop("We only support Poisson and Gaussian data") } - data_distribution_lower <- tolower(data_distribution) - got_coords <- has_coords(spatial_data = spatial_data) if (!got_coords) { stop("Please make sure you have set coordinates on spatial_data using sp::coordinates.") @@ -33,10 +31,10 @@ priors_shiny <- function(spatial_data, spatial_crs <- sp::proj4string(spatial_data) mesh_crs <- mesh$crs$input - if (is.na(mesh_crs) && is.na(spatial_crs)) { + if ((is.null(mesh_crs) || is.na(mesh_crs)) && (is.na(spatial_crs) || is.null(spatial_crs))) { warning("Cannot read CRS from mesh or spatial_data, using default CRS = +proj=longlat +datum=WGS84") crs <- "+proj=longlat +datum=WGS84" - } else if (is.na(mesh_crs)) { + } else if (is.na(mesh_crs) || is.null(mesh_crs)) { crs <- spatial_crs } else { crs <- mesh_crs @@ -159,7 +157,7 @@ priors_shiny <- function(spatial_data, ), shiny::column( 6, - shiny::selectInput(inputId = "data_dist", label = "Data distribution", choices = c("Poisson", "Gaussian")), + shiny::selectInput(inputId = "data_dist", label = "Data distribution", choices = c("Poisson", "Gaussian"), selected = data_distribution), ) ), shiny::checkboxGroupInput(inputId = "features", label = "Features", choices = features), @@ -329,7 +327,7 @@ priors_shiny <- function(spatial_data, formula_str() }) - data_distribution <- shiny::reactive({ + data_distribution_internal <- shiny::reactive({ tolower(input$data_dist) }) @@ -338,7 +336,7 @@ priors_shiny <- function(spatial_data, formula_local <- inla_formula() measurement_data_local <- measurement_data - data_dist_local <- data_distribution() + data_dist_local <- data_distribution_internal() family_control <- NULL if (data_dist_local == "poisson") { family_control <- list(link = "log") @@ -443,7 +441,7 @@ priors_shiny <- function(spatial_data, create_prediction_field( mesh = mesh, plot_type = "predicted_mean_fields", - data_dist = data_distribution(), + data_dist = data_distribution_internal(), var_a = data[["mean_post"]], var_b = data[["fixed_mean"]] ) @@ -451,7 +449,7 @@ priors_shiny <- function(spatial_data, create_prediction_field( mesh = mesh, plot_type = "random_effect_fields", - data_dist = data_distribution(), + data_dist = data_distribution_internal(), var_a = data[["mean_post"]] ) } @@ -529,7 +527,7 @@ priors_shiny <- function(spatial_data, params <- model_vals$run_params[[input$select_run_code]] family_control_str <- "NULL" - if (data_distribution() == "poisson") { + if (data_distribution_internal() == "poisson") { family_control_str <- "list(link = 'log')," } @@ -546,7 +544,7 @@ priors_shiny <- function(spatial_data, )", "\n\n", paste0("model_output <- inlabru::bru(formula, data = measurement_data, - family = '", data_distribution(), "', + family = '", data_distribution_internal(), "', E = measurement_data[[", input$exposure_param, "]], control.family = ", family_control_str, " options = list( @@ -571,6 +569,6 @@ priors_shiny <- function(spatial_data, #' #' @return shiny::app #' @export -interactive_priors <- function(spatial_data, measurement_data, time_variable, mesh, log_folder = NULL) { - shiny::runApp(priors_shiny(spatial_data = spatial_data, measurement_data = measurement_data, time_variable = time_variable, mesh = mesh, log_folder = log_folder)) +interactive_priors <- function(spatial_data, measurement_data, time_variable, mesh, data_distribution = "Poisson", log_folder = NULL) { + shiny::runApp(priors_shiny(spatial_data = spatial_data, measurement_data = measurement_data, time_variable = time_variable, mesh = mesh, data_distribution = data_distribution, log_folder = log_folder)) } diff --git a/vignettes/hydro.Rmd b/vignettes/hydro.Rmd index 6dcdd889..f0e7e20d 100644 --- a/vignettes/hydro.Rmd +++ b/vignettes/hydro.Rmd @@ -327,9 +327,11 @@ inlabru_model2 <- inlabru::bru(formula2, Now it's time to see its output, we can use the `fdmr::model_viewer` and pass in our model output, mesh and measurement data. This Shiny app helps analyse model output by automatically creating plots and a map of predictions for your data. Make sure you have a CRS set on your mesh. If no CRS can be read a default of `"+proj=longlat +datum=WGS84"` will be used which may lead to errors with data using UTM as we've done here. ```{r eval=FALSE} -fdmr::model_viewer(model_output = inlabru_model2, mesh = mesh, measurement_data = inla_data, data_type = "Gaussian") +fdmr::model_viewer(model_output = inlabru_model2, mesh = mesh, measurement_data = inla_data, data_distribution = "Gaussian") ``` +We can also view the summary of the model outputs + ```{r error=TRUE} summary(inlabru_model1) ``` From 18c3ccc5fb43a744335ea8bc141cf6cb7c105c37 Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Mon, 9 Oct 2023 17:20:17 +0100 Subject: [PATCH 20/21] Add rename of priors app back --- NAMESPACE | 2 +- R/shiny_priors.R | 6 +++--- man/create_prediction_field.Rd | 4 ++-- man/model_builder.Rd | 1 + man/model_builder_shiny.Rd | 3 +++ man/model_viewer.Rd | 9 ++++++++- man/model_viewer_shiny.Rd | 4 +++- 7 files changed, 21 insertions(+), 8 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 4235daea..2d2f89f6 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -8,8 +8,8 @@ export(load_tutorial_data) export(mesh_builder) export(mesh_checker) export(mesh_to_spatial) -export(model_viewer) export(model_builder) +export(model_viewer) export(numbers_only) export(parse_model_output) export(plot_barchart) diff --git a/R/shiny_priors.R b/R/shiny_priors.R index b36a63a8..fe374223 100644 --- a/R/shiny_priors.R +++ b/R/shiny_priors.R @@ -11,7 +11,7 @@ #' #' @return shiny::app #' @keywords internal -priors_shiny <- function(spatial_data, +model_builder_shiny <- function(spatial_data, measurement_data, time_variable, mesh, @@ -568,6 +568,6 @@ priors_shiny <- function(spatial_data, #' #' @return shiny::app #' @export -interactive_priors <- function(spatial_data, measurement_data, time_variable, mesh, data_distribution = "Poisson", log_folder = NULL) { - shiny::runApp(priors_shiny(spatial_data = spatial_data, measurement_data = measurement_data, time_variable = time_variable, mesh = mesh, data_distribution = data_distribution, log_folder = log_folder)) +model_builder <- function(spatial_data, measurement_data, time_variable, mesh, data_distribution = "Poisson", log_folder = NULL) { + shiny::runApp(model_builder_shiny(spatial_data = spatial_data, measurement_data = measurement_data, time_variable = time_variable, mesh = mesh, data_distribution = data_distribution, log_folder = log_folder)) } diff --git a/man/create_prediction_field.Rd b/man/create_prediction_field.Rd index 5ad702f1..aa3a115a 100644 --- a/man/create_prediction_field.Rd +++ b/man/create_prediction_field.Rd @@ -7,7 +7,7 @@ create_prediction_field( mesh, plot_type = "predicted_mean_fields", - data_type = "poisson", + data_dist = "poisson", var_a = NULL, var_b = NULL ) @@ -17,7 +17,7 @@ create_prediction_field( \item{plot_type}{Type of plot to create, "predicted_mean_fields" etc} -\item{data_type}{Type of data, "poisson" etc} +\item{data_dist}{Type of data, "poisson" etc} \item{var_a}{Data for variable a, required for "predicted_mean_fields" and "random_effect_fields"} diff --git a/man/model_builder.Rd b/man/model_builder.Rd index 720f11f8..40bfeed7 100644 --- a/man/model_builder.Rd +++ b/man/model_builder.Rd @@ -9,6 +9,7 @@ model_builder( measurement_data, time_variable, mesh, + data_distribution = "Poisson", log_folder = NULL ) } diff --git a/man/model_builder_shiny.Rd b/man/model_builder_shiny.Rd index dabf61a3..4918cabf 100644 --- a/man/model_builder_shiny.Rd +++ b/man/model_builder_shiny.Rd @@ -9,6 +9,7 @@ model_builder_shiny( measurement_data, time_variable, mesh, + data_distribution = "Poisson", log_folder = NULL ) } @@ -21,6 +22,8 @@ model_builder_shiny( \item{mesh}{INLA mesh} +\item{data_distribution}{Data distribution, Poisson or Gaussian} + \item{log_folder}{Folder to write out logs} } \value{ diff --git a/man/model_viewer.Rd b/man/model_viewer.Rd index 5947add3..8ea96f7e 100644 --- a/man/model_viewer.Rd +++ b/man/model_viewer.Rd @@ -4,7 +4,12 @@ \alias{model_viewer} \title{Mesh building shiny app. Creates and visualises a mesh from some spatial data.} \usage{ -model_viewer(model_output, mesh, measurement_data) +model_viewer( + model_output, + mesh, + measurement_data, + data_distribution = "Poisson" +) } \arguments{ \item{model_output}{INLA model output} @@ -12,6 +17,8 @@ model_viewer(model_output, mesh, measurement_data) \item{mesh}{INLA mesh} \item{measurement_data}{Measurement data} + +\item{data_distribution}{Type of data, Poisson or Gaussian} } \value{ shiny::app diff --git a/man/model_viewer_shiny.Rd b/man/model_viewer_shiny.Rd index 7bd6353b..1e0fbcb9 100644 --- a/man/model_viewer_shiny.Rd +++ b/man/model_viewer_shiny.Rd @@ -4,7 +4,7 @@ \alias{model_viewer_shiny} \title{Parse inlabru model output} \usage{ -model_viewer_shiny(model_output, mesh, measurement_data) +model_viewer_shiny(model_output, mesh, measurement_data, data_distribution) } \arguments{ \item{model_output}{INLA model output} @@ -12,6 +12,8 @@ model_viewer_shiny(model_output, mesh, measurement_data) \item{mesh}{INLA mesh} \item{measurement_data}{Measurement data} + +\item{data_distribution}{Type of data, Poisson or Gaussian} } \value{ shiny::app From 27572088bd7bb767e73307e7d431ea23b54e735b Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Mon, 9 Oct 2023 17:30:02 +0100 Subject: [PATCH 21/21] Formatting fix [skip ci] --- R/shiny_priors.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/R/shiny_priors.R b/R/shiny_priors.R index fe374223..169e8630 100644 --- a/R/shiny_priors.R +++ b/R/shiny_priors.R @@ -12,11 +12,11 @@ #' @return shiny::app #' @keywords internal model_builder_shiny <- function(spatial_data, - measurement_data, - time_variable, - mesh, - data_distribution = "Poisson", - log_folder = NULL) { + measurement_data, + time_variable, + mesh, + data_distribution = "Poisson", + log_folder = NULL) { future::plan(future::multisession()) if (!(data_distribution %in% c("Poisson", "Gaussian"))) {