reslr.R
In reslr: Modelling Relative Sea Level Data

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"#,
  #fig.path = "vig/"
)
options(rmarkdown.html_vignette.check_title = FALSE)

## ----eval = FALSE,message=FALSE-----------------------------------------------
#  # install.packages("reslr")
#  # library(devtools)
#  # devtools::install()
#  #devtools::install_github("maeveupton/reslr")
#  install_github("maeveupton/reslr")

## ---- readpkg,eval = TRUE, message=FALSE--------------------------------------
library(reslr)

## ---- example1site,eval = TRUE------------------------------------------------
# For 1 site
data_1site <- reslr::NAACproxydata %>% dplyr::filter(Site == "Cedar Island")
# For multiple sites
data_multisite <- reslr::NAACproxydata %>% dplyr::filter(Site %in% c(
  "Snipe Key", "Cheesequake",
  "Placentia", "Leeds Point"
))

## ---- example1dataagain, eval = TRUE------------------------------------------
# For 1 site
CedarIslandNC <- NAACproxydata %>% dplyr::filter(Site == "Cedar Island")

## ---- loadslr, eval = TRUE----------------------------------------------------
CedarIslandNC_input <- reslr_load(
  data = CedarIslandNC,
  include_tide_gauge = FALSE,
  include_linear_rate = FALSE,
  TG_minimum_dist_proxy = FALSE,
  list_preferred_TGs = NULL,
  all_TG_1deg = FALSE,
  prediction_grid_res = 50,
  input_age_type = "CE",
  sediment_average_TG = 10
)

## ---- slrdata,eval=TRUE-------------------------------------------------------
data <- CedarIslandNC_input$data

## ---- datagridslr, eval = TRUE------------------------------------------------
data_grid <- CedarIslandNC_input$data_grid

## ---- printdata, eval=TRUE----------------------------------------------------
print(CedarIslandNC_input)

## ---- plotdata,fig.align = 'center',fig.width = 7,fig.height = 5,eval = TRUE----
plot(
  x = CedarIslandNC_input,
  title = "Plot of the raw data",
  xlab = "Year (CE)",
  ylab = "Relative Sea Level (m)",
  plot_tide_gauges = FALSE,
  plot_proxy_records = TRUE,
  plot_caption = TRUE
)

## ---- runslr,eval = TRUE------------------------------------------------------
res_eiv_slr_t <- reslr_mcmc(
  input_data = CedarIslandNC_input,
  model_type = "eiv_slr_t",
  CI = 0.95
)

## ---- printrunslr,eval=TRUE---------------------------------------------------
print(res_eiv_slr_t)

## ---- summaryslr, eval = TRUE-------------------------------------------------
summary(res_eiv_slr_t)

## ---- runslrmoreit,eval = FALSE-----------------------------------------------
#  res_eiv_slr_t <- reslr_mcmc(
#    input_data = CedarIslandNC_input,
#    model_type = "eiv_slr_t",
#    # Update these values
#    n_iterations = 6000, # Number of iterations
#    n_burnin = 1000, # Number of iterations to discard at the beginning
#    n_thin = 4, # Reduces number of output samples to save memory and computation time
#    n_chains = 3 # Number of Markov chains
#  )

## ---- plotslrres, fig.align = 'center',fig.width = 7,fig.height = 5,eval = TRUE----
plot(res_eiv_slr_t,
  xlab = "Year (CE)",
  ylab = "Relative Sea Level (m)"
)

## ---- dataframeslrres, eval = TRUE--------------------------------------------
output_dataframes <- res_eiv_slr_t$output_dataframes
head(output_dataframes)

## ----exampledata1, eval = TRUE------------------------------------------------
# For 1 site
CedarIslandNC <- reslr::NAACproxydata %>% dplyr::filter(Site == "Cedar Island")

## ----loaddatacp, eval = TRUE--------------------------------------------------
CedarIslandNC_input <- reslr_load(
  data = CedarIslandNC,
  include_tide_gauge = FALSE,
  include_linear_rate = FALSE,
  TG_minimum_dist_proxy = FALSE,
  list_preferred_TGs = NULL,
  all_TG_1deg = FALSE,
  prediction_grid_res = 50,
  sediment_average_TG = 10
)

## ----datacp,eval=TRUE---------------------------------------------------------
data <- CedarIslandNC_input$data
head(data)

## ----datagridcp, eval = TRUE--------------------------------------------------
data_grid <- CedarIslandNC_input$data_grid
head(data_grid)

## ----printcp, eval=TRUE-------------------------------------------------------
print(CedarIslandNC_input)

## ----plotdatacp,fig.align = 'center',fig.width = 7,fig.height = 5,eval = TRUE----
plot(
  x = CedarIslandNC_input,
  title = "Plot of the raw data",
  xlab = "Year (CE)",
  ylab = "Relative Sea Level (m)",
  plot_proxy_records = TRUE,
  plot_tide_gauges = FALSE
)

## ----runcp1, eval = TRUE------------------------------------------------------
res_eiv_cp1_t <- reslr_mcmc(
  input_data = CedarIslandNC_input,
  model_type = "eiv_cp_t",
  n_cp = 1,
  CI = 0.95
)

## ----runcp2, eval = FALSE-----------------------------------------------------
#  res_eiv_cp2_t <- reslr_mcmc(
#    input_data = CedarIslandNC_input,
#    model_type = "eiv_cp_t",
#    n_cp = 2, # Updating the default setting to include an additional change point.
#    CI = 0.95
#  )

## ----printcpout, eval=TRUE----------------------------------------------------
print(res_eiv_cp1_t)

## ----summarycp, eval = TRUE---------------------------------------------------
summary(res_eiv_cp1_t)

## ----runcpmoreit,eval = FALSE-------------------------------------------------
#  res_eiv_cp1_t <- reslr_mcmc(
#    input_data = CedarIslandNC_input,
#    model_type = "eiv_cp_t",
#    # Update these values
#    n_iterations = 6000, # Number of iterations
#    n_burnin = 1000, # Number of iterations to discard at the beginning
#    n_thin = 4, # Reduces number of output samples to save memory and computation time
#    n_chains = 3 # Number of Markov chains
#  )

## ---- plotcpres, fig.align = 'center',fig.width = 7,fig.height = 5,eval = TRUE----
plot(res_eiv_cp1_t,
  xlab = "Year (CE)",
  ylab = "Relative Sea Level (m)",
)

## ----cpdataframe, eval = TRUE-------------------------------------------------
output_dataframes <- res_eiv_cp1_t$output_dataframes
head(output_dataframes)

## ----exampledataset2 ,eval = TRUE---------------------------------------------
# For 1 site
CedarIslandNC <- reslr::NAACproxydata %>% dplyr::filter(Site == "Cedar Island")

## ---- loadigp ,eval = TRUE----------------------------------------------------
CedarIslandNC_input <- reslr_load(
  data = CedarIslandNC,
  include_tide_gauge = FALSE,
  include_linear_rate = FALSE,
  TG_minimum_dist_proxy = FALSE,
  list_preferred_TGs = NULL,
  all_TG_1deg = FALSE,
  prediction_grid_res = 50,
  sediment_average_TG = 10
)

## ----dataigp,eval=FALSE-------------------------------------------------------
#  data <- CedarIslandNC_input$data

## ----datagridigp, eval = FALSE------------------------------------------------
#  data_grid <- CedarIslandNC_input$data_grid

## ----printigp, eval=TRUE------------------------------------------------------
print(CedarIslandNC_input)

## ---- plotigpdata,fig.align = 'center',fig.width = 7,fig.height = 5,eval = TRUE----
plot(
  x = CedarIslandNC_input,
  title = "Plot of the raw data",
  xlab = "Year (CE)",
  ylab = "Relative Sea Level (m)",
  plot_proxy_records = TRUE,
  plot_tide_gauges = FALSE
)

## ----runigp, eval = FALSE-----------------------------------------------------
#  res_eiv_igp_t <- reslr_mcmc(
#    input_data = CedarIslandNC_input,
#    model_type = "eiv_igp_t",
#    CI = 0.95,
#  
#  )

## ----printigpout, eval=FALSE--------------------------------------------------
#  print(res_eiv_igp_t)

## ----summaryigp, eval = FALSE-------------------------------------------------
#  summary(res_eiv_igp_t)

## ----runigpmore, eval = FALSE-------------------------------------------------
#  res_eiv_igp_t <- reslr_mcmc(
#    input_data = CedarIslandNC_input,
#    model_type = "eiv_igp_t",
#    # Update these values
#    n_iterations = 6000, # Number of iterations
#    n_burnin = 1000, # Number of iterations to discard at the beginning
#    n_thin = 4, # Reduces number of output samples to save memory and computation time
#    n_chains = 3 # Number of Markov chains
#  )

## ----plotigpres, fig.align = 'center',fig.width = 7,fig.height = 5,eval = FALSE----
#  plot(res_eiv_igp_t,
#    plot_type = "model_fit_plot",
#    xlab = "Year (CE)",
#    ylab = "Relative Sea Level (m)",
#    plot_proxy_records = TRUE,
#    plot_tide_gauges = FALSE
#  )

## ----plotigpres_load, eval = TRUE,echo=FALSE, fig.align = 'center',out.width="100%"----
url <- "https://raw.githubusercontent.com/maeveupton/reslr/main/reslrvigplots/plotigpres-1.png"
knitr::include_graphics(url)

## ----plotigpresrate, fig.align = 'center',fig.width = 7,fig.height = 5,eval = FALSE----
#  plot(res_eiv_igp_t,
#    plot_type = "rate_plot",
#    xlab = "Year (CE)",
#    y_rate_lab = "Rate of Change (mm per year)"
#  )

## ----plotigpresrate_load, eval = TRUE,echo=FALSE, fig.align = 'center',out.width="100%"----
url <- "https://raw.githubusercontent.com/maeveupton/reslr/main/reslrvigplots/plotigpresrate-1.png"
knitr::include_graphics(url)

## ----igpdataout, eval = FALSE-------------------------------------------------
#  output_dataframes <- res_eiv_igp_t$output_dataframes

## ---- dataexample3,eval = TRUE------------------------------------------------
# For 1 site
CedarIslandNC <- reslr::NAACproxydata %>% dplyr::filter(Site == "Cedar Island")

## ----loadspt, eval = TRUE-----------------------------------------------------
CedarIslandNC_input <- reslr_load(
  data = CedarIslandNC,
  include_tide_gauge = FALSE,
  include_linear_rate = FALSE,
  TG_minimum_dist_proxy = FALSE,
  list_preferred_TGs = NULL,
  all_TG_1deg = FALSE,
  prediction_grid_res = 50,
  sediment_average_TG = 10
)

## ----dataspt, eval=TRUE-------------------------------------------------------
data <- CedarIslandNC_input$data

## ----datagridspt, eval = TRUE-------------------------------------------------
data_grid <- CedarIslandNC_input$data_grid

## ----printspt, eval=TRUE------------------------------------------------------
print(CedarIslandNC_input)

## ----plotspt,fig.align = 'center',fig.width = 7,fig.height = 5,eval = TRUE----
plot(
  x = CedarIslandNC_input,
  title = "Plot of the raw data",
  xlab = "Year (CE)",
  ylab = "Relative Sea Level (m)",
  plot_proxy_records = TRUE,
  plot_tide_gauges = FALSE
)

## ----runspt,eval = TRUE,message=FALSE,results='hide'--------------------------
res_ni_spline_t <- reslr_mcmc(
  input_data = CedarIslandNC_input,
  model_type = "ni_spline_t",
  CI = 0.95
)

## ----printresspt, eval=TRUE---------------------------------------------------
print(res_ni_spline_t)

## ----summaryspt, eval = TRUE--------------------------------------------------
summary(res_ni_spline_t)

## ---- runsptmore,eval = FALSE-------------------------------------------------
#  res_ni_spline_t <- reslr_mcmc(
#    input_data = CedarIslandNC,
#    model_type = "ni_spline_t",
#    # Update these values
#    n_iterations = 6000, # Number of iterations
#    n_burnin = 1000, # Number of iterations to discard at the beginning
#    n_thin = 4, # Reduces number of output samples to save memory and computation time
#    n_chains = 3 # Number of Markov chains
#  )

## ---- plotsptres,fig.align = 'center',fig.width = 7,fig.height = 5,eval = TRUE----
plot(res_ni_spline_t,
  plot_type = "model_fit_plot",
  xlab = "Year (CE)",
  ylab = "Relative Sea Level (m)"
)

## ----plotsptresrate, fig.align = 'center',fig.width = 7,fig.height = 5,eval = TRUE----
plot(res_ni_spline_t,
  plot_type = "rate_plot",
  xlab = "Year (CE)",
  y_rate_lab = "Rate of Change (mm per year)"
)

## ----outdfspt, eval = TRUE----------------------------------------------------
output_dataframes <- res_ni_spline_t$output_dataframes
head(output_dataframes)

## ----data2sites, eval = TRUE--------------------------------------------------
# For 2 site
multi_site <- reslr::NAACproxydata %>%
  dplyr::filter(Site %in% c("Cedar Island", "Nassau"))

## ----loadspst, eval = TRUE----------------------------------------------------
multi_site_input <- reslr_load(
  data = multi_site,
  include_tide_gauge = FALSE,
  include_linear_rate = FALSE,
  TG_minimum_dist_proxy = FALSE,
  list_preferred_TGs = NULL,
  all_TG_1deg = FALSE,
  prediction_grid_res = 50,
  sediment_average_TG = 10
)

## ----dataspst,eval=TRUE-------------------------------------------------------
data <- multi_site_input$data
head(data)

## ----datagridspst, eval = TRUE------------------------------------------------
data_grid <- multi_site_input$data_grid
head(data_grid)

## ----printspst, eval=TRUE-----------------------------------------------------
print(multi_site_input)

## ----plotspst,fig.align = 'center',fig.width = 7,fig.height = 10,eval = TRUE----
plot(
  x = multi_site_input,
  title = "Plot of the raw data",
  xlab = "Year (CE)",
  ylab = "Relative Sea Level (m)",
  plot_proxy_records = TRUE,
  plot_tide_gauges = FALSE
)

## ----runspst,eval = FALSE,results='hide',message=FALSE------------------------
#  res_ni_spline_st <- reslr_mcmc(
#    input_data = multi_site_input,
#    model_type = "ni_spline_st",
#    CI = 0.95
#  )

## ----printspstout,eval=FALSE--------------------------------------------------
#  print(res_ni_spline_st)

## ----summaryspst, eval = FALSE------------------------------------------------
#  summary(res_ni_spline_st)

## ---- runspstmore,eval = FALSE------------------------------------------------
#  res_ni_spline_st <- reslr::reslr_mcmc(
#    input_data = multi_site_input,
#    model_type = "ni_spline_st",
#    # Update these values
#    n_iterations = 6000, # Number of iterations
#    n_burnin = 1000, # Number of iterations to discard at the beginning
#    n_thin = 4, # Reduces number of output samples to save memory and computation time
#    n_chains = 3 # Number of Markov chains
#  )

## ----plotspstres, eval = FALSE,fig.align = 'center',fig.width = 7,fig.height = 5----
#  plot(res_ni_spline_st,
#    plot_type = "model_fit_plot",
#    xlab = "Year (CE)",
#    ylab = "Relative Sea Level (m)"
#  )

## ----plotspstres_load, eval = TRUE,echo=FALSE, fig.align = 'center',out.width="100%"----
url <- "https://raw.githubusercontent.com/maeveupton/reslr/main/reslrvigplots/plotspstres-1.png"
knitr::include_graphics(url)

## ----plotspstresrate, fig.align = 'center',fig.width = 7,fig.height = 5,eval =FALSE,results='hide',message=FALSE----
#  plot(res_ni_spline_st,
#    plot_type = "rate_plot",
#    xlab = "Year (CE)",
#    y_rate_lab = "Rate of Change (mm per year)"
#  )

## ----plotspstresrate_load, eval = TRUE,echo=FALSE, fig.align = 'center',out.width="100%"----
url <- "https://raw.githubusercontent.com/maeveupton/reslr/main/reslrvigplots/plotspstresrate-1.png"
knitr::include_graphics(url)

## ----outdfspst, eval = FALSE--------------------------------------------------
#  output_dataframes <- res_ni_spline_st$output_dataframes

## ----data2sitesmore, eval = TRUE----------------------------------------------
# For 9 site
multi_9_sites <- reslr::NAACproxydata %>%
  dplyr::filter(Site %in% c(
    "Cedar Island", "Nassau", "Snipe Key",
    "Placentia", "Cape May Courthouse", "East River Marsh",
    "Fox Hill Marsh", "Swan Key", "Big River Marsh"
  ))

## ----loadnigam, eval = TRUE---------------------------------------------------
multi_9_sites_input <- reslr_load(
  data = multi_9_sites,
  include_tide_gauge = TRUE,
  include_linear_rate = TRUE,
  TG_minimum_dist_proxy = FALSE,
  list_preferred_TGs = NULL,
  all_TG_1deg = TRUE,
  prediction_grid_res = 50,
  sediment_average_TG = 10
)

## ----datanigam,eval=FALSE-----------------------------------------------------
#  data <- multi_9_sites_input$data

## ----datagridnigam, eval = FALSE----------------------------------------------
#  data_grid <- multi_9_sites_input$data_grid

## ----printnigam, eval=TRUE----------------------------------------------------
print(multi_9_sites_input)

## ---- plotnigam,fig.align = 'center',fig.width = 7,fig.height = 10,eval = TRUE----
plot(
  x = multi_9_sites_input,
  title = "Plot of the raw data",
  xlab = "Year (CE)",
  ylab = "Relative Sea Level (m)",
  plot_proxy_records = TRUE,
  plot_tide_gauges = TRUE
)

## ----runnigam,eval = FALSE,results='hide',message=FALSE-----------------------
#  res_ni_gam_decomp <- reslr_mcmc(
#    input_data = multi_9_sites_input,
#    model_type = "ni_gam_decomp",
#    CI = 0.95
#  )

## ----printnigamout, eval=FALSE------------------------------------------------
#  print(res_ni_gam_decomp)

## ----summarynigam, eval = FALSE-----------------------------------------------
#  summary(res_ni_gam_decomp)

## ---- runnigammore,eval = FALSE-----------------------------------------------
#  res_ni_gam_decomp <- reslr_mcmc(
#    input_data = multi_9_sites_input,
#    model_type = "ni_gam_decomp",
#    # Update these values
#    n_iterations = 6000, # Number of iterations
#    n_burnin = 1000, # Number of iterations to discard at the beginning
#    n_thin = 4, # Reduces number of output samples to save memory and computation time
#    n_chains = 3 # Number of Markov chains
#  )

## ----plotnigamres, eval = FALSE,fig.align = 'center',fig.width = 7,fig.height = 5----
#  plot(res_ni_gam_decomp,
#    plot_type = "model_fit_plot",
#    plot_tide_gauge = FALSE
#  )

## ----plotnigamres_load, eval = TRUE,echo=FALSE, fig.align = 'center',out.width="100%"----
knitr::include_graphics("https://raw.githubusercontent.com/maeveupton/reslr/main/reslrvigplots/plotnigamres-1.png")

## ----plotnigamresrate, fig.align = 'center',fig.width = 7,fig.height = 5,eval = FALSE----
#  plot(res_ni_gam_decomp,
#    plot_type = "rate_plot"
#  )

## ----plotnigamresrate_load, eval = TRUE,echo=FALSE, fig.align = 'center',out.width="100%"----
knitr::include_graphics("https://raw.githubusercontent.com/maeveupton/reslr/main/docs/reslrvigplots/plotnigamresrate-1.png")

## ----totaldf, eval = FALSE----------------------------------------------------
#  total_model_fit_df <- res_ni_gam_decomp$output_dataframes$total_model_fit_df

## ----plotnigamregres, fig.align = 'center',fig.width = 7,fig.height = 5,eval = FALSE----
#  plot(res_ni_gam_decomp, plot_type = "regional_plot")

## ----plotnigamregres_load, eval = TRUE,echo=FALSE, fig.align = 'center',out.width="100%"----
knitr::include_graphics("https://raw.githubusercontent.com/maeveupton/reslr/main/reslrvigplots/plotnigamregres-1.png")

## ----regdf, eval = FALSE------------------------------------------------------
#  regional_component_df <- res_ni_gam_decomp$output_dataframes$regional_component_df

## ----plotnigamregresrate, fig.align = 'center',eval = FALSE-------------------
#  plot(res_ni_gam_decomp, plot_type = "regional_rate_plot")

## ----plotnigamregresrate_load, eval = TRUE,echo=FALSE, fig.align = 'center',out.width="100%"----
knitr::include_graphics("https://raw.githubusercontent.com/maeveupton/reslr/main/reslrvigplots/plotnigamregresrate-1.png")

## ----plotnigamlinres, fig.align = 'center',eval = FALSE-----------------------
#  plot(res_ni_gam_decomp, plot_type = "linear_local_plot")

## ----plotnigamlinres_load, eval = TRUE,echo=FALSE, fig.align = 'center',out.width="100%"----
url <- "https://raw.githubusercontent.com/maeveupton/reslr/main/reslrvigplots/plotnigamlinres-1.png"
knitr::include_graphics(url)

## ----linlocdf, eval = FALSE---------------------------------------------------
#  lin_loc_component_df <- res_ni_gam_decomp$output_dataframes$lin_loc_component_df

## ----linlocrate, eval = FALSE-------------------------------------------------
#  lin_loc_component_rates <- lin_loc_component_df %>%
#    dplyr::group_by(SiteName) %>%
#    dplyr::summarise(
#      linear_rate = unique(linear_rate),
#      linear_rate_err = unique(linear_rate_err)
#    )

## ----plotnigamnonlinres, fig.align = 'center',fig.width = 7,fig.height = 5,eval = FALSE----
#  plot(res_ni_gam_decomp, plot_type = "non_linear_local_plot")

## ----plotnigamnonlinres_load, eval = TRUE,echo=FALSE, fig.align = 'center',out.width="100%"----
url <- "https://raw.githubusercontent.com/maeveupton/reslr/main/reslrvigplots/plotnigamnonlinres-1.png"
knitr::include_graphics(url)

## ----non_lindf, eval = FALSE--------------------------------------------------
#  non_lin_loc_component_df <- res_ni_gam_decomp$output_dataframes$non_lin_loc_component_df

## ----plotnigamnonlinresrate, fig.align = 'center',fig.width = 7,fig.height = 5,eval = FALSE----
#  plot(res_ni_gam_decomp, plot_type = "non_linear_local_rate_plot")

## ----plotnigamnonlinresrate_load, eval = TRUE,echo=FALSE, fig.align = 'center',out.width="100%"----
url <- "https://raw.githubusercontent.com/maeveupton/reslr/main/reslrvigplots/plotnigamnonlinresrate-1.png"
knitr::include_graphics(url)

## ----plotnigamall, fig.align = 'center',fig.width = 7,fig.height = 5,eval = FALSE,message=FALSE,results='hide'----
#  plot(res_ni_gam_decomp, plot_type = "nigam_component_plot")

## ----plotnigamall_load, eval = TRUE,echo=FALSE, fig.align = 'center',out.width="100%"----
url <- "https://raw.githubusercontent.com/maeveupton/reslr/main/reslrvigplots/plotnigamall-1.png"
knitr::include_graphics(url)