README.qmd

---
#title: "README"
#author: "Bhan, Lam"
format: gfm
editor: visual
---

# Replication Code for "Automating Urban Soundscape Enhancements with AI: In-situ Assessment of Quality and Restorativeness in Traffic-Exposed Residential Areas"

The GitHub repository contains the code to replicate the analysis, figures and tables for the paper titled: "Automating Urban Soundscape Enhancements with AI: In-situ Assessment of Quality and Restorativeness in Traffic-Exposed Residential Areas". The initial stable release v1.0.0 has been archived on Zenodo: [![](images/clipboard-1276697786.png)](https://doi.org/10.5281/zenodo.11141691 "Link to archived stable release on Zenodo")

The data that support the findings of this study are openly available in NTU research\
data repository DR-NTU (Data) at <https://doi.org/10.21979/N9/NEH5TR>.

The subheadings in this repository follows the headings in the paper (after the [Data Loading] section) for consistency.

The following figures are produced by this replication code:

-   [`Table 2`](#tab_2) in [2.3. Stimuli and automatic masker selection](#sec2.3)

-   [`Table 3`](#tab_3) in [2.4. Non-acoustic environmental conditions for in-situ validation study](#sec2.4)

-   [`Table 4`](#tab_4) in [2.5. Participants](#sec2.5)

-   [`Table 5`](#tab_5) in [3. Results](#sec3)

-   [`Table B.1`](#tab_b1) in [3. Results](#sec3)

-   [`Figure 3`](#fig_3) in [3. Results](#sec3)

-   [`Table 6`](#tab_6) in [3.3. Correlation between subjective metrics](#sec3.3)

-   [`Table B.2`](#tab_b2) in [3.4. Effect of order, group size and initial conditions](#sec3.4)

-   [`Table 7`](#tab_7) in [4. Results: Objective binaural measurements](#sec4)

-   [`Table B.3`](#tab_b3) in [4. Results: Objective binaural measurements](#sec4)

-   [`Figure 4`](#fig_4) in [4. Results: Objective binaural measurements](#sec4)

## Initialisation

```{r init, echo=FALSE, message=FALSE}
knitr::opts_chunk$set(warning = FALSE, message = FALSE)

if(!requireNamespace("pacman", quietly=TRUE)) install.packages("pacman")
pacman::p_load(
        #data wrangling
        dataverse, tidyverse, 
        #visualisation
        gt, gtsummary, gtools, paletteer, ggpubr, pals,
        #stats
        effectsize, emmeans, psych, lme4
)

#load script
source("scripts/supportScripts.R")

metadata<-data.frame(
        category=c(
                "categorical","categorical","categorical",
                "interval","interval",
                "categorical","categorical","categorical",
                "interval","interval",
                "interval","interval","interval","interval","interval"
                ),
        variable=c(
                "dom_noise","dom_natural","dom_human",
                "PosAff","NegAff",
                "overall","appropriate","loudness",
                "ISOPL","ISOEV",
                "PRSSFas","PRSSBA","PRSSCom","PRSSEC","PRSSES"
                )
)

```

## Data Loading

First, check if the RData file (`fullData.Rdata`) exists. If it does not exist, download the RData file containing all the data tables required to replicate all the figures, tables and analyses in this paper from the Dataverse repository (<https://doi.org/10.21979/N9/NEH5TR>).

```{r dataloading, echo=TRUE, message=FALSE}

#check if RData file exists then download from the DOI if it doesn't
if (!file.exists("data/fullData.RData")){
        as_binary<-dataverse::get_file(
                file = "fullData.RData",
                dataset = "doi:10.21979/N9/NEH5TR",
                server = "https://researchdata.ntu.edu.sg"
        )
        
        #write binary file
        writeBin(as_binary, con = "data/fullData.RData")
}

load("data/fullData.RData")
load("data/splvstime.RData")
```

## 2. Method

### 2.3. Stimuli and automatic masker selection {#sec2.3}

#### `Table 2:` Frequency distribution of the maskers chosen by the AMSS during the 10-min listening period across all "AMSS" group participants. Description and availability of the corresponding maskers as detailed by Ooi et al. \[48\] in the ARAUS dataset. {#tab_2}

```{r table2, message=FALSE, echo=TRUE}

tbl2_predictions <- predict_session_data |> 
        tbl_summary(
                include = c(predictions),
                #only show percentage
                statistic = list(all_categorical() ~ "{p}%")
                ) |> 
        remove_row_type(predictions, type = "header") |>
        as_gt() |>
        cols_add(
                 Description = c(
                        "Bahama Mockingbird",
                        "Baltimore Oriole",
                        "Northern Cardinal",
                        "Veery",
                        "Common Redshank"
                )
        ) |>
        cols_label(
                label = md("**Maskers**"),
                stat_0 = md("**Frequency (%)**"),
                Description = md("**Description**")
        ) |>
        rm_footnotes()
tbl2_predictions

```

### 2.4. Non-acoustic environmental conditions for in-situ validation study {#sec2.4}

#### `Table 3`: Summary statistics of environmental parameters captured at ROOF during the 10-min listening period across all participants. {#tab_3}

```{r table3, message=FALSE, echo=TRUE}

#list of metrics to compute mean
metric_list <- c("temperature","humidity","lux","wind_speed","24h_psi","pm25")

enviro_summary_tbl <- enviro_session_data |>
        tbl_summary(
                include = metric_list,
                by = condition,
                digits = all_continuous() ~ 2,
                type = list(everything() ~ 'continuous'),
                statistic = list(all_continuous() ~ "{mean} ({sd})"),
                missing = "no"
                ) |>
        #add_n() %>% # add column with total number of non-missing observations
        modify_header(label = "**Environmental Parameter**") |>
        add_p(
                everything() ~ "wilcox.test",
                pvalue_fun = function(x) style_pvalue(x, digits = 3)
                ) |>
        as_gt()
        
enviro_summary_tbl
```

### 2.5. Participants {#sec2.5}

#### `Table 4:` Summary of participant demographics and non-acoustic factors (PSS-10, WNSS, WHO-5, baseline annoyance) across each condition (AMSS and AMB). {#tab_4}

```{r table4, message=FALSE, echo=TRUE}

#custom function for ks-test between ambient and amss groups across {}
ks_test <- function(data, variable, by, ...) {
  data <- data[c(variable, by)] %>% dplyr::filter(complete.cases(.))
  ks.test(data[[variable]] ~ factor(data[[by]])) %>%
  broom::tidy()
}

amss_insitu_participant_data |>
        dplyr::select(c(condition,gender,age,pss,wnss,wbi), 
                              starts_with("annoy-")) |>
        drop_na() |>
        tbl_summary(
                by = condition,
                type = list(!c(gender) ~ 'continuous'),
                statistic = list(all_continuous() ~ "{mean} ({sd})"),
                digits = list(all_continuous() ~ c(2, 2)),
                label = list(
                        pss ~ "PSS-10",
                        wnss ~ "INS",
                        wbi ~ "WHO-5",
                        `annoy-aircraft` ~ "BA@aircraft~",
                        `annoy-mrt` ~ "BA@mrt~",
                        `annoy-consite` ~ "BA@consite~",
                        `annoy-reno` ~ "BA@reno~",
                        `annoy-traffic` ~ "BA@traffic~",
                        `annoy-animals` ~ "BA@animals~",
                        `annoy-children` ~ "BA@children~",
                        `annoy-people` ~ "BA@people~",
                        `annoy-others` ~ "BA@others~"
                        )
        ) |>
        add_p(
                test = list(
                        gender ~ "prop.test",
                        c(age,pss,wnss,wbi,
                          starts_with("annoy-")) ~ c("ks_test")
                        ),
                pvalue_fun = function(x) style_pvalue(x, digits = 2)
                ) |>
        add_overall() |>
        as_gt() |>
        text_transform(
                locations = cells_body(),
                fn = function(x) {
                        str_replace_all(x,pattern = "@",
                                        replacement = "<sub>") |>
                                str_replace_all("~","</sub>") 
                        }
                )
```

## 3. Results: Site evaluation questionnaire {#sec3}

#### `Table 5:` Mean responses $\mu$ (standard deviation $\sigma$) of perceptual attributes in the site evaluation questionnaire investigated for the validation study, organized by *site* and *condition*. The scales for all attributes are normalised to the range $[-1,1]$. Percentage changes are computed between the AMB and AMSS for *site*, and between ROOF and GND for *condition* as scale changes on the $[-1,1]$ range with respect to the former. For instance, a change from $-0.25$ in the AMB condition to $0.75$ in the AMSS condition would be reported as a $50$% change. Significant changes as determined by posthoc tests are indicated in bold. {#tab_5}

```{r table5, message=FALSE, echo=TRUE}

#compute statistical tests for all variables
stat.results<-twoWLMERMANOVA(metadata,amss_insitu_participant_data)

#retrieve significant results
signif_posthoc <- stat.results |>
        dplyr::filter(
                p.value<0.05 & 
                        !grepl("MP",term) &
                        grepl("Contrasts",test)
                      ) 

signif_posthoc_bycond <- signif_posthoc |>
        filter(grepl("Condition",test)) |>
        select(variable,term)

signif_posthoc_bysite <- signif_posthoc |>
        filter(grepl("Site",test)) |>
        select(variable,term)

#summarise by condition at each site
summary_bycond_gt <- amss_insitu_participant_data |>
        dplyr::filter(!site=="MP") |>
        dplyr::select(!c(pID,order,partGrp,p:m,pss:who5)) |>
        tbl_strata(
            strata = site,
            .tbl_fun =
              ~ .x %>%
                tbl_summary(
                        by = condition, 
                        missing = "no",
                        type = everything() ~ "continuous",
                        statistic = all_continuous() ~ "{mean} ({sd})",
                        label = list(
                                dom_noise ~ "!DOM#@Noi~",
                                dom_human ~ "!DOM#@Hum~",
                                dom_natural ~ "!DOM#@Nat~",
                                overall ~ "!OSQ#",
                                appropriate ~ "!APPR#",
                                loudness ~ "!PLN#",
                                ISOPL ~ "!ISOPL#",
                                ISOEV ~ "!ISOEV#",
                                PRSSFas ~ "!PRSS#@Fas~",
                                PRSSBA ~ "!PRSS#@BA~",
                                PRSSCom ~ "!PRSS#@Com~",
                                PRSSEC ~ "!PRSS#@EC~",
                                PRSSES ~ "!PRSS#@ES~",
                                PosAff ~ "!PA#",
                                NegAff ~ "!NA#"
                                )
                        ) %>%
                add_difference(
                        estimate_fun = everything() ~ 
                                function(x) paste0(style_sigfig((-x/2) * 100), "%")
                ),
            .header = "**{strata}**"
            ) |>
        modify_column_hide(columns = c(p.value_1,ci_1,p.value_2,ci_2)) |>
        # remove difference footnote
        modify_footnote(update = everything() ~ NA) |>
        modify_header(all_stat_cols() ~ "**{level}**") 

#summarise by site for each cond
summary_bysite_gt <- amss_insitu_participant_data |>
        dplyr::filter(!site=="MP") |>
        dplyr::select(!c(pID,order,partGrp,p:m,pss:who5)) |>
        dplyr::mutate(site=factor(site, levels=c("GND","ROOF"))) |>
        tbl_strata(
            strata = condition,
            .tbl_fun =
              ~ .x %>%
                tbl_summary(
                        by = site, 
                        missing = "no",
                        type = everything() ~ "continuous",
                        statistic = all_continuous() ~ "{mean} ({sd})",
                        label = list(
                                dom_noise ~ "!DOM#@Noi~",
                                dom_human ~ "!DOM#@Hum~",
                                dom_natural ~ "!DOM#@Nat~",
                                overall ~ "!OSQ#",
                                appropriate ~ "!APPR#",
                                loudness ~ "!PLN#",
                                ISOPL ~ "!ISOPL#",
                                ISOEV ~ "!ISOEV#",
                                PRSSFas ~ "!PRSS#@Fas~",
                                PRSSBA ~ "!PRSS#@BA~",
                                PRSSCom ~ "!PRSS#@Com~",
                                PRSSEC ~ "!PRSS#@EC~",
                                PRSSES ~ "!PRSS#@ES~",
                                PosAff ~ "!PA#",
                                NegAff ~ "!NA#"
                                )
                        ) %>%
                add_difference(
                        estimate_fun = everything() ~ 
                                function(x) paste0(style_sigfig((-x/2) * 100), "%")
                        ),
            .header = "**{strata}**"
            ) |>
        modify_column_hide(columns = c(p.value_1,ci_1,p.value_2,ci_2)) |>
        # remove difference footnote
        modify_footnote(update = everything() ~ NA) |>
        modify_header(all_stat_cols() ~ "**{level}**") 

#merge tables columnwise        
tbl_merge_cond_site <-
        tbl_merge(
                tbls = list(summary_bycond_gt, summary_bysite_gt),
                tab_spanner = c(
                        "**Contrasts by condition at ...**",
                        "**Contrasts by site under ...**"
                        )
                ) |>
        as_gt() |>
        text_transform(
                locations = cells_body(),
                fn = function(x) {
                        str_replace_all(
                                x,pattern = "@",
                                replacement = "<sub>"
                                ) |>
                                str_replace_all("~","</sub>") |>
                                str_replace_all("!","<i>") |>
                                str_replace_all("#","</i>")
                        }
                ) |>
        #highlight significant results at GND site
        tab_style(
                style = list(
                        cell_fill(color = "#FDE992"),
                        cell_text(weight = "bold")
                        ),
                locations = cells_body(
                        columns = c(estimate_1_1),
                        rows = variable %in% 
                                signif_posthoc_bycond[
                                        str_detect(signif_posthoc_bycond$term,
                                                   "GND"),
                                        ]$variable
                        )
        ) |>
        #highlight significant results at ROOF site
        tab_style(
                style = list(
                        cell_fill(color = "#FDE992"),
                        cell_text(weight = "bold")
                        ),
                locations = cells_body(
                        columns = c(estimate_2_1),
                        rows = variable %in% 
                                signif_posthoc_bycond[
                                        str_detect(signif_posthoc_bycond$term,
                                                   "ROOF"),
                                        ]$variable
                        )
        ) |>
        #highlight significant results under AMB condition
        tab_style(
                style = list(
                        cell_fill(color = "#FDE992"),
                        cell_text(weight = "bold")
                        ),
                locations = cells_body(
                        columns = c(estimate_1_2),
                        rows = variable %in% 
                                signif_posthoc_bysite[
                                        str_detect(signif_posthoc_bysite$term,
                                                   "AMB"),
                                        ]$variable
                        )
        ) |>
        #highlight significant results under AMSS condition
        tab_style(
                style = list(
                        cell_fill(color = "#FDE992"),
                        cell_text(weight = "bold")
                        ),
                locations = cells_body(
                        columns = c(estimate_2_2),
                        rows = variable %in% 
                                signif_posthoc_bysite[
                                        str_detect(signif_posthoc_bysite$term,
                                                   "AMSS"),
                                        ]$variable
                        )
        )
        
tbl_merge_cond_site

```

#### `Table B.1:` Summary of statistical tests for attributes in soundscape evaluation questionnaire (sound source dominance, overall quality, appropriateness, loudness, *ISOPL*, *ISOEV*, and PRSS dimensions) across site (GND and ROOF), condition (AMSS and AMB), and their interaction (*site*:*condition*). Test abbreviations and symbols for significance levels and effect sizes are defined in the footnote. {#tab_b1}

```{r tableB1, message=FALSE, echo=TRUE}

#plot statistical test reults in a table
stat.results |>
        dplyr::mutate(
                p.value=paste0(
                        gtools::stars.pval(p.value),
                        formatC(p.value, format = "f", digits = 4)
                        ),
                eff.size=case_when(
                        eff.size > 0.14 ~ paste0(
                                "(L)",
                                formatC(eff.size,format = "f",digits = 4)
                                ),
                        eff.size > 0.06 ~ paste0(
                                "(M)",
                                formatC(eff.size,format = "f",digits = 4)
                                ),
                        eff.size > 0.01 ~ paste0(
                                "(S)",formatC(eff.size,format = "f",digits = 4)
                                ),
                        .default = formatC(eff.size,format = "f",digits = 4))) |>
        dplyr::group_by(variable) |>
        gt::gt() |>
        cols_label(
                term ~ html("<b>Term</b>"),
                test ~ html("<b>Test</b>"),
                p.value ~ html("<i>p</i>-<b>value</b>"),
                eff.size ~ html("<b>Effect Size</b>")
        )
```

#### `Figure 3:` Simple contrast of means across all perceptual attributes organized by *condition* and *site* Contrasts by *condition* are between groups at each *site*, whereas contrasts by *site* are within group for each *condition*. The scales for all attributes are normalised to the range \[-1,1\]. Significant differences as determined by posthoc contrast tests are accentuated {#fig_3}

```{r figure3, message=FALSE}

#prepare dataframe for plotting
plot_df<- amss_insitu_participant_data |>
        dplyr::filter(!site=="MP") |> #remove the meeting point
        pivot_longer(
                cols = metadata$variable,
                values_to = "score",
                names_to = "Attribute"
                ) |>
        dplyr::select(!c(p:who5)) |>
        dplyr::mutate(
                Attribute=factor(Attribute,levels=metadata$variable)
        )

#Contrasts by condition under GND and ROOF site conditions
plot_site_df<-plot_df |>
        group_by(site,condition,Attribute) |>
        summarise(
                mean=mean(score,na.rm=TRUE),
                sd=sd(score,na.rm=TRUE)
                ) |>
        ungroup() |>
        dplyr::mutate(
                significant=ifelse(
                        site=="ROOF" & 
                                Attribute %in% c("dom_natural","PosAff",
                                                 "ISOPL","overall",
                                                 "PRSSFas","PRSSBA",
                                                 "PRSSCom"),
                                   TRUE,FALSE)
                )

#Contrasts by site under AMSS and AMB conditions
plot_cond_df<-plot_df |>
        group_by(condition,site,Attribute) |>
        summarise(mean=mean(score,na.rm=TRUE),
                  sd=sd(score,na.rm=TRUE)) |>
        ungroup() |>
        dplyr::mutate(
                significant=ifelse((
                        condition=="AMB" & 
                                Attribute %in% c("dom_natural","ISOPL",
                                              "overall","loudness",
                                              "PRSSFas","PRSSBA","PRSSCom")) |
                                (condition=="AMSS" & 
                                         Attribute %in% c("loudness")),
                                   TRUE,FALSE)
                )

#prepare legend labels
plot_legends <- c(
                bquote(~DOM[Noi]),
                bquote(~DOM[Nat]),
                bquote(~DOM[Hum]),
                "PA","NA","OQ","APPR",
                "PLN","ISOPL","ISOEV",
                bquote(~PRSS[Fas]),
                bquote(~PRSS[BA]),
                bquote(~PRSS[Com]),
                bquote(~PRSS[EC]),
                bquote(~PRSS[ES])
                )

site_plot<-ggplot(data = plot_site_df,
       aes(x=condition,y=mean,group=Attribute,
           alpha=significant,color=Attribute)) +
        geom_line() + 
        geom_point() +
        scale_alpha_discrete(
                range=c(0.2, 1),
                guide = 'none' #turn off legend
                ) +
        scale_color_paletteer_d(
                palette = "awtools::bpalette",
                labels=plot_legends
                ) +
        facet_wrap(vars(site)) +
        theme_minimal() +
        theme(
                panel.grid.major.x = element_blank(),
                panel.grid.minor.x = element_blank()
        ) +
        geom_hline(
                yintercept = 0,
                color="darkgrey",
                size=1
                ) +
        ylim(c(-1,1)) +
        ylab("Normalised mean") +
        ggpubr::labs_pubr()
        

cond_plot<-ggplot(data = plot_cond_df,
       aes(x=site,y=mean,group=Attribute,
           alpha=significant,color=Attribute)) +
        geom_line() + 
        geom_point() +
        scale_alpha_discrete(
                range=c(0.2, 1),
                guide = 'none' #turn off legend
                ) +
        scale_color_paletteer_d(
                palette = "awtools::bpalette",
                labels=plot_legends) +
        facet_wrap(vars(condition)) +
        theme_minimal() +
        theme(
                panel.grid.major.x = element_blank(),
                panel.grid.minor.x = element_blank()
        ) +
        geom_hline(
                yintercept = 0,
                color="darkgrey",
                size=1
                ) +
        ylim(c(-1,1)) +
        ylab("Normalised mean") +
        ggpubr::labs_pubr()


comb_plot<-ggarrange(
        site_plot, cond_plot, #labels = c("Site", "Condition"),
        common.legend = TRUE, legend = "right"
)
comb_plot

```

### 3.3. Correlation between subjective metrics {#sec3.3}

#### `Table 6:` Kendall correlation matrix between all attributes in the site evaluation questionnaire where the significance of each entry in the upper triangle is denoted with a Holm-adjusted $p$-value and each entry in the lower triangle is denoted with an unadjusted $p$-value. Asterisks indicate \*$p<0.05$; \*\*$p<0.01$; \*\*\*$p<0.001$; \*\*\*\*$p<0.0001$. The unit diagonal has been removed for clarity. {#tab_6}

```{r table6, message=FALSE, echo=FALSE}

#compute Kendall correlation using the psych package
corr.all<-psych::corr.test(
        amss_insitu_participant_data |> 
            dplyr::select(dom_noise,dom_human,dom_natural,
                          PosAff,NegAff,overall,appropriate,loudness,
                          ISOPL,ISOEV,PRSSFas,PRSSBA,PRSSCom,PRSSEC,PRSSES
                          ),
    method = "kendall")

corr.r <- as.data.frame(corr.all$r) |>
        dplyr::mutate_all(.funs = list(~formatC(.,digits = 2,
                                                format = "f"))) 

corr.p <- as.data.frame(corr.all$p) |>
        dplyr::mutate_all(.funs=list(~gtools::stars.pval(.)))

corr.mat <- matrix(
        paste0(as.matrix(corr.p),as.matrix(corr.r)),
        nrow=nrow(corr.p), 
        dimnames=dimnames(corr.p)
        )

#remove diagonal values
diag(corr.mat)=NA

corr.mat <- as.data.frame(corr.mat) |>
        `rownames<-`(colnames(corr.p)) 

corr.mat |>
        gt(rownames_to_stub = TRUE) |>
         sub_missing(
            columns = everything(),
            missing_text = ""
            )

```

### 3.4. Effect of order, group size and initial conditions {#sec3.4}

#### `Table B.2:` Summary of exact two-sample Kolmogorov-Smirnov tests to examine effect of order (GND--ROOF or ROOF--GND) and group size (1 or $>1$) on each soundscape evaluation attribute (sound source dominance, overall quality, appropriateness, loudness, ISOPL, ISOEV, and PRSS dimensions) across each condition (AMSS and AMB). All the $p$-values were adjusted for multiple comparisons within conditions with the Benjamini-Hochberg (BH) method. {#tab_b2}

```{r tableB2, message=FALSE}

ks.df <- amss_insitu_participant_data |>
        dplyr::select(c(condition,partGrp,order,metadata$variable)) |>
        pivot_longer(names_to = "variable",
                     values_to = "score",
                     cols = metadata$variable)

ks.order.grpsize <- rbind(
        #KS Test by order
        ks.df |> 
                dplyr::group_by(condition,variable) |>
                dplyr::summarise(
                        ks_test = list(ks.test(score[order == 18],
                                               score[order == 81],
                                               exact = NULL,
                                               alternative = "two.sided")),
                        ks.pvalue = ks_test[[1]]$p.value
                ) |>
                dplyr::ungroup() |>
                dplyr::group_by(condition) |>
                dplyr::mutate(ks.padj = p.adjust(ks.pvalue, method="BH"),
                              confvar = "order"),
        #KS Test by group size
        ks.df |> 
                dplyr::group_by(condition,variable) |>
                dplyr::summarise(
                        ks_test = list(ks.test(score[partGrp == "single"],
                                               score[partGrp == "multi"],
                                               exact = NULL,
                                               alternative = "two.sided")),
                        ks.pvalue = ks_test[[1]]$p.value
                ) |>
                dplyr::ungroup() |>
                dplyr::group_by(condition) |>
                dplyr::mutate(ks.padj = p.adjust(ks.pvalue, method="BH"),
                              confvar = "group size")
        ) |>
        #add significance stars
        dplyr::mutate(
                        ks.padj=paste0(gtools::stars.pval(ks.padj),
                               formatC(ks.padj,
                                       format = "f",
                                       digits = 2))
        ) |>
        dplyr::select(!ks_test) |>
        pivot_wider(values_from = ks.padj,
                    names_from = variable,
                    id_cols = c(condition,confvar)) |>
        #reorder columns
        dplyr::select(
                dom_noise, dom_human, dom_natural, 
                PosAff, NegAff, 
                overall, appropriate, loudness, 
                ISOPL, ISOEV,
                PRSSFas, PRSSBA, PRSSCom, PRSSEC, PRSSES, confvar) |>
        dplyr::group_by(confvar) |>
        gt::gt() |>
        cols_label(
                ISOEV ~ html("<i>ISOEV</i>"),
                ISOPL ~ html("<i>ISOPL</i>"),
                NegAff ~ html("<i>NA</i>"),
                PosAff ~ html("<i>PA</i>"),
                PRSSFas ~ html("<i>PRSS</i><sub>Fas</sub>"),
                PRSSBA ~ html("<i>PRSS</i><sub>BA</sub>"),
                PRSSCom ~ html("<i>PRSS</i><sub>Com</sub>"),
                PRSSEC ~ html("<i>PRSS</i><sub>EC</sub>"),
                PRSSES ~ html("<i>PRSS</i><sub>ES</sub>"),
                appropriate ~ html("<i>APPR</i>"),
                dom_natural ~ html("<i>DOM</i><sub>Nat</sub>"),
                dom_human ~ html("<i>DOM</i><sub>Hum</sub>"),
                dom_noise ~ html("<i>DOM</i><sub>Noi</sub>"),
                loudness ~ html("<i>PLN</i>"),
                overall ~ html("<i>OSQ</i>"),
        ) 
ks.order.grpsize
```

## 4. Results: Objective binaural measurements {#sec4}

#### `Table 7:` Summary of mean *L*<sub>Aeq</sub>, *L*<sub>Ceq</sub>, *N*<sub>95</sub>, *ISOPL*, *OSQ*, *PRSS*<sub>Fas</sub>, *PRSS*<sub>BA</sub>, and *PRSS*<sub>Com</sub> values across 20 AMSS and AMB sessions in each of the GND and RTGP sites. Supplemented mean values for the AMSS sessions excluding aircraft flyby (3 in GND; 1 in ROOF) are included. {#tab_7}

```{r table7, message=FALSE, echo=TRUE}

#summary including aircraft flyby
isopl_obj_tbl <- comb_obj_insitu_session_data |>
        ungroup() |>
        pivot_wider(
                names_from = "attribute",
                values_from = "score"
        ) |>
        tbl_strata2(
                strata = condition,
                .tbl_fun =
                        ~ .x %>%
                        tbl_summary(include = c(`L[Aeq]`,`L[Ceq]`,
                                                `N[95]`,`ISOPL`,
                                                `OSQ`,`PRSS[Fas]`,
                                                `PRSS[BA]`,`PRSS[Com]`),
                            by = site,
                            type = list(everything() ~ 'continuous'),
                            statistic = list(everything() ~ "{mean} ({sd})"),
                            missing = "no",
                            digits = list(all_continuous() ~ c(2, 2)))
                )

#remove aircraft noise
isopl_obj_noaircraft_tbl<-comb_obj_insitu_session_data |> 
        ungroup() |>
        pivot_wider(
                names_from = "attribute",
                values_from = "score"
        ) |>
        dplyr::filter(
                !((date=="20230914" & sessionTime=="9" & site=="GND") |
                        (date=="20230915" & sessionTime=="14" 
                         & site %in% c("GND","ROOF")) |
                                 (date=="20230915" & sessionTime=="16" 
                                  & site=="GND"))
        ) |>
        #dplyr::group_by(condition,site) |>
        tbl_strata2(
                strata = condition,
                .tbl_fun =
                        ~ .x %>%
                        tbl_summary(include = c(`L[Aeq]`,`L[Ceq]`,
                                                `N[95]`,`ISOPL`,
                                                `OSQ`,`PRSS[Fas]`,
                                                `PRSS[BA]`,`PRSS[Com]`),
                            by = site,
                            type = list(everything() ~ 'continuous'),
                            statistic = list(everything() ~ "{mean} ({sd})"),
                            missing = "no",
                            digits = list(all_continuous() ~ c(2, 2))),
                .header = "**{strata}** (without aircraft flyby)"
                ) |>
        modify_header(label="") |>
        modify_column_hide(columns = c(stat_1_1,stat_2_1))

#merge tables columnwise        
tbl_merge_cond_site <-
        tbl_merge(
                tbls = list(isopl_obj_tbl, isopl_obj_noaircraft_tbl)
        ) |>
        modify_spanning_header(
                c(stat_1_1_1,stat_2_1_1) ~ "**AMB**"
        ) |>
        modify_spanning_header(
                c(stat_1_2_1,stat_2_2_1) ~ "**AMSS**"
        ) |>
        modify_spanning_header(
                c(stat_1_2_2,stat_2_2_2) ~ "**AMSS (without aircraft flyby)**"
        ) |>
        as_gt()
tbl_merge_cond_site
```

### `Table B.3:` Kendall correlation matrix between all objective acoustic measures and perceptual attributes in the site evaluation questionnaire where the significance of each entry in the upper triangle is denoted with a Holm-adjusted $p$-value and each entry in the lower triangle is denoted with an unadjusted $p$-value. Asterisks indicate \*$p<0.05$; \*\*$p<0.01$; \*\*\*$p<0.001$; \*\*\*\*$p<0.0001$. The unit diagonal has been removed for clarity. {#tab_b3}

```{r tableB3, echo=TRUE, message=FALSE}

corr_obj<-psych::corr.test(
        comb_obj_insitu_session_data |>
                ungroup() |>
                pivot_wider(
                        names_from = attribute,
                        values_from = score,
                ) |>
                dplyr::select(ISOPL,
                              OSQ,
                              `PRSS[Fas]`,`PRSS[BA]`,`PRSS[Com]`,
                              `L[Aeq]`,`L[Ceq]`,`N[95]`) |>
                dplyr::mutate_all(.funs = as.numeric),
    method = "kendall")

corr_obj_r <- as.data.frame(corr_obj$r) |>
        dplyr::mutate_all(
                .funs = list(~formatC(.,digits = 2,format = "f"))
                ) 

corr_obj_p <- as.data.frame(corr_obj$p) |>
        dplyr::mutate_all(
                .funs=list(~gtools::stars.pval(.))
                )

corr_obj_mat <- matrix(
        paste0(as.matrix(corr_obj_p),as.matrix(corr_obj_r)),
        nrow=nrow(corr_obj_p), 
        dimnames=dimnames(corr_obj_p)
        ) 

#remove diagonal values
diag(corr_obj_mat)=NA

#convert to data frame
corr_obj_mat <- as.data.frame(corr_obj_mat) |>
        `rownames<-`(colnames(corr_obj_p)) 


corr_obj_mat |>
        gt(rownames_to_stub = TRUE) |>
        sub_missing(
            columns = everything(),
            missing_text = ""
            )

```

### `Figure 4:` Mean perceptual *ISOPL*, *OSQ*, *PRSS*<sub>Fas</sub>, *PRSS*<sub>BA</sub>, and *PRSS*<sub>Com</sub> scores across all participants per session (y-axis) as a function of normalized objective *L*<sub>Aeq</sub>, *L*<sub>Ceq</sub>, *N*<sub>95</sub>, scores of each session (x-axis). Fifty percent of the sessions lie within the median contours computed for AMB--GND, AMB--ROOF, AMSS--GND, AMSS--ROOF contrast subgroups. The left to right columns represent *L*<sub>Aeq</sub>, *L*<sub>Ceq</sub>, and *N*<sub>95</sub>, and each row represents each of the perceptual metrics, respectively. {#fig_4}

```{r densityplot, message=FALSE, echo=TRUE}

#plot isopl vs decibel score

#density colors
densityClr<-pals::stevens.pinkgreen()[c(9,5,7,3)]

#ISOPL vs LAeq
isopl_LA_plot <- ggplot(comb_obj_insitu_session_data) +
  facet_wrap(~attribute,
    labeller = label_parsed,
    scales = "free", 
    ncol = 1
  ) +
  # add mean points of each session
  geom_point(aes(y = `L[Aeq]`, x = score, shape = site, color = pair)) +
  # add kde contours
  stat_density_2d(
    bins = 3, contour_var = "ndensity", breaks = c(0.5),
    geom = "density_2d",
    aes(y = `L[Aeq]`, x = score, color = pair)
  ) +
  scale_fill_manual(values = densityClr) +
  scale_color_manual(values = densityClr, name = "subgroup") +
  ylim(c(50, 70)) +
  xlim(c(-1, 1)) +
  geom_vline(
    xintercept = 0,
    color = "darkgrey",
    size = 1
  ) +
  #xlab("ISOPL") +
  ylab(bquote(paste(L[Aeq], ", dB(A)"))) +
  ggpubr::labs_pubr() +
  theme(
    legend.position = "bottom",
    strip.background = element_blank(),
    strip.text.x = element_blank(),
    axis.title.x = element_blank()
  )

isopl_LC_plot <- ggplot(comb_obj_insitu_session_data) +
  # add mean points of each session
  geom_point(aes(y = `L[Ceq]`, x = score, shape = site, color = pair)) +
  facet_wrap(~attribute,
    labeller = label_parsed,
    scales = "free",
    ncol = 1
  ) +
  # add kde contours
  stat_density_2d(
    bins = 3, contour_var = "ndensity", breaks = c(0.5),
    geom = "density_2d",
    aes(y = `L[Ceq]`, x = score, color = pair)
  ) +
  scale_fill_manual(values = densityClr) +
  scale_color_manual(values = densityClr, name = "subgroup") +
  ylim(c(60, 80)) +
  xlim(c(-1, 1)) +
  geom_vline(
    xintercept = 0,
    color = "darkgrey",
    size = 1
  ) +
  ylab(bquote(paste(L[Ceq], ", dB(C)"))) +
  ggpubr::labs_pubr() +
  theme(
    legend.position = "bottom",
    strip.background = element_blank(),
    strip.text.x = element_blank(),
    axis.title.x = element_blank(),
  )

isopl_N95_plot <- ggplot(comb_obj_insitu_session_data) +
  facet_wrap(~attribute,
    labeller = label_parsed,
    scales = "free",
    ncol = 1,
    strip.position = "right"
  ) +
  # add mean points of each session
  geom_point(aes(y = `N[95]`, x = score, shape = site, color = pair)) +
  # add kde contours
  stat_density_2d(
    bins = 3, contour_var = "ndensity", breaks = c(0.5),
    geom = "density_2d",
    aes(y = `N[95]`, x = score, color = pair)
  ) +
  scale_fill_manual(values = densityClr) +
  scale_color_manual(values = densityClr, name = "subgroup") +
  ylim(c(5, 20)) +
  xlim(c(-1, 1)) +
  geom_vline(
    xintercept = 0,
    color = "darkgrey",
    size = 1
  ) +
  ylab(bquote(paste(N[95], ", soneGF"))) +
  ggpubr::labs_pubr() +
  theme(
    legend.position = "bottom",
    axis.title.x = element_blank()
  )

comb_isopl_obj_plot <- ggarrange(
  isopl_LA_plot, 
  isopl_LC_plot,
  isopl_N95_plot,
  common.legend = TRUE, legend = "bottom",
  nrow = 1
)
comb_isopl_obj_plot
```

## SPL vs Time

```{r}

#aggregate SPL data by second to create smoother graph
spl_vs_time_data_agg <- spl_vs_time_data |> 
        #remove aircraft flyby sessions
        dplyr::filter(
                !((date=="20230914" & sessionTime=="9" & site=="GND") |
                        (date=="20230915" & sessionTime=="14" 
                         & site %in% c("GND","ROOF")) |
                                 (date=="20230915" & sessionTime=="16" 
                                  & site=="GND"))
                ) |>
        #remove the GND site
        dplyr::filter(site=="ROOF") |>
        #remove night sessions
        dplyr::filter(!sessionTime==20) |>
        #choose louder channel except for 20230818 sessions
        dplyr::mutate(
                binLoud = ifelse(
                        date == "20230818",
                        binL,
                        ifelse(
                                binL>binR,
                                binL,
                                binR
                        )
                )
        ) |>
        group_by(date,condition,sessionTime) |>
        dplyr::mutate(time=as_datetime(time)) |>
        timetk::summarise_by_time(
                .date_var = time,
                .by = "0.125 second",
                binL = seewave::meandB(binL),
                binR = seewave::meandB(binR),
                binLoud = seewave::meandB(binLoud)
        ) |>
        ungroup() |>
        distinct(date, condition, sessionTime, time, .keep_all = TRUE)

#plot SPL vs Time over 10 min listening period

p_spltime <- ggplot(
        data = spl_vs_time_data_agg,
        aes(
                x = time,
                y = binLoud,
                color = condition,
                fill = condition
        )
        ) +
        #geom_smooth(alpha=0.5, method = "gam") +
        stat_summary(
                fun = seewave::meandB,
                geom = "line",
                linewidth = 0.5
        ) +
        stat_summary(
                fun = seewave::meandB,
                geom = "ribbon",
                alpha = .3,
                linetype = 0,
                fun.max = function(x) seewave::meandB(x) + sd(x) / sqrt(length(x)),
                fun.min = function(x) seewave::meandB(x) - sd(x) / sqrt(length(x))
        ) 
        
p_spltime_save <- p_spltime + 
        scale_color_paletteer_d("Redmonder::qPBI", direction = -1) +
        scale_fill_paletteer_d("Redmonder::qPBI", direction = -1) +
        xlim(0,600) + 
        ylim(60,70) +
        labs(y=expression(italic(L)[AF], dB(A))) +
        xlab("Time (s)")

ggsave(path = "output/",
       plot = p_spltime_save,
       filename = "spltime.pdf",
       device = "pdf",
       units = "px",
       width = 1200,
       height = 600,
       scale = 3.5,
       dpi = 600) 
```