Combined Meta Analysis

1 Meta Analyses

1.1 Data Prepping

# Load packages ----
library(metafor)
library(meta)
library(metaforest)
library(dplyr)
library(tidyverse)
library(readxl)
library(grid)
library(gridExtra)
library(ggplot2)

# PFT Prepping -----------
## Load the PFT data
spiro_data <- read.csv("../data/spiro_data_01.04.25.csv")

# Clean data: select required columns and remove missing values
spiro_data_clean <- spiro_data %>% 
  select(studyID, ESID, author, year, PFT, n_exp, mn_exp, std_exp, n_ctl, mn_ctl, std_ctl) %>% 
  filter(!is.na(std_exp) & !is.na(std_ctl) &
         !is.na(n_exp) & !is.na(mn_exp) &
         !is.na(n_ctl) & !is.na(mn_ctl))

# Convert ESID to factor and sort by ESID
spiro_data_clean$ESID <- as.factor(spiro_data_clean$ESID)
spiro_data_clean <- spiro_data_clean %>% arrange(ESID)

# MIP/MEP Prepping --------------
# Load the MIP/MEP data
data1 <- read.csv("../data/WI_MA_data_25.02.25.csv")

## 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")

1.2 MIP MA

# Calculations ----
## 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")

## 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 = TRUE,
  method.random.ci = TRUE,  # Hartung-Knapp adjustment
  prediction = TRUE,  # Include prediction interval
  sm = "MD"
)

# Forest Plots ----
## Set minimal margins (bottom, left, top, right)
png("MIP_forest_plot.png", width = 8, height = 6, units = "in", res = 300)

par(mar = c(2, 4, 1, 2))  

forest(MIP_meta,
       leftcols = c("studlab", "n.e", "mean.e", "sd.e", "n.c", "mean.c", "sd.c"),
       leftlabs = c("Author", "n", "Mean", "SD", "n", "Mean", "SD"),
       rightcols = c("effect", "ci"),
       rightlabs = c("MD", "95% CI"),
       comb.fixed = TRUE,
       comb.random = TRUE,
       prediction = TRUE,
       print.tau2 = TRUE,
       print.I2 = TRUE,
       print.H = TRUE,  # Add H^2 statistic
       col.predict = "red",  # Prediction interval in red
       col.diamond = "blue",  # Confidence interval in blue
       hetstat = TRUE,
       overall = TRUE,
       overall.hetstat = TRUE,
       test.overall.common = TRUE,
       test.overall.random = TRUE,
       main = "Maximum Inspiratory Pressure (MIP) generation in wind instrumentalists vs. controls",
       fontsize = 8,     # Reduced from 10 to 8
       cex = 0.8,        # Added to control element size
       xlim = c(-50, 50), # Adjust as needed for your data range
       header.height = 0.5) # Reduce header height

dev.off()

png 
  2 

## Save image
knitr::include_graphics("MIP_forest_plot.png", dpi = 300)

# Summary Statistics ----
## Print Analyses
MIP_meta

Number of studies: k = 8
Number of observations: o = 358 (o.e = 201, o.c = 157)

                          MD              95%-CI  z|t  p-value
Common effect model  14.7971 [  8.9443; 20.6500] 4.96 < 0.0001
Random effects model 11.3964 [ -2.2482; 25.0409] 1.98   0.0888
Prediction interval          [-26.3817; 49.1744]              

Quantifying heterogeneity (with 95%-CIs):
 tau^2 = 217.1740 [52.8257; 920.6605]; tau = 14.7368 [7.2681; 30.3424]
 I^2 = 79.4% [59.8%; 89.4%]; H = 2.20 [1.58; 3.07]

Test of heterogeneity:
     Q d.f.  p-value
 33.94    7 < 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
- Hartung-Knapp adjustment for random effects model (df = 7)
- Prediction interval based on t-distribution (df = 7)

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

Summary for Maximum Inspiratory Pressure (MIP):

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

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

print(MIP_rma)

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

tau^2 (estimated amount of total heterogeneity): 217.1740 (SE = 161.6789)
tau (square root of estimated tau^2 value):      14.7368
I^2 (total heterogeneity / total variability):   74.56%
H^2 (total variability / sampling variability):  3.93

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

Model Results:

estimate      se    tval  df    pval    ci.lb    ci.ub    
 11.3964  5.7703  1.9750   7  0.0888  -2.2482  25.0409  . 

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

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

Fixed-effect model:

print(MIP_fe)

Fixed-Effects Model (k = 8)

I^2 (total heterogeneity / total variability):   79.38%
H^2 (total variability / sampling variability):  4.85

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

Model Results:

estimate      se    zval    pval   ci.lb    ci.ub      
 14.7971  2.9862  4.9552  <.0001  8.9443  20.6500  *** 

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

cat("\nHeterogeneity statistics:\n")

Heterogeneity statistics:

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

I² = 74.6 %

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

H² = 3.93 

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

τ² = 217.1740 

1.3 MEP MA

#Calculations -------------

## 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")

## 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 = TRUE,
  method.random.ci = TRUE,  # Hartung-Knapp adjustment
  prediction = TRUE,  # Include prediction interval
  sm = "MD"
)

# Forest Plots ----
## Set smaller margins (bottom, left, top, right)
png("MEP_forest_plot.png", width = 8, height = 6, units = "in", res = 300)

par(mar = c(2, 4, 1, 2))  

forest(MEP_meta,
       leftcols = c("studlab", "n.e", "mean.e", "sd.e", "n.c", "mean.c", "sd.c"),
       leftlabs = c("Author", "n", "Mean", "SD", "n", "Mean", "SD"),
       rightcols = c("effect", "ci"),
       rightlabs = c("MD", "95% CI"),
       comb.fixed = TRUE,
       comb.random = TRUE,
       prediction = TRUE,
       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 = TRUE,
       main = "Maximum Expiratory Pressure (MEP) generation in wind instrumentalists vs. controls",
       fontsize = 8,
       cex = 0.8,
       xlim = c(-50, 50),
       header.height = 0.5)  # Reduce header height

dev.off()

png 
  2 

## Save image
knitr::include_graphics("MIP_forest_plot.png", dpi = 300)

# Summary Statistics ----
## Print Analyses
MEP_meta

Number of studies: k = 8
Number of observations: o = 392 (o.e = 218, o.c = 174)

                          MD              95%-CI  z|t  p-value
Common effect model  12.2426 [  6.2282; 18.2569] 3.99 < 0.0001
Random effects model 12.2735 [ -0.9235; 25.4706] 2.20   0.0638
Prediction interval          [-14.6029; 39.1499]              

Quantifying heterogeneity (with 95%-CIs):
 tau^2 = 102.4046 [0.0000; >1024.0461]; tau = 10.1195 [0.0000; >32.0007]
 I^2 = 50.2% [0.0%; 77.8%]; H = 1.42 [1.00; 2.12]

Test of heterogeneity:
     Q d.f. p-value
 14.07    7  0.0500

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
- Hartung-Knapp adjustment for random effects model (df = 7)
- Prediction interval based on t-distribution (df = 7)

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

Summary for Maximum Expiratory Pressure (MEP):

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

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

print(MEP_rma)

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

tau^2 (estimated amount of total heterogeneity): 102.4046 (SE = 110.7729)
tau (square root of estimated tau^2 value):      10.1195
I^2 (total heterogeneity / total variability):   53.62%
H^2 (total variability / sampling variability):  2.16

Test for Heterogeneity:
Q(df = 7) = 14.0676, p-val = 0.0500

Model Results:

estimate      se    tval  df    pval    ci.lb    ci.ub    
 12.2735  5.5810  2.1991   7  0.0638  -0.9235  25.4706  . 

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

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

Fixed-effect model:

print(MEP_fe)

Fixed-Effects Model (k = 8)

I^2 (total heterogeneity / total variability):   50.24%
H^2 (total variability / sampling variability):  2.01

Test for Heterogeneity:
Q(df = 7) = 14.0676, p-val = 0.0500

Model Results:

estimate      se    zval    pval   ci.lb    ci.ub      
 12.2426  3.0686  3.9896  <.0001  6.2282  18.2569  *** 

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

cat("\nHeterogeneity statistics:\n")

Heterogeneity statistics:

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

I² = 53.6 %

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

H² = 2.16 

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

τ² = 102.4046 

1.4 PFT MA

# Calculations -------------
## Calculate effect sizes (Hedges' g)
es <- escalc(measure = "SMD", 
             n1i = spiro_data_clean$n_exp, 
             m1i = spiro_data_clean$mn_exp, 
             sd1i = spiro_data_clean$std_exp,
             n2i = spiro_data_clean$n_ctl, 
             m2i = spiro_data_clean$mn_ctl, 
             sd2i = spiro_data_clean$std_ctl,
             data = spiro_data_clean)

## Create custom labels (adding outcome from PFT)
es$slab <- paste(spiro_data_clean$author, spiro_data_clean$year, "\
Outcome:", spiro_data_clean$PFT)

## Run the multi-level meta-analysis with ESID nested within studyID 
res_mv <- rma.mv(yi, vi, random = ~ 1 | studyID/ESID, data = es, method = "REML")

## Print summary of meta-analysis
summary_result <- summary(res_mv)
print(summary_result)

Multivariate Meta-Analysis Model (k = 111; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-389.3753   778.7505   784.7505   792.8520   784.9769   

Variance Components:

            estim    sqrt  nlvls  fixed        factor 
sigma^2.1  0.2744  0.5238     20     no       studyID 
sigma^2.2  0.4620  0.6797     35     no  studyID/ESID 

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

Model Results:

estimate      se    zval    pval    ci.lb   ci.ub    
  0.2116  0.1715  1.2340  0.2172  -0.1245  0.5477    

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

## Determine custom row positions with gaps between ESID groups
## And add space for subtitles and summary statistics
grp_counts <- table(spiro_data_clean$ESID)
grp_names <- names(grp_counts)
n_groups <- length(grp_counts)

### - 2 rows at the top for padding (extra space under column labels)
### - 2 rows for each group (title + summary)
### - 1 gap row between each group
extra_top_padding <- 2  # Increased from 1 to 2
extra_rows_per_group <- 2  # 1 for title + 1 for summary
total_extra_rows <- extra_top_padding + (extra_rows_per_group * n_groups) + (n_groups - 1)
total_rows <- sum(grp_counts) + total_extra_rows

## Define ESID group titles
esid_titles <- c(
  "1" = "Expiratory Flow Outcomes",
  "2" = "Lung Capacity Outcomes",
  "3" = "Cardiopulmonary Outcomes",
  "4" = "Ventilatory Outcomes"
)

## Calculate row positions for each study, title, and summary
row_positions <- numeric(nrow(es))
title_positions <- numeric(n_groups)
summary_positions <- numeric(n_groups)
gap_positions <- numeric(n_groups - 1)

current_row <- extra_top_padding + 1  # Start after top padding

for (i in 1:n_groups) {
  # Position for group title
  title_positions[i] <- current_row
  current_row <- current_row + 1
  
  # Position for group summary (right after the title)
  summary_positions[i] <- current_row
  current_row <- current_row + 1
  
  # Positions for studies in this group
  group_indices <- which(spiro_data_clean$ESID == grp_names[i])
  group_size <- length(group_indices)
  row_positions[group_indices] <- current_row:(current_row + group_size - 1)
  current_row <- current_row + group_size
  
  # Add gap after each group except the last
  if (i < n_groups) {
    gap_positions[i] <- current_row - 0.5
    current_row <- current_row + 1
  }
}

## Calculate dynamic height based on number of rows
row_height_inches <- 0.4
dynamic_height <- total_rows * row_height_inches + 3  # Add 3 inches for margins and title

## Calculate group-specific meta-analysis results
group_results <- list()
for (i in 1:n_groups) {
  group_indices <- which(spiro_data_clean$ESID == grp_names[i])
  if (length(group_indices) > 1) {
    # Only run meta-analysis if there's more than one study in the group
    group_data <- es[group_indices, ]
    group_res <- try(rma.mv(yi, vi, random = ~ 1 | studyID, data = group_data, method = "REML"), silent = TRUE)
    
    if (!inherits(group_res, "try-error")) {
      group_results[[i]] <- group_res
    } else {
      # Fallback to simpler model if the multilevel model fails
      group_res <- try(rma(yi, vi, data = group_data, method = "REML"), silent = TRUE)
      if (!inherits(group_res, "try-error")) {
        group_results[[i]] <- group_res
      } else {
        group_results[[i]] <- NULL
      }
    }
  } else if (length(group_indices) == 1) {
    # For single study, just store the effect size and variance
    group_results[[i]] <- list(
      b = es$yi[group_indices],
      se = sqrt(es$vi[group_indices]),
      ci.lb = es$yi[group_indices] - 1.96 * sqrt(es$vi[group_indices]),
      ci.ub = es$yi[group_indices] + 1.96 * sqrt(es$vi[group_indices])
    )
  } else {
    group_results[[i]] <- NULL
  }
}

## Print group-specific results
for (i in 1:n_groups) {
  if (!is.null(group_results[[i]])) {
    if (is.list(group_results[[i]]) && !is.null(group_results[[i]]$b)) {
      cat("ESID", grp_names[i], ":", esid_titles[grp_names[i]], "\
")
      est <- sprintf("%.3f", group_results[[i]]$b)
      ci_l <- sprintf("%.3f", group_results[[i]]$ci.lb)
      ci_u <- sprintf("%.3f", group_results[[i]]$ci.ub)
      cat("Group Summary: ", est, " [", ci_l, ", ", ci_u, "]\
\
")
    }
  }
}

ESID 1 : Expiratory Flow Outcomes 
Group Summary:  0.059  [ -0.198 ,  0.316 ]

ESID 2 : Lung Capacity Outcomes 
Group Summary:  0.017  [ -0.658 ,  0.691 ]

ESID 3 : Cardiopulmonary Outcomes 
Group Summary:  0.251  [ -0.370 ,  0.871 ]

ESID 4 : Ventilatory Outcomes 
Group Summary:  0.769  [ -0.289 ,  1.827 ]

## Extract meta-analysis statistics for annotation
mv_sum <- summary(res_mv)
overall_est <- sprintf("%.3f", mv_sum$b)
overall_se <- sprintf("%.3f", mv_sum$se)
ci_lb <- sprintf("%.3f", mv_sum$ci.lb)
ci_ub <- sprintf("%.3f", mv_sum$ci.ub)
tau2 <- sprintf("%.3f", sum(mv_sum$sigma2))
## I² calculation - this is an approximation for multilevel models
I2 <- round(100 * sum(mv_sum$sigma2) / (sum(mv_sum$sigma2) + mean(mv_sum$vi)), 1)

## Compose the meta-analysis text annotation (centered)
meta_text <- paste0("Overall Estimate: ", overall_est, " (SE=", overall_se, 
                    ", 95% CI: [", ci_lb, ", ", ci_ub, "])     ",
                    "Tau² = ", tau2, ", I² = ", I2, "%")

## Print overall meta-analysis statistics
cat("Overall Meta-analysis Statistics:\
")

Overall Meta-analysis Statistics:

cat(meta_text, "\
\
")

Overall Estimate: 0.212 (SE=0.171, 95% CI: [-0.124, 0.548])     Tau² = 0.736, I² = 91.2% 

# Forest and Funnel Plots -------------------
## Set up the plotting area with adjusted margins
oldpar <- par(mar = c(6, 4, 2, 1))  # Increase bottom margin for meta text

## Create the forest plot with dynamic sizing
## First save as PDF
pdf("forest_plot_final.pdf", width = 10, height = dynamic_height)

## Create the basic forest plot
forest(res_mv,
       slab = es$slab,
       xlab = "Effect Size (Hedges' g)",
       main = "Forest Plot for Multi-level Meta-analysis\
(with ESID grouping and PFT outcome)",
       cex = 0.8,           # Text size for study labels
       cex.lab = 1.0,       # Label size
       cex.axis = 0.9,      # Axis text size
       rows = row_positions,
       ylim = c(0, total_rows + 1)
)

## Add horizontal lines to separate ESID groups
if(length(gap_positions) > 0) {
  abline(h = gap_positions, lty = 2, col = "gray")
}

## Add ESID group titles and summaries
for (i in 1:n_groups) {
  # Add group title
  text(-6, title_positions[i], paste("ESID", grp_names[i], ":", esid_titles[grp_names[i]]), 
       cex = 1.1, font = 2, adj = 0)
  
  # Add group summary statistics in red under the title
  if (!is.null(group_results[[i]])) {
    if (is.list(group_results[[i]]) && !is.null(group_results[[i]]$b)) {
      est <- sprintf("%.3f", group_results[[i]]$b)
      ci_l <- sprintf("%.3f", group_results[[i]]$ci.lb)
      ci_u <- sprintf("%.3f", group_results[[i]]$ci.ub)
      summary_text <- paste0("Group Summary: ", est, " [", ci_l, ", ", ci_u, "]")
      text(-6, summary_positions[i], summary_text, col = "red", cex = 0.9, adj = 0)
    }
  }
}

## Add meta-analysis statistics text at the bottom (centered)
mtext(meta_text, side = 1, line = 4, cex = 0.8, adj = 0.5)

dev.off()

png 
  2 

## Also save as PNG for easier viewing
png("forest_plot_final2.png", width = 10*150, height = dynamic_height*150, res = 150)

## Create the basic forest plot
forest(res_mv,
       slab = es$slab,
       xlab = "Effect Size (Hedges' g)",
       main = "Forest Plot for Multi-level Meta-analysis\
(with ESID grouping and PFT outcome)",
       cex = 0.8,           # Text size for study labels
       cex.lab = 1.0,       # Label size
       cex.axis = 0.9,      # Axis text size
       rows = row_positions,
       ylim = c(0, total_rows + 1)
)

## Add horizontal lines to separate ESID groups
if(length(gap_positions) > 0) {
  abline(h = gap_positions, lty = 2, col = "gray")
}

## Add ESID group titles and summaries
for (i in 1:n_groups) {
  # Add group title
  text(-6, title_positions[i], paste("ESID", grp_names[i], ":", esid_titles[grp_names[i]]), 
       cex = 1.1, font = 2, adj = 0)
  
  # Add group summary statistics in red under the title
  if (!is.null(group_results[[i]])) {
    if (is.list(group_results[[i]]) && !is.null(group_results[[i]]$b)) {
      est <- sprintf("%.3f", group_results[[i]]$b)
      ci_l <- sprintf("%.3f", group_results[[i]]$ci.lb)
      ci_u <- sprintf("%.3f", group_results[[i]]$ci.ub)
      summary_text <- paste0("Group Summary: ", est, " [", ci_l, ", ", ci_u, "]")
      text(-6, summary_positions[i], summary_text, col = "red", cex = 0.9, adj = 0)
    }
  }
}

## Add meta-analysis statistics text at the bottom (centered)
mtext(meta_text, side = 1, line = 4, cex = 0.8, adj = 0.5)

dev.off()

png 
  2 

## Create the funnel plot
pdf("funnel_plot_final.pdf", width = 8, height = 8)
funnel(res_mv, 
       main = "Funnel Plot for Publication Bias Assessment", 
       cex = 1.2, 
       cex.axis = 1.2)
dev.off()

png 
  2 

png("funnel_plot_final.png", width = 8*150, height = 8*150, res = 150)
funnel(res_mv, 
       main = "Funnel Plot for Publication Bias Assessment", 
       cex = 1.2, 
       cex.axis = 1.2)
dev.off()

png 
  2 

## Display the forest plot in the R graphics device
## Create the basic forest plot
forest(res_mv,
       slab = es$slab,
       xlab = "Effect Size (Hedges' g)",
       main = "Forest Plot for Multi-level Meta-analysis\
(with ESID grouping and PFT outcome)",
       cex = 0.8,           # Text size for study labels
       cex.lab = 1.0,       # Label size
       cex.axis = 0.9,      # Axis text size
       rows = row_positions,
       ylim = c(0, total_rows + 1)
)

## Add horizontal lines to separate ESID groups
if(length(gap_positions) > 0) {
  abline(h = gap_positions, lty = 2, col = "gray")
}

## Add ESID group titles and summaries
for (i in 1:n_groups) {
  # Add group title
  text(-6, title_positions[i], paste("ESID", grp_names[i], ":", esid_titles[grp_names[i]]), 
       cex = 1.1, font = 2, adj = 0)
  
  # Add group summary statistics in red under the title
  if (!is.null(group_results[[i]])) {
    if (is.list(group_results[[i]]) && !is.null(group_results[[i]]$b)) {
      est <- sprintf("%.3f", group_results[[i]]$b)
      ci_l <- sprintf("%.3f", group_results[[i]]$ci.lb)
      ci_u <- sprintf("%.3f", group_results[[i]]$ci.ub)
      summary_text <- paste0("Group Summary: ", est, " [", ci_l, ", ", ci_u, "]")
      text(-6, summary_positions[i], summary_text, col = "red", cex = 0.9, adj = 0)
    }
  }
}

## Add meta-analysis statistics text at the bottom (centered)
mtext(meta_text, side = 1, line = 4, cex = 0.8, adj = 0.5)

## Display the funnel plot in the R graphics device
funnel(res_mv, 
       main = "Funnel Plot for Publication Bias Assessment", 
       cex = 1.2, 
       cex.axis = 1.2)

## Restore original graphical parameters
par(oldpar)

## Print confirmation of saved files
print("Final files saved:")

[1] "Final files saved:"

print("1. forest_plot_final.pdf")

[1] "1. forest_plot_final.pdf"

print("2. forest_plot_final.png")

[1] "2. forest_plot_final.png"

print("3. funnel_plot_final.pdf")

[1] "3. funnel_plot_final.pdf"

print("4. funnel_plot_final.png")

[1] "4. funnel_plot_final.png"

2 Methods

This report synthesizes data from multiple meta-analyses examining respiratory function in wind instrumentalists compared to non-instrumentalist controls. The analytical methods employed included:

1. Random-effects and fixed-effects models for Maximum Inspiratory Pressure (MIP) and Maximum Expiratory Pressure (MEP) meta-analyses, using the inverse variance method.
2. Restricted maximum-likelihood estimator (REML) for tau² calculation, which estimates the amount of total heterogeneity between studies.
3. Hartung-Knapp adjustment applied to random effects models to provide more conservative estimates for small sample sizes.
4. Q-Profile method for confidence interval calculations of tau² and tau.
5. Multivariate meta-analysis model for Pulmonary Function Tests (PFT), facilitating analysis of multiple outcomes simultaneously.
6. Forest plots and funnel plots for visualization of effect sizes and assessment of publication bias.
7. Heterogeneity assessment using I² and H² statistics, which quantify the proportion of observed variance reflecting real differences in effect size rather than sampling error.

3 Results

Maximum Inspiratory Pressure (MIP) - Random-effects model estimate: 11.40 cmH₂O (95% CI: -2.25 to 25.04, p = 0.0888) - **Fixed-effects model estimate: 14.80 cmH₂O (95% CI: 8.94 to 20.65, p < 0.0001) - Heterogeneity statistics*: - I² = 74.6% (indicating substantial heterogeneity) - H² = 3.93 - τ² = 217.17 - Q(df = 7) = 33.94, p < 0.0001

Maximum Expiratory Pressure (MEP) - Random-effects model estimate: 12.27 cmH₂O (95% CI: -0.92 to 25.47, p = 0.0638) - Fixed-effects model estimate: 12.24 cmH₂O (95% CI: 6.23 to 18.26, p < 0.0001) - Heterogeneity statistics: - I² = 53.6% (indicating moderate heterogeneity) - H² = 2.16 - τ² = 102.40 - Q(df = 7) = 14.07, p = 0.0500

Pulmonary Function Tests (PFT) - Overall effect size estimate: 0.212 (95% CI: -0.124 to 0.548, p = 0.2172) - Heterogeneity statistics: - τ² = 0.736 - I² = 91.2% (indicating very high heterogeneity) - Q(df = 110) = 1916.97, p < 0.0001

PFT Subgroup Results (Standardized Mean Differences) - Expiratory Flow Outcomes: 0.059 (95% CI: -0.198 to 0.316) - Lung Capacity Outcomes: 0.017 (95% CI: -0.658 to 0.691) - Cardiopulmonary Outcomes: 0.251 (95% CI: -0.370 to 0.871) - Ventilatory Outcomes: 0.769 (95% CI: -0.289 to 1.827)

4 Discussion

The meta-analyses reveal a pattern wherein respiratory muscle strength (MIP and MEP) shows greater differences between wind instrumentalists and controls compared to standard pulmonary function measures.

Respiratory Muscle Strength (MIP and MEP)

The fixed-effects models for both MIP and MEP show statistically significant differences favoring wind instrumentalists, with mean differences of 14.80 cmH₂O and 12.24 cmH₂O, respectively. However, the random-effects models yielded marginally non-significant results, though with appreciable effect sizes in the same direction. This pattern aligns with previous research by Jiang et al. (2016), who reported that wind instrumentalists develop specific adaptations in respiratory muscles due to the controlled breathing patterns required during performance.

These findings are physiologically plausible given that wind instrumentalists regularly engage in activities requiring sustained respiratory pressure control, effectively serving as a form of respiratory muscle training (RMT). As noted by Sapienza et al. (2011), extended periods of controlled airflow against resistance can lead to adaptations similar to those seen in targeted RMT programs, including increased muscle fiber cross-sectional area and enhanced neuromuscular efficiency.

Pulmonary Function Tests

The smaller and non-significant differences in standard PFT measures (overall SMD = 0.212, p = 0.2172) suggest that wind instrument playing may not substantially alter overall lung capacity or airflow characteristics. This finding is consistent with Bouhuys (1964), who proposed that while respiratory muscles may adapt to playing demands, the structural properties of the lungs themselves may be less responsive to this specific form of training.

The trend toward larger effects in ventilatory outcomes (SMD = 0.769) compared to lung capacity outcomes (SMD = 0.017) further supports the hypothesis that wind instrument playing primarily affects the control of breathing rather than pulmonary structure. As Deniz et al. (2006) noted, the technical requirements of wind instrument performance focus more on precise respiratory control than on maximizing lung volumes.

Shifting the focus onto the respiratory muscles

These findings provide several compelling reasons to focus research and clinical attention on respiratory muscle function rather than standard pulmonary function measures when examining the physiological impacts of wind instrument playing:

1. Magnitude of Effects: The mean differences in MIP and MEP (11.40-14.80 cmH₂O) represent clinically meaningful differences in respiratory muscle strength, whereas the smaller effects on PFT measures suggest minimal structural pulmonary adaptations.

2. Training Specificity: Wind instrument playing involves controlled breathing against resistance, directly training the respiratory muscles rather than challenging maximal ventilatory capacity. This aligns with training specificity principles outlined by McConnell (2013), suggesting that adaptations would be expected primarily in the muscles being trained rather than in lung structure.

3. Mechanistic Plausibility: The respiratory demands of wind instrument playing—maintaining specific pressures, controlling airflow rate, and sustaining exhalation—involve precise coordination of inspiratory and expiratory muscles, similar to targeted RMT exercises (Illi et al., 2012).

4. Clinical Applications: The observed differences in respiratory muscle strength have potential therapeutic implications. Bausewein et al. (2008) noted that targeted respiratory muscle training might benefit patients with respiratory weakness, suggesting that modified wind instrument playing could potentially serve as an engaging form of RMT.

5. Consistency with Broader Literature: The results align with research by Sakakibara and Onodera (2014), who found that respiratory muscle performance improvements can occur independently of changes in standard pulmonary function measures in other forms of respiratory training.

5 Limitations

Several important limitations must be considered when interpreting these results:

1. Heterogeneity: The substantial heterogeneity observed (I² = 74.6% for MIP, 53.6% for MEP, and 91.2% for PFT) indicates considerable variability among studies, potentially due to differences in study design, measurement techniques, participant characteristics, or types of wind instruments studied.

2. Confidence Intervals for Random Effects Models: The 95% confidence intervals for the random effects models of both MIP and MEP include zero, suggesting uncertainty about the true effect direction despite the point estimates favoring wind instrumentalists.

3. Small Number of Studies: The MIP and MEP meta-analyses each included only 8 studies, which limits the precision of the pooled estimates and the ability to explore potential moderating factors.

4. Lack of Longitudinal Data: Most included studies were likely cross-sectional, making it difficult to distinguish between training-induced adaptations and self-selection effects (people with naturally stronger respiratory muscles may be more likely to continue wind instrument playing).

5. Potential Publication Bias: Though funnel plots were generated, specific results regarding publication bias were not provided in the data. The possibility remains that unpublished studies with null findings may exist.

6. Variation in Wind Instrument Types: Different wind instruments place varying demands on respiratory muscles, but subgroup analyses by instrument type were not included in these analyses. This was primarily due to a large amount of the included studies not specifying the individual instruments played, i.e., general ‘wind instruments’.

6 Conclusions

This meta-analysis provides evidence that wind instrumentalists demonstrate greater respiratory muscle strength (as measured by MIP and MEP) compared to non-playing controls, while differences in standard pulmonary function tests are less pronounced. The findings suggest that wind instrument playing may serve as a form of specific respiratory muscle training, with adaptations primarily occurring in the muscles responsible for generating and controlling respiratory pressures rather than in lung structure or capacity.

The larger effects observed for respiratory muscle strength measures justify a shift in research focus toward respiratory muscle function when studying wind instrumentalist health and performance. Future studies should employ longitudinal designs to establish causal relationships between wind instrument practice and respiratory adaptations, control for potential confounding factors, and investigate differences between instrument types and playing techniques.

These findings also suggest potential therapeutic applications, as the respiratory muscle strengthening associated with wind instrument playing might be beneficial for individuals with respiratory weakness or conditions requiring respiratory rehabilitation. Further research exploring wind instrument playing as a form of respiratory muscle training in clinical populations could yield valuable insights for rehabilitation purposes.

Overall, this meta-analysis provides useful insights into the respiratory physiology of wind instruments compared to controls, and elucidates a pressing need for further research into respiratory muscle performance for wind instrumentalists.