Adding the fit_model function

NAMESPACE

@@ -2,6 +2,7 @@
 
 export(clear_caches)
 export(create_prediction_field)
+export(fit_model)
 export(get_tutorial_datapath)
 export(latlong_to_utm)
 export(load_tutorial_data)

R/fit_model.R

@@ -0,0 +1,72 @@
+#' Fit a model using inlabru::bru
+#'
+#' @param y outcome variable name
+#' @param data Spatial data
+#' @param locations list of coordinates associated with process
+#' @param space_variables list of spatial variables
+#' @param spacetime_variables list of spacetime variables
+#' @param family gaussian or poisson
+#' @param latitude_col name of latitude column
+#' @param longitude_col name of longitude column
+#' @param time_variable name of time variable in data
+#' @param initial_range_mesh Initial range in mesh creation
+#' @param max_edge_mesh Max edge value used in mesh creation
+#'
+#' @return list of model output and mesh
+#' @export
+fit_model <- function(
+    y,
+    data,
+    locations,
+    space_variables,
+    spacetime_variables,
+    time_variable = "time",
+    family = "gaussian",
+    latitude_col = "LAT",
+    longitude_col = "LONG",
+    initial_range_mesh = 0.1,
+    max_edge_mesh = 0.05,
+    verbose = FALSE) {
+    library("inlabru")
+    library("INLA")
+
+    print("Creating mesh...")
+    mesh <- fmesher::fm_mesh_2d_inla(
+        loc = locations
+    )
+
+    print("Creating SPDE...")
+    spde <- INLA::inla.spde2.pcmatern(
+        mesh = mesh,
+        prior.range = c(initial_range_mesh, 0.5),
+        prior.sigma = c(1, 0.01)
+    )
+
+    space_variables <- append(space_variables, y)
+    spacetime_variables <- append(spacetime_variables, y)
+
+    rhoprior <- base::list(theta = list(prior = "pccor1", param = c(0, 0.9)))
+
+    print("Running model...")
+    model_output <- inlabru::bru(
+        components = ~ space(coordinates, model = spde) +
+            spacetime(
+                coordinates,
+                model = spde,
+                group = time,
+                control.group = list(
+                    model = "ar1",
+                    hyper = rhoprior
+                )
+            ) +
+            beta_u(1, prec.linear = 1),
+        inlabru::like(formula = data[[y]] ~ space, family = "gaussian", data = data[space_variables]),
+        inlabru::like(formula = data[[y]] ~ beta_u * space + spacetime, family = "gaussian", data = data[spacetime_variables]),
+        options = list(
+            verbose = verbose
+        )
+    )
+
+    print("Model run complete.")
+    list(model_output = model_output, mesh = mesh)
+}

vignettes/TEST_fitmodel.Rmd

@@ -0,0 +1,74 @@
+---
+title: "TESTING FITMODEL"
+output: 
+  bookdown::html_document2:
+    base_format: rmarkdown::html_vignette
+    fig_caption: yes
+bibliography: "covid.bib"
+link-citations: yes
+vignette: >
+  %\VignetteIndexEntry{Fitting a model}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+this file tests all of the different datasets using `fit_model` function
+
+
+```{r}
+devtools::load_all()
+```
+
+# COVID-19
+```{r}
+# dont run because it takes forever
+# fdmr::retrieve_tutorial_data(dataset = "covid")
+# covid19_data <- fdmr::load_tutorial_data(dataset = "covid", filename = "covid19_data.rds")
+# covid19_data = covid19_data %>% dplyr::filter(week < 5)
+# sp::coordinates(covid19_data) <- c("LONG", "LAT")
+# fdmr::fit_model(y='cases'
+#                 , locations = covid19_data
+#                 ,time_variable = 'week'
+#                 ,latitude_col = 'LAT'
+#                 ,longitude_col = 'LONG'
+#                 ,data = covid19_data
+#                 ,space_variables = c('Population','carebeds.ratio')
+#                 ,spacetime_variables = c('pm25', 'no2'))
+```
+
+# COVID19 Bristol
+```{r}
+fdmr::retrieve_tutorial_data(dataset = "priors")
+covid19_data <- fdmr::load_tutorial_data(dataset = "priors", filename = "covid19_dataBris.rds")
+sp::coordinates(covid19_data) <- c("LONG", "LAT")
+run_output <- fdmr::fit_model(y = 'cases'
+                ,locations=covid19_data
+                ,time_variable = 'week'
+                ,latitude_col = 'LAT'
+                ,longitude_col = 'LONG'
+                ,data = covid19_data
+                ,space_variables = c('Population')
+                ,spacetime_variables = c('prevalence'))
+summary(run_output$model_output)
+```
+
+
+# hydro
+```{r}
+fdmr::retrieve_tutorial_data(dataset = "tibetan_lakes")
+lakes_data_filepath <- fdmr::get_tutorial_datapath(dataset = "tibetan_lakes", filename = "example_lake_data.csv")
+lakes_data <- read.csv(lakes_data_filepath)
+head(lakes_data)
+sp::coordinates(lakes_data) <- c("centroid_lon", "centroid_lat")
+run_output <- fdmr::fit_model(y = 'water_balance_m3'
+                ,locations=lakes_data
+                ,time_variable = 'year'
+                ,latitude_col = 'centroid_lat'
+                ,longitude_col = 'centroid_lon'
+                ,data = lakes_data
+                ,space_variables = c('SUB_AREA')
+                ,spacetime_variables = c('precip'))
+summary(run_output$model_output)
+```
+
+

vignettes/fit_model.Rmd

@@ -0,0 +1,55 @@
+---
+title: "Fitting a model"
+output: 
+  bookdown::html_document2:
+    base_format: rmarkdown::html_vignette
+    fig_caption: yes
+bibliography: "covid.bib"
+link-citations: yes
+vignette: >
+  %\VignetteIndexEntry{Fitting a model}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+```{r}
+library(magrittr)
+```
+
+In order to show off the `fit_model` function of 4d-modeller we will try and predict when it rains during the year 1999, the peak year of human civilization. First we will import the data.
+
+```{r error=true}
+raindata <- read.csv('fitmodel_data/rain_data.csv')
+# TODO : gareth make this work
+# fdmr::plot_map(raster_data = raindata)
+df = raindata %>% dplyr::filter(Rain > 5) %>% head(15) %>% dplyr::arrange("Time")
+df
+```
+```{r}
+sp::coordinates(df) <- c("Longitude", "Latitude")
+```
+
+## Fit the model
+
+Now we're ready to fit the model
+
+
+# ```{r}
+# fdmr::fit_model
+# ```
+
+```{r}
+run_output <- fdmr::fit_model(y="Rain", data=df, process_coords = df, time_variable="Time_Index", latitude_col = "Latitude", longitude_col = "Longitude", verbose=TRUE)
+# barbara streisand
+```
+
+Here `run_output` is a list containing the mesh created and the model output.
+
+```{r}
+summary(run_output$model_output)
+```
+
+We can also view the output using the `model_viewer` Shiny app.
+
+```{r}
+fdmr::model_viewer(model_output = run_output$model_output, mesh  = run_output$mesh, measurement_data = data, data_distribution = "Gaussian")
+```

vignettes/fitmodel_data/datasets.md

@@ -0,0 +1,5 @@
+
+# Data sets for fitmodel functions
+
+# NCAR data
+https://www.kaggle.com/datasets/shiratorizawa/ncarcsv2

vignettes/fitmodel_data/make_final_data.py

@@ -0,0 +1,21 @@
+import numpy as np # linear algebra
+import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
+import matplotlib.pyplot as plt
+
+def main():
+    df = pd.read_csv('NCAR_raw.csv')
+    final_data = df[['lat', 'lon', 'time', 'rain', 'slp', 'tmp']].replace(-9999.0, 0.0).fillna(0).copy()
+    final_data.columns = ['Latitude', 'Longitude', 'Time', 'Rain', 'Sea_Level_Pressure', 'Temperature']
+    final_data.sort_values(by='Time', inplace=True)
+    final_data.reset_index(inplace=True)
+    final_data['Time'] = pd.to_datetime(final_data['Time'])
+    final_data = final_data[final_data.Time.between('1999-01-01', '1999-12-31')]
+    # convert timestamps to time index for inlabru
+    time_index = final_data.Time.apply(lambda t: t.timestamp())
+    final_data['Time_Index'] = ((time_index - time_index.min())/86400).astype(int)
+    final_data['Sea_Level_Pressure'] = final_data.Sea_Level_Pressure.replace(0, np.nan).fillna(method='ffill')
+    final_data.fillna(final_data['Sea_Level_Pressure'].mean(), inplace=True)
+    final_data.to_csv('rain_data.csv', index=False)
+
+if __name__=='__main__':
+    main()