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:

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 in 2.3. Stimuli and automatic masker selection

• Table 3 in 2.4. Non-acoustic environmental conditions for in-situ validation study

• Table 4 in 2.5. Participants

• Table 5 in 3. Results

• Table B.1 in 3. Results

• Figure 3 in 3. Results

• Table 6 in 3.3. Correlation between subjective metrics

• Table B.2 in 3.4. Effect of order, group size and initial conditions

• Table 7 in 4. Results: Objective binaural measurements

• Table B.3 in 4. Results: Objective binaural measurements

• Figure 4 in 4. Results: Objective binaural measurements

Initialisation

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

check if RData file exists then download from the DOI if it doesn't

if (!file.exists("data/fullData.RData")){ asbinary<-dataverse::getfile( 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") ```

2. Method

2.3. Stimuli and automatic masker selection

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.

r 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

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

``` r

list of metrics to compute mean

metriclist <- c("temperature","humidity","lux","windspeed","24h_psi","pm25")

envirosummarytbl <- envirosessiondata |> tblsummary( include = metriclist, by = condition, digits = allcontinuous() ~ 2, type = list(everything() ~ 'continuous'), statistic = list(allcontinuous() ~ "{mean} ({sd})"), missing = "no" ) |> #addn() %>% # add column with total number of non-missing observations modifyheader(label = "Environmental Parameter") |> addp( everything() ~ "wilcox.test", pvaluefun = function(x) stylepvalue(x, digits = 3) ) |> asgt()

envirosummarytbl ```

2.5. Participants

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

``` r

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() }

amssinsituparticipantdata |> dplyr::select(c(condition,gender,age,pss,wnss,wbi), startswith("annoy-")) |> dropna() |> tblsummary( by = condition, type = list(!c(gender) ~ 'continuous'), statistic = list(allcontinuous() ~ "{mean} ({sd})"), digits = list(allcontinuous() ~ 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~" ) ) |> addp( test = list( gender ~ "prop.test", c(age,pss,wnss,wbi, startswith("annoy-")) ~ c("kstest") ), pvaluefun = function(x) stylepvalue(x, digits = 2) ) |> addoverall() |> asgt() |> texttransform( locations = cellsbody(), fn = function(x) { strreplaceall(x,pattern = "@", replacement = "<sub>") |> str</sub>replace_all("~","") } ) ```

3. Results: Site evaluation questionnaire

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.

``` r

compute statistical tests for all variables

stat.results<-twoWLMERMANOVA(metadata,amssinsituparticipant_data) ```

Category is: categorical; Variable is: dom_noise
Two-Way Mixed Effects Repeated Measures Formula: 
rank(dom_noise) ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Category is: categorical; Variable is: dom_natural
Two-Way Mixed Effects Repeated Measures Formula: 
rank(dom_natural) ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Interaction Effect is significant: p=0 (p<0.05) 
Category is: categorical; Variable is: dom_human
Two-Way Mixed Effects Repeated Measures Formula: 
rank(dom_human) ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Category is: interval; Variable is: PosAff
Two-Way Mixed Effects Repeated Measures Formula: 
PosAff ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Interaction Effect is significant: p=0.021 (p<0.05) 
Category is: interval; Variable is: NegAff
Two-Way Mixed Effects Repeated Measures Formula: 
NegAff ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
NegAff is non-normal: p<0.05 (p=3.40444416059014e-14)Two-Way Mixed Effects Repeated Measures Formula: 
rank(NegAff) ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Category is: categorical; Variable is: overall
Two-Way Mixed Effects Repeated Measures Formula: 
rank(overall) ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Interaction Effect is significant: p=0.027 (p<0.05) 
Category is: categorical; Variable is: appropriate
Two-Way Mixed Effects Repeated Measures Formula: 
rank(appropriate) ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Category is: categorical; Variable is: loudness
Two-Way Mixed Effects Repeated Measures Formula: 
rank(loudness) ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Interaction Effect is significant: p=0.022 (p<0.05) 
Category is: interval; Variable is: ISOPL
Two-Way Mixed Effects Repeated Measures Formula: 
ISOPL ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Interaction Effect is significant: p=0.008 (p<0.05) 
Category is: interval; Variable is: ISOEV
Two-Way Mixed Effects Repeated Measures Formula: 
ISOEV ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Category is: interval; Variable is: PRSSFas
Two-Way Mixed Effects Repeated Measures Formula: 
PRSSFas ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Interaction Effect is significant: p=0.008 (p<0.05) 
Category is: interval; Variable is: PRSSBA
Two-Way Mixed Effects Repeated Measures Formula: 
PRSSBA ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Interaction Effect is significant: p=0.019 (p<0.05) 
Category is: interval; Variable is: PRSSCom
Two-Way Mixed Effects Repeated Measures Formula: 
PRSSCom ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Interaction Effect is significant: p=0.046 (p<0.05) 
Category is: interval; Variable is: PRSSEC
Two-Way Mixed Effects Repeated Measures Formula: 
PRSSEC ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>
Category is: interval; Variable is: PRSSES
Two-Way Mixed Effects Repeated Measures Formula: 
PRSSES ~ (1 | pID) + site * condition
<environment: 0x0000018a1081b358>

``` r

retrieve significant results

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

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

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

summarise by condition at each site

summarybycondgt <- amssinsituparticipantdata |> dplyr::filter(!site=="MP") |> dplyr::select(!c(pID,order,partGrp,p:m,pss:who5)) |> tblstrata( strata = site, .tblfun = ~ .x %>% tblsummary( by = condition, missing = "no", type = everything() ~ "continuous", statistic = allcontinuous() ~ "{mean} ({sd})", label = list( domnoise ~ "!DOM#@Noi~", domhuman ~ "!DOM#@Hum~", domnatural ~ "!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#" ) ) %>% adddifference( estimatefun = everything() ~ function(x) paste0(stylesigfig((-x/2) * 100), "%") ), .header = "{strata}" ) |> modifycolumnhide(columns = c(p.value1,ci1,p.value2,ci2)) |> # remove difference footnote modifyfootnote(update = everything() ~ NA) |> modifyheader(allstat_cols() ~ "{level}")

summarise by site for each cond

summarybysitegt <- amssinsituparticipantdata |> dplyr::filter(!site=="MP") |> dplyr::select(!c(pID,order,partGrp,p:m,pss:who5)) |> dplyr::mutate(site=factor(site, levels=c("GND","ROOF"))) |> tblstrata( strata = condition, .tblfun = ~ .x %>% tblsummary( by = site, missing = "no", type = everything() ~ "continuous", statistic = allcontinuous() ~ "{mean} ({sd})", label = list( domnoise ~ "!DOM#@Noi~", domhuman ~ "!DOM#@Hum~", domnatural ~ "!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#" ) ) %>% adddifference( estimatefun = everything() ~ function(x) paste0(stylesigfig((-x/2) * 100), "%") ), .header = "{strata}" ) |> modifycolumnhide(columns = c(p.value1,ci1,p.value2,ci2)) |> # remove difference footnote modifyfootnote(update = everything() ~ NA) |> modifyheader(allstat_cols() ~ "{level}")

merge tables columnwise

tblmergecondsite <- tblmerge( tbls = list(summarybycondgt, summarybysitegt), tabspanner = c( "Contrasts by condition at ...", "Contrasts by site under ..." ) ) |> asgt() |> texttransform( locations = cellsbody(), fn = function(x) { strreplaceall( x,pattern = "@", replacement = "<sub>" ) |> strreplaceall("~","</sub>") |> strreplaceall("!","") |> strreplaceall("#","") } ) |> #highlight significant results at GND site tabstyle( style = list( cellfill(color = "#FDE992"), celltext(weight = "bold") ), locations = cellsbody( columns = c(estimate11), rows = variable %in% signifposthocbycond[ strdetect(signifposthocbycond$term, "GND"), ]$variable ) ) |> #highlight significant results at ROOF site tabstyle( style = list( cellfill(color = "#FDE992"), celltext(weight = "bold") ), locations = cellsbody( columns = c(estimate21), rows = variable %in% signifposthocbycond[ strdetect(signifposthocbycond$term, "ROOF"), ]$variable ) ) |> #highlight significant results under AMB condition tabstyle( style = list( cellfill(color = "#FDE992"), celltext(weight = "bold") ), locations = cellsbody( columns = c(estimate12), rows = variable %in% signifposthocbysite[ strdetect(signifposthocbysite$term, "AMB"), ]$variable ) ) |> #highlight significant results under AMSS condition tabstyle( style = list( cellfill(color = "#FDE992"), celltext(weight = "bold") ), locations = cellsbody( columns = c(estimate22), rows = variable %in% signifposthocbysite[ strdetect(signifposthocbysite$term, "AMSS"), ]$variable ) )

tblmergecond_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.

``` r

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=casewhen( 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::groupby(variable) |> gt::gt() |> cols_label( term ~ html("Term"), test ~ html("Test"), p.value ~ html("p-value"), eff.size ~ html("Effect Size") ) ```

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

``` r

prepare dataframe for plotting

plotdf<- amssinsituparticipantdata |> dplyr::filter(!site=="MP") |> #remove the meeting point pivotlonger( cols = metadata$variable, valuesto = "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

plotsitedf<-plotdf |> groupby(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

plotconddf<-plotdf |> groupby(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]) )

siteplot<-ggplot(data = plotsitedf, aes(x=condition,y=mean,group=Attribute, alpha=significant,color=Attribute)) + geomline() + geompoint() + scalealphadiscrete( range=c(0.2, 1), guide = 'none' #turn off legend ) + scalecolorpaletteerd( palette = "awtools::bpalette", labels=plotlegends ) + facetwrap(vars(site)) + thememinimal() + theme( panel.grid.major.x = elementblank(), panel.grid.minor.x = elementblank() ) + geomhline( yintercept = 0, color="darkgrey", size=1 ) + ylim(c(-1,1)) + ylab("Normalised mean") + ggpubr::labs_pubr()

condplot<-ggplot(data = plotconddf, aes(x=site,y=mean,group=Attribute, alpha=significant,color=Attribute)) + geomline() + geompoint() + scalealphadiscrete( range=c(0.2, 1), guide = 'none' #turn off legend ) + scalecolorpaletteerd( palette = "awtools::bpalette", labels=plotlegends) + facetwrap(vars(condition)) + thememinimal() + theme( panel.grid.major.x = elementblank(), panel.grid.minor.x = elementblank() ) + geomhline( yintercept = 0, color="darkgrey", size=1 ) + ylim(c(-1,1)) + ylab("Normalised mean") + ggpubr::labs_pubr()

combplot<-ggarrange( siteplot, condplot, #labels = c("Site", "Condition"), common.legend = TRUE, legend = "right" ) combplot ```

3.3. Correlation between subjective metrics

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.

3.4. Effect of order, group size and initial conditions

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.

``` r ks.df <- amssinsituparticipantdata |> dplyr::select(c(condition,partGrp,order,metadata$variable)) |> pivotlonger(namesto = "variable", valuesto = "score", cols = metadata$variable)

ks.order.grpsize <- rbind( #KS Test by order ks.df |> dplyr::groupby(condition,variable) |> dplyr::summarise( kstest = list(ks.test(score[order == 18], score[order == 81], exact = NULL, alternative = "two.sided")), ks.pvalue = kstest[[1]]$p.value ) |> dplyr::ungroup() |> dplyr::groupby(condition) |> dplyr::mutate(ks.padj = p.adjust(ks.pvalue, method="BH"), confvar = "order"), #KS Test by group size ks.df |> dplyr::groupby(condition,variable) |> dplyr::summarise( kstest = list(ks.test(score[partGrp == "single"], score[partGrp == "multi"], exact = NULL, alternative = "two.sided")), ks.pvalue = kstest[[1]]$p.value ) |> dplyr::ungroup() |> dplyr::groupby(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(!kstest) |> pivotwider(valuesfrom = ks.padj, namesfrom = variable, idcols = c(condition,confvar)) |> #reorder columns dplyr::select( domnoise, domhuman, domnatural, PosAff, NegAff, overall, appropriate, loudness, ISOPL, ISOEV, PRSSFas, PRSSBA, PRSSCom, PRSSEC, PRSSES, confvar) |> dplyr::groupby(confvar) |> gt::gt() |> colslabel( ISOEV ~ html("ISOEV"), ISOPL ~ html("ISOPL"), NegAff ~ html("NA"), PosAff ~ html("PA"), PRSSFas ~ html("PRSS<sub>Fas</sub>"), PRSSBA ~ html("PRSS<sub>BA</sub>"), PRSSCom ~ html("PRSS<sub>Com</sub>"), PRSSEC ~ html("PRSS<sub>EC</sub>"), PRSSES ~ html("PRSS<sub>ES</sub>"), appropriate ~ html("APPR"), domnatural ~ html("DOM<sub>Nat</sub>"), domhuman ~ html("DOM<sub>Hum</sub>"), dom_noise ~ html("DOM<sub>Noi</sub>"), loudness ~ html("PLN"), overall ~ html("OSQ"), ) ks.order.grpsize ```

4. Results: Objective binaural measurements

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.

``` r

summary including aircraft flyby

isoplobjtbl <- combobjinsitusessiondata |> ungroup() |> pivotwider( namesfrom = "attribute", valuesfrom = "score" ) |> tblstrata2( strata = condition, .tblfun = ~ .x %>% tblsummary(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

isoplobjnoaircrafttbl<-combobjinsitusessiondata |> ungroup() |> pivotwider( namesfrom = "attribute", valuesfrom = "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::groupby(condition,site) |> tblstrata2( strata = condition, .tblfun = ~ .x %>% tblsummary(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(allcontinuous() ~ c(2, 2))), .header = "{strata} (without aircraft flyby)" ) |> modifyheader(label="") |> modifycolumnhide(columns = c(stat11,stat21))

merge tables columnwise

tblmergecondsite <- tblmerge( tbls = list(isoplobjtbl, isoplobjnoaircrafttbl) ) |> modifyspanningheader( c(stat111,stat211) ~ "AMB" ) |> modifyspanningheader( c(stat121,stat221) ~ "AMSS" ) |> modifyspanningheader( c(stat122,stat222) ~ "AMSS (without aircraft flyby)" ) |> asgt() tblmergecond_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.

``r 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")

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

corrobjp <- as.data.frame(corrobj$p) |> dplyr::mutateall( .funs=list(~gtools::stars.pval(.)) )

corrobjmat <- matrix( paste0(as.matrix(corrobjp),as.matrix(corrobjr)), nrow=nrow(corrobjp), dimnames=dimnames(corrobjp) )

remove diagonal values

diag(corrobjmat)=NA

convert to data frame

corrobjmat <- as.data.frame(corrobjmat) |> rownames<-(colnames(corrobjp))

corrobjmat |> gt(rownamestostub = TRUE) |> submissing( columns = everything(), missingtext = "" ) ```

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.

``` r

plot isopl vs decibel score

density colors

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

ISOPL vs LAeq

isoplLAplot <- ggplot(combobjinsitusessiondata) + facetwrap(~attribute, labeller = labelparsed, scales = "free", ncol = 1 ) + # add mean points of each session geompoint(aes(y = L[Aeq], x = score, shape = site, color = pair)) + # add kde contours statdensity2d( bins = 3, contourvar = "ndensity", breaks = c(0.5), geom = "density2d", aes(y = L[Aeq], x = score, color = pair) ) + scalefillmanual(values = densityClr) + scalecolormanual(values = densityClr, name = "subgroup") + ylim(c(50, 70)) + xlim(c(-1, 1)) + geomvline( xintercept = 0, color = "darkgrey", size = 1 ) + #xlab("ISOPL") + ylab(bquote(paste(L[Aeq], ", dB(A)"))) + ggpubr::labspubr() + theme( legend.position = "bottom", strip.background = elementblank(), strip.text.x = elementblank(), axis.title.x = elementblank() )

isoplLCplot <- ggplot(combobjinsitusessiondata) + # add mean points of each session geompoint(aes(y = L[Ceq], x = score, shape = site, color = pair)) + facetwrap(~attribute, labeller = labelparsed, scales = "free", ncol = 1 ) + # add kde contours statdensity2d( bins = 3, contourvar = "ndensity", breaks = c(0.5), geom = "density2d", aes(y = L[Ceq], x = score, color = pair) ) + scalefillmanual(values = densityClr) + scalecolormanual(values = densityClr, name = "subgroup") + ylim(c(60, 80)) + xlim(c(-1, 1)) + geomvline( xintercept = 0, color = "darkgrey", size = 1 ) + ylab(bquote(paste(L[Ceq], ", dB(C)"))) + ggpubr::labspubr() + theme( legend.position = "bottom", strip.background = elementblank(), strip.text.x = elementblank(), axis.title.x = elementblank(), )

isoplN95plot <- ggplot(combobjinsitusessiondata) + facetwrap(~attribute, labeller = labelparsed, scales = "free", ncol = 1, strip.position = "right" ) + # add mean points of each session geompoint(aes(y = N[95], x = score, shape = site, color = pair)) + # add kde contours statdensity2d( bins = 3, contourvar = "ndensity", breaks = c(0.5), geom = "density2d", aes(y = N[95], x = score, color = pair) ) + scalefillmanual(values = densityClr) + scalecolormanual(values = densityClr, name = "subgroup") + ylim(c(5, 20)) + xlim(c(-1, 1)) + geomvline( xintercept = 0, color = "darkgrey", size = 1 ) + ylab(bquote(paste(N[95], ", soneGF"))) + ggpubr::labspubr() + theme( legend.position = "bottom", axis.title.x = elementblank() )

combisoplobjplot <- ggarrange( isoplLAplot, isoplLCplot, isoplN95plot, common.legend = TRUE, legend = "bottom", nrow = 1 ) combisoplobjplot ```