Section 8 Statistical Analysis: Stand Basal Area
In this section, we’ll evaluate the influence of the processing parameters on UAS-determined stand basal area. The accuracy UAS-derived stand basal area is impacted by the detection of individual trees (modeled in this section) and the estimated DBH of those detected trees (modeled in this section).
Even though there are concerns with the under-detection of smaller trees in the understory (Tinkham & Swayze 2021), the use of tree lists derived from remote sensing data (i.e. UAS SfM and Lidar) can still be useful for monitoring aspects of ecosystem function that are driven by large trees (Jeronimo et al. 2018).
Jeronimo et al. 2018 summarize this idea:
Large trees dominate basal area and biomass (Lutz et al. 2012; Bastin et al. 2015), carbon accumulation and growth (Sillett et al. 2010; Stephenson et al. 2014), and stand spatial heterogeneity (Lutz et al. 2013), and they provide unique structures that are keystone elements of many vertebrate species’ habitat requirements (Tews et al. 2004; North et al. 2017). Analysis of ITD results can capture the majority of variation in these processes even while omitting many smaller trees. For example, several studies have predicted biomass and volume directly from ITD results using allometric equations, explaining 67–93 percent of variation in even-aged pine plantations (Bortolot and Wynne 2005; Popescu 2007) and achieving error rates less than about 30 percent for Douglas-fir forests in various successional stages (Edson and Wing 2011). (p. 341)
We will utilize our “full” model presented here
8.1 Setup
load the data if needed
# load data if needed
if(ls()[ls() %in% "ptcld_validation_data"] %>% length()==0){
ptcld_validation_data = readr::read_csv("../data/ptcld_full_analysis_data.csv") %>%
dplyr::mutate(
depth_maps_generation_quality = factor(
depth_maps_generation_quality %>%
tolower() %>%
stringr::str_replace_all("ultrahigh", "ultra high")
, ordered = TRUE
, levels = c(
"lowest"
, "low"
, "medium"
, "high"
, "ultra high"
)
) %>% forcats::fct_rev()
, depth_maps_generation_filtering_mode = factor(
depth_maps_generation_filtering_mode %>% tolower()
, ordered = TRUE
, levels = c(
"disabled"
, "mild"
, "moderate"
, "aggressive"
)
) %>% forcats::fct_rev()
)
}load our plotting functions if needed (not showing these functions here but see the prior section for function definitions)
8.2 BA Percentage Error (bias)
8.2.1 Summary Statistics
# summarize data
dta_temp = ptcld_validation_data %>%
dplyr::group_by(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode) %>%
# collapse across study site
dplyr::summarise(
basal_area_pct_error = mean(basal_area_pct_error, na.rm = T)
, n = dplyr::n()
)
# set limits for color scale
lmt_basal_area_pct_error = ceiling(10*max(abs(range(dta_temp$basal_area_pct_error, na.rm = T))))/10
# scales::show_col(scales::pal_dichromat("BluetoOrange.10")(10))
# scales::show_col(scales::pal_div_gradient()(seq(0, 1, length.out = 7)))
# plot it
dta_temp %>%
ggplot(mapping = aes(
y = depth_maps_generation_quality
, x = depth_maps_generation_filtering_mode
, fill = basal_area_pct_error
, label = paste0(scales::percent(basal_area_pct_error,accuracy = 0.1), "\n(n=", n,")")
)) +
geom_tile(color = "white") +
geom_text(color = "white", size = 3) +
facet_grid(cols = vars(software)) +
scale_x_discrete(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0)) +
scale_fill_stepsn(
n.breaks = 7
, colors = scales::pal_div_gradient()(seq(0, 1, length.out = 7))
, limits = c(-lmt_basal_area_pct_error,lmt_basal_area_pct_error)
, labels = scales::percent_format(accuracy = 1)
, show.limits = T
) +
labs(
x = "filtering mode"
, y = "quality"
, fill = "Basal Area % Error"
, title = "mean basal area percentage error and # of study sites"
, subtitle = paste(
"negative values = UAS basal area < field basal area"
, " || "
, "positive values = UAS basal area > field basal area"
)
) +
theme_light() +
theme(
# legend.position = "none"
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
, panel.background = element_blank()
, panel.grid = element_blank()
# , plot.title = element_text(hjust = 0.5)
# , plot.subtitle = element_text(hjust = 0.5)
, strip.text = element_text(color = "black", face = "bold")
)
let’s check the distribution (of our dependent or \(y\) variable)
# distribution
ptcld_validation_data %>%
ggplot(mapping = aes(x = basal_area_pct_error)) +
geom_hline(yintercept = 0, color = "gray77") +
geom_vline(xintercept = 0, color = "gray77") +
# geom_vline(xintercept = c(0,1)) +
geom_density(fill = "lightblue", alpha = 0.7, color = NA) +
labs(y="",x="Basal Area % Error") +
scale_y_continuous(breaks = c(0)) +
scale_x_continuous(breaks = scales::extended_breaks(10), labels = scales::percent_format(accuracy = 1)) +
theme_light() +
theme(panel.grid = element_blank())
and the summary statistics
ptcld_validation_data %>%
dplyr::ungroup() %>%
dplyr::select(basal_area_pct_error) %>%
dplyr::summarise(
dplyr::across(
dplyr::everything()
, .fns = list(
mean = ~ mean(.x,na.rm=T), median = ~ median(.x,na.rm=T), sd = ~ sd(.x,na.rm=T)
, min = ~ min(.x,na.rm=T), max = ~ max(.x,na.rm=T)
, q25 = ~ quantile(.x, 0.25, na.rm = T)
, q75 = ~ quantile(.x, 0.75, na.rm = T)
)
, .names = "{.fn}"
)
) %>%
tidyr::pivot_longer(everything()) %>%
kableExtra::kbl(caption = "summary: `basal_area_pct_error`", digits = 3, col.names = NULL) %>%
kableExtra::kable_styling()| mean | -0.131 |
| median | -0.104 |
| sd | 0.338 |
| min | -1.000 |
| max | 0.519 |
| q25 | -0.275 |
| q75 | 0.066 |
8.2.2 Model: BA percentage error
# from Kurz:
gamma_a_b_from_omega_sigma = function(mode, sd) {
if (mode <= 0) stop("mode must be > 0")
if (sd <= 0) stop("sd must be > 0")
rate = (mode + sqrt(mode^2 + 4 * sd^2)) / (2 * sd^2)
shape = 1 + mode * rate
return(list(shape = shape, rate = rate))
}
mean_y_temp = mean(ptcld_validation_data$basal_area_pct_error, na.rm = T)
sd_y_temp = sd(ptcld_validation_data$basal_area_pct_error, na.rm = T)
omega_temp = sd_y_temp / 2
sigma_temp = 2 * sd_y_temp
s_r_temp = gamma_a_b_from_omega_sigma(mode = omega_temp, sd = sigma_temp)
stanvars_temp =
brms::stanvar(mean_y_temp, name = "mean_y") +
brms::stanvar(sd_y_temp, name = "sd_y") +
brms::stanvar(s_r_temp$shape, name = "alpha") +
brms::stanvar(s_r_temp$rate, name = "beta")8.2.2.1 Prior distributions
#### setting priors
# required libraries: tidyverse, tidybayes, brms, palettetown, latex2exp
priors_temp <- c(
brms::prior(normal(mean_y_temp, sd_y_temp * 5), class = "Intercept")
, brms::prior(gamma(s_r_temp$shape, s_r_temp$rate), class = "sd")
, brms::prior(cauchy(0, sd_y_temp), class = "sigma")
)
# plot
plt_prior11 =
priors_temp %>%
tidybayes::parse_dist() %>%
tidybayes::marginalize_lkjcorr(K = 2) %>%
tidyr::separate(
.args
, sep = ","
, into = c("a","b")
, remove = F
) %>%
dplyr::mutate(
distrib = dplyr::case_when(
!is.na(b) ~ paste0(
class, " ~ "
, .dist
, "("
, a %>% readr::parse_number() %>% round(2)
, ","
, b %>% readr::parse_number() %>% round(2)
, ")"
)
, T ~ paste0(
class, " ~ "
, .dist
, "("
, a %>% readr::parse_number() %>% round(2)
, ")"
)
)
) %>%
ggplot(., aes(dist = .dist, args = .args)) +
facet_grid(cols = vars(distrib), scales = "free") +
ggdist::stat_halfeye(
aes(fill = prior),
n = 10e2,
show.legend = F
, fill = "slategray"
) +
coord_flip() +
theme_light() +
theme(
strip.text = element_text(face = "bold", color = "black"),
axis.text.y = element_blank(),
axis.ticks = element_blank()
, axis.text = element_text(size = 6)
)+
labs(
x = ""
, title = "Priors: Basal Area Percentage Error"
, y = ""
)
plt_prior11
8.2.2.2 Fit the model
Now fit the model.
brms_ba_pe_mod = brms::brm(
formula = basal_area_pct_error ~
# baseline
1 +
# main effects
(1 | depth_maps_generation_quality) +
(1 | depth_maps_generation_filtering_mode) +
(1 | software) +
(1 | study_site) + # only fitting main effects of site and not interactions
# two-way interactions
(1 | depth_maps_generation_quality:depth_maps_generation_filtering_mode) +
(1 | depth_maps_generation_quality:software) +
(1 | depth_maps_generation_filtering_mode:software) +
# three-way interactions
(1 | depth_maps_generation_quality:depth_maps_generation_filtering_mode:software)
, data = ptcld_validation_data
, family = brms::brmsfamily(family = "gaussian")
, iter = 20000, warmup = 10000, chains = 4
, control = list(adapt_delta = 0.999, max_treedepth = 13)
, cores = round(parallel::detectCores()/2)
, prior = c(
brms::prior(normal(mean_y, sd_y * 5), class = "Intercept")
, brms::prior(gamma(alpha, beta), class = "sd")
, brms::prior(cauchy(0, sd_y), class = "sigma")
)
, stanvars = stanvars_temp
, file = paste0(rootdir, "/fits/brms_ba_pe_mod")
)8.2.2.3 Quality:Filtering - interaction
draws_temp =
ptcld_validation_data %>%
dplyr::distinct(depth_maps_generation_quality, depth_maps_generation_filtering_mode) %>%
tidybayes::add_epred_draws(
brms_ba_pe_mod, allow_new_levels = T
# this part is crucial
, re_formula = ~ (1 | depth_maps_generation_quality) +
(1 | depth_maps_generation_filtering_mode) +
(1 | depth_maps_generation_quality:depth_maps_generation_filtering_mode)
) %>%
dplyr::rename(value = .epred) %>%
dplyr::mutate(med = tidybayes::median_hdci(value)$y)
# plot
draws_temp %>%
dplyr::mutate(depth_maps_generation_quality = depth_maps_generation_quality %>% forcats::fct_rev()) %>%
# plot
ggplot(
mapping = aes(
y = value, x = depth_maps_generation_filtering_mode
, fill = med
)
) +
geom_hline(yintercept = 0, color = "gray33") +
tidybayes::stat_eye(
point_interval = median_hdi, .width = 0.95
, slab_alpha = 0.98
, interval_color = "black", linewidth = 1
, point_color = "black", point_fill = "black", point_size = 1
) +
scale_fill_stepsn(
n.breaks = 7
, colors = scales::pal_div_gradient()(seq(0, 1, length.out = 7))
, limits = c(-lmt_basal_area_pct_error*.5,lmt_basal_area_pct_error*.5)
, labels = scales::percent_format(accuracy = 1)
, show.limits = T
) +
facet_grid(cols = vars(depth_maps_generation_quality)) +
scale_y_continuous(
labels = scales::percent_format(accuracy = 1)
, limits = c(-lmt_basal_area_pct_error*1.2,lmt_basal_area_pct_error*1.2)
) +
labs(
x = "filtering mode", y = "Basal Area % Error"
, subtitle = "posterior predictive distribution of basal area percentage error with 95% HDI\nby dense cloud quality"
) +
theme_light() +
theme(
legend.position = "none"
, legend.direction = "horizontal"
, axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
, strip.text = element_text(color = "black", face = "bold")
) 
and a table of these 95% HDI values
table_temp =
draws_temp %>%
tidybayes::median_hdi(value) %>%
dplyr::select(-c(.point,.interval, .width,.row)) %>%
dplyr::arrange(depth_maps_generation_quality,depth_maps_generation_filtering_mode) %>%
dplyr::mutate(dplyr::across(.cols = is.numeric, .fns = ~ scales::percent(.x,accuracy = 0.1)))
table_temp %>%
dplyr::select(-c(depth_maps_generation_quality)) %>%
kableExtra::kbl(
digits = 2
, caption = "Basal Area % Error<br>95% HDI of the posterior predictive distribution"
, col.names = c(
"filtering mode"
, "Basal Area % Error<br>median"
, "HDI low", "HDI high"
)
, escape = F
) %>%
kableExtra::kable_styling() %>%
kableExtra::pack_rows(index = table(forcats::fct_inorder(table_temp$depth_maps_generation_quality))) %>%
kableExtra::scroll_box(height = "8in")| filtering mode |
Basal Area % Error median |
HDI low | HDI high |
|---|---|---|---|
| ultra high | |||
| aggressive | -1.8% | -54.6% | 52.8% |
| moderate | 0.9% | -51.2% | 54.3% |
| mild | 3.3% | -49.2% | 56.1% |
| disabled | 3.1% | -49.6% | 56.3% |
| high | |||
| aggressive | -5.9% | -57.4% | 48.5% |
| moderate | -3.2% | -54.5% | 49.6% |
| mild | -0.9% | -52.9% | 51.2% |
| disabled | -0.1% | -51.5% | 52.3% |
| medium | |||
| aggressive | -16.4% | -71.7% | 35.3% |
| moderate | -12.1% | -63.7% | 41.5% |
| mild | -9.5% | -61.4% | 42.7% |
| disabled | -8.5% | -62.4% | 42.5% |
| low | |||
| aggressive | -25.3% | -79.0% | 27.1% |
| moderate | -20.7% | -72.6% | 32.2% |
| mild | -16.9% | -66.7% | 37.6% |
| disabled | -15.5% | -66.9% | 37.2% |
| lowest | |||
| aggressive | -37.2% | -89.6% | 20.2% |
| moderate | -34.1% | -86.5% | 22.4% |
| mild | -30.3% | -83.9% | 25.0% |
| disabled | -29.5% | -81.9% | 27.5% |
8.2.2.4 Software:Quality - interaction
draws_temp =
ptcld_validation_data %>%
dplyr::distinct(depth_maps_generation_quality, software) %>%
tidybayes::add_epred_draws(
brms_ba_pe_mod, allow_new_levels = T
# this part is crucial
, re_formula = ~ (1 | depth_maps_generation_quality) +
(1 | software) +
(1 | depth_maps_generation_quality:software)
) %>%
dplyr::rename(value = .epred) %>%
dplyr::mutate(med = tidybayes::median_hdci(value)$y)
# plot
draws_temp %>%
dplyr::mutate(depth_maps_generation_quality = depth_maps_generation_quality %>% forcats::fct_rev()) %>%
# plot
ggplot(
mapping = aes(
y = value, x = software
, fill = med
)
) +
geom_hline(yintercept = 0, color = "gray33") +
tidybayes::stat_eye(
point_interval = median_hdi, .width = .95
, slab_alpha = 0.98
, interval_color = "black", linewidth = 1
, point_color = "black", point_fill = "black", point_size = 1
) +
scale_fill_stepsn(
n.breaks = 7
, colors = scales::pal_div_gradient()(seq(0, 1, length.out = 7))
, limits = c(-lmt_basal_area_pct_error*.5,lmt_basal_area_pct_error*.5)
, labels = scales::percent_format(accuracy = 1)
, show.limits = T
) +
facet_grid(cols = vars(depth_maps_generation_quality)) +
scale_y_continuous(
labels = scales::percent_format(accuracy = 1)
, limits = c(-lmt_basal_area_pct_error*1.2,lmt_basal_area_pct_error*1.2)
) +
labs(
x = "software", y = "Basal Area % Error"
, subtitle = "posterior predictive distribution of basal area percentage error with 95% HDI\nby dense cloud quality"
) +
theme_light() +
theme(
legend.position = "none"
, legend.direction = "horizontal"
, axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
, strip.text = element_text(color = "black", face = "bold")
) 
and a table of these 95% HDI values
table_temp =
draws_temp %>%
tidybayes::median_hdi(value) %>%
# remove out-of-sample obs
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(depth_maps_generation_quality, software)
, by = dplyr::join_by(depth_maps_generation_quality, software)
) %>%
dplyr::select(-c(.point,.interval, .width,.row)) %>%
dplyr::arrange(software,depth_maps_generation_quality) %>%
dplyr::mutate(dplyr::across(.cols = is.numeric, .fns = ~ scales::percent(.x,accuracy = 0.1)))
table_temp %>%
dplyr::select(-c(software)) %>%
kableExtra::kbl(
digits = 2
, caption = "Basal Area % Error<br>95% HDI of the posterior predictive distribution"
, col.names = c(
"quality"
, "Basal Area % Error<br>median"
, "HDI low", "HDI high"
)
, escape = F
) %>%
kableExtra::kable_styling() %>%
kableExtra::pack_rows(index = table(forcats::fct_inorder(table_temp$software))) %>%
kableExtra::scroll_box(height = "8in")| quality |
Basal Area % Error median |
HDI low | HDI high |
|---|---|---|---|
| METASHAPE | |||
| ultra high | 7.5% | -30.6% | 46.2% |
| high | 3.7% | -35.5% | 41.8% |
| medium | -5.1% | -43.7% | 33.1% |
| low | -23.4% | -62.0% | 15.1% |
| lowest | -50.7% | -89.0% | -12.0% |
| OPENDRONEMAP | |||
| ultra high | -1.7% | -39.9% | 37.9% |
| high | -10.4% | -48.9% | 29.0% |
| medium | -31.7% | -70.2% | 7.6% |
| low | -32.4% | -71.1% | 6.8% |
| lowest | -33.8% | -71.8% | 5.8% |
| PIX4D | |||
| ultra high | 7.4% | -31.9% | 47.1% |
| high | 5.9% | -33.8% | 45.1% |
| medium | 2.6% | -36.0% | 43.4% |
| low | -6.9% | -46.6% | 33.4% |
8.2.2.5 Software:Filtering - interaction
draws_temp =
ptcld_validation_data %>%
dplyr::distinct(depth_maps_generation_filtering_mode, software) %>%
tidybayes::add_epred_draws(
brms_ba_pe_mod, allow_new_levels = T
# this part is crucial
, re_formula = ~ (1 | depth_maps_generation_filtering_mode) +
(1 | software) +
(1 | depth_maps_generation_filtering_mode:software)
) %>%
dplyr::rename(value = .epred) %>%
dplyr::mutate(med = tidybayes::median_hdci(value)$y)
# plot
draws_temp %>%
# plot
ggplot(
mapping = aes(
y = value, x = software
, fill = med
)
) +
geom_hline(yintercept = 0, color = "gray33") +
tidybayes::stat_eye(
point_interval = median_hdi, .width = .95
, slab_alpha = 0.98
, interval_color = "black", linewidth = 1
, point_color = "black", point_fill = "black", point_size = 1
) +
scale_fill_stepsn(
n.breaks = 7
, colors = scales::pal_div_gradient()(seq(0, 1, length.out = 7))
, limits = c(-lmt_basal_area_pct_error*.5,lmt_basal_area_pct_error*.5)
, labels = scales::percent_format(accuracy = 1)
, show.limits = T
) +
facet_grid(cols = vars(depth_maps_generation_filtering_mode)) +
scale_y_continuous(
labels = scales::percent_format(accuracy = 1)
, limits = c(-lmt_basal_area_pct_error*1.2,lmt_basal_area_pct_error*1.2)
) +
labs(
x = "software", y = "Basal Area % Error"
, subtitle = "posterior predictive distribution of basal area percentage error with 95% HDI\nby filtering mode"
) +
theme_light() +
theme(
legend.position = "none"
, legend.direction = "horizontal"
, axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
, strip.text = element_text(color = "black", face = "bold")
) 
and a table of these 95% HDI values
table_temp =
draws_temp %>%
tidybayes::median_hdi(value) %>%
dplyr::select(-c(.point,.interval, .width,.row)) %>%
dplyr::arrange(software,depth_maps_generation_filtering_mode) %>%
dplyr::mutate(dplyr::across(.cols = is.numeric, .fns = ~ scales::percent(.x,accuracy = 0.1)))
table_temp %>%
dplyr::select(-c(software)) %>%
kableExtra::kbl(
digits = 2
, caption = "Basal Area % Error<br>95% HDI of the posterior predictive distribution"
, col.names = c(
"filtering mode"
, "Basal Area % Error<br>median"
, "HDI low", "HDI high"
)
, escape = F
) %>%
kableExtra::kable_styling() %>%
kableExtra::pack_rows(index = table(forcats::fct_inorder(table_temp$software))) %>%
kableExtra::scroll_box(height = "8in")| filtering mode |
Basal Area % Error median |
HDI low | HDI high |
|---|---|---|---|
| METASHAPE | |||
| aggressive | -31.0% | -73.2% | 11.4% |
| moderate | -15.4% | -57.2% | 27.5% |
| mild | -5.2% | -46.6% | 37.7% |
| disabled | -2.9% | -44.4% | 40.1% |
| OPENDRONEMAP | |||
| aggressive | -19.3% | -61.5% | 22.8% |
| moderate | -22.1% | -65.7% | 18.6% |
| mild | -22.8% | -65.7% | 19.2% |
| disabled | -22.8% | -68.1% | 16.6% |
| PIX4D | |||
| moderate | -5.3% | -49.5% | 36.7% |
| mild | -0.5% | -43.9% | 42.1% |
| disabled | 1.5% | -42.6% | 43.4% |
8.2.2.6 Software:Quality:Filtering - interaction
# get draws
fltr_sftwr_draws_temp =
tidyr::crossing(
depth_maps_generation_quality = unique(ptcld_validation_data$depth_maps_generation_quality)
, depth_maps_generation_filtering_mode = unique(ptcld_validation_data$depth_maps_generation_filtering_mode)
, software = unique(ptcld_validation_data$software)
) %>%
tidybayes::add_epred_draws(
brms_ba_pe_mod, allow_new_levels = T
# this part is crucial
, re_formula = ~ (1 | depth_maps_generation_quality) + # main effects
(1 | depth_maps_generation_quality) +
(1 | depth_maps_generation_filtering_mode) +
(1 | software) +
# two-way interactions
(1 | depth_maps_generation_quality:depth_maps_generation_filtering_mode) +
(1 | depth_maps_generation_quality:software) +
(1 | depth_maps_generation_filtering_mode:software) +
# three-way interactions
(1 | depth_maps_generation_quality:depth_maps_generation_filtering_mode:software)
) %>%
dplyr::rename(value = .epred) %>%
dplyr::mutate(med = tidybayes::median_hdci(value)$y)
# plot
qlty_fltr_sftwr_ba_pe =
fltr_sftwr_draws_temp %>%
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode)
, by = dplyr::join_by(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode)
) %>%
dplyr::mutate(depth_maps_generation_quality = depth_maps_generation_quality %>% forcats::fct_rev()) %>%
ggplot(
mapping = aes(
y = value
, x = depth_maps_generation_filtering_mode
, fill = med
)
) +
geom_hline(yintercept = 0, color = "gray33") +
tidybayes::stat_eye(
point_interval = median_hdi, .width = .95
, slab_alpha = 0.98
, interval_color = "black", linewidth = 1
, point_color = "black", point_fill = "black", point_size = 1
) +
scale_fill_stepsn(
n.breaks = 7
, colors = scales::pal_div_gradient()(seq(0, 1, length.out = 7))
, limits = c(-lmt_basal_area_pct_error*.5,lmt_basal_area_pct_error*.5)
, labels = scales::percent_format(accuracy = 1)
, show.limits = T
) +
facet_grid(
rows = vars(software)
, cols = vars(depth_maps_generation_quality)
# , switch = "y"
) +
scale_y_continuous(
labels = scales::percent_format(accuracy = 1)
, limits = c(-lmt_basal_area_pct_error*1.5,lmt_basal_area_pct_error*1.2)
) +
labs(
x = "filtering mode", y = "Basal Area % Error"
, subtitle = "posterior predictive distribution of basal area percentage error with 95% HDI\nby dense cloud quality and software"
# , caption = form_temp
) +
theme_light() +
theme(
legend.position = "none"
, legend.direction = "horizontal"
, axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
, strip.text = element_text(color = "black", face = "bold")
, panel.grid = element_blank()
# , strip.placement = "outside"
) +
guides(
fill = guide_legend(override.aes = list(shape = NA, size = 6, alpha = 0.9, lwd = NA))
)
# print it
qlty_fltr_sftwr_ba_pe
and a table of these 95% HDI values
table_temp =
fltr_sftwr_draws_temp %>%
tidybayes::median_hdi(value) %>%
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode)
, by = dplyr::join_by(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode)
) %>%
dplyr::mutate(depth_maps_generation_quality = depth_maps_generation_quality %>% forcats::fct_rev()) %>%
dplyr::select(c(
software, depth_maps_generation_quality, depth_maps_generation_filtering_mode
, value, .lower, .upper
)) %>%
dplyr::ungroup() %>%
dplyr::arrange(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode) %>%
dplyr::mutate(dplyr::across(.cols = is.numeric, .fns = ~ scales::percent(.x,accuracy = 0.1)))
table_temp %>%
# dplyr::select(-c(software)) %>%
kableExtra::kbl(
digits = 2
, caption = "Basal Area % Error<br>95% HDI of the posterior predictive distribution"
, col.names = c(
"software", "quality", "filtering mode"
, "Basal Area % Error<br>median"
, "HDI low", "HDI high"
)
, escape = F
) %>%
kableExtra::kable_styling() %>%
# kableExtra::pack_rows(index = table(forcats::fct_inorder(table_temp$software))) %>%
kableExtra::collapse_rows(columns = 1:2, valign = "top") %>%
kableExtra::scroll_box(height = "8in")| software | quality | filtering mode |
Basal Area % Error median |
HDI low | HDI high |
|---|---|---|---|---|---|
| METASHAPE | lowest | aggressive | -69.1% | -103.6% | -34.5% |
| moderate | -53.9% | -88.2% | -19.4% | ||
| mild | -41.9% | -78.0% | -8.9% | ||
| disabled | -39.6% | -74.8% | -5.8% | ||
| low | aggressive | -43.5% | -78.5% | -9.5% | |
| moderate | -25.9% | -61.0% | 7.8% | ||
| mild | -14.2% | -49.6% | 19.5% | ||
| disabled | -11.1% | -44.7% | 24.3% | ||
| medium | aggressive | -23.5% | -58.7% | 10.1% | |
| moderate | -6.5% | -40.5% | 28.5% | ||
| mild | 3.4% | -30.6% | 38.4% | ||
| disabled | 6.0% | -27.3% | 41.6% | ||
| high | aggressive | -12.7% | -47.1% | 21.9% | |
| moderate | 2.1% | -30.5% | 38.5% | ||
| mild | 11.5% | -22.7% | 45.9% | ||
| disabled | 13.9% | -21.0% | 48.3% | ||
| ultra high | aggressive | -8.3% | -42.9% | 26.1% | |
| moderate | 6.0% | -28.4% | 40.8% | ||
| mild | 15.6% | -19.1% | 49.7% | ||
| disabled | 16.7% | -18.3% | 50.5% | ||
| OPENDRONEMAP | lowest | aggressive | -33.0% | -67.4% | 1.2% |
| moderate | -35.7% | -70.3% | -1.7% | ||
| mild | -36.1% | -70.5% | -1.6% | ||
| disabled | -35.8% | -70.7% | -1.9% | ||
| low | aggressive | -32.9% | -67.4% | 2.1% | |
| moderate | -34.4% | -69.1% | -0.2% | ||
| mild | -34.3% | -68.0% | 0.8% | ||
| disabled | -33.9% | -68.5% | 0.5% | ||
| medium | aggressive | -31.3% | -65.9% | 3.1% | |
| moderate | -33.0% | -67.3% | 1.7% | ||
| mild | -34.1% | -69.5% | -0.3% | ||
| disabled | -34.0% | -68.6% | 0.4% | ||
| high | aggressive | -8.2% | -43.7% | 25.2% | |
| moderate | -12.0% | -47.1% | 21.6% | ||
| mild | -13.1% | -46.3% | 22.3% | ||
| disabled | -13.0% | -48.0% | 21.2% | ||
| ultra high | aggressive | 0.9% | -33.5% | 35.4% | |
| moderate | -2.8% | -36.7% | 32.5% | ||
| mild | -3.7% | -39.5% | 29.3% | ||
| disabled | -5.4% | -40.0% | 29.0% | ||
| PIX4D | low | moderate | -8.5% | -43.6% | 25.7% |
| mild | -2.3% | -38.1% | 31.6% | ||
| disabled | 0.7% | -33.6% | 36.0% | ||
| medium | moderate | 2.3% | -32.1% | 37.2% | |
| mild | 6.9% | -28.0% | 41.4% | ||
| disabled | 8.9% | -24.5% | 44.8% | ||
| high | moderate | 5.3% | -29.5% | 39.2% | |
| mild | 9.4% | -24.9% | 44.6% | ||
| disabled | 11.4% | -22.3% | 46.8% | ||
| ultra high | moderate | 6.8% | -26.9% | 42.2% | |
| mild | 11.4% | -23.7% | 45.7% | ||
| disabled | 12.1% | -22.8% | 46.3% |
8.3 BA Absolute Percentage Error (precision)
8.3.1 Summary Statistics
# summarize data
dta_temp = ptcld_validation_data %>%
dplyr::group_by(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode) %>%
# collapse across study site
dplyr::summarise(
basal_area_abs_pct_error = mean(basal_area_abs_pct_error, na.rm = T)
, n = dplyr::n()
)
# set limits for color scale
lmt_basal_area_abs_pct_error = ceiling(1.02*10*max(abs(range(dta_temp$basal_area_abs_pct_error, na.rm = T))))/10
# scales::show_col(viridis::mako(n = 10, begin = 0.2, end = 0.9, direction = -1))
# scales::show_col(scales::pal_div_gradient()(seq(0, 1, length.out = 7)))
# plot it
dta_temp %>%
ggplot(mapping = aes(
y = depth_maps_generation_quality
, x = depth_maps_generation_filtering_mode
, fill = basal_area_abs_pct_error
, label = paste0(scales::percent(basal_area_abs_pct_error,accuracy = 0.1), "\n(n=", n,")")
)) +
geom_tile(color = "white") +
geom_text(color = "white", size = 3) +
facet_grid(cols = vars(software)) +
scale_x_discrete(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0)) +
scale_fill_stepsn(
n.breaks = 5
, colors = viridis::mako(n = 5, begin = 0.2, end = 0.9, direction = -1)
, limits = c(0,lmt_basal_area_abs_pct_error)
, labels = scales::percent_format(accuracy = 1)
, show.limits = T
) +
labs(
x = "filtering mode"
, y = "quality"
, fill = "Basal Area Abs. % Error"
, title = "mean basal area absolute percentage error and # of study sites"
) +
theme_light() +
theme(
legend.position = "none"
, axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
, panel.background = element_blank()
, panel.grid = element_blank()
# , plot.title = element_text(hjust = 0.5)
# , plot.subtitle = element_text(hjust = 0.5)
, strip.text = element_text(color = "black", face = "bold")
)
let’s check the distribution (of our dependent or \(y\) variable)
# distribution
ptcld_validation_data %>%
ggplot(mapping = aes(x = basal_area_abs_pct_error)) +
geom_hline(yintercept = 0, color = "gray77") +
geom_vline(xintercept = 0, color = "gray77") +
# geom_vline(xintercept = c(0,1)) +
geom_density(fill = "lightblue", alpha = 0.7, color = NA) +
labs(y="",x="Basal Area Abs. % Error") +
scale_y_continuous(breaks = c(0)) +
scale_x_continuous(breaks = scales::extended_breaks(10), labels = scales::percent_format(accuracy = 1)) +
theme_light() +
theme(panel.grid = element_blank())
and the summary statistics
ptcld_validation_data %>%
dplyr::ungroup() %>%
dplyr::select(basal_area_abs_pct_error) %>%
dplyr::summarise(
dplyr::across(
dplyr::everything()
, .fns = list(
mean = ~ mean(.x,na.rm=T), median = ~ median(.x,na.rm=T), sd = ~ sd(.x,na.rm=T)
, min = ~ min(.x,na.rm=T), max = ~ max(.x,na.rm=T)
, q25 = ~ quantile(.x, 0.25, na.rm = T)
, q75 = ~ quantile(.x, 0.75, na.rm = T)
)
, .names = "{.fn}"
)
) %>%
tidyr::pivot_longer(everything()) %>%
kableExtra::kbl(caption = "summary: `basal_area_abs_pct_error`", digits = 3, col.names = NULL) %>%
kableExtra::kable_styling()| mean | 0.263 |
| median | 0.190 |
| sd | 0.248 |
| min | 0.003 |
| max | 1.000 |
| q25 | 0.088 |
| q75 | 0.359 |
8.3.2 Model: BA Absolute percentage error
We’ll model the BA APE using the gamma likelihood to accurately represent the dependent variable which is continuous and strictly positive (i.e. it is impossible to have a negative APE). The gamma distribution is a great alternative that accounts for data with a zero lower limit and any right skew.
We borrow here from the excellent series on causal inference by A. Solomon Kurz
8.3.2.1 Prior distributions
First, let’s set the priors
# desired priors
m_temp <- mean(ptcld_validation_data$basal_area_abs_pct_error) # desired mean
s_temp <- sd(ptcld_validation_data$basal_area_abs_pct_error) # desired SD
# use the equations to get the lognormal parameter values
mu_temp <- log(m_temp / sqrt((s_temp*5)^2 / m_temp^2 + 1))
sigma_temp <- sqrt(log((s_temp*5)^2 / m_temp^2 + 1))
# what are the lognormal parameter values?
# mu_temp
# sigma_temp
# exp(mu_temp)
stanvars_temp =
brms::stanvar(mu_temp, name = "mu_temp") +
brms::stanvar(sigma_temp, name = "sigma_temp")
#### setting priors
# required libraries: tidyverse, tidybayes, brms, palettetown, latex2exp
brms_ba_ape_mod_priors_temp <- c(
brms::prior(normal(mu_temp, sigma_temp), class = "Intercept")
, brms::prior(student_t(3, 0, 2.5), class = "sd")
, brms::prior(gamma(0.1, 0.1), class = shape)
# , brms::prior(exponential(0.5), class = shape)
)plot the priors
# plot
plt_prior12 =
brms_ba_ape_mod_priors_temp %>%
tidybayes::parse_dist() %>%
tidybayes::marginalize_lkjcorr(K = 2) %>%
tidyr::separate(
.args
, sep = ","
, into = c("a","b")
, remove = F
) %>%
dplyr::mutate(
distrib = dplyr::case_when(
!is.na(b) ~ paste0(
class, " ~ "
, .dist
, "("
, a %>% readr::parse_number() %>% round(2)
, ","
, b %>% readr::parse_number() %>% round(2)
, ")"
)
, T ~ paste0(
class, " ~ "
, .dist
, "("
, a %>% readr::parse_number() %>% round(2)
, ")"
)
)
) %>%
ggplot(., aes(dist = .dist, args = .args)) +
facet_grid(cols = vars(distrib), scales = "free") +
ggdist::stat_halfeye(
aes(fill = prior),
n = 10e2,
show.legend = F
, fill = "slategray"
) +
coord_flip() +
theme_light() +
theme(
strip.text = element_text(face = "bold", color = "black"),
axis.text.y = element_blank(),
axis.ticks = element_blank()
, axis.text = element_text(size = 6)
)+
labs(
x = ""
, title = "Priors: Basal Area Absolute Percentage Error"
, y = ""
)
plt_prior12
8.3.2.2 Fit the model
Now fit the model.
#### Fit the model
brms_ba_ape_mod = brms::brm(
formula = basal_area_abs_pct_error ~
# baseline
1 +
# main effects
(1 | depth_maps_generation_quality) +
(1 | depth_maps_generation_filtering_mode) +
(1 | software) +
(1 | study_site) + # only fitting main effects of site and not interactions
# two-way interactions
(1 | depth_maps_generation_quality:depth_maps_generation_filtering_mode) +
(1 | depth_maps_generation_quality:software) +
(1 | depth_maps_generation_filtering_mode:software) +
# three-way interactions
(1 | depth_maps_generation_quality:depth_maps_generation_filtering_mode:software)
, data = ptcld_validation_data
# , family = brms::brmsfamily("Gamma", link = "log")
, family = brms::brmsfamily("weibull", link = "log")
# priors
, prior = brms_ba_ape_mod_priors_temp
, stanvars = stanvars_temp
# mcmc
, iter = 20000, warmup = 10000, chains = 4
, control = list(adapt_delta = 0.999, max_treedepth = 13)
, cores = round(parallel::detectCores()/2)
, file = paste0(rootdir, "/fits/brms_ba_ape_mod")
)8.3.2.3 Quality:Filtering - interaction
draws_temp =
ptcld_validation_data %>%
dplyr::distinct(depth_maps_generation_quality, depth_maps_generation_filtering_mode) %>%
tidybayes::add_epred_draws(
brms_ba_ape_mod, allow_new_levels = T
# this part is crucial
, re_formula = ~ (1 | depth_maps_generation_quality) +
(1 | depth_maps_generation_filtering_mode) +
(1 | depth_maps_generation_quality:depth_maps_generation_filtering_mode)
) %>%
dplyr::rename(value = .epred) %>%
dplyr::mutate(med = tidybayes::median_hdci(value)$y)
# plot
draws_temp %>%
dplyr::mutate(depth_maps_generation_quality = depth_maps_generation_quality %>% forcats::fct_rev()) %>%
# plot
ggplot(
mapping = aes(
y = value, x = depth_maps_generation_filtering_mode
, fill = med
)
) +
geom_hline(yintercept = 0, color = "gray33") +
tidybayes::stat_eye(
point_interval = median_hdi, .width = .95
, slab_alpha = 0.98
, interval_color = "black", linewidth = 1
, point_color = "black", point_fill = "black", point_size = 1
) +
scale_fill_stepsn(
n.breaks = 5
, colors = viridis::mako(n = 5, begin = 0.2, end = 0.9, direction = -1)
, limits = c(0,lmt_basal_area_abs_pct_error*.6)
, labels = scales::percent_format(accuracy = 0.1)
, show.limits = T
) +
facet_grid(cols = vars(depth_maps_generation_quality)) +
scale_y_continuous(
limits = c(-0.005,lmt_basal_area_abs_pct_error*1.1)
, labels = scales::percent_format(accuracy = 1)
) +
labs(
x = "filtering mode", y = "BA Abs. % Error"
, subtitle = "posterior predictive distribution of BA Abs. % Error with 95% HDI\nby dense cloud quality"
) +
theme_light() +
theme(
legend.position = "none"
, legend.direction = "horizontal"
, axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
, strip.text = element_text(color = "black", face = "bold")
) 
and a table of these 95% HDI values
table_temp =
draws_temp %>%
tidybayes::median_hdi(value) %>%
dplyr::select(-c(.point,.interval, .width,.row)) %>%
dplyr::arrange(depth_maps_generation_quality,depth_maps_generation_filtering_mode) %>%
dplyr::mutate(dplyr::across(.cols = is.numeric, .fns = ~ scales::percent(.x,accuracy = 0.1)))
table_temp %>%
dplyr::select(-c(depth_maps_generation_quality)) %>%
kableExtra::kbl(
digits = 2
, caption = "BA Abs. % Error<br>95% HDI of the posterior predictive distribution"
, col.names = c(
"filtering mode"
, "BA Abs. % Error<br>median"
, "HDI low", "HDI high"
)
, escape = F
) %>%
kableExtra::kable_styling() %>%
kableExtra::pack_rows(index = table(forcats::fct_inorder(table_temp$depth_maps_generation_quality))) %>%
kableExtra::scroll_box(height = "8in")| filtering mode |
BA Abs. % Error median |
HDI low | HDI high |
|---|---|---|---|
| ultra high | |||
| aggressive | 15.8% | 1.9% | 36.2% |
| moderate | 14.8% | 1.8% | 33.1% |
| mild | 14.4% | 1.6% | 31.9% |
| disabled | 13.3% | 1.7% | 29.6% |
| high | |||
| aggressive | 17.9% | 2.0% | 39.8% |
| moderate | 16.6% | 2.2% | 36.3% |
| mild | 16.1% | 2.1% | 35.1% |
| disabled | 15.3% | 2.0% | 33.9% |
| medium | |||
| aggressive | 23.5% | 2.9% | 52.3% |
| moderate | 20.7% | 2.4% | 44.4% |
| mild | 19.8% | 2.8% | 43.0% |
| disabled | 18.7% | 2.5% | 40.9% |
| low | |||
| aggressive | 29.3% | 3.7% | 65.2% |
| moderate | 26.1% | 3.4% | 56.6% |
| mild | 23.9% | 3.2% | 52.1% |
| disabled | 22.2% | 3.2% | 49.0% |
| lowest | |||
| aggressive | 37.5% | 3.9% | 86.5% |
| moderate | 34.6% | 3.6% | 78.9% |
| mild | 32.2% | 3.6% | 73.1% |
| disabled | 30.6% | 2.7% | 69.0% |
we can also make pairwise comparisons
# first we need to define the contrasts to make
contrast_list =
tidyr::crossing(
x1 = unique(ptcld_validation_data$depth_maps_generation_quality)
, x2 = unique(ptcld_validation_data$depth_maps_generation_quality)
) %>%
dplyr::mutate(
dplyr::across(
dplyr::everything()
, .fns = function(x){factor(
x, ordered = T
, levels = levels(ptcld_validation_data$depth_maps_generation_quality)
)}
)
) %>%
dplyr::filter(x1<x2) %>%
dplyr::arrange(x1,x2) %>%
dplyr::mutate(dplyr::across(dplyr::everything(), as.character)) %>%
purrr::transpose()
# make the contrasts using compare_levels
brms_contrast_temp = draws_temp %>%
tidybayes::compare_levels(
value
, by = depth_maps_generation_quality
, comparison =
contrast_list
# tidybayes::emmeans_comparison("revpairwise")
#"pairwise"
) %>%
dplyr::rename(contrast = depth_maps_generation_quality)
# separate contrast
brms_contrast_temp = brms_contrast_temp %>%
dplyr::ungroup() %>%
tidyr::separate_wider_delim(
cols = contrast
, delim = " - "
, names = paste0(
"sorter"
, 1:(max(stringr::str_count(brms_contrast_temp$contrast, "-"))+1)
)
, too_few = "align_start"
, cols_remove = F
) %>%
dplyr::filter(sorter1!=sorter2) %>%
dplyr::mutate(
dplyr::across(
tidyselect::starts_with("sorter")
, .fns = function(x){factor(
x, ordered = T
, levels = levels(ptcld_validation_data$depth_maps_generation_quality)
)}
)
, contrast = contrast %>%
forcats::fct_reorder(
paste0(as.numeric(sorter1), as.numeric(sorter2)) %>%
as.numeric()
)
, depth_maps_generation_filtering_mode = depth_maps_generation_filtering_mode %>%
factor(
levels = levels(ptcld_validation_data$depth_maps_generation_filtering_mode)
, ordered = T
)
) %>%
# median_hdi summary for coloring
dplyr::group_by(contrast,depth_maps_generation_filtering_mode) %>%
make_contrast_vars()plot it
plt_contrast(
brms_contrast_temp
# , caption_text = form_temp
, y_axis_title = "quality"
, facet = "depth_maps_generation_filtering_mode"
, label_size = 1.55
) +
scale_x_continuous(
labels = scales::percent_format(accuracy = 1)
, limits = c(-1.2,0.3)
, breaks = seq(-1,1,0.5)
) +
labs(
subtitle = "posterior predictive distribution of group constrasts with 95% & 50% HDI\nby filtering mode"
, x = "constrast BA Abs. % Error"
) +
theme(
axis.text.x = element_text(size = 7)
)
ggplot2::ggsave(
"../data/qlty_fltr_comp_ba_ape.png"
, plot = ggplot2::last_plot() + labs(subtitle = "")
, height = 7, width = 10.5
)and summarize these contrasts
brms_contrast_temp %>%
dplyr::group_by(contrast, depth_maps_generation_filtering_mode, pr_lt_zero) %>%
tidybayes::median_hdi(value) %>%
dplyr::arrange(contrast, depth_maps_generation_filtering_mode) %>%
dplyr::select(contrast, depth_maps_generation_filtering_mode, value, .lower, .upper, pr_lt_zero) %>%
dplyr::mutate(dplyr::across(.cols = is.numeric, .fns = ~ scales::percent(.x,accuracy = 0.1))) %>%
kableExtra::kbl(
digits = 2
, caption = "BA Abs. % Error<br>95% HDI of the posterior predictive distribution of group contrasts"
, col.names = c(
"quality contrast"
, "filtering mode"
, "median difference<br>BA Abs. % Error"
, "HDI low", "HDI high"
, "Pr(diff<0)"
)
, escape = F
) %>%
kableExtra::kable_styling() %>%
kableExtra::scroll_box(height = "8in")| quality contrast | filtering mode |
median difference BA Abs. % Error |
HDI low | HDI high | Pr(diff<0) |
|---|---|---|---|---|---|
| ultra high - high | aggressive | -1.8% | -14.4% | 8.5% | 69% |
| ultra high - high | moderate | -1.4% | -12.3% | 8.4% | 67% |
| ultra high - high | mild | -1.4% | -12.2% | 8.0% | 67% |
| ultra high - high | disabled | -1.8% | -12.1% | 7.0% | 72% |
| ultra high - medium | aggressive | -7.0% | -25.6% | 5.3% | 92% |
| ultra high - medium | moderate | -5.3% | -20.5% | 5.5% | 88% |
| ultra high - medium | mild | -4.9% | -19.6% | 5.2% | 88% |
| ultra high - medium | disabled | -4.9% | -18.9% | 4.6% | 89% |
| ultra high - low | aggressive | -12.7% | -38.1% | 2.5% | 97% |
| ultra high - low | moderate | -10.7% | -31.8% | 2.9% | 96% |
| ultra high - low | mild | -8.9% | -28.7% | 3.8% | 94% |
| ultra high - low | disabled | -8.4% | -26.9% | 3.4% | 94% |
| ultra high - lowest | aggressive | -21.1% | -59.9% | 2.7% | 98% |
| ultra high - lowest | moderate | -19.1% | -53.7% | 2.7% | 98% |
| ultra high - lowest | mild | -17.3% | -49.5% | 3.1% | 97% |
| ultra high - lowest | disabled | -16.8% | -47.3% | 2.7% | 97% |
| high - medium | aggressive | -5.0% | -22.6% | 7.3% | 84% |
| high - medium | moderate | -3.6% | -17.9% | 7.4% | 80% |
| high - medium | mild | -3.3% | -17.3% | 7.0% | 79% |
| high - medium | disabled | -2.9% | -15.9% | 7.3% | 78% |
| high - low | aggressive | -10.6% | -34.6% | 4.3% | 95% |
| high - low | moderate | -8.9% | -28.7% | 4.5% | 94% |
| high - low | mild | -7.3% | -25.6% | 5.3% | 91% |
| high - low | disabled | -6.4% | -23.4% | 5.8% | 89% |
| high - lowest | aggressive | -18.9% | -56.7% | 3.5% | 97% |
| high - lowest | moderate | -17.3% | -50.9% | 3.3% | 97% |
| high - lowest | mild | -15.5% | -46.3% | 4.0% | 96% |
| high - lowest | disabled | -14.7% | -44.7% | 3.4% | 96% |
| medium - low | aggressive | -5.1% | -27.1% | 10.0% | 81% |
| medium - low | moderate | -4.8% | -24.8% | -24.5% | 82% |
| medium - low | moderate | -4.8% | -23.3% | 8.5% | 82% |
| medium - low | mild | -3.6% | -20.6% | 9.3% | 77% |
| medium - low | disabled | -3.1% | -19.1% | 9.1% | 75% |
| medium - lowest | aggressive | -13.1% | -47.8% | 8.1% | 92% |
| medium - lowest | moderate | -13.1% | -45.4% | 5.4% | 94% |
| medium - lowest | mild | -11.7% | -41.7% | -41.5% | 93% |
| medium - lowest | mild | -11.7% | -40.1% | 6.7% | 93% |
| medium - lowest | disabled | -11.1% | -38.4% | 6.0% | 93% |
| low - lowest | aggressive | -7.2% | -38.9% | 14.8% | 80% |
| low - lowest | moderate | -7.5% | -36.9% | 11.2% | 83% |
| low - lowest | mild | -7.4% | -34.7% | 9.8% | 84% |
| low - lowest | disabled | -7.5% | -33.4% | 8.9% | 86% |
8.3.2.4 Software:Quality - interaction
draws_temp =
tidyr::crossing(
depth_maps_generation_quality = unique(ptcld_validation_data$depth_maps_generation_quality)
, software = unique(ptcld_validation_data$software)
) %>%
tidybayes::add_epred_draws(
brms_ba_ape_mod, allow_new_levels = T
# this part is crucial
, re_formula = ~ (1 | depth_maps_generation_quality) +
(1 | software) +
(1 | depth_maps_generation_quality:software)
) %>%
dplyr::rename(value = .epred) %>%
dplyr::mutate(med = tidybayes::median_hdci(value)$y)
# plot
draws_temp %>%
# remove out-of-sample obs
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(depth_maps_generation_quality, software)
, by = dplyr::join_by(depth_maps_generation_quality, software)
) %>%
dplyr::mutate(depth_maps_generation_quality = depth_maps_generation_quality %>% forcats::fct_rev()) %>%
# plot
ggplot(
mapping = aes(
y = value, x = software
, fill = med
)
) +
geom_hline(yintercept = 0, color = "gray33") +
tidybayes::stat_eye(
point_interval = median_hdi, .width = .95
, slab_alpha = 0.98
, interval_color = "black", linewidth = 1
, point_color = "black", point_fill = "black", point_size = 1
) +
scale_fill_stepsn(
n.breaks = 5
, colors = viridis::mako(n = 5, begin = 0.2, end = 0.9, direction = -1)
, limits = c(0,lmt_basal_area_abs_pct_error*.6)
, labels = scales::percent_format(accuracy = 0.1)
, show.limits = T
) +
facet_grid(cols = vars(depth_maps_generation_quality)) +
scale_y_continuous(
limits = c(-0.005,lmt_basal_area_abs_pct_error*1.1)
, labels = scales::percent_format(accuracy = 1)
) +
labs(
x = "software", y = "BA Abs. % Error"
, subtitle = "posterior predictive distribution of BA Abs. % Error with 95% HDI\nby dense cloud quality"
) +
theme_light() +
theme(
legend.position = "none"
, legend.direction = "horizontal"
, axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
, strip.text = element_text(color = "black", face = "bold")
) 
and a table of these 95% HDI values
table_temp =
draws_temp %>%
tidybayes::median_hdi(value) %>%
# remove out-of-sample obs
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(depth_maps_generation_quality, software)
, by = dplyr::join_by(depth_maps_generation_quality, software)
) %>%
dplyr::select(-c(.point,.interval, .width,.row)) %>%
dplyr::arrange(software,depth_maps_generation_quality) %>%
dplyr::mutate(dplyr::across(.cols = is.numeric, .fns = ~ scales::percent(.x,accuracy = 0.1)))
table_temp %>%
dplyr::select(-c(software)) %>%
kableExtra::kbl(
digits = 2
, caption = "BA Abs. % Error<br>95% HDI of the posterior predictive distribution"
, col.names = c(
"quality"
, "BA Abs. % Error<br>median"
, "HDI low", "HDI high"
)
, escape = F
) %>%
kableExtra::kable_styling() %>%
kableExtra::pack_rows(index = table(forcats::fct_inorder(table_temp$software))) %>%
kableExtra::scroll_box(height = "8in")| quality |
BA Abs. % Error median |
HDI low | HDI high |
|---|---|---|---|
| METASHAPE | |||
| ultra high | 14.9% | 4.9% | 28.8% |
| high | 16.3% | 5.4% | 31.3% |
| medium | 18.5% | 6.1% | 35.6% |
| low | 26.6% | 9.3% | 51.2% |
| lowest | 45.7% | 16.1% | 88.8% |
| OPENDRONEMAP | |||
| ultra high | 14.0% | 4.6% | 27.1% |
| high | 18.9% | 6.5% | 36.6% |
| medium | 32.1% | 10.9% | 62.4% |
| low | 34.4% | 10.7% | 65.4% |
| lowest | 38.0% | 13.4% | 74.8% |
| PIX4D | |||
| ultra high | 13.5% | 4.3% | 26.6% |
| high | 14.0% | 4.3% | 27.5% |
| medium | 15.9% | 5.3% | 31.8% |
| low | 21.3% | 7.1% | 41.8% |
we can also make pairwise comparisons
# calculate contrast
brms_contrast_temp = draws_temp %>%
tidybayes::compare_levels(
value
, by = depth_maps_generation_quality
, comparison = contrast_list
) %>%
dplyr::rename(contrast = depth_maps_generation_quality)
# separate contrast
brms_contrast_temp = brms_contrast_temp %>%
dplyr::ungroup() %>%
tidyr::separate_wider_delim(
cols = contrast
, delim = " - "
, names = paste0(
"sorter"
, 1:(max(stringr::str_count(brms_contrast_temp$contrast, "-"))+1)
)
, too_few = "align_start"
, cols_remove = F
) %>%
dplyr::filter(sorter1!=sorter2) %>%
dplyr::mutate(
dplyr::across(
tidyselect::starts_with("sorter")
, .fns = function(x){factor(
x, ordered = T
, levels = levels(ptcld_validation_data$depth_maps_generation_quality)
)}
)
, contrast = contrast %>%
forcats::fct_reorder(
paste0(as.numeric(sorter1), as.numeric(sorter2)) %>%
as.numeric()
)
) %>%
# median_hdi summary for coloring
dplyr::group_by(contrast, software) %>%
make_contrast_vars()
# remove out-of-sample obs
brms_contrast_temp = brms_contrast_temp %>%
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(software, depth_maps_generation_quality)
, by = dplyr::join_by(software, sorter1 == depth_maps_generation_quality)
) %>%
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(software, depth_maps_generation_quality)
, by = dplyr::join_by(software, sorter2 == depth_maps_generation_quality)
)plot it
plt_contrast(
brms_contrast_temp
# , caption_text = form_temp
, y_axis_title = "quality"
, facet = "software"
, label_size = 1.7
) +
scale_x_continuous(
labels = scales::percent_format(accuracy = 1)
, limits = c(-1.1,0.3)
, breaks = seq(-1,1,0.5)
) +
labs(
subtitle = "posterior predictive distribution of group constrasts with 95% & 50% HDI\nby software"
, x = "constrast BA Abs. % Error"
)
ggplot2::ggsave(
"../data/qlty_sftwr_comp_ba_ape.png"
, plot = ggplot2::last_plot() + labs(subtitle = "")
, height = 7, width = 10.5
)create plot for combining with other APE contrasts for publication
ptchwrk_qlty_sftwr_comp_ba_ape =
plt_contrast(
brms_contrast_temp
, y_axis_title = "quality contrast"
, facet = "software"
, label_size = 1.5
, label = "pr_diff_lab_sm"
, annotate_size = 1.6
) +
labs(
subtitle = "" # "constrast BA Abs. % Error"
, x = "BA Abs. % Error constrast"
) +
theme(
legend.position="none"
, axis.title.y = element_text(size = 10, face = "bold")
, axis.title.x = element_text(size = 8)
)
# ptchwrk_qlty_sftwr_comp_ba_apeand summarize these contrasts
table_temp = brms_contrast_temp %>%
dplyr::group_by(contrast, software, pr_lt_zero) %>%
tidybayes::median_hdi(value) %>%
dplyr::arrange(contrast, software) %>%
dplyr::select(contrast, software, value, .lower, .upper, pr_lt_zero) %>%
dplyr::arrange(software, contrast) %>%
dplyr::mutate(dplyr::across(.cols = is.numeric, .fns = ~ scales::percent(.x,accuracy = 0.1)))
table_temp %>%
dplyr::select(-c(software)) %>%
kableExtra::kbl(
digits = 2
, caption = "BA Abs. % Error<br>95% HDI of the posterior predictive distribution of group contrasts"
, col.names = c(
"quality contrast"
, "median difference<br>BA Abs. % Error"
, "HDI low", "HDI high"
, "Pr(diff<0)"
)
, escape = F
) %>%
kableExtra::kable_styling() %>%
kableExtra::pack_rows(index = table(forcats::fct_inorder(table_temp$software))) %>%
kableExtra::scroll_box(height = "8in")| quality contrast |
median difference BA Abs. % Error |
HDI low | HDI high | Pr(diff<0) |
|---|---|---|---|---|
| METASHAPE | ||||
| ultra high - high | -1.3% | -9.8% | 6.5% | 66% |
| ultra high - medium | -3.3% | -13.8% | 5.4% | 82% |
| ultra high - low | -11.3% | -27.0% | -0.4% | 99% |
| ultra high - lowest | -30.3% | -62.2% | -7.8% | 100% |
| high - medium | -2.0% | -12.1% | 7.0% | 71% |
| high - low | -9.9% | -25.0% | 0.9% | 98% |
| high - lowest | -28.8% | -60.3% | -7.5% | 100% |
| medium - low | -7.7% | -22.2% | 3.3% | 94% |
| medium - lowest | -26.6% | -57.3% | -5.9% | 100% |
| low - lowest | -18.3% | -45.2% | -0.8% | 99% |
| OPENDRONEMAP | ||||
| ultra high - high | -4.6% | -15.5% | 3.2% | 90% |
| ultra high - medium | -17.6% | -39.3% | -3.2% | 100% |
| ultra high - low | -19.9% | -42.2% | -4.0% | 100% |
| ultra high - lowest | -23.5% | -50.0% | -5.6% | 100% |
| high - medium | -12.7% | -32.2% | -0.3% | 99% |
| high - low | -14.9% | -34.8% | -1.2% | 100% |
| high - lowest | -18.5% | -42.4% | -1.7% | 100% |
| medium - low | -2.2% | -19.6% | 14.6% | 63% |
| medium - lowest | -5.4% | -28.8% | 13.7% | 74% |
| low - lowest | -3.4% | -24.0% | 16.7% | 67% |
| PIX4D | ||||
| ultra high - high | -0.5% | -8.7% | 7.5% | 56% |
| ultra high - medium | -2.2% | -12.1% | 6.6% | 73% |
| ultra high - low | -7.5% | -21.1% | 2.5% | 96% |
| high - medium | -1.8% | -11.4% | 6.9% | 69% |
| high - low | -7.0% | -20.5% | 2.5% | 95% |
| medium - low | -5.1% | -18.4% | 4.5% | 88% |
The contrasts above address the question “are there differences in APE based on dense point cloud generation quality within each software?”.
To address the different question of “are there differences in APE based on the processing software used at a given dense point cloud generation quality?” we need to utilize a different formulation of the comparison parameter within our call to the tidybayes::compare_levels function and calculate the contrast by software instead
# calculate contrast
brms_contrast_temp = draws_temp %>%
tidybayes::compare_levels(
value
, by = software
, comparison = "pairwise"
) %>%
dplyr::rename(contrast = software) %>%
# median_hdi summary for coloring
dplyr::group_by(contrast, depth_maps_generation_quality) %>%
make_contrast_vars()
# remove out-of-sample obs
brms_contrast_temp = brms_contrast_temp %>%
tidyr::separate_wider_delim(
cols = contrast
, delim = " - "
, names = paste0(
"sorter"
, 1:(max(stringr::str_count(brms_contrast_temp$contrast, "-"))+1)
)
, too_few = "align_start"
, cols_remove = F
) %>%
dplyr::filter(sorter1!=sorter2) %>%
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(software, depth_maps_generation_quality)
, by = dplyr::join_by(sorter1 == software, depth_maps_generation_quality)
) %>%
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(software, depth_maps_generation_quality)
, by = dplyr::join_by(sorter2 == software, depth_maps_generation_quality)
)plot it
plt_contrast(
brms_contrast_temp
# , caption_text = form_temp
, y_axis_title = "software"
, facet = "depth_maps_generation_quality"
, label_size = 1.7
) +
scale_x_continuous(
labels = scales::percent_format(accuracy = 1)
, limits = c(-0.7,0.6)
, breaks = seq(-0.6,0.6,0.3)
) +
facet_wrap(facets = vars(depth_maps_generation_quality), ncol = 2) +
labs(
subtitle = "posterior predictive distribution of group constrasts with 95% & 50% HDI\nby dense cloud quality"
, x = "constrast BA Abs. % Error"
) +
theme(
legend.position = c(.75, .13)
) +
guides(fill = guide_colorbar(theme = theme(
legend.key.width = unit(1, "lines"),
legend.key.height = unit(7, "lines")
)))
ggplot2::ggsave(
"../data/sftwr_qlty_comp_ba_ape.png"
, plot = ggplot2::last_plot() + labs(subtitle = "")
, height = 7, width = 10.5
)create plot for combining with other APE contrasts for publication
ptchwrk_sftwr_qlty_comp_ba_ape =
plt_contrast(
brms_contrast_temp
, y_axis_title = "software contrast"
, facet = "depth_maps_generation_quality"
, label_size = 1.7
, label = "pr_diff_lab_sm"
, annotate_size = 1.8
) +
facet_wrap(facets = vars(depth_maps_generation_quality), ncol = 3) +
labs(
subtitle = ""
, x = "BA Abs. % Error constrast"
) +
theme(
legend.position = "inside"
, legend.position.inside = c(.8, .10)
, axis.title.y = element_text(size = 10, face = "bold")
, axis.title.x = element_text(size = 8)
) +
guides(fill = guide_colorbar(theme = theme(
legend.key.width = unit(1, "lines"),
legend.key.height = unit(6.5, "lines")
)))
# ptchwrk_sftwr_qlty_comp_ba_apeand summarize these contrasts
table_temp = brms_contrast_temp %>%
dplyr::group_by(contrast, depth_maps_generation_quality, pr_lt_zero) %>%
tidybayes::median_hdi(value) %>%
dplyr::arrange(contrast, depth_maps_generation_quality) %>%
dplyr::select(contrast, depth_maps_generation_quality, value, .lower, .upper, pr_lt_zero) %>%
dplyr::mutate(dplyr::across(.cols = is.numeric, .fns = ~ scales::percent(.x,accuracy = 0.1)))
table_temp %>%
dplyr::select(-c(contrast)) %>%
kableExtra::kbl(
digits = 2
, caption = "BA Abs. % Error<br>95% HDI of the posterior predictive distribution of group contrasts"
, col.names = c(
"quality"
, "median difference<br>BA Abs. % Error"
, "HDI low", "HDI high"
, "Pr(diff<0)"
)
, escape = F
) %>%
kableExtra::kable_styling() %>%
kableExtra::pack_rows(index = table(forcats::fct_inorder(table_temp$contrast))) %>%
kableExtra::scroll_box(height = "8in")| quality |
median difference BA Abs. % Error |
HDI low | HDI high | Pr(diff<0) |
|---|---|---|---|---|
| OPENDRONEMAP - METASHAPE | ||||
| ultra high | -0.8% | -9.1% | 7.5% | 60% |
| high | 2.3% | -6.4% | 13.0% | 26% |
| medium | 13.0% | -0.9% | 34.0% | 2% |
| low | 7.2% | -7.8% | 27.0% | 13% |
| lowest | -7.0% | -34.1% | 15.1% | 77% |
| PIX4D - METASHAPE | ||||
| ultra high | -1.3% | -10.1% | 7.1% | 66% |
| high | -2.1% | -12.0% | 6.2% | 73% |
| medium | -2.3% | -12.9% | 8.1% | 72% |
| low | -4.9% | -20.6% | 8.7% | 81% |
| PIX4D - OPENDRONEMAP | ||||
| ultra high | -0.4% | -9.2% | 8.8% | 55% |
| high | -4.6% | -16.0% | 4.7% | 87% |
| medium | -15.6% | -37.2% | -0.6% | 99% |
| low | -12.3% | -33.2% | 3.4% | 96% |
8.3.2.5 Software:Filtering - interaction
draws_temp =
tidyr::crossing(
depth_maps_generation_filtering_mode = unique(ptcld_validation_data$depth_maps_generation_filtering_mode)
, software = unique(ptcld_validation_data$software)
) %>%
tidybayes::add_epred_draws(
brms_ba_ape_mod, allow_new_levels = T
# this part is crucial
, re_formula = ~ (1 | depth_maps_generation_filtering_mode) +
(1 | software) +
(1 | depth_maps_generation_filtering_mode:software)
) %>%
dplyr::rename(value = .epred) %>%
dplyr::mutate(med = tidybayes::median_hdci(value)$y)
# plot
draws_temp %>%
# remove out-of-sample obs
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(depth_maps_generation_filtering_mode, software)
, by = dplyr::join_by(depth_maps_generation_filtering_mode, software)
) %>%
# plot
ggplot(
mapping = aes(
y = value, x = software
, fill = med
)
) +
geom_hline(yintercept = 0, color = "gray33") +
tidybayes::stat_eye(
point_interval = median_hdi, .width = .95
, slab_alpha = 0.98
, interval_color = "black", linewidth = 1
, point_color = "black", point_fill = "black", point_size = 1
) +
scale_fill_stepsn(
n.breaks = 5
, colors = viridis::mako(n = 5, begin = 0.2, end = 0.9, direction = -1)
, limits = c(0,lmt_basal_area_abs_pct_error*.6)
, labels = scales::percent_format(accuracy = 0.1)
, show.limits = T
) +
facet_grid(cols = vars(depth_maps_generation_filtering_mode)) +
scale_y_continuous(
limits = c(-0.005,lmt_basal_area_abs_pct_error*1.1)
, labels = scales::percent_format(accuracy = 1)
) +
labs(
x = "software", y = "BA Abs. % Error"
, subtitle = "posterior predictive distribution of BA Abs. % Error with 95% HDI\nby filtering mode"
) +
theme_light() +
theme(
legend.position = "none"
, legend.direction = "horizontal"
, axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
, strip.text = element_text(color = "black", face = "bold")
) 
and a table of these 95% HDI values
table_temp =
draws_temp %>%
tidybayes::median_hdi(value) %>%
dplyr::select(-c(.point,.interval, .width,.row)) %>%
dplyr::arrange(software,depth_maps_generation_filtering_mode) %>%
dplyr::mutate(dplyr::across(.cols = is.numeric, .fns = ~ scales::percent(.x,accuracy = 0.1)))
table_temp %>%
dplyr::select(-c(software)) %>%
kableExtra::kbl(
digits = 2
, caption = "BA Abs. % Error<br>95% HDI of the posterior predictive distribution"
, col.names = c(
"filtering mode"
, "BA Abs. % Error<br>median"
, "HDI low", "HDI high"
)
, escape = F
) %>%
kableExtra::kable_styling() %>%
kableExtra::pack_rows(index = table(forcats::fct_inorder(table_temp$software))) %>%
kableExtra::scroll_box(height = "8in")| filtering mode |
BA Abs. % Error median |
HDI low | HDI high |
|---|---|---|---|
| METASHAPE | |||
| aggressive | 28.2% | 7.5% | 59.1% |
| moderate | 21.5% | 6.5% | 44.9% |
| mild | 20.0% | 5.3% | 40.9% |
| disabled | 18.6% | 5.4% | 38.7% |
| OPENDRONEMAP | |||
| aggressive | 25.5% | 6.5% | 53.3% |
| moderate | 25.9% | 6.9% | 53.8% |
| mild | 25.1% | 6.9% | 52.3% |
| disabled | 22.5% | 5.9% | 46.7% |
| PIX4D | |||
| aggressive | 20.3% | 4.7% | 45.8% |
| moderate | 18.9% | 5.1% | 39.5% |
| mild | 17.1% | 4.3% | 35.9% |
| disabled | 16.2% | 4.0% | 34.0% |
we can also make pairwise comparisons
# calculate contrast
brms_contrast_temp = draws_temp %>%
tidybayes::compare_levels(
value
, by = depth_maps_generation_filtering_mode
, comparison = "pairwise"
) %>%
dplyr::rename(contrast = depth_maps_generation_filtering_mode)
# separate contrast
brms_contrast_temp = brms_contrast_temp %>%
dplyr::ungroup() %>%
tidyr::separate_wider_delim(
cols = contrast
, delim = " - "
, names = paste0(
"sorter"
, 1:(max(stringr::str_count(brms_contrast_temp$contrast, "-"))+1)
)
, too_few = "align_start"
, cols_remove = F
) %>%
dplyr::filter(sorter1!=sorter2) %>%
dplyr::mutate(
dplyr::across(
tidyselect::starts_with("sorter")
, .fns = function(x){factor(
x, ordered = T
, levels = levels(ptcld_validation_data$depth_maps_generation_filtering_mode)
)}
)
, contrast = contrast %>%
forcats::fct_reorder(
paste0(as.numeric(sorter1), as.numeric(sorter2)) %>%
as.numeric()
) %>%
# re order for filtering mode
forcats::fct_rev()
) %>%
# median_hdi summary for coloring
dplyr::group_by(contrast, software) %>%
make_contrast_vars()
# remove out-of-sample obs
brms_contrast_temp = brms_contrast_temp %>%
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(software, depth_maps_generation_filtering_mode)
, by = dplyr::join_by(software, sorter1 == depth_maps_generation_filtering_mode)
) %>%
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(software, depth_maps_generation_filtering_mode)
, by = dplyr::join_by(software, sorter2 == depth_maps_generation_filtering_mode)
)plot it
plt_contrast(
brms_contrast_temp
# , caption_text = form_temp
, y_axis_title = "filtering mode"
, facet = "software"
, label_size = 1.6
) +
scale_x_continuous(
labels = scales::percent_format(accuracy = 1)
, limits = c(-0.8,0.65)
, breaks = seq(-0.8,0.8,0.4)
) +
labs(
subtitle = "posterior predictive distribution of group constrasts with 95% & 50% HDI\nby software"
, x = "constrast BA Abs. % Error"
)
ggplot2::ggsave(
"../data/fltr_sftwr_comp_ba_ape.png"
, plot = ggplot2::last_plot() + labs(subtitle = "")
, height = 7, width = 10.5
)create plot for combining with other APE contrasts for publication
ptchwrk_fltr_sftwr_comp_ba_ape =
plt_contrast(
brms_contrast_temp
, y_axis_title = "filtering mode contrast"
, facet = "software"
, label_size = 1.7
, label = "pr_diff_lab_sm"
, annotate_size = 1.6
) +
labs(
subtitle = "" # "constrast BA Abs. % Error"
, x = "BA Abs. % Error constrast"
) +
theme(
legend.position="none"
, axis.title.y = element_text(size = 10, face = "bold")
, axis.title.x = element_text(size = 8)
)
# ptchwrk_fltr_sftwr_comp_ba_apeand summarize these contrasts
table_temp = brms_contrast_temp %>%
dplyr::group_by(contrast, software, pr_lt_zero) %>%
tidybayes::median_hdi(value) %>%
dplyr::arrange(contrast, software) %>%
dplyr::select(contrast, software, value, .lower, .upper, pr_lt_zero) %>%
dplyr::arrange(software, contrast) %>%
dplyr::mutate(dplyr::across(.cols = is.numeric, .fns = ~ scales::percent(.x,accuracy = 0.1)))
table_temp %>%
dplyr::select(-c(software)) %>%
kableExtra::kbl(
digits = 2
, caption = "BA Abs. % Error<br>95% HDI of the posterior predictive distribution of group contrasts"
, col.names = c(
"filtering contrast"
, "median difference<br>BA Abs. % Error"
, "HDI low", "HDI high"
, "Pr(diff<0)"
)
, escape = F
) %>%
kableExtra::kable_styling() %>%
kableExtra::pack_rows(index = table(forcats::fct_inorder(table_temp$software))) %>%
kableExtra::scroll_box(height = "8in")| filtering contrast |
median difference BA Abs. % Error |
HDI low | HDI high | Pr(diff<0) |
|---|---|---|---|---|
| METASHAPE | ||||
| disabled - mild | -1.2% | -9.7% | 6.3% | 67% |
| disabled - moderate | -2.6% | -12.6% | 4.7% | 80% |
| disabled - aggressive | -8.9% | -25.9% | 1.1% | 98% |
| mild - moderate | -1.3% | -10.7% | 6.6% | 67% |
| mild - aggressive | -7.5% | -24.0% | 2.0% | 96% |
| moderate - aggressive | -6.0% | -22.0% | 3.0% | 93% |
| OPENDRONEMAP | ||||
| disabled - mild | -2.2% | -13.8% | 6.4% | 75% |
| disabled - moderate | -3.1% | -14.9% | 5.8% | 80% |
| disabled - aggressive | -2.7% | -16.3% | 7.6% | 74% |
| mild - moderate | -0.8% | -11.8% | 9.5% | 59% |
| mild - aggressive | -0.4% | -13.3% | 11.9% | 54% |
| moderate - aggressive | 0.2% | -11.2% | 13.0% | 48% |
| PIX4D | ||||
| disabled - mild | -0.8% | -8.6% | 6.4% | 62% |
| disabled - moderate | -2.4% | -11.8% | 4.4% | 80% |
| mild - moderate | -1.5% | -10.5% | 5.7% | 71% |
8.3.2.6 Software:Quality:Filtering - interaction
# get draws
fltr_sftwr_draws_temp =
tidyr::crossing(
depth_maps_generation_quality = unique(ptcld_validation_data$depth_maps_generation_quality)
, depth_maps_generation_filtering_mode = unique(ptcld_validation_data$depth_maps_generation_filtering_mode)
, software = unique(ptcld_validation_data$software)
) %>%
tidybayes::add_epred_draws(
brms_ba_ape_mod, allow_new_levels = T
# this part is crucial
, re_formula = ~ (1 | depth_maps_generation_quality) + # main effects
(1 | depth_maps_generation_quality) +
(1 | depth_maps_generation_filtering_mode) +
(1 | software) +
# two-way interactions
(1 | depth_maps_generation_quality:depth_maps_generation_filtering_mode) +
(1 | depth_maps_generation_quality:software) +
(1 | depth_maps_generation_filtering_mode:software) +
# three-way interactions
(1 | depth_maps_generation_quality:depth_maps_generation_filtering_mode:software)
) %>%
dplyr::rename(value = .epred) %>%
dplyr::mutate(med = tidybayes::median_hdci(value)$y)
# plot
qlty_fltr_sftwr_ba_ape =
fltr_sftwr_draws_temp %>%
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode)
, by = dplyr::join_by(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode)
) %>%
dplyr::mutate(depth_maps_generation_quality = depth_maps_generation_quality %>% forcats::fct_rev()) %>%
ggplot(
mapping = aes(
y = value
, x = depth_maps_generation_filtering_mode
, fill = med
)
) +
geom_hline(yintercept = 0, color = "gray33") +
tidybayes::stat_eye(
point_interval = median_hdi, .width = .95
, slab_alpha = 0.98
, interval_color = "black", linewidth = 1
, point_color = "black", point_fill = "black", point_size = 1
) +
scale_fill_stepsn(
n.breaks = 5
, colors = viridis::mako(n = 5, begin = 0.2, end = 0.9, direction = -1)
, limits = c(0,lmt_basal_area_abs_pct_error*.6)
, labels = scales::percent_format(accuracy = 0.1)
, show.limits = T
) +
facet_grid(
rows = vars(software)
, cols = vars(depth_maps_generation_quality)
# , switch = "y"
) +
scale_y_continuous(
limits = c(-0.005,lmt_basal_area_abs_pct_error*1.45)
, labels = scales::percent_format(accuracy = 1)
) +
labs(
x = "filtering mode", y = "BA Abs. % Error"
, subtitle = "posterior predictive distribution of BA Abs. % Error with 95% HDI\nby dense cloud quality and software"
# , caption = form_temp
) +
theme_light() +
theme(
legend.position = "none"
, legend.direction = "horizontal"
, axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
, strip.text = element_text(color = "black", face = "bold")
, panel.grid = element_blank()
# , strip.placement = "outside"
) +
guides(
fill = guide_legend(override.aes = list(shape = NA, size = 6, alpha = 0.9, lwd = NA))
)
# print it
qlty_fltr_sftwr_ba_ape
and a table of these 95% HDI values
table_temp =
fltr_sftwr_draws_temp %>%
tidybayes::median_hdi(value) %>%
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode)
, by = dplyr::join_by(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode)
) %>%
dplyr::mutate(depth_maps_generation_quality = depth_maps_generation_quality %>% forcats::fct_rev()) %>%
dplyr::select(c(
software, depth_maps_generation_quality, depth_maps_generation_filtering_mode
, value, .lower, .upper
)) %>%
dplyr::ungroup() %>%
dplyr::arrange(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode) %>%
dplyr::mutate(dplyr::across(.cols = is.numeric, .fns = ~ scales::percent(.x,accuracy = 0.1)))
table_temp %>%
# dplyr::select(-c(software)) %>%
kableExtra::kbl(
digits = 2
, caption = "BA Abs. % Error<br>95% HDI of the posterior predictive distribution"
, col.names = c(
"software", "quality", "filtering mode"
, "BA Abs. % Error<br>median"
, "HDI low", "HDI high"
)
, escape = F
) %>%
kableExtra::kable_styling() %>%
# kableExtra::pack_rows(index = table(forcats::fct_inorder(table_temp$software))) %>%
kableExtra::collapse_rows(columns = 1:2, valign = "top") %>%
kableExtra::scroll_box(height = "8in")| software | quality | filtering mode |
BA Abs. % Error median |
HDI low | HDI high |
|---|---|---|---|---|---|
| METASHAPE | lowest | aggressive | 61.6% | 23.2% | 116.5% |
| moderate | 47.5% | 18.3% | 88.1% | ||
| mild | 41.8% | 15.4% | 77.7% | ||
| disabled | 39.3% | 14.6% | 73.3% | ||
| low | aggressive | 38.0% | 14.5% | 71.0% | |
| moderate | 27.0% | 10.0% | 49.9% | ||
| mild | 23.8% | 8.7% | 44.1% | ||
| disabled | 21.6% | 7.6% | 40.6% | ||
| medium | aggressive | 25.3% | 9.4% | 47.9% | |
| moderate | 17.9% | 6.4% | 33.7% | ||
| mild | 16.7% | 6.2% | 31.3% | ||
| disabled | 15.6% | 6.0% | 29.4% | ||
| high | aggressive | 20.7% | 7.5% | 39.4% | |
| moderate | 15.8% | 6.2% | 29.8% | ||
| mild | 15.2% | 5.9% | 28.3% | ||
| disabled | 14.4% | 5.5% | 26.5% | ||
| ultra high | aggressive | 18.7% | 6.5% | 35.8% | |
| moderate | 14.7% | 5.7% | 27.7% | ||
| mild | 14.1% | 5.5% | 26.2% | ||
| disabled | 12.9% | 5.0% | 23.8% | ||
| OPENDRONEMAP | lowest | aggressive | 40.4% | 14.5% | 77.1% |
| moderate | 41.2% | 15.4% | 77.7% | ||
| mild | 39.3% | 14.5% | 74.0% | ||
| disabled | 35.3% | 13.6% | 67.0% | ||
| low | aggressive | 38.3% | 14.8% | 73.5% | |
| moderate | 37.7% | 14.4% | 70.8% | ||
| mild | 35.2% | 13.3% | 66.0% | ||
| disabled | 31.2% | 11.5% | 57.9% | ||
| medium | aggressive | 35.2% | 13.1% | 67.1% | |
| moderate | 34.2% | 12.6% | 64.2% | ||
| mild | 33.5% | 12.0% | 62.3% | ||
| disabled | 30.0% | 11.5% | 56.0% | ||
| high | aggressive | 19.3% | 7.2% | 36.6% | |
| moderate | 20.1% | 7.8% | 37.2% | ||
| mild | 19.8% | 7.3% | 36.7% | ||
| disabled | 17.9% | 6.5% | 33.1% | ||
| ultra high | aggressive | 14.2% | 5.3% | 27.2% | |
| moderate | 15.1% | 5.7% | 28.2% | ||
| mild | 14.7% | 5.7% | 27.8% | ||
| disabled | 12.6% | 4.6% | 23.7% | ||
| PIX4D | low | moderate | 23.0% | 8.1% | 42.4% |
| mild | 19.5% | 7.6% | 37.1% | ||
| disabled | 17.9% | 6.5% | 34.1% | ||
| medium | moderate | 16.2% | 5.8% | 30.4% | |
| mild | 14.7% | 5.3% | 27.8% | ||
| disabled | 13.9% | 5.1% | 26.4% | ||
| high | moderate | 14.2% | 5.3% | 26.9% | |
| mild | 13.1% | 4.8% | 24.8% | ||
| disabled | 12.7% | 4.6% | 23.7% | ||
| ultra high | moderate | 14.0% | 5.0% | 26.2% | |
| mild | 13.0% | 4.7% | 24.1% | ||
| disabled | 12.1% | 4.6% | 22.6% |
we can also make pairwise comparisons
# calculate contrast
brms_contrast_temp = fltr_sftwr_draws_temp %>%
tidybayes::compare_levels(
value
, by = depth_maps_generation_filtering_mode
, comparison = "pairwise"
) %>%
dplyr::rename(contrast = depth_maps_generation_filtering_mode)
# separate contrast
brms_contrast_temp = brms_contrast_temp %>%
dplyr::ungroup() %>%
tidyr::separate_wider_delim(
cols = contrast
, delim = " - "
, names = paste0(
"sorter"
, 1:(max(stringr::str_count(brms_contrast_temp$contrast, "-"))+1)
)
, too_few = "align_start"
, cols_remove = F
) %>%
dplyr::filter(sorter1!=sorter2) %>%
dplyr::mutate(
dplyr::across(
tidyselect::starts_with("sorter")
, .fns = function(x){factor(
x, ordered = T
, levels = levels(ptcld_validation_data$depth_maps_generation_filtering_mode)
)}
)
, contrast = contrast %>%
forcats::fct_reorder(
paste0(as.numeric(sorter1), as.numeric(sorter2)) %>%
as.numeric()
) %>%
# re order for filtering mode
forcats::fct_rev()
) %>%
# median_hdi summary for coloring
dplyr::group_by(contrast, software, depth_maps_generation_quality) %>%
make_contrast_vars()
# remove out-of-sample obs
brms_contrast_temp = brms_contrast_temp %>%
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode)
, by = dplyr::join_by(software, depth_maps_generation_quality, sorter1==depth_maps_generation_filtering_mode)
) %>%
dplyr::inner_join(
ptcld_validation_data %>% dplyr::distinct(software, depth_maps_generation_quality, depth_maps_generation_filtering_mode)
, by = dplyr::join_by(software, depth_maps_generation_quality, sorter2==depth_maps_generation_filtering_mode)
) %>%
dplyr::mutate(depth_maps_generation_quality = depth_maps_generation_quality %>% forcats::fct_rev())plot it
brms_contrast_temp %>%
plt_contrast(
facet = c("depth_maps_generation_quality", "software")
, y_axis_title = "filtering mode"
, label_size = 0
) +
facet_grid(
rows = vars(software)
, cols = vars(depth_maps_generation_quality)
) +
scale_x_continuous(
labels = scales::percent_format(accuracy = 1)
, limits = c(-0.8,0.4)
, breaks = seq(-0.8,0.8,0.4)
) +
labs(
subtitle = "posterior predictive distribution of group constrasts with 95% & 50% HDI\nby dense cloud quality and software"
, x = "constrast BA Abs. % Error"
)
ggplot2::ggsave(
"../data/qlty_fltr_sftwr_comp_ba_ape.png"
, plot = ggplot2::last_plot() + labs(subtitle = "")
, height = 7, width = 10.5
)Export some final images for publication
p1_temp = qlty_fltr_sftwr_ba_pe + labs(subtitle = "A: Basal Area % Error") + theme(plot.subtitle = element_text(face="bold"))
p2_temp = qlty_fltr_sftwr_ba_ape + labs(subtitle = "B: Basal Area Abs. % Error") + theme(plot.subtitle = element_text(face="bold"))
# export
p1_temp / p2_temp
ggplot2::ggsave(
filename = paste0("../data/qlty_fltr_sftwr_ba_comb.jpeg")
, plot = ggplot2::last_plot()
, width = 8.5
, height = 11
, units = "in"
, dpi = "print"
)patchwork of Height APE contrasts
layout_temp = c(
# area(t, l, b, r)
patchwork::area(2, 1, 2, 1)
, patchwork::area(2, 3, 2, 3)
, patchwork::area(4, 1, 4, 3)
)
# check the layout
# plot(layout_temp)
############################
# patchwork for height
############################
ptchwrk_qlty_sftwr_comp_ba_ape +
labs(subtitle = "A: Quality Contrast by Software") +
theme(plot.background = element_rect(colour = "gray88", fill=NA, size=3)) +
ptchwrk_fltr_sftwr_comp_ba_ape +
labs(subtitle = "B: Filtering Mode Contrast by Software") +
theme(plot.background = element_rect(colour = "gray88", fill=NA, size=3)) +
patchwork::free(
ptchwrk_sftwr_qlty_comp_ba_ape +
labs(subtitle = "C: Software Contrast by Quality") +
theme(plot.background = element_rect(colour = "gray88", fill=NA, size=3)) ) +
# plot_annotation(tag_levels = list(c('#', '&'), '1')) +
patchwork::plot_layout(
design = layout_temp
, widths = c(1,0.01,1)
, heights = c(0.01,1,0.01,1,0.01)
) &
scale_x_continuous(
labels = scales::percent_format(accuracy = 1)
, limits = c(-1.1,0.52)
, breaks = seq(-0.8,0.8,0.4)
) &
theme(
axis.title.y = element_blank()
, plot.subtitle = element_text(face = "bold", hjust = 0.0)
# , plot.background = element_rect(colour = "gray88", fill=NA, size=3)
)
8.4 Posterior Predictive Checks
Markov chain Monte Carlo (MCMC) simulations were conducted using the brms package (Bürkner 2017) to estimate posterior predictive distributions of the parameters of interest. We ran three chains of 100,000 iterations with the first 50,000 discarded as burn-in. Trace-plots were utilized to visually assess model convergence and sufficient convergence was checked with \(\hat{R}\) values near 1 (Brooks & Gelman, 1998). Posterior predictive checks were used to evaluate model goodness-of-fit by comparing data simulated from the model with the observed data used to estimate the model parameters (Hobbs & Hooten, 2015). Calculating the proportion of MCMC iterations in which the test statistic (i.e., mean and sum of squares) from the simulated data and observed data are more extreme than one another provides the Bayesian P-value. Lack of fit is indicated by a value close to 0 or 1 while a value of 0.5 indicates perfect fit (Hobbs & Hooten, 2015).
To learn more about this approach to posterior predictive checks, check out Gabry’s (2022) vignette, Graphical posterior predictive checks using the bayesplot package.
8.4.3 \(\hat{R}\) values
plt_rhat_temp <- function(my_mod) {
my_mod %>%
brms::rhat() %>%
as.data.frame() %>%
tibble::rownames_to_column(var = "parameter") %>%
dplyr::rename_with(tolower) %>%
dplyr::rename(rhat = 2) %>%
dplyr::filter(
stringr::str_starts(parameter, "b_")
| stringr::str_starts(parameter, "r_")
| stringr::str_starts(parameter, "sd_")
| parameter == "phi"
| parameter == "sigma"
) %>%
dplyr::mutate(
parameter = parameter %>%
stringr::str_replace_all("depth_maps_generation_quality", "quality") %>%
stringr::str_replace_all("depth_maps_generation_filtering_mode", "filtering")
, chk = (rhat <= 1*0.998 | rhat >= 1*1.002)
) %>%
ggplot(aes(x = rhat, y = parameter, color = chk, fill = chk)) +
geom_vline(xintercept = 1, linetype = "dashed", color = "gray44", lwd = 1.2) +
geom_vline(xintercept = 1*0.998, lwd = 1.5) +
geom_vline(xintercept = 1*1.002, lwd = 1.5) +
geom_vline(xintercept = 1*0.999, lwd = 1.2, color = "gray33") +
geom_vline(xintercept = 1*1.001, lwd = 1.2, color = "gray33") +
geom_point() +
scale_fill_manual(values = c("navy", "firebrick")) +
scale_color_manual(values = c("navy", "firebrick")) +
scale_y_discrete(NULL, breaks = NULL) +
labs(
x = latex2exp::TeX("$\\hat{R}$")
, subtitle = latex2exp::TeX("MCMC chain convergence check for $\\hat{R}$ values")
) +
theme_light() +
theme(
legend.position = "none"
, axis.text.y = element_text(size = 3)
, panel.grid.major.x = element_blank()
, panel.grid.minor.x = element_blank()
, plot.title = element_text(size = 9)
, plot.subtitle = element_text(size = 8)
, axis.title.y = element_blank()
, axis.title.x = element_text(size = 9)
, panel.grid = element_blank()
, axis.text.x = element_text(size = 6.5)
)
}
# plot together with patchwork
brms_ba_pe_mod %>% plt_rhat_temp() + labs(title = "Basal Area Percentage Error") +
brms_ba_ape_mod %>% plt_rhat_temp() + labs(title = "Basal Area Absolute Percentage Error") &
theme(plot.subtitle = element_blank())
8.4.4 ESS values
The effective length of an MCMC chain is indicated by the effective sample size (ESS), which refers to the sample size of the MCMC chain not to the sample size of the data Kruschke (2015) notes:
One simple guideline is this: For reasonably accurate and stable estimates of the limits of the 95% HDI, an ESS of 10,000 is recommended. This is merely a heuristic based on experience with practical applications, it is not a requirement. If accuracy of the HDI limits is not crucial for your application, then a smaller ESS may be sufficient. (p.184)
plt_ess_temp <- function(my_mod) {
# get ess values from model summary
dplyr::bind_rows(
summary(my_mod) %>%
purrr::pluck("random") %>%
purrr::flatten() %>%
purrr::keep_at(~ .x == "Bulk_ESS") %>%
unlist() %>%
dplyr::as_tibble()
, summary(my_mod) %>%
purrr::pluck("fixed") %>%
purrr::flatten() %>%
purrr::keep_at(~ .x == "Bulk_ESS") %>%
unlist() %>%
dplyr::as_tibble()
) %>%
dplyr::rename(ess = 1) %>%
dplyr::mutate(parameter = dplyr::row_number(), chk = ess<10000) %>%
ggplot(aes(x = ess, y = parameter, color = chk, fill = chk)) +
geom_vline(xintercept = 10000, linetype = "dashed", color = "gray44", lwd = 1.2) +
geom_segment( aes(x = 0, xend=ess, yend=parameter), color="black") +
geom_point() +
scale_fill_manual(values = c("blue4", "blue3")) +
scale_color_manual(values = c("blue4", "blue3")) +
scale_y_continuous(NULL, breaks = NULL) +
scale_x_continuous(labels = scales::comma) +
labs(
x = "ESS"
, subtitle = "MCMC chain resolution check for effective sample size (ESS) values"
, y = ""
) +
theme_light() +
theme(
legend.position = "none"
, axis.text.y = element_text(size = 4)
, panel.grid.major.x = element_blank()
, panel.grid.minor.x = element_blank()
, plot.title = element_text(size = 9)
, plot.subtitle = element_text(size = 8)
, axis.title.y = element_blank()
, axis.title.x = element_text(size = 9)
, axis.text.x = element_text(size = 6.5)
)
}
# plot together with patchwork
brms_ba_pe_mod %>% plt_ess_temp() + labs(title = "Basal Area Percentage Error") +
brms_ba_ape_mod %>% plt_ess_temp() + labs(title = "Basal Area Absolute Percentage Error") &
theme(plot.subtitle = element_blank())
8.4.5 Mean and SD
posterior predictive check for the overall model combining mean and sd
plt_pp_temp <- function(my_mod) {
my_mod %>%
brms::pp_check(type = "stat_2d", ndraws = 5000) +
theme_light() +
theme(
legend.position = "top", legend.direction = "horizontal"
, legend.text = element_text(size = 8)
, plot.title = element_text(size = 9)
)
}
# plot together with patchwork
brms_ba_pe_mod %>% plt_pp_temp() + labs(title = "Basal Area Percentage Error") +
brms_ba_ape_mod %>% plt_pp_temp() + labs(title = "Basal Area Absolute Percentage Error") &
theme(plot.subtitle = element_blank())
8.4.5.1 Bayesian p-value
The Bayesian p-value is the probability that a test statistic in the reference distribution exceeds its value in the data. The Bayesian p-value is calculated from the posterior predictive distribution of the new data and the distribution of the observed data. We estimate the probability that the test statistic calculated from “new” data arising from our model (\(y_{new}\)) is more extreme than the test statistic calculated from the observed data (\(y\)): \(\text{P-value}(y) = Pr(T(y_{new}) > T(y))\) where the test statistic \(T(y)\) describes the distribution of the data as a summary of the data; it could be the mean, variance, the coefficient of variation, the kurtosis, the maximum, or the minimum of the observed data set, or it might be an “omnibus” statistic like a squared discrepancy or a chi-square value Hobbs and Hooten (2015, p. 188)
Bayesian P values for mean and standard deviation test statistics The P values for the mean (P mean) give the probability that the mean of the data of new, out-of-sample observations simulated from the model exceeds the mean of the observed data. The P values for the standard deviation (P SD) give the probability that the standard deviation of new, out-of-sample observations simulated from the model exceeds the standard deviation of the observed data. Large (\(\gtrapprox 0.90\)) or small (\(\lessapprox 0.10\)) values indicate lack of fit. Hobbs and Hooten (2015);Hobbs et al. (2024, Appendix S2 p. 8)
Check the Bayesian p-values between the models
# get the model p-values
set.seed(16) # replicate results
dplyr::bind_rows(
get_mod_p_val(brms_ba_pe_mod, ndraws = 5000) %>%
dplyr::mutate(model = "brms_ba_pe_mod")
, get_mod_p_val(brms_ba_ape_mod, ndraws = 5000) %>%
dplyr::mutate(model = "brms_ba_ape_mod")
) %>%
dplyr::relocate(model) %>%
kableExtra::kbl(digits = 3) %>%
kableExtra::kable_styling()| model | P.mean | P.sd |
|---|---|---|
| brms_ba_pe_mod | 0.493 | 0.564 |
| brms_ba_ape_mod | 0.671 | 0.894 |
8.4.5.2 \(\sigma\) posteriors
Finally, we can quantify the variation in our dependent variable by comparing the \(\sigma\) (sd) posteriors
plt_sigma_temp <- function(my_mod) {
# extract the posterior draws
brms::as_draws_df(my_mod) %>%
dplyr::select(c(tidyselect::starts_with("sd_"))) %>%
tidyr::pivot_longer(dplyr::everything()) %>%
# dplyr::group_by(name) %>%
# tidybayes::median_hdi(value) %>%
dplyr::mutate(
name = name %>%
stringr::str_replace_all("depth_maps_generation_quality", "quality") %>%
stringr::str_replace_all("depth_maps_generation_filtering_mode", "filtering") %>%
stringr::str_remove_all("sd_") %>%
stringr::str_remove_all("__Intercept") %>%
stringr::str_replace_all("_", " ") %>%
forcats::fct_reorder(value)
) %>%
# plot
ggplot(aes(x = value, y = name)) +
tidybayes::stat_dotsinterval(
point_interval = median_hdi, .width = .95
, justification = -0.04
, shape = 21 #, point_size = 3
, quantiles = 100
) +
scale_x_continuous(breaks = NULL) +
labs(x = "", y = ""
, subtitle = latex2exp::TeX("$\\sigma_\\beta$ posterior distributions")
) +
theme_light() +
theme(
plot.subtitle = element_text(size = 8)
, plot.title = element_text(size = 9)
)
}
# plot together with patchwork
brms_ba_pe_mod %>% plt_sigma_temp() + labs(title = "Basal Area Percentage Error") +
brms_ba_ape_mod %>% plt_sigma_temp() + labs(title = "Basal Area Absolute Percentage Error")
Variance of study site is slightly stronger than variance of depth map generation quality, but the posterior predictive distributions overlap a good deal. The study site (the “subjects” in our study) seems to have the overall strongest effect, but this comes with high uncertainty. Taken alone, the influence of quality and software comes with huge uncertainty. This makes sense as the influence of software largely depends on the depth map generation quality. Likewise, the influence of quality depends on the software. We are fairly certain of this conditional influence based on the uncertainty associated with the interaction of these two terms. Filtering mode by itself has the weakest effect on tree detection and this comes with relatively high certainty. As with the influence of software, we are certain that the effect of filtering mode is conditional on the depth map generation quality.