1 Data

library(aqp)
library(soilDB)
library(sharpshootR)
library(cluster)
library(reshape2)
library(latticeExtra)
library(Gmedian)
library(MASS)
library(tactile)

# get some interesting data
x <- fetchKSSL('plainfield')

# remove partial-profile records
x <- subset(x, pedon_completeness_index > 3)
x <- subset(x, pedon_key != '6678')

# normalize taxonname for use as a grouping variable
x$taxonname <- tolower(x$taxonname)

# compute weighted-mean particle diameter (mm) for later
x$wmpd <- with(horizons(x), ((vcs * 1.5) + (cs * 0.75) + (ms * 0.375) + (fs * 0.175) + (vfs *0.075) + (silt * 0.026) + (clay * 0.0015)) / (vcs + cs + ms + fs + vfs + silt + clay))

2 Aggregation

# aggregation formula
fm <- taxonname ~ sand + silt + clay + estimated_ph_h2o + bs82 + wmpd

# L1 medians
z <- L1_profiles(x, fm = fm, method = 'constant', maxDepthConstant = 150)

# marginal medians
a <- slab(x, fm = fm)

2.1 Sanity Check

# check that sand + silt + clay sum to 100
(z$ssc_total <- z$sand + z$silt + z$clay)
##   [1] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
##  [24] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
##  [47] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
##  [70] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
##  [93] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
## [116] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
## [139] 100 100 100 100 100 100 100 100 100 100 100 100 100

2.2 Merge Marginal + L1

## TODO: this involves a lot of dumb re-naming
vars <- c("sand", "silt", "clay", "estimated_ph_h2o", "bs82", "wmpd")
z.long <- reshape2::melt(as(z, 'data.frame'), id.vars = c(idname(z)