Bhan, Lam
Replication code for “Do neonates hear what we measure? Assessing neonatal ward soundscapes at the neonates’ ears”
The GitHub repository contains the code to replicate the analysis, figures and tables for the paper titled: “Do neonates hear what we measure? Assessing neonatal ward soundscapes at the neonates’ ears”.
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/8GHNGX
The subheadings in this repository follows the headings in the paper (after the Data Preparation section) for consistency.
The following figures and tables are produced by this replication code:
Table 4
in 3.1. HD wardFigure 3
in 4.1. A- and C- weighted metricsTable 5
in 4.2. Acoustic guidelinesFigure 4
in 4.3. Occurence ratesTable A.1
in Appendix A Statistical test resultsTable A.2
in Appendix A Statistical test resultsTable A.3
in Appendix A Statistical test results
Download the dataset from Dataverse if it does not exist. Note that this code repository is configure to ignore the data folder during git commits due to the large size (3.4 GB) of the dataset csv file.
[1] "The dataset exists"
Load the “timeSeries1sec.csv” dataset file containing the measurement data
The bed positions in the wards were shifted according to the schedule below.
From NICU:
- Start to 2022-03-17 11:20: NICU-A
- 2022-03-17 11:20 to 2022-03-21 12:45: NICU-B
- 2022-03-21 12:45 onwards: NICU-A
From HD:
- Start to 2022-04-03 09:00: HD-A
- 2022-04-03 09:00 to 2022-04-03 10:00: HD-B
- 2022-04-03 10:00 onwards: HD-A
Changeover dates to mark out:
- start: 20220303 16:00:00 /
- change: 20220318 16:00:00 to 17:00:00
- change: 20220324 16:00:00 to 17:00:00
timeseries.dt <- timeseries.dt |>
#account for transition period of \pm 15 mins
dplyr::mutate(
sub_location = case_when(
#trim start for both
datetime < ymd_hms(
"2022-03-03 17:00:00", tz = "Singapore"
) ~ "transition",
#NICU-A
datetime < ymd_hms(
"2022-03-17 11:05:00", tz = "Singapore"
) & location == "NICU" ~ "NICU-A",
#30 min gap; NICU-B
datetime > ymd_hms(
"2022-03-17 11:35:00", tz = "Singapore"
) &
datetime < ymd_hms(
"2022-03-21 12:45:00", tz = "Singapore"
) &
location == "NICU" ~ "NICU-B",
#30 min gap; NICU-A and trim end for nicu
datetime > ymd_hms(
"2022-03-21 13:15:00", tz = "Singapore"
) &
datetime < ymd_hms(
"2022-03-24 21:00:00", tz = "Singapore"
) &
location == "NICU" ~ "NICU-A",
#HD-A
datetime < ymd_hms(
"2022-04-03 08:45:00", tz = "Singapore"
) & location == "HD" ~ "HD-A",
#30 min gap; HD-B
datetime > ymd_hms(
"2022-04-03 09:15:00", tz = "Singapore"
) &
datetime < ymd_hms(
"2022-04-03 09:45:00", tz = "Singapore"
) & location == "HD" ~ "HD-B",
#30 min gap: HD-A and trim end for HD
datetime > ymd_hms(
"2022-04-03 10:15:00", tz = "Singapore"
) &
datetime < ymd_hms(
"2022-04-13 15:00:00", tz = "Singapore"
) &
location == "HD" ~ "HD-A",
.default = "transition"
)
) |>
dplyr::filter(!sub_location == "transition")
#metrics to be summarised
params.art.df <- data.frame(
aweight=c("L[AS]","L[ASmax]","L[AS10]","L[AS50]","L[AS90]"),
cweight=c("L[CS]","L[CSmax]","L[CS10]","L[CS50]","L[CS90]"),
tuhms=c("T","T[max]","T[10]","T[50]","T[90]")
)
# Function to calculate metrics
calculate_metrics <- function(data, prefix) {
data |> summarise(
!!sym(prefix) := if(prefix %in% c("LAS", "LCS")) meandB(score_1min) else mean(score_1min),
!!sym(paste0(prefix, "max")) := max(score_1min),
!!sym(paste0(prefix, "10")) := quantile(score_1min, 0.90),
!!sym(paste0(prefix, "50")) := quantile(score_1min, 0.50),
!!sym(paste0(prefix, "90")) := quantile(score_1min, 0.10)
)
}
# Process and combine all metrics
timeseries_1hour_metrics <- timeseries.dt |>
dplyr::mutate(hour = floor_date(datetime,unit = "hour"))
timeseries_1hour_metrics <-
bind_cols(
#a-weighted
timeseries_1hour_metrics |>
dplyr::filter(acoUnit == "dBA") |>
group_by(sub_location, microphone, hour) |>
group_modify(~ calculate_metrics(., "LAS")) |>
ungroup(),
#c-weighted
timeseries_1hour_metrics |>
dplyr::filter(acoUnit == "dBC") |>
group_by(sub_location, microphone, hour) |>
group_modify(~ calculate_metrics(., "LCS")) |>
ungroup() |>
dplyr::select(c("LCS":"LCS90")),
#tonality
timeseries_1hour_metrics |>
dplyr::filter(acoUnit == "tuHMS") |>
group_by(sub_location, microphone, hour) |>
group_modify(~ calculate_metrics(., "T")) |>
ungroup() |>
dplyr::select(c("T":"T90"))
) |>
dplyr::mutate(
microphone = as.factor(microphone),
sub_location = as.factor(sub_location),
`LCS-LAS` = LCS - LAS,
`LAS10-LAS90` = LAS10 - LAS90,
location = ifelse(
grepl("NICU",sub_location),
"NICU",
"HD"
), .before="sub_location"
)
# add week number for plotting
firstHour<-min(timeseries_1hour_metrics$hour)
timeseries_1hour_metrics_wk <- timeseries_1hour_metrics |>
dplyr::mutate(
week=case_when(
difftime(hour, firstHour, units="days")<=7 ~ "Week 1",
difftime(hour, firstHour, units="days")<=14 ~ "Week 2",
difftime(hour, firstHour, units="days")<=21 ~ "Week 3",
difftime(hour, firstHour, units="days")<=28 ~ "Week 4",
difftime(hour, firstHour, units="days")<=35 ~ "Week 5",
difftime(hour, firstHour, units="days")<=42 ~ "Week 6"
)) |>
group_by(location,microphone) |>
arrange(hour) |>
ungroup()
# convert to dataframe to long form
timeseries_1hour_metrics_wk_long <- timeseries_1hour_metrics_wk |>
dplyr::filter(microphone %in% c("binL","binR")) |>
pivot_longer(cols = c(LAS:T90,`LCS-LAS`,`LAS10-LAS90`),
names_to = "acoUnit_stat",
values_to = "score") |>
pivot_wider(
names_from = microphone,
values_from = score
) |>
#find binaural average
dplyr::mutate(
bin_LR = ifelse(
grepl("^L", acoUnit_stat),
10*log10((10^(binL/10)+10^(binR/10))/2),
(binL+binR)/2
),
hour=factor(hour(hour), levels=c(0:23))
)
#1hour aggregated dBA metrics for both locations
timeseries_1hour_metrics_fixed <- timeseries_1hour_metrics |>
dplyr::filter(
hour >= ymd_hms("20220303 17:00:00", tz ="Singapore") &
hour <= ymd_hms("20220317 10:30:00",
tz ="Singapore")
) |>
dplyr::mutate(
location = as.factor(location)
)
Table 4
: Mean and standard deviation of the differences in 1-h slow time-weighted metrics at the HD-A and NICU-A bed positions measured from 03/03/2022 17:00:00 to 17/03/2022 10:30:00. Difference pairs with significant ART contrasts are indicated in bold.
#compute differences separately for each location due to the different
#number of microphones
#hd
diff.summary.1h.hd <- timeseries_1hour_metrics_fixed |>
dplyr::select(!c(LAS90,LCS90:T90)) |>
dplyr::filter(location=="HD") |>
pivot_longer(
names_to = "acoUnit",
values_to = "score",
cols = c(
LAS,LAS10,LAS50,LASmax,
LCS,LCS10,LCS50,LCSmax,
`LCS-LAS`,`LAS10-LAS90`)
) |>
pivot_wider(names_from = microphone, values_from = score) |>
# mutate(`binL-binR`=abs(binL-binR),
# `binL-146AE`=abs(binL-`146AE`),
# `binR-146AE`=abs(binR-`146AE`)) %>%
mutate(`binL-binR`=binL-binR,
`binL-146AE`=binL-`146AE`,
`binR-146AE`=binR-`146AE`) |>
pivot_longer(names_to = "microphone",
values_to = "score",
cols = c(`146AE`:`binR-146AE`)) |>
#summary by acounit and microphone
dplyr::group_by(microphone,acoUnit) |>
summarise(
mean=round(mean(score),3),
sd=sd(score),
score=list(score)
) |>
dplyr::mutate(location="HD",.before = microphone)
#nicu
diff.summary.1h.nicu <- timeseries_1hour_metrics_fixed |>
dplyr::filter(location=="NICU") |>
dplyr::select(!c(LAS90,LCS90:T90)) |>
pivot_longer(
names_to = "acoUnit",
values_to = "score",
cols = c(
LAS,LAS10,LAS50,LASmax,
LCS,LCS10,LCS50,LCSmax,
`LCS-LAS`,`LAS10-LAS90`)
) |>
pivot_wider(names_from = microphone, values_from = score) |>
mutate(
`binL-binR`=binL-binR,
`binL-146AEIn`=binL-`146AEIn`,
`binL-146AEOut`=binL-`146AEOut`,
`binR-146AEIn`=binR-`146AEIn`,
`binR-146AEOut`=binR-`146AEOut`,
`146AEIn-146AEOut`=`146AEIn`-`146AEOut`
) |>
pivot_longer(
names_to = "microphone",
values_to = "score",
cols = c(`146AEIn`:`146AEIn-146AEOut`)
) |>
#summary by acounit and microphone
dplyr::group_by(microphone,acoUnit) |>
summarise(mean=round(mean(score),3),
sd=sd(score),
score=list(score)) |>
dplyr::mutate(location="NICU",.before = microphone)
# Define a function to apply md() to each column label
md_label <- function(label) {
md(label)
}
#plot table of differences across all locations
gt.diff.summary.1h <- rbind(diff.summary.1h.hd,diff.summary.1h.nicu) |>
dplyr::filter(!microphone %in% c("146AE","binL","binR",
"146AEIn","146AEOut")) |>
dplyr::mutate(mean_sd=paste0(formatC(mean,
format = "f",
digits = 2),
" (",
formatC(sd,
format = "f",
digits = 2),
")")) |>
pivot_wider(
id_cols = acoUnit,
names_from = c(location,microphone),
values_from = mean_sd
) |>
gt::gt() |>
tab_spanner(
label = "HD",
columns = starts_with("HD")
) |>
tab_spanner(
label = "NICU",
columns = starts_with("NICU")
) |>
cols_label_with(
fn = ~ gsub("HD_|NICU_","", .)
) |>
cols_label_with(
fn = ~ gsub("146AE ","mout ",.)
) |>
cols_label_with(
fn = ~ gsub("bin|146AE","m",.)
) |>
cols_label_with(
fn = ~ gsub("(In|Out)", "\\L\\1", ., perl = TRUE)
) |>
cols_label_with(
fn = ~ gsub("(R|L|out|in)", "<sub>\\1</sub>", .,
perl = TRUE)
) |>
cols_label_with (
fn = ~ md_label(.)
)
gt.diff.summary.1h |> gtsave(filename = "diff_1h_meansd.tex",
path = "output/")
gt.diff.summary.1h
acoUnit |
HD
|
NICU
|
|||||||
---|---|---|---|---|---|---|---|---|---|
mL-m | mL-mR | mR-m | min-mout | mL-min | mL-mout | mL-mR | mR-min | mR-mout | |
LAS | -1.35 (1.18) | -0.07 (0.82) | -1.27 (0.79) | -0.38 (0.40) | 3.02 (0.54) | 2.64 (0.50) | 0.04 (0.55) | 2.98 (0.38) | 2.60 (0.51) |
LAS10 | -1.62 (0.97) | -0.17 (0.50) | -1.45 (0.87) | -0.47 (0.46) | 2.90 (0.56) | 2.44 (0.54) | -0.00 (0.65) | 2.91 (0.46) | 2.44 (0.61) |
LAS10-LAS90 | -0.76 (0.77) | -0.03 (0.51) | -0.73 (0.69) | -0.18 (0.44) | -0.35 (0.47) | -0.52 (0.44) | -0.11 (0.64) | -0.23 (0.53) | -0.41 (0.63) |
LAS50 | -1.24 (0.71) | -0.20 (0.24) | -1.04 (0.66) | -0.36 (0.42) | 3.29 (0.58) | 2.93 (0.57) | 0.11 (0.70) | 3.17 (0.43) | 2.81 (0.56) |
LASmax | -1.27 (2.77) | 0.04 (2.13) | -1.31 (1.83) | -0.60 (1.20) | 3.04 (1.62) | 2.44 (1.61) | 0.55 (1.46) | 2.49 (1.39) | 1.89 (1.61) |
LCS | -0.29 (0.64) | -0.96 (0.55) | 0.67 (0.36) | -0.35 (0.07) | 1.09 (0.19) | 0.73 (0.20) | 0.48 (0.16) | 0.60 (0.25) | 0.25 (0.27) |
LCS-LAS | 1.05 (0.78) | -0.89 (0.43) | 1.94 (0.65) | 0.03 (0.35) | -1.93 (0.53) | -1.90 (0.55) | 0.44 (0.42) | -2.37 (0.45) | -2.35 (0.57) |
LCS10 | -0.52 (0.44) | -0.92 (0.16) | 0.40 (0.54) | -0.33 (0.11) | 1.28 (0.35) | 0.95 (0.38) | 0.41 (0.24) | 0.87 (0.46) | 0.54 (0.50) |
LCS50 | -0.16 (0.16) | -1.08 (0.06) | 0.92 (0.19) | -0.36 (0.06) | 0.94 (0.13) | 0.57 (0.14) | 0.54 (0.15) | 0.40 (0.17) | 0.04 (0.19) |
LCSmax | -1.16 (2.94) | -0.28 (1.76) | -0.88 (1.93) | -0.39 (0.79) | 2.19 (0.84) | 1.81 (0.98) | 0.31 (0.71) | 1.88 (0.87) | 1.49 (1.14) |
Figure 3
: A- and C-weighted decibel metrics averaged by hour of the day across the entire measurement duration at NICU-A, NICU-B and HD-A measurement points.
# define labels for metrics
acoUnit_stat_labs <- c("L[AS]","L[ASmax]","L[AS10]","L[AS50]","L[AS90]",
"L[CS]","L[CSmax]","L[CS10]","L[CS50]","L[CS90]",
"T","T[max]","T[10]","T[50]","T[90]",
"L[CS]-L[AS]","L[AS10]-L[AS90]")
names(acoUnit_stat_labs) <- unique(timeseries_1hour_metrics_wk_long$acoUnit_stat)
# define plot customisation parameters
acoUnit.plotlist<-list()
acoMetric.units<-c("dB(A)","dB(C)","tu[HMS]")
ylim.list<-list(c(65,85),
c(65,85),
c(0,2.5))
acoMetric_list<-c("LA","LC","T")
# generate plot
idx<-0
for (acoMetric in acoMetric_list) {
idx<-idx+1
acoUnit.plotlist[[idx]] <- ggplot(
data = timeseries_1hour_metrics_wk_long |>
dplyr::filter(
grepl(paste0("^",acoMetric), acoUnit_stat) &
!acoUnit_stat %in%
c(
"LCS-LAS",
"LAS10-LAS90",
"LAS90","LCS90","T90"
) &
!sub_location == "HD-B"
),
aes(
x = as.numeric(hour),
y = bin_LR,
color = sub_location,
fill = sub_location),
) +
stat_summary(
fun = mean,
geom = "line",
linewidth = 1
) +
stat_summary(
fun = mean,
geom = "ribbon",
alpha = .3,
#fill = "#EB5286",
fun.max = function(x) mean(x) + sd(x) / sqrt(length(x)),
fun.min = function(x) mean(x) - sd(x) / sqrt(length(x))
) +
scale_color_paletteer_d("Redmonder::qPBI", direction = -1,
name = "Location"
) +
scale_fill_paletteer_d("Redmonder::qPBI", direction = -1,
name = "Location"
) +
facet_wrap(
~acoUnit_stat, ncol = 6,
labeller = labeller(
#acoUnit=acoUnit.labs,
acoUnit_stat=as_labeller(
acoUnit_stat_labs,label_parsed))) +
ggthemes::theme_hc() +
scale_x_continuous(
breaks = seq(0, 24, by = 6),
minor_breaks = seq(0, 24, 2)
) +
theme(panel.grid.minor.y = element_line(color = 1,
size = 0.1,
linetype = 3)) +
theme(panel.grid.major.x = element_line(color = 1,
size = 0.1,
linetype = 1)) +
theme(panel.grid.minor.x = element_line(color = 1,
size = 0.1,
linetype = 3)) +
#xlim(0,23) +
#ylim(ylim.list[[idx]])+
xlab("Hour") +
labs(
y=parse(text = acoMetric.units[idx])
)
if (acoMetric == "LA" ) {
acoUnit.plotlist[[idx]] <- acoUnit.plotlist[[idx]] +
#65 dB limit for LAS10, CC 9th ed
geom_hline(
data = data.frame(yint=65, acoUnit_stat = "LAS10"),
aes(yintercept = yint),
color = "maroon",
linetype = "dashed"
) +
#50 dB limit for LAS50, CC 9th ed
geom_hline(
data = data.frame(yint=50, acoUnit_stat = "LAS50"),
aes(yintercept = yint),
color = "maroon",
linetype = "dashed"
) +
#45 dB limit for LAS, 8th ed
geom_hline(
data = data.frame(yint=45, acoUnit_stat = "LAS"),
aes(yintercept = yint),
color = "maroon",
linetype = "dashed"
) +
#65 dB limit for LASmax, 8th ed
geom_hline(
data = data.frame(yint=65, acoUnit_stat = "LASmax"),
aes(yintercept = yint),
color = "maroon",
linetype = "dashed"
)
}
}
#Combine plot
p.all <- (acoUnit.plotlist[[1]] +
theme(legend.position = "none") +
xlab("")) /
# (acoUnit.plotlist[[2]] + theme(legend.position = "none") +
# xlab("")) /
acoUnit.plotlist[[2]]
#display plot
p.all
#save plot
ggsave(path = "output/",
plot = p.all,
filename = "fig3_acmetrics_24h.pdf",
device = "pdf",
units = "px",
width = 1600,
height = 1500,
scale = 3.5,
dpi = 600)
Table 5:
Summary of mean A-weighted metrics and percentage of time where the metrics were within guidelines over $N=981$ , $N=392$ , and $N=98$ 1-h periods at HD-A, NICU-A and NICU-B bed positions, respectively.
timeseries_1hour_CC_regulations <- timeseries_1hour_metrics |>
dplyr::mutate(CC9_LAS10=ifelse(LAS10<=65,"Yes","No"),
CC9_LAS50=ifelse(LAS50<=50,"Yes","No"),
CC8_LASmax=ifelse(LASmax<=65,"Yes","No"),
CC8_LAS=ifelse(LAS<=45,"Yes","No")) |>
dplyr::mutate(across(c(CC9_LAS10:CC8_LAS),
~factor(.,levels=c("Yes","No"))))
gt_subloc_reg_summary <- timeseries_1hour_CC_regulations |>
dplyr::select(
sub_location,microphone,LAS,LAS10,LAS50,LASmax,
CC9_LAS10,CC9_LAS50,CC8_LASmax,CC8_LAS
) |>
dplyr::filter(!sub_location=="HD-B") |>
gtsummary::tbl_strata(
strata = sub_location,
.tbl_fun =
~ .x |>
dplyr::mutate(microphone = droplevels(microphone)) |>
gtsummary::tbl_summary(
by = microphone,
type = list(c(CC9_LAS10:CC8_LAS) ~ "categorical",
c(LAS:LASmax) ~ "continuous"),
statistic = list(c(LAS:LASmax) ~ "{mean} ({sd})"),
digits = c(LAS:LASmax) ~ 2
),
.header = "**{strata}**"
)
gt_subloc_reg_summary
Characteristic |
HD-A
|
NICU-A
|
NICU-B
|
||||||||
---|---|---|---|---|---|---|---|---|---|---|---|
146AE N = 9811 |
binL N = 9811 |
binR N = 9811 |
146AEIn N = 3921 |
146AEOut N = 3921 |
binL N = 3921 |
binR N = 3921 |
146AEIn N = 981 |
146AEOut N = 981 |
binL N = 981 |
binR N = 981 |
|
LAS | 57.18 (3.64) | 56.05 (3.87) | 56.09 (3.48) | 56.45 (1.93) | 56.78 (2.09) | 59.33 (1.90) | 59.38 (1.90) | 55.05 (1.88) | 55.87 (1.78) | 57.47 (2.01) | 58.29 (1.81) |
LAS10 | 60.25 (4.09) | 58.74 (3.85) | 58.98 (3.87) | 58.54 (2.73) | 58.96 (2.86) | 61.33 (2.65) | 61.40 (2.70) | 56.73 (2.54) | 57.55 (2.33) | 59.15 (2.74) | 59.89 (2.27) |
LAS50 | 52.48 (3.71) | 51.35 (3.30) | 51.59 (3.36) | 54.53 (1.43) | 54.79 (1.55) | 57.61 (1.43) | 57.63 (1.38) | 53.47 (1.10) | 54.60 (1.27) | 55.97 (1.44) | 56.80 (1.15) |
LASmax | 74.85 (4.14) | 73.97 (5.28) | 73.76 (4.20) | 70.96 (4.00) | 71.55 (3.99) | 73.94 (4.21) | 73.40 (4.12) | 68.44 (4.85) | 68.78 (4.92) | 70.71 (5.05) | 71.47 (4.81) |
CC9_LAS10 | |||||||||||
Yes | 877 (89%) | 953 (97%) | 942 (96%) | 387 (99%) | 386 (98%) | 360 (92%) | 353 (90%) | 98 (100%) | 98 (100%) | 97 (99%) | 97 (99%) |
No | 104 (11%) | 28 (2.9%) | 39 (4.0%) | 5 (1.3%) | 6 (1.5%) | 32 (8.2%) | 39 (9.9%) | 0 (0%) | 0 (0%) | 1 (1.0%) | 1 (1.0%) |
CC9_LAS50 | |||||||||||
Yes | 309 (31%) | 389 (40%) | 368 (38%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) |
No | 672 (69%) | 592 (60%) | 613 (62%) | 392 (100%) | 392 (100%) | 392 (100%) | 392 (100%) | 98 (100%) | 98 (100%) | 98 (100%) | 98 (100%) |
CC8_LASmax | |||||||||||
Yes | 8 (0.8%) | 21 (2.1%) | 19 (1.9%) | 24 (6.1%) | 14 (3.6%) | 6 (1.5%) | 4 (1.0%) | 24 (24%) | 22 (22%) | 12 (12%) | 10 (10%) |
No | 973 (99%) | 960 (98%) | 962 (98%) | 368 (94%) | 378 (96%) | 386 (98%) | 388 (99%) | 74 (76%) | 76 (78%) | 86 (88%) | 88 (90%) |
CC8_LAS | |||||||||||
Yes | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) |
No | 981 (100%) | 981 (100%) | 981 (100%) | 392 (100%) | 392 (100%) | 392 (100%) | 392 (100%) | 98 (100%) | 98 (100%) | 98 (100%) | 98 (100%) |
1 Mean (SD); n (%) |
gt_subloc_reg_summary |>
gtsummary::as_gt() |>
gt::gtsave(
filename = paste0(
"regulation.summary.tbl.subloc.tex"
),
path = "./output/"
)
Figure 4
: Occurrence rate of $\textit{OR}^h_\text{SNR}(5)$ and $\textit{OR}_{T}^h(0.4)$ averaged over the same daily 1-h period throughout the entire measurement campaign. A-weighted decibel metrics and tonality metrics averaged by hour of the day across the entire measurement duration at NICU-A, NICU-B and HD-A measurement points.
SNR = 5 #SNR threshold for dBA occurence rate
#compute occurrence rate
OR_dBA_tuHMS <- timeseries.dt |>
#filter only binaural and dBA and tuHMS metrics
dplyr::filter(
acoUnit %in% c("dBA","tuHMS") &
microphone %in% c("binL","binR") #&
# datetime < ymd_hms("2022-03-05 23:00:00",
# tz = "Singapore") #for testing
) |>
#prepare data frame
pivot_wider(
names_from = microphone,
values_from = score_1min,
values_fn = ~ mean(.x, na.rm = TRUE) #to handle duplicates
) |>
#aggregate both binL and binR into binLR
dplyr::mutate(
binLR = ifelse(
acoUnit == "dBA",
10*log10((10^(binL/10)+10^(binR/10))/2),
(binL+binR)/2
)
) |>
#drop binL and binR columns
dplyr::select(!c(binL,binR)) |>
#filter only required sub locations
dplyr::filter(sub_location %in% c("NICU-A","NICU-B","HD-A")) |>
#factorise sub_location, create hour and minute column
dplyr::mutate(
sub_location = as.factor(sub_location),
minute = lubridate::floor_date(datetime, unit = "minute"),
hour = lubridate::floor_date(datetime, unit = "hour")
) |>
#aggregate by minute
group_by(sub_location, acoUnit, minute) |>
dplyr::mutate(
#max
binLR_1min_max = max(binLR)
) |>
ungroup() |>
distinct() |>
#find 1h average as the noise floor
dplyr::group_by(sub_location,acoUnit,hour) |>
dplyr::mutate(
# 50% exceedance level
binLR_50 = quantile(binLR, 0.5)
) |>
ungroup() |>
distinct () |>
summarise(
.by = c(sub_location,acoUnit,hour),
count = n(), #number of events in an hour
OR = ifelse(
acoUnit == "dBA", #if dBA
#find the number of events > SNR of 5 dBA
sum(
binLR_1min_max > (binLR_50 + SNR),
na.rm = TRUE
)*(100/count), #1min max SPL > 1h SPL
sum(binLR_1min_max > 0.4)*(100/count) #tonality > 0.4
)
) |>
distinct() |>
#add column to aggregate hour and convert to 0 to 23 range
dplyr::mutate(hour_only = as.numeric(hour(hour)))
or_plot <-
ggplot(
OR_dBA_tuHMS,
aes(x = as.numeric(hour_only), #convert to numeric for plot
y = OR,
color = sub_location,
fill = sub_location)
) +
facet_wrap(
~acoUnit,
nrow = 2,
scales = "free",
labeller = labeller(
acoUnit = as_labeller(
c("dBA" = "italic(L)[AS]", "tuHMS" = "italic(T)"),
label_parsed)
)
) +
stat_summary(
fun = mean,
geom = "line",
linewidth = 1
) +
stat_summary(
fun = mean,
geom = "ribbon",
alpha = .3,
linetype = 0,
fun.max = function(x) mean(x) + sd(x) / sqrt(length(x)),
fun.min = function(x) mean(x) - sd(x) / sqrt(length(x))
) +
scale_color_paletteer_d("Redmonder::qPBI", direction = -1) +
scale_fill_paletteer_d("Redmonder::qPBI", direction = -1) +
xlim(0,24) + ylim(0,100) +
ylab("% of time over threshold") +
xlab("Hour") +
scale_x_continuous(
breaks = seq(0, 24, by = 6),
minor_breaks = seq(0, 24, 2)
) +
ggthemes::theme_hc() +
theme(panel.grid.minor.y = element_line(color = 1,
size = 0.1,
linetype = 3)) +
theme(panel.grid.major.x = element_line(color = 1,
size = 0.1,
linetype = 1)) +
theme(panel.grid.minor.x = element_line(color = 1,
size = 0.1,
linetype = 3))
or_plot
ggsave(path = "output/",
plot = or_plot,
filename = "fig4_or_dBA_tuHMS.pdf",
device = "pdf",
units = "px",
width = 1600,
height = 1500,
#scale = 3.5,
scale = 2.0,
dpi = 600)
Table A.1
: Summary of LME-ART-ANOVA and posthoc contrast tests with microphone type as the fixed effect, and 1-h time periods as the random effect for each acoustic metric at the HD ward.
#ART ANOVA repeated measures between microphones at each site
#initialise data frame
#data frame to store main effects of 1-W RM ART ANOVA
art_main_mic_df <- data.frame(
term = as.character(),
test = as.character(),
location = as.character(),
acoUnit = as.character(),
p.value = as.numeric(),
eff.size = as.numeric()
)
metric_list <- c("LAS","LAS10","LAS50","LASmax",
"LCS","LCS10","LCS50","LCSmax",
"LCS-LAS","LAS10-LAS90")
#iterate through all the metrics for 1W RM ART ANOVA main effects
for (metricIdx in metric_list) {
#for each location (each location has different number of microphones)
for (loc in unique(timeseries_1hour_metrics_fixed$location)) {
#ART ANOVA
model = art(
as.formula(
paste0("`",metricIdx,"` ~ microphone + (1|hour)")
),
data = timeseries_1hour_metrics_fixed |>
dplyr::filter(location==loc) #|>
#dplyr::select(!c(LAS90,LCS:LCS90))
)
#effect size partial eta squared
Result.model = anova(model)
model_eff<-omega_squared(model)
art_main_mic_df <- rbind(
art_main_mic_df,
data.frame(
term=Result.model$Term,
test="LME-ART-ANOVA",
location=loc,
acoUnit=metricIdx,
p.value=Result.model$`Pr(>F)`,
eff.size=model_eff$Omega2_partial)
)
#contrast tests on microphones
marginal <- art.con(model, "microphone")
art_main_mic_df <- rbind(
art_main_mic_df,
data.frame(
term=data.frame(marginal)$contrast,
test="ART Contrasts",
location=loc,
acoUnit=metricIdx,
p.value=data.frame(marginal)$p.value,
eff.size=NA
)
)
}
}
#}
#add significance symbols
art_main_mic_df<-art_main_mic_df |>
dplyr::mutate(
#add signif symbol
sig = symnum(
p.value, corr = FALSE, na = FALSE,
cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", " ")
),
#add effect size symbol
eff = symnum(
abs(eff.size), corr = FALSE, na = FALSE,
cutpoints = c(0, 0.01, 0.06, 0.14, 1),
symbols = c(" ", "S", "M", "L")
),
pvalue_sig = paste0(
sig,
formatC(p.value, format = "f", digits=4)
),
eff_sym = ifelse(
is.na(eff.size),
"",
paste0(
"(", eff, ")",
formatC(eff.size, format = "f", digits=2)
)
)
)
Plot summary of LME-ART-ANOVA and postboc contrast tests for each metric at HD ward
gt_art_mic_HD <- art_main_mic_df |>
dplyr::filter(location == "HD") |>
dplyr::select(!c("p.value","eff.size","sig","eff","location")) |>
gt(
groupname_col = "acoUnit"
)
gt_art_mic_HD |>
gtsave(filename = "art_mic_HD.tex",
path = "output/")
gt_art_mic_HD
term | test | pvalue_sig | eff_sym |
---|---|---|---|
LAS | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.57 |
146AE - binL | ART Contrasts | ****0.0000 | |
146AE - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | 0.1118 | |
LAS10 | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.67 |
146AE - binL | ART Contrasts | ****0.0000 | |
146AE - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | ***0.0002 | |
LAS50 | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.74 |
146AE - binL | ART Contrasts | ****0.0000 | |
146AE - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | ****0.0000 | |
LASmax | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.25 |
146AE - binL | ART Contrasts | ****0.0000 | |
146AE - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | 0.7208 | |
LCS | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.83 |
146AE - binL | ART Contrasts | ****0.0000 | |
146AE - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | ****0.0000 | |
LCS10 | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.75 |
146AE - binL | ART Contrasts | ****0.0000 | |
146AE - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | ****0.0000 | |
LCS50 | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.92 |
146AE - binL | ART Contrasts | ****0.0000 | |
146AE - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | ****0.0000 | |
LCSmax | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.24 |
146AE - binL | ART Contrasts | ****0.0000 | |
146AE - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | ****0.0000 | |
LCS-LAS | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.75 |
146AE - binL | ART Contrasts | ****0.0000 | |
146AE - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | ****0.0000 | |
LAS10-LAS90 | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.48 |
146AE - binL | ART Contrasts | ****0.0000 | |
146AE - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | 0.2912 |
Table A.2
: Summary of LME-ART-ANOVA and posthoc contrast tests with microphone type as the fixed effect, and 1-h time periods as the random effect for each acoustic metric at the NICU ward
gt_art_mic_NICU <- art_main_mic_df |>
dplyr::filter(location == "NICU") |>
dplyr::select(!c("p.value","eff.size","sig","eff","location")) |>
gt(
groupname_col = "acoUnit"
)
gt_art_mic_NICU |>
gtsave(filename = "art_mic_NICU.tex",
path = "output/")
gt_art_mic_NICU
term | test | pvalue_sig | eff_sym |
---|---|---|---|
LAS | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.91 |
146AEIn - 146AEOut | ART Contrasts | ****0.0000 | |
146AEIn - binL | ART Contrasts | ****0.0000 | |
146AEIn - binR | ART Contrasts | ****0.0000 | |
146AEOut - binL | ART Contrasts | ****0.0000 | |
146AEOut - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | 0.3211 | |
LAS10 | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.89 |
146AEIn - 146AEOut | ART Contrasts | ****0.0000 | |
146AEIn - binL | ART Contrasts | ****0.0000 | |
146AEIn - binR | ART Contrasts | ****0.0000 | |
146AEOut - binL | ART Contrasts | ****0.0000 | |
146AEOut - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | 0.8876 | |
LAS50 | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.91 |
146AEIn - 146AEOut | ART Contrasts | ****0.0000 | |
146AEIn - binL | ART Contrasts | ****0.0000 | |
146AEIn - binR | ART Contrasts | ****0.0000 | |
146AEOut - binL | ART Contrasts | ****0.0000 | |
146AEOut - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | **0.0022 | |
LASmax | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.64 |
146AEIn - 146AEOut | ART Contrasts | ****0.0000 | |
146AEIn - binL | ART Contrasts | ****0.0000 | |
146AEIn - binR | ART Contrasts | ****0.0000 | |
146AEOut - binL | ART Contrasts | ****0.0000 | |
146AEOut - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | ****0.0000 | |
LCS | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.89 |
146AEIn - 146AEOut | ART Contrasts | ****0.0000 | |
146AEIn - binL | ART Contrasts | ****0.0000 | |
146AEIn - binR | ART Contrasts | ****0.0000 | |
146AEOut - binL | ART Contrasts | ****0.0000 | |
146AEOut - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | ****0.0000 | |
LCS10 | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.86 |
146AEIn - 146AEOut | ART Contrasts | ****0.0000 | |
146AEIn - binL | ART Contrasts | ****0.0000 | |
146AEIn - binR | ART Contrasts | ****0.0000 | |
146AEOut - binL | ART Contrasts | ****0.0000 | |
146AEOut - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | ****0.0000 | |
LCS50 | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.90 |
146AEIn - 146AEOut | ART Contrasts | ****0.0000 | |
146AEIn - binL | ART Contrasts | ****0.0000 | |
146AEIn - binR | ART Contrasts | ****0.0000 | |
146AEOut - binL | ART Contrasts | ****0.0000 | |
146AEOut - binR | ART Contrasts | 0.9282 | |
binL - binR | ART Contrasts | ****0.0000 | |
LCSmax | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.70 |
146AEIn - 146AEOut | ART Contrasts | ****0.0000 | |
146AEIn - binL | ART Contrasts | ****0.0000 | |
146AEIn - binR | ART Contrasts | ****0.0000 | |
146AEOut - binL | ART Contrasts | ****0.0000 | |
146AEOut - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | ****0.0000 | |
LCS-LAS | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.91 |
146AEIn - 146AEOut | ART Contrasts | 0.3148 | |
146AEIn - binL | ART Contrasts | ****0.0000 | |
146AEIn - binR | ART Contrasts | ****0.0000 | |
146AEOut - binL | ART Contrasts | ****0.0000 | |
146AEOut - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | ****0.0000 | |
LAS10-LAS90 | |||
microphone | LME-ART-ANOVA | ****0.0000 | (L)0.27 |
146AEIn - 146AEOut | ART Contrasts | ****0.0000 | |
146AEIn - binL | ART Contrasts | ****0.0000 | |
146AEIn - binR | ART Contrasts | ****0.0000 | |
146AEOut - binL | ART Contrasts | ****0.0000 | |
146AEOut - binR | ART Contrasts | ****0.0000 | |
binL - binR | ART Contrasts | **0.0045 |
Table A.3
: Summary of LME-ART-ANOVA and posthoc contrast tests with bed position as the fixed effect, and 1-h time periods and microphone type as the random effects for each acoustic metric.
The difference between sites are investigated independently at each of the NICU and HD wards. In the NICU, differences were investigated between NICU-A and NICU-B; and between HD-A and HD-B in the HD ward.
metric_list <- c(
"LAS","LAS10","LAS50","LASmax",
"LCS","LCS10","LCS50","LCSmax"
)
#data frame to store main effects of 1-W RM ART ANOVA
art_main_subloc_df <- data.frame(
acoUnit=as.character(),
term=as.character(),
test=as.character(),
p.value=as.numeric(),
eff.size=as.numeric()
)
# Create a progress bar
pb <- progress_bar$new(
format = " Processing metrics [:bar] :percent (:current/:total) :eta",
total = length(metric_list),
clear = FALSE,
width = 60
)
#iterate through all the metrics for 1W RM ART ANOVA main effects
#linear mixed-effects model where the site is a fixed
#effect, and time and microphone positions are random effects
for (metricIdx in metric_list) {
# update the progress bar
pb$tick()
data_filtered <- timeseries_1hour_metrics |>
dplyr::filter(
sub_location %in% c(
"HD-A","NICU-A","NICU-B"
) &
microphone %in% c("binL","binR")
)
m.art = art(
as.formula(paste0(
"`",metricIdx,
"` ~ sub_location + (1|hour) + (1|microphone)"
)),
data = data_filtered
)
#effect size partial eta squared
Result.model = anova(m.art)
model_eff<-omega_squared(m.art)
#add main effects test results
art_main_subloc_df <- rbind(
art_main_subloc_df,
data.frame(
#location=loc,
acoUnit=metricIdx,
test="LME-ART-ANOVA",
term=Result.model$Term,
p.value=Result.model$`Pr(>F)`,
eff.size=model_eff$Omega2_partial)
)
#if main effect p < 0.05 run contrast tests
if (Result.model$`Pr(>F)` < 0.05) {
con_test_sum <- summary(
art.con(
m.art,
"sub_location",
adjust = "bonferroni"
)
)
#compute cohen's d effect siz with art model
m.art.subloc <- artlm(m.art, "sub_location")
con_test_sum$d = con_test_sum$estimate / car::sigmaHat(m.art.subloc)
art_main_subloc_df <- rbind(
art_main_subloc_df,
data.frame(
acoUnit=metricIdx,
test="ART Contrasts",
term=con_test_sum$contrast,
p.value=con_test_sum$p.value,
eff.size=con_test_sum$d)
)
}
}
#add significance stars
art_main_subloc_df <- art_main_subloc_df |>
mutate(
#add signif symbol
sig = symnum(
p.value, corr = FALSE, na = FALSE,
cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", " ")
),
#add effect size symbol
eff = symnum(
abs(eff.size), corr = FALSE, na = FALSE,
cutpoints = c(0, 0.01, 0.06, 0.14, 10000),
symbols = c(" ", "S", "M", "L")
),
pvalue_sig=paste0(
sig,
formatC(p.value, format = "f", digits=4)
),
eff_sym=ifelse(
eff==" ",
formatC(eff.size, format = "f", digits=2),
paste0(
"(", eff, ")",
formatC(eff.size, format = "f", digits=2)
)
)
) |>
dplyr::select(!c("p.value","eff.size","sig","eff"))
art_main_subloc_df |>
gt() |>
gtsave(filename = "art_subloc.tex",
path = "output/")
art_main_subloc_df |>
gt()
acoUnit | test | term | pvalue_sig | eff_sym |
---|---|---|---|---|
LAS | LME-ART-ANOVA | sub_location | ****0.0000 | (L)0.47 |
LAS | ART Contrasts | (HD-A) - (NICU-A) | ****0.0000 | (L)-2.16 |
LAS | ART Contrasts | (HD-A) - (NICU-B) | ****0.0000 | (L)-0.92 |
LAS | ART Contrasts | (NICU-A) - (NICU-B) | ****0.0000 | (L)1.24 |
LAS10 | LME-ART-ANOVA | sub_location | ****0.0000 | (L)0.25 |
LAS10 | ART Contrasts | (HD-A) - (NICU-A) | ****0.0000 | (L)-1.34 |
LAS10 | ART Contrasts | (HD-A) - (NICU-B) | 1.0000 | (S)0.02 |
LAS10 | ART Contrasts | (NICU-A) - (NICU-B) | ****0.0000 | (L)1.36 |
LAS50 | LME-ART-ANOVA | sub_location | ****0.0000 | (L)0.81 |
LAS50 | ART Contrasts | (HD-A) - (NICU-A) | ****0.0000 | (L)-4.50 |
LAS50 | ART Contrasts | (HD-A) - (NICU-B) | ****0.0000 | (L)-3.27 |
LAS50 | ART Contrasts | (NICU-A) - (NICU-B) | ****0.0000 | (L)1.23 |
LASmax | LME-ART-ANOVA | sub_location | ****0.0000 | (S)0.04 |
LASmax | ART Contrasts | (HD-A) - (NICU-A) | 1.0000 | (S)0.03 |
LASmax | ART Contrasts | (HD-A) - (NICU-B) | ****0.0000 | (L)0.91 |
LASmax | ART Contrasts | (NICU-A) - (NICU-B) | ****0.0000 | (L)0.88 |
LCS | LME-ART-ANOVA | sub_location | ****0.0000 | (L)0.48 |
LCS | ART Contrasts | (HD-A) - (NICU-A) | ****0.0000 | (L)-2.27 |
LCS | ART Contrasts | (HD-A) - (NICU-B) | 0.5197 | (M)0.13 |
LCS | ART Contrasts | (NICU-A) - (NICU-B) | ****0.0000 | (L)2.40 |
LCS10 | LME-ART-ANOVA | sub_location | ****0.0000 | (L)0.34 |
LCS10 | ART Contrasts | (HD-A) - (NICU-A) | ****0.0000 | (L)-1.64 |
LCS10 | ART Contrasts | (HD-A) - (NICU-B) | ****0.0000 | (L)0.56 |
LCS10 | ART Contrasts | (NICU-A) - (NICU-B) | ****0.0000 | (L)2.20 |
LCS50 | LME-ART-ANOVA | sub_location | ****0.0000 | (L)0.54 |
LCS50 | ART Contrasts | (HD-A) - (NICU-A) | ****0.0000 | (L)-2.62 |
LCS50 | ART Contrasts | (HD-A) - (NICU-B) | *0.0290 | (L)-0.23 |
LCS50 | ART Contrasts | (NICU-A) - (NICU-B) | ****0.0000 | (L)2.40 |
LCSmax | LME-ART-ANOVA | sub_location | ****0.0000 | (M)0.06 |
LCSmax | ART Contrasts | (HD-A) - (NICU-A) | ****0.0000 | (L)-0.28 |
LCSmax | ART Contrasts | (HD-A) - (NICU-B) | ****0.0000 | (L)1.05 |
LCSmax | ART Contrasts | (NICU-A) - (NICU-B) | ****0.0000 | (L)1.33 |