In lvaudor/riverbed: Surfaces and volumes based on profiles.

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  echo = TRUE, message=FALSE, warning=FALSE
)

library(riverbed)
library(dplyr)
library(ggplot2)

Example datasets

data(s1)
data(s2)

Sources of uncertainty

Now suppose that the series $s_1$ and $s_2$ are measured with a certain level of uncertainty, due to

imprecision in the measure of height $z$ ($\sigma_z$) and/or
imprecision in the measure of longitudinal coordinate $l$ ($\sigma_l$).

These imprecisions might be different for the two series (change in sampling protocol, improvement in sampling gear between two dates, etc.).

The uncertainty in measures results in a certain amount of uncertainty in the estimate of area between the two curves. This is provided by the function area_uncertainty().

Here, with errors in the measures of height (0.1 and 0.3) higher than in the measures of longitudinal coordinates (0.05 and 0.2), and errors more important for series $s_2$ (0.3 and 0.2) than for $s_2$ (0.1 and 0.05).

Estimation of the uncertainty with area_between()

result_area <- area_between(s1,s2,
                            sigma_z=c(0.1,0.3),
                            sigma_l=c(0.05,0.2))

The plot_area() function can provide a visual hint of the uncertainties in $l$ and $z$ measures with horizontal and vertical error bars:

plot_area(result_area, show_uncertainty=TRUE)

A_est <- result_area$area
sigma_area <- round(result_area$sigma_area,2)
A_liminf <- round(A_est-1.96*A_est,2)
A_limsup <- round(A_est+1.96*A_est,2)

Here, for instance,

The estimate of area is r round(A_est,2).
The uncertainty in measures results in an error in the estimate of area of $\sigma_{area}$=r sigma_area. Hence we have a 95% confidence interval for the estimate of area of $[A-1.96\cdot\sigma_{area}, A+1.96\cdot\sigma_{area}]$=[r A_liminf,r A_limsup].

Check uncertainty through simulation

We can check the calculation of error through a simple simulation with 100 series $s_{1tmp}$ and $s_{2tmp}$ varying around $s_1$ and $s_2$ respectively, with variations corresponding to estimation errors $\sigma_z$ and $\sigma_l$:

set.seed(33)
sigma_z=c(0.1,0.3)
sigma_l=c(0.05,0.2)
f=function(i){
  s1_tmp <- tibble(l=s1$l+rnorm(nrow(s1),0,sigma_l[1]),
                   z=s1$z+rnorm(nrow(s1),0,sigma_z[1]))
  s2_tmp <- tibble(l=s2$l+rnorm(nrow(s2),0,sigma_l[2]),
                   z=s2$z+rnorm(nrow(s2),0,sigma_z[2]))
  return(area_between(s1_tmp,s2_tmp)$area)
}
area_vals=purrr::map_dbl(1:100,f)
sd(area_vals)

res=area_between(s1,s2,
                 sigma_z=sigma_z, sigma_l=sigma_l)
res$sigma_area