Respiratory Muslce Performance in Wind Instrumentalists: Outputs from a Systematic Review and Meta-Analysis

Author

Sarah Morris

Published

March 9, 2027

1 Meta-Analyses

Code

#### 1. Load packages ----
library(svglite)
library(metafor)
library(meta)
library(metaforest)
library(dplyr)
library(tidyverse)
library(readxl)
library(grid)
library(gridExtra)

#### 2. Load the data ----
data1 <- read.csv("data/WI_MA_data_30.05.25.csv")

#### 3. Process the data ----
# Calculate effect sizes using escalc() for all data
data1 <- escalc(measure="MD", m1i=mn_exp, sd1i=std_exp, n1i=n_exp, 
                m2i=mn_ctl, sd2i=std_ctl, n2i=n_ctl, data=data1, 
                slab=paste(author, year, sep=", "))

# Split data into MIP and MEP datasets
data_MIP <- data1 %>% filter(RMP == "MIP")
data_MEP <- data1 %>% filter(RMP == "MEP")

#### 4. Meta-analysis for MIP ----
# Random-effects model with Hartung-Knapp adjustment
MIP_rma <- rma(yi, vi, data = data_MIP, method = "REML", test = "knha")

# Fixed-effect model
MIP_fe <- rma(yi, vi, data = data_MIP, method = "FE")

# Sort data by effect size (MD) from largest to smallest
data_MIP <- data_MIP[order(-data_MIP$yi), ] 

# Create meta object for forest plot
MIP_meta <- metacont(
  n.e = data_MIP$n_exp, mean.e = data_MIP$mn_exp, sd.e = data_MIP$std_exp, 
  n.c = data_MIP$n_ctl, mean.c = data_MIP$mn_ctl, sd.c = data_MIP$std_ctl, 
  data = data_MIP, studlab = data_MIP$author,
  common = TRUE, random = FALSE,
  method.random.ci = FALSE,  # Hartung-Knapp adjustment
  prediction = FALSE,  # Include prediction interval
  sm = "MD",
  title = "Maximum Inspiratory Pressure (MIP) generation in wind instrumentalists vs. controls"
)

#### 5. Meta-analysis for MEP ----
# Random-effects model with Hartung-Knapp adjustment
MEP_rma <- rma(yi, vi, data = data_MEP, method = "REML", test = "knha")

# Fixed-effect model
MEP_fe <- rma(yi, vi, data = data_MEP, method = "FE")

# Sort data by effect size (MD) from largest to smallest
data_MEP <- data_MEP[order(-data_MEP$yi), ] 

# Create meta object for forest plot
MEP_meta <- metacont(
  n.e = data_MEP$n_exp, mean.e = data_MEP$mn_exp, sd.e = data_MEP$std_exp, 
  n.c = data_MEP$n_ctl, mean.c = data_MEP$mn_ctl, sd.c = data_MEP$std_ctl, 
  data = data_MEP, studlab = data_MEP$author,
  common = TRUE, random = FALSE,
  method.random.ci = FALSE,  # Hartung-Knapp adjustment
  prediction = FALSE,  # Include prediction interval
  sm = "MD",
    title = "Maximum Expiratory Pressure (MEP) generation in wind instrumentalists vs. controls"
)

1.1 MIP MA Results

Code

#### 6. Create Forest Plots ----

#### MIP Forest Plot ----
# Define output directory
output_dir <- "C:/Users/sjmor/OneDrive - The University of Sydney (Staff)/PhD Studies/Project 1. SysRev/New Review/Base/Figures"

# Create directory if it doesn't exist
if (!dir.exists(output_dir)) {
  dir.create(output_dir, recursive = TRUE)
}

# Function to create the MIP forest plot
create_mip_plot <- function() {
  par(mar = c(2, 4, 1, 2))
  
  forest(MIP_meta,
         leftcols = c("studlab", "year", "n.e", "mean.e", "sd.e", "n.c", "mean.c", "sd.c"),
         leftlabs = c("Author", "Year", "n", "Mean", "SD", "n", "Mean", "SD"),
         rightcols = c("effect", "ci"),
         rightlabs = c("MD", "95% CI"),
         digits = 1,
         digits.mean = 1,
         digits.sd = 1,
         comb.fixed = TRUE,
         comb.random = FALSE,
         prediction = FALSE,
         print.tau2 = TRUE,
         print.I2 = TRUE,
         print.H = TRUE,
         col.predict = "red",
         col.diamond = "blue",
         hetstat = TRUE,
         overall = TRUE,
         overall.hetstat = TRUE,
         test.overall.common = TRUE,
         test.overall.random = FALSE,
         main = "Maximum Inspiratory Pressure (MIP) generation in wind instrumentalists vs. controls",
         fontsize = 8,
         cex = 0.8,
         xlim = c(-30, 60),
         header.height = 0.5)
}

# Save as PNG
png(file.path(output_dir, "Figure 2. MIP forest plot.png"), width = 8, height = 6, units = "in", res = 300)
create_mip_plot()
dev.off()

png 
  2

Code

# Save as JPEG
jpeg(file.path(output_dir, "Figure 2. MIP forest plot.jpeg"), width = 8, height = 6, units = "in", res = 300)
create_mip_plot()
dev.off()

png 
  2

Code

# Save as TIFF
tiff(file.path(output_dir, "Figure 2. MIP forest plot.tiff"), width = 8, height = 6, units = "in", res = 300)
create_mip_plot()
dev.off()

png 
  2

Code

# Save as PDF
pdf(file.path(output_dir, "Figure 2. MIP forest plot.pdf"), width = 8, height = 6)
create_mip_plot()
dev.off()

png 
  2

Code

# Save as SVG
svg(file.path(output_dir, "Figure 2. MIP forest plot.svg"), width = 8, height = 6)
create_mip_plot()
dev.off()

png 
  2

Code

# Display in R Markdown
knitr::include_graphics(file.path(output_dir, "Figure 2. MIP forest plot.png"), dpi = 300)

Code

#### 7. Summary Statistics ----
# Print Analyses
MIP_meta

Review:     Maximum Inspiratory Pressure (MIP) generation in wind instrument ...

Number of studies: k = 4
Number of observations: o = 211 (o.e = 120, o.c = 91)

                         MD             95% CI    z  p-value
Common effect model 20.8060 [13.7368; 27.8753] 5.77 < 0.0001

Quantifying heterogeneity (with 95% CIs):
 tau^2 = 322.5526 [68.0641; >3225.5258]; tau = 17.9597 [8.2501; >56.7937]
 I^2 = 86.9% [68.4%; 94.6%]; H = 2.76 [1.78; 4.29]

Test of heterogeneity:
     Q d.f.  p-value
 22.87    3 < 0.0001

Details of meta-analysis methods:
- Inverse variance method
- Restricted maximum-likelihood estimator for tau^2
- Q-Profile method for confidence interval of tau^2 and tau
- Calculation of I^2 based on Q

Code

# Display summary statistics for MIP
cat("\nSummary for Maximum Inspiratory Pressure (MIP):\n")


Summary for Maximum Inspiratory Pressure (MIP):

Code

cat("Random-effects model (with Hartung-Knapp adjustment):\n")

Random-effects model (with Hartung-Knapp adjustment):

Code

print(MIP_rma)


Random-Effects Model (k = 4; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 322.5526 (SE = 309.7010)
tau (square root of estimated tau^2 value):      17.9597
I^2 (total heterogeneity / total variability):   85.74%
H^2 (total variability / sampling variability):  7.01

Test for Heterogeneity:
Q(df = 3) = 22.8665, p-val < .0001

Model Results:

estimate      se    tval  df    pval     ci.lb    ci.ub    
 17.5094  9.6336  1.8175   3  0.1667  -13.1490  48.1678    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

cat("\nFixed-effect model:\n")


Fixed-effect model:

Code

print(MIP_fe)


Fixed-Effects Model (k = 4)

I^2 (total heterogeneity / total variability):   86.88%
H^2 (total variability / sampling variability):  7.62

Test for Heterogeneity:
Q(df = 3) = 22.8665, p-val < .0001

Model Results:

estimate      se    zval    pval    ci.lb    ci.ub      
 20.8060  3.6068  5.7685  <.0001  13.7368  27.8753  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

cat("\nHeterogeneity statistics:\n")


Heterogeneity statistics:

Code

cat("I² =", formatC(MIP_rma$I2, digits=1, format="f"), "%\n")

I² = 85.7 %

Code

cat("H² =", formatC(MIP_rma$H2, digits=2, format="f"), "\n")

H² = 7.01

Code

cat("τ² =", formatC(MIP_rma$tau2, digits=4, format="f"), "\n\n")

τ² = 322.5526

1.2 MEP MA Results

Code

#### 6. Create Forest Plots ----
# MEP Forest Plot

# Function to create the MEP forest plot
create_mep_plot <- function() {
  par(mar = c(2, 4, 3, 2))
  
  forest(MEP_meta,
         leftcols = c("studlab", "year", "n.e", "mean.e", "sd.e", "n.c", "mean.c", "sd.c"),
         leftlabs = c("Author", "Year", "n", "Mean", "SD", "n", "Mean", "SD"),
         rightcols = c("effect", "ci"),
         rightlabs = c("MD", "95% CI"),
         digits = 1,
         digits.mean = 1,
         digits.sd = 1,
         comb.fixed = TRUE,
         comb.random = FALSE,
         prediction = FALSE,
         print.tau2 = TRUE,
         print.I2 = TRUE,
         print.H = TRUE,
         col.predict = "red",
         col.diamond = "blue",
         hetstat = TRUE,
         overall = TRUE,
         overall.hetstat = TRUE,
         test.overall.common = TRUE,
         test.overall.random = FALSE,
         main = "Maximum Expiratory Pressure (MEP) generation in wind instrumentalists vs. controls",
         fontsize = 7,
         cex = 0.7,
         xlim = c(-20, 90),
         header.height = 0.5)
}

# Save as PNG
png(file.path(output_dir, "Figure 3. MEP forest plot.png"), width = 8, height = 6, units = "in", res = 300)
create_mep_plot()
dev.off()

png 
  2

Code

# Save as JPEG
jpeg(file.path(output_dir, "Figure 3. MEP forest plot.jpeg"), width = 8, height = 6, units = "in", res = 300)
create_mep_plot()
dev.off()

png 
  2

Code

# Save as TIFF
tiff(file.path(output_dir, "Figure 3. MEP forest plot.tiff"), width = 8, height = 6, units = "in", res = 300)
create_mep_plot()
dev.off()

png 
  2

Code

# Save as PDF
pdf(file.path(output_dir, "Figure 3. MEP forest plot.pdf"), width = 8, height = 6)
create_mep_plot()
dev.off()

png 
  2

Code

# Save as SVG
svg(file.path(output_dir, "Figure 3. MEP forest plot.svg"), width = 8, height = 6)
create_mep_plot()
dev.off()

png 
  2

Code

# Display in R Markdown
knitr::include_graphics(file.path(output_dir, "Figure 3. MEP forest plot.png"), dpi = 300)

Code

#### 7. Summary Statistics ----
# Print Analyses
MEP_meta

Review:     Maximum Expiratory Pressure (MEP) generation in wind instrumenta ...

Number of studies: k = 3
Number of observations: o = 185 (o.e = 107, o.c = 78)

                         MD            95% CI    z p-value
Common effect model 15.3436 [7.0355; 23.6518] 3.62  0.0003

Quantifying heterogeneity (with 95% CIs):
 tau^2 = 122.8165 [0.0000; >1228.1648]; tau = 11.0823 [0.0000; >35.0452]
 I^2 = 50.2% [0.0%; 85.6%]; H = 1.42 [1.00; 2.64]

Test of heterogeneity:
    Q d.f. p-value
 4.02    2  0.1341

Details of meta-analysis methods:
- Inverse variance method
- Restricted maximum-likelihood estimator for tau^2
- Q-Profile method for confidence interval of tau^2 and tau
- Calculation of I^2 based on Q

Code

# Display summary statistics for MEP
cat("\nSummary for Maximum Expiratory Pressure (MEP):\n")


Summary for Maximum Expiratory Pressure (MEP):

Code

cat("Random-effects model (with Hartung-Knapp adjustment):\n")

Random-effects model (with Hartung-Knapp adjustment):

Code

print(MEP_rma)


Random-Effects Model (k = 3; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 122.8165 (SE = 251.1952)
tau (square root of estimated tau^2 value):      11.0823
I^2 (total heterogeneity / total variability):   50.13%
H^2 (total variability / sampling variability):  2.01

Test for Heterogeneity:
Q(df = 2) = 4.0189, p-val = 0.1341

Model Results:

estimate      se    tval  df    pval     ci.lb    ci.ub    
 21.1392  8.9177  2.3705   2  0.1412  -17.2307  59.5091    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

cat("\nFixed-effect model:\n")


Fixed-effect model:

Code

print(MEP_meta)

Review:     Maximum Expiratory Pressure (MEP) generation in wind instrumenta ...

Number of studies: k = 3
Number of observations: o = 185 (o.e = 107, o.c = 78)

                         MD            95% CI    z p-value
Common effect model 15.3436 [7.0355; 23.6518] 3.62  0.0003

Quantifying heterogeneity (with 95% CIs):
 tau^2 = 122.8165 [0.0000; >1228.1648]; tau = 11.0823 [0.0000; >35.0452]
 I^2 = 50.2% [0.0%; 85.6%]; H = 1.42 [1.00; 2.64]

Test of heterogeneity:
    Q d.f. p-value
 4.02    2  0.1341

Details of meta-analysis methods:
- Inverse variance method
- Restricted maximum-likelihood estimator for tau^2
- Q-Profile method for confidence interval of tau^2 and tau
- Calculation of I^2 based on Q

Code

cat("\nHeterogeneity statistics:\n")


Heterogeneity statistics:

Code

cat("I² =", formatC(MEP_rma$I2, digits=1, format="f"), "%\n")

I² = 50.1 %

Code

cat("H² =", formatC(MEP_rma$H2, digits=2, format="f"), "\n")

H² = 2.01

Code

cat("τ² =", formatC(MEP_rma$tau2, digits=4, format="f"), "\n")

τ² = 122.8165

2 Risk of Bias

Code

# Risk of Bias Visualization with Symbols ----
library(ggplot2)
library(tidyr)
library(dplyr)

# Read the CSV file
data1 <- read.csv("C:/Users/sjmor/OneDrive - The University of Sydney (Staff)/PhD Studies/Project 1. SysRev/New Review/Base/Figures/ROB_data.csv", header = TRUE, check.names = FALSE)

# Define output directory
output_dir <- "C:/Users/sjmor/OneDrive - The University of Sydney (Staff)/PhD Studies/Project 1. SysRev/New Review/Base/Figures"

# Get all RoB columns (columns 2 through 19, excluding "Overall")
rob_columns <- names(data1)[2:18]  # All columns except Study and Overall

# Create a subset with just the RoB data
rob_subset <- data1[1:4, c("Study", rob_columns)]

# Prepare data for plotting
rob_long <- rob_subset %>%
  pivot_longer(cols = -Study, 
               names_to = "Domain", 
               values_to = "Risk") %>%
  mutate(
    Risk_category = factor(Risk, levels = c("Low", "Unclear", "High")),
    Symbol = case_when(
      Risk == "Low" ~ "+",      # Green plus
      Risk == "Unclear" ~ "—",  # Orange minus
      Risk == "High" ~ "X",     # Red X
      TRUE ~ "?"                # Question mark for missing
    ),
    # Preserve original column order
    Domain = factor(Domain, levels = rob_columns)
  )

# Create color mapping for risk levels
risk_colors <- c("Low" = "#02C100",      # Green
                 "Unclear" = "#FFA500",  # Orange
                 "High" = "#BF0000")     # Red

# Create the risk of bias plot with symbols
rob_plot <- ggplot(rob_long, aes(x = Domain, y = Study)) +
  geom_tile(aes(fill = Risk_category), color = "white", size = 1) +
  geom_text(aes(label = Symbol), size = 5, fontface = "bold") +
  scale_fill_manual(values = risk_colors, 
                    na.value = "gray90",
                    name = "Risk of Bias",
                    labels = c("Low (+)", "Unclear (-)", "High (X)")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 14),
        axis.text.y = element_text(size = 14),
        axis.title = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        panel.grid = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 18, face = "bold")) +
  labs(title = "Risk of Bias Summary")

# Display the plot
print(rob_plot)

Code

# Save the plot in multiple formats
ggsave(file.path(output_dir, "Figure 4. ROB Plot.svg"),
       plot = rob_plot,
       width = 14, height = 5, units = "in", dpi = 300)

ggsave(file.path(output_dir, "Figure 4. ROB Plot.tiff"),
       plot = rob_plot,
       width = 14, height = 5, units = "in", dpi = 300)

ggsave(file.path(output_dir, "Figure 4. ROB Plot.png"),
       plot = rob_plot,
       width = 14, height = 5, units = "in", dpi = 300)

ggsave(file.path(output_dir, "Figure 4. ROB Plot.jpeg"),
       plot = rob_plot,
       width = 14, height = 5, units = "in", dpi = 300, quality = 100)

cat("Risk of bias plots saved to:", output_dir, "\n")

Risk of bias plots saved to: C:/Users/sjmor/OneDrive - The University of Sydney (Staff)/PhD Studies/Project 1. SysRev/New Review/Base/Figures

Code

# Optional: Create a summary bar chart showing proportion of each risk level
rob_summary <- rob_long %>%
  filter(!is.na(Risk)) %>%
  group_by(Domain, Risk_category) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(Domain) %>%
  mutate(proportion = count / sum(count) * 100)

rob_summary_plot <- ggplot(rob_summary, aes(x = Domain, y = proportion, fill = Risk_category)) +
  geom_bar(stat = "identity", color = "white") +
  scale_fill_manual(values = risk_colors,
                    name = "Risk of Bias",
                    labels = c("Low (+)", "Unclear (-)", "High (X)")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 14),
        axis.text.y = element_text(size = 14),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 14),
        legend.position = "bottom",
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        plot.title = element_text(hjust = 0.5, size = 18, face = "bold")) +
  labs(title = "Risk of Bias Summary - Proportion by Domain",
       y = "Percentage (%)") +
  coord_flip()

# Display summary plot
print(rob_summary_plot)

Code

# Save summary plot in multiple formats
ggsave(file.path(output_dir, "Figure 5. ROB Summary.svg"),
       plot = rob_summary_plot,
       width = 10, height = 8, units = "in", dpi = 300)

ggsave(file.path(output_dir, "Figure 5. ROB Summary.tiff"),
       plot = rob_summary_plot,
       width = 10, height = 8, units = "in", dpi = 300)

ggsave(file.path(output_dir, "Figure 5. ROB Summary.png"),
       plot = rob_summary_plot,
       width = 10, height = 8, units = "in", dpi = 300)

ggsave(file.path(output_dir, "Figure 5. ROB Summary.jpeg"),
       plot = rob_summary_plot,
       width = 10, height = 8, units = "in", dpi = 300, quality = 100)

cat("Risk of bias summary plots saved to:", output_dir, "\n")

Risk of bias summary plots saved to: C:/Users/sjmor/OneDrive - The University of Sydney (Staff)/PhD Studies/Project 1. SysRev/New Review/Base/Figures