integrate examples, demonstration, or even support for {errors} package

[dylanbeaudette]

Nearly automatic propagation of errors.

https://cran.r-project.org/web/packages/errors/vignettes/rjournal.pdf

Related to #20.

Some tinkering with SSURGO data by incorporating low, rv, and high values.

## probably need to do this via error prorogation / multivariate simulation

# https://cran.r-project.org/web/packages/errors/vignettes/rjournal.pdf

library(aqp)
library(soilDB)
library(lattice)
library(errors)

y <- fetchSDA(WHERE = "compname = 'clear lake' AND areasymbol != 'US'")

idx <- unique(y, vars = c(horizonDepths(y)))
y <- y[idx, ]

x <- trunc(y, z1 = 0, z2 = 30)



vars <- c('cokey', 'hzdept_r', 'hzdepb_r', 'om_l', 'om_r', 'om_h', 'total_frags_pct', 'dbthirdbar_l', 'dbthirdbar_r', 'dbthirdbar_h')

h <- horizons(x)[, vars]

# may not be required
# NULL -> 0
h$total_frags_pct <- ifelse(is.na(h$total_frags_pct), 0, h$total_frags_pct)

# remove NA
idx <- which(complete.cases(h))
h <- h[idx, ]

id <- h$cokey

thick <- h$hzdepb_r - h$hzdept_r

soc_l <- (h$om_l / 100) * 0.58
soc_h <- (h$om_h / 100) * 0.58
soc <- (h$om_r / 100) * 0.58

fraction <- (1 - (h$total_frags_pct / 100))

db_l <- h$dbthirdbar_l
db_h <- h$dbthirdbar_h
db <- h$dbthirdbar_r

# constant errors applied to thickness and fine earth fraction
errors(thick) <- 5
errors(fraction) <- 0.05

# constant error term: use boxplot notch concepts
# +/- 1.58 IQR/sqrt(n) ...  but that isn't quite right

# use 1/2 IQR for now
errors(soc) <- IQR(soc_h - soc_l, na.rm = TRUE) / 2
errors(db) <- IQR(db_h - db_l, na.rm = TRUE) / 2

# variable error term: average of (high - rv) and (rv - low)
errors(soc) <- ((soc_h - soc) + (soc - soc_l)) / 2
errors(db) <- ((db_h - db) + (db - db_l)) / 2

const <- 10
errors(const) <- 0


plot(soc, db)

stock <- thick * soc * fraction * db * 10

s <- tapply(stock, id, sum, simplify = FALSE)
s <- do.call('c', s)

# kg / m^2
s

format(median(s), digits = 2, notation = 'plus-minus')