#### 1. Load packages ----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_25.02.25.csv")#### 3. Process the data ----# Calculate effect sizes using escalc() for all datadata1 <-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 datasetsdata_MIP <- data1 %>%filter(RMP =="MIP")data_MEP <- data1 %>%filter(RMP =="MEP")
1 Process the data
The escalc() function comes from the metafor package and is used for computing effect sizes (e.g., mean differences, odds ratios, etc.) and their sampling variances from supplied raw data.
measure=“MD”: The measure argument specifies the type of effect size to calculate.
Individual Parameters:
m1i=mn_exp: The mean of the experimental group. sd1i=std_exp: The standard deviation of the experimental group. n1i=n_exp: The sample size of the experimental group. m2i=mn_ctl: The mean of the control group. sd2i=std_ctl: The standard deviation of the control group. n2i=n_ctl: The sample size of the control group. These columns should already exist in data1.
slab argument:
slab=paste(author, year, sep=“,”): This slab argument creates study labels that combine the author and year columns with a comma separator. These labels are useful for plotting or referencing studies in meta-analyses. For example, if author=“Smith” and year=2020, the slab would be “Smith, 2020”.
What happens in this step?
The escalc() function:
Calculates the effect size (mean difference) and its variance for each study. Appends these new columns (typically labeled as yi and vi) to the data1 dataframe: yi: The calculated standardized effect size (in this case, MD). vi: The variance of the calculated effect size.
Separating Data Based on RMP
This step uses the piping functionality (%>%) from the dplyr package to filter the rows of data1 into two separate datasets (data_MIP and data_MEP) based on the value of the RMP column.
filter() function: The filter() function is used to select rows in a dataframe based on logical conditions.
Here:
RMP == “MIP” selects rows where the value of RMP is “MIP”, creating the dataset data_MIP. RMP == “MEP” selects rows where the value of RMP is “MEP”, creating the dataset data_MEP. What happens?
The rows of data1 are divided into two new datasets: data_MIP: Contains only rows where RMP is “MIP”. data_MEP: Contains only rows where RMP is “MEP”.
Summary of Purpose:
Effect Size Calculation: The escalc() function calculates the mean difference (MD) effect size and its variance for each row in the dataset (data1), using the provided means, standard deviations, and sample sizes of experimental and control groups.
Dataset Splitting: The data1 dataframe is split into two smaller datasets, data_MIP and data_MEP, based on the grouping variable RMP, which takes on two values: “MIP” and “MEP”. This allows the analysis of subsets of the data based on the RMP grouping.
Code
#### 4. Meta-analysis for MIP ----# Random-effects model with Hartung-Knapp adjustmentMIP_rma <-rma(yi, vi, data = data_MIP, method ="REML", test ="knha")# Fixed-effect modelMIP_fe <-rma(yi, vi, data = data_MIP, method ="FE")# Create meta object for forest plotMIP_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 adjustmentprediction =TRUE, # Include prediction intervalsm ="MD")
2 Meta-analysis for MIP
1. Random-Effects Model with Hartung-Knapp Adjustment
Function: rma()
This function from the metafor package fits a meta-analysis model.
method = “REML”:
Method for fitting the random-effects model. “REML” indicates restricted maximum-likelihood estimation, a common and robust method for estimating between-study variance (τ²).
In random-effects meta-analysis, it’s assumed that studies are drawn from a population of studies with some true heterogeneity (study-to-study variance). The model estimates a single overall effect while accounting for this heterogeneity.
test = “knha”:
This specifies how to calculate confidence intervals and p-values for the model. “knha” stands for the Hartung-Knapp adjustment.
This adjustment is a more robust way of testing the significance of the random-effects model, particularly when there are few studies or high heterogeneity. It adjusts the confidence intervals to be wider and more conservative than traditional methods.
What happens here?
The function fits a random-effects meta-analysis model for the subset data_MIP with mean differences (yi) as the effect sizes.
It accounts for between-study variability using the REML method and adjusts confidence intervals using the Hartung-Knapp method.
2. Fixed-Effect Model
In a fixed-effect meta-analysis, it is assumed that:
All studies estimate a common true effect size (no between-study heterogeneity).
Differences between studies are due solely to sampling error.
Details:
method = “FE”:
This selects the fixed-effect meta-analysis model.
Other parameters:
yi, vi, and data = data_MIP do the same job as in the previous random-effects model, i.e., providing the input data (effect sizes, variances, and dataset).
What happens here?
A fixed-effect model is fitted to data_MIP. No between-study variability (τ²) is considered. Essentially, all the weight in this analysis comes purely from within-study sample sizes (large studies contribute more weight to the pooled estimate).
3. Create a Meta-Analysis Object for Forest Plot
Function: metacont()
This function comes from the meta package and is used to organize meta-analysis data for generating plots (e.g., forest plots) or performing further meta-analytic tasks.
Inputs (n.e, mean.e, sd.e, etc.):
n.e: Sample size of the experimental group for each study. mean.e: Mean of the experimental group (data_MIP\(mn_exp).
sd.e: Standard deviation of the experimental group (data_MIP\)std_exp). n.c: Sample size of the control group for each study. mean.c: Mean of the control group (data_MIP\(mn_ctl).
sd.c: Standard deviation of the control group (data_MIP\)std_ctl). These correspond to the respective columns from data_MIP.
studlab = data_MIP$author: Labels the studies in the output (e.g., in a forest plot) using the study author names from the author column in data_MIP.
Model Parameters:
common = TRUE: Fit a fixed-effect model. random = TRUE: Fit a random-effects model.
Both models are fitted, allowing comparisons between fixed- and random-effects estimates.
method.random.ci = TRUE:
This specifies that the random-effects model should use the Hartung-Knapp adjustment for calculating confidence intervals, similar to what was used in rma() for the random-effects model.
prediction = TRUE:
Includes prediction intervals (not to be confused with confidence intervals).
Prediction intervals provide the range where future study effect sizes are likely to fall, accounting for between-study heterogeneity.
TALK TO TROY ABOUT ABOVE!!!
sm = “MD”:
Specifies the summary measure as Mean Difference.
What happens here?
A meta-analysis object (MIP_meta) is created for data_MIP that combines both fixed- and random-effects models.
This object can be used to create a forest plot, showing: Study-level effect sizes. Overall effect size estimates (fixed and random). Confidence intervals. Prediction intervals for the random-effects model (if random = TRUE).
Code
#### 5. Meta-analysis for MEP ----# Random-effects model with Hartung-Knapp adjustmentMEP_rma <-rma(yi, vi, data = data_MEP, method ="REML", test ="knha")# Fixed-effect modelMEP_fe <-rma(yi, vi, data = data_MEP, method ="FE")# Create meta object for forest plotMEP_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 adjustmentprediction =TRUE, # Include prediction intervalsm ="MD")
3 Meta-analysis for MEP —-
Random-effects model with Hartung-Knapp adjustment
rma() function (from metafor package) fits a random-effects model.
Arguments:
Same as above
Key Output:
The MEP_rma object stores the results of the random-effects model meta-analysis, including pooled effect size, confidence intervals (adjusted using Hartung-Knapp), heterogeneity statistics (e.g., τ², I²), and other relevant measures.
par(mar): Adjusts the margins of the plotting area. The mar argument specifies the size of the margins in the order: bottom, left, top, and right (here c(4, 4, 4, 2)).
This ensures there is enough space around the plot for labels and other elements.
4.1 Forest Plot for MEP_meta
Parameters in the forest() Function:
MIP_meta
Input: This is the meta-analysis object generated using metacont(), which contains the meta-analytic data (effect sizes, confidence intervals, heterogeneity statistics, etc.).
The plot will visualize the results stored in this meta-analysis object for MIP.
Left-Side Columns and Labels
leftcols: Specifies the data to be displayed in the left side of the forest plot. Here: studlab: Study labels (e.g., authors’ names or study identifiers). n.e, mean.e, sd.e: Sample size (n), mean, and standard deviation (SD) for the experimental group. n.c, mean.c, sd.c: Sample size, mean, and SD for the control group. leftlabs: Customizes the labels for these columns in the plot: “Author” is the header for studlab. “n”, “Mean”, and “SD” are headers for the corresponding columns for both groups.
Right-Side Columns and Labels
rightcols: Specifies the data to be displayed on the right side of the forest plot.
effect: The effect size (mean difference, MD, in this case). ci: The confidence interval of the effect size. rightlabs: Customizes the labels for these columns in the plot: “MD” is the label for the effect size. “95% CI” is the label for the confidence interval.
Results to Include
R comb.fixed = TRUE, comb.random = TRUE,
comb.fixed and comb.random: Whether to show the pooled effect size (summary/combined effect size) for fixed- and random-effects models. Both are set to TRUE, so results for both models are displayed.
R prediction = TRUE,
prediction: Displays the prediction interval (the range in which the effect size of a new study is expected to lie). This is useful for assessing the generalizability of the results.
Heterogeneity Statistics
R print.tau2 = TRUE, print.I2 = TRUE, print.H = TRUE,
print.tau2: Displays the τ² (tau-squared) statistic, which quantifies the variance of the true effect sizes across studies (i.e., heterogeneity). print.I2: Displays the I² statistic, which represents the proportion of total variation in observed effects due to heterogeneity. print.H: Displays the H² statistic, which is another measure of heterogeneity (a ratio comparing total variability to sampling variability).
Colors and Customization
R col.predict = “red”,
col.predict: Color of the prediction interval. It is set to red, making it easily distinguishable from the confidence intervals.
R col.diamond = “blue”,
col.diamond: Color of the diamonds that represent the pooled effect sizes for the fixed- and random-effects models. It is set to blue.
Options for Overall Statistics
R hetstat = TRUE, overall = TRUE, overall.hetstat = TRUE,
hetstat: If TRUE, shows the heterogeneity statistics (I², τ², H²) below the forest plot. overall: If TRUE, includes the overall summary effect in the forest plot. overall.hetstat: If TRUE, includes heterogeneity statistics for the overall summary effect.
Hypothesis Tests
R test.overall.common = TRUE, test.overall.random = TRUE,
test.overall.common: If TRUE, displays the p-value testing the null hypothesis that the pooled effect size (fixed-effects model) is zero. test.overall.random: If TRUE, displays the p-value testing the null hypothesis that the pooled effect size (random-effects model) is zero.
Title and Font Size
R main = “Maximum Inspiratory Pressure (MIP) generation in wind instrumentalists vs. controls”, fontsize = 10
main: The main title of the forest plot. fontsize: Adjusts the font size of text in the forest plot.
4.2 Forest Plot for MEP_meta
R forest(MEP_meta, … main = “Maximum Expiratory Pressure (MEP) generation in wind instrumentalists vs. controls”, …)
The plot for MEP_meta is identical to the MIP_meta plot, except for the title (main). The process and parameters are the same as explained above, tailored for the maximum expiratory pressure dataset.
##Summary
These forest plots give a comprehensive visualization of the meta-analysis results for MIP and MEP:
Individual studies with their effect sizes and confidence intervals are displayed.
Pooled results for fixed- and random-effects models are included. Heterogeneity statistics (τ², I², H²) and prediction intervals are shown.
Colors and labels are customized for clarity.
The plots allow easy comparison of results across studies while assessing overall patterns and heterogeneity.
Code
#### 7. Summary Statistics ----# Display summary statistics for MIPcat("\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
cat(): Displays custom text in the console without quotes. adds a newline at the beginning for better readability in the output. This prints the header for the MIP summary.
Random-Effects Model Statistics for MIP
R cat(“Random-effects model (with Hartung-Knapp adjustment):”) print(MIP_rma)
Hartung-Knapp adjustment: This is a modification of the random-effects model to provide more accurate confidence intervals, especially when the number of studies is small. The random-effects model accounts for both within-study and between-study variations.
print(MIP_rma): Outputs the summary of the random-effects meta-analysis object (MIP_rma).
This typically includes: Pooled effect size (e.g., standardized mean difference or mean difference). Confidence intervals for the pooled effect size. Test statistics (e.g., z-value or p-value for the null hypothesis that the effect size is 0). Measures of heterogeneity (e.g., I², H², τ²). These results are generated by the metafor package in R.
Fixed-Effects Model Statistics for MIP
R cat(“-effect model:”) print(MIP_fe)
Fixed-effects models assume that all studies estimate the same underlying effect size (i.e., no between-study variance), unlike random-effects models.
print(MIP_fe): Shows the summary of the fixed-effects meta-analysis object (MIP_fe). This includes the pooled effect size, confidence intervals, and test statistics based on the fixed-effects assumption.
Heterogeneity statistics provide information on how much the effects vary between studies: I²: Percentage of total variability in effect sizes due to between-study variation (heterogeneity). Values range from 0% (no heterogeneity) to 100% (high heterogeneity). MIP_rma\(I2: Extracts the I² value from the MIP_rma object.
formatC(MIP_rma\)I2, digits=1, format=“f”): Formats the I² value to one decimal place as a percentage. H²: The ratio of total variability to within-study variability. Larger values of H² indicate more heterogeneity. Extracted using MIP_rma$H2 and formatted to 2 decimal places. τ²: The estimated between-study variance. Indicates the extent of variability in effect sizes across the studies.
Extracted using MIP_rma$tau2 and formatted to 4 decimal places. These statistics help assess whether a random-effects model is appropriate.
5.2 Section for MEP
The code here mirrors the structure for MIP, but the operations are applied to the Maximum Expiratory Pressure (MEP) data.
R cat(“for Maximum Expiratory Pressure (MEP):”)
Prints the header for the MEP summary.
Random-Effects Model Statistics for MEP
R cat(“Random-effects model (with Hartung-Knapp adjustment):”) print(MEP_rma)
Displays the results of the random-effects meta-analysis for MEP, including pooled effect size, confidence intervals, and heterogeneity measures.
Fixed-Effects Model Statistics for MEP
R cat(“-effect model:”) print(MEP_fe)
Displays the results of the fixed-effects meta-analysis for MEP.
Outputs heterogeneity statistics for the MEP random-effects meta-analysis: I²: Percentage of heterogeneity among studies. H²: Variance ratio for total and within-study variability. τ²: Between-study variance (absolute measure of heterogeneity).
The values are extracted from the MEP_rma object and formatted similarly to those for MIP.
Key Functions and Outputs
cat(): Used to print custom text and formatted results. It ensures good readability for the console output. print(): Displays meta-analysis object summaries (created via metafor or similar packages). These include effect sizes, confidence intervals, and heterogeneity. formatC(): Formats numerical values with specific formatting options (digits, format): digits: Number of decimal places. format: Specifies the format of the output (e.g., “f” for fixed-point notation). Objects (MIP_rma, MIP_fe, MEP_rma, MEP_fe): Contain the results of meta-analyses (random- and fixed-effects models) and heterogeneity statistics.
Source Code
---title: "Meta Analysis Rerun 21.02.25 EXPLANATIONS"author: "Sarah Morris"date: "2025-01-04"format: html: toc: true toc-depth: 2 toc-title: "Table of Contents" toc-location: right number-sections: true theme: cosmo code-fold: true code-tools: true highlight-style: githubexecute: echo: true warning: false error: false---```{r}#### 1. Load packages ----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_25.02.25.csv")#### 3. Process the data ----# Calculate effect sizes using escalc() for all datadata1 <-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 datasetsdata_MIP <- data1 %>%filter(RMP =="MIP")data_MEP <- data1 %>%filter(RMP =="MEP")```# Process the dataThe escalc() function comes from the metafor package and is used for computing effect sizes (e.g., mean differences, odds ratios, etc.) and their sampling variances from supplied raw data.measure="MD": The measure argument specifies the type of effect size to calculate.*Individual Parameters:*m1i=mn_exp: The mean of the experimental group. sd1i=std_exp: The standard deviation of the experimental group. n1i=n_exp: The sample size of the experimental group. m2i=mn_ctl: The mean of the control group. sd2i=std_ctl: The standard deviation of the control group. n2i=n_ctl: The sample size of the control group. These columns should already exist in data1.*slab argument:*slab=paste(author, year, sep=", "): This slab argument creates study labels that combine the author and year columns with a comma separator. These labels are useful for plotting or referencing studies in meta-analyses. For example, if author="Smith" and year=2020, the slab would be "Smith, 2020".*What happens in this step?*The escalc() function:Calculates the effect size (mean difference) and its variance for each study. Appends these new columns (typically labeled as yi and vi) to the data1 dataframe: yi: The calculated standardized effect size (in this case, MD). vi: The variance of the calculated effect size.**Separating Data Based on RMP**This step uses the piping functionality (%\>%) from the dplyr package to filter the rows of data1 into two separate datasets (data_MIP and data_MEP) based on the value of the RMP column.filter() function: The filter() function is used to select rows in a dataframe based on logical conditions.Here:RMP == "MIP" selects rows where the value of RMP is "MIP", creating the dataset data_MIP. RMP == "MEP" selects rows where the value of RMP is "MEP", creating the dataset data_MEP. What happens?The rows of data1 are divided into two new datasets: data_MIP: Contains only rows where RMP is "MIP". data_MEP: Contains only rows where RMP is "MEP".**Summary of Purpose:**Effect Size Calculation: The escalc() function calculates the mean difference (MD) effect size and its variance for each row in the dataset (data1), using the provided means, standard deviations, and sample sizes of experimental and control groups.Dataset Splitting: The data1 dataframe is split into two smaller datasets, data_MIP and data_MEP, based on the grouping variable RMP, which takes on two values: "MIP" and "MEP". This allows the analysis of subsets of the data based on the RMP grouping.```{r}#### 4. Meta-analysis for MIP ----# Random-effects model with Hartung-Knapp adjustmentMIP_rma <-rma(yi, vi, data = data_MIP, method ="REML", test ="knha")# Fixed-effect modelMIP_fe <-rma(yi, vi, data = data_MIP, method ="FE")# Create meta object for forest plotMIP_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 adjustmentprediction =TRUE, # Include prediction intervalsm ="MD")```# Meta-analysis for MIP**1. Random-Effects Model with Hartung-Knapp Adjustment**Function: rma()This function from the metafor package fits a meta-analysis model.*method = "REML":*Method for fitting the random-effects model. "REML" indicates restricted maximum-likelihood estimation, a common and robust method for estimating between-study variance (τ²).In random-effects meta-analysis, it’s assumed that studies are drawn from a population of studies with some true heterogeneity (study-to-study variance). The model estimates a single overall effect while accounting for this heterogeneity.*test = "knha":*This specifies how to calculate confidence intervals and p-values for the model. "knha" stands for the Hartung-Knapp adjustment.This adjustment is a more robust way of testing the significance of the random-effects model, particularly when there are few studies or high heterogeneity. It adjusts the confidence intervals to be wider and more conservative than traditional methods.*What happens here?*The function fits a random-effects meta-analysis model for the subset data_MIP with mean differences (yi) as the effect sizes.It accounts for between-study variability using the REML method and adjusts confidence intervals using the Hartung-Knapp method.**2. Fixed-Effect Model**In a fixed-effect meta-analysis, it is assumed that:All studies estimate a common true effect size (no between-study heterogeneity).Differences between studies are due solely to sampling error.*Details:*method = "FE":This selects the fixed-effect meta-analysis model.Other parameters:yi, vi, and data = data_MIP do the same job as in the previous random-effects model, i.e., providing the input data (effect sizes, variances, and dataset).*What happens here?*A fixed-effect model is fitted to data_MIP. No between-study variability (τ²) is considered. Essentially, all the weight in this analysis comes purely from within-study sample sizes (large studies contribute more weight to the pooled estimate).**3. Create a Meta-Analysis Object for Forest Plot**Function: metacont()This function comes from the meta package and is used to organize meta-analysis data for generating plots (e.g., forest plots) or performing further meta-analytic tasks.*Inputs (n.e, mean.e, sd.e, etc.):*n.e: Sample size of the experimental group for each study. mean.e: Mean of the experimental group (data_MIP$mn_exp).sd.e: Standard deviation of the experimental group (data_MIP$std_exp). n.c: Sample size of the control group for each study. mean.c: Mean of the control group (data_MIP$mn_ctl).sd.c: Standard deviation of the control group (data_MIP$std_ctl). These correspond to the respective columns from data_MIP.studlab = data_MIP\$author: Labels the studies in the output (e.g., in a forest plot) using the study author names from the author column in data_MIP.*Model Parameters:*common = TRUE: Fit a fixed-effect model. random = TRUE: Fit a random-effects model.Both models are fitted, allowing comparisons between fixed- and random-effects estimates.method.random.ci = TRUE:This specifies that the random-effects model should use the Hartung-Knapp adjustment for calculating confidence intervals, similar to what was used in rma() for the random-effects model.prediction = TRUE:Includes prediction intervals (not to be confused with confidence intervals).Prediction intervals provide the range where future study effect sizes are likely to fall, accounting for between-study heterogeneity.TALK TO TROY ABOUT ABOVE!!!sm = "MD":Specifies the summary measure as Mean Difference.*What happens here?*A meta-analysis object (MIP_meta) is created for data_MIP that combines both fixed- and random-effects models.This object can be used to create a forest plot, showing: Study-level effect sizes. Overall effect size estimates (fixed and random). Confidence intervals. Prediction intervals for the random-effects model (if random = TRUE).```{r}#### 5. Meta-analysis for MEP ----# Random-effects model with Hartung-Knapp adjustmentMEP_rma <-rma(yi, vi, data = data_MEP, method ="REML", test ="knha")# Fixed-effect modelMEP_fe <-rma(yi, vi, data = data_MEP, method ="FE")# Create meta object for forest plotMEP_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 adjustmentprediction =TRUE, # Include prediction intervalsm ="MD")```# Meta-analysis for MEP ----1. Random-effects model with Hartung-Knapp adjustmentrma() function (from metafor package) fits a random-effects model.*Arguments:*Same as above*Key Output:*The MEP_rma object stores the results of the random-effects model meta-analysis, including pooled effect size, confidence intervals (adjusted using Hartung-Knapp), heterogeneity statistics (e.g., τ², I²), and other relevant measures.```{r}#### 6. Create Forest Plots ----# MIP Forest Plotpar(mar =c(4, 4, 4, 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 statisticcol.predict ="red", # Prediction interval in redcol.diamond ="blue", # Confidence interval in bluehetstat =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 =10)# MEP Forest Plotforest(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, # Add H^2 statisticcol.predict ="red", # Prediction interval in redcol.diamond ="blue", # Confidence interval in bluehetstat =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 =10)```# Creating Forest Plots**General Setup for MIP Forest Plot**par(mar): Adjusts the margins of the plotting area. The mar argument specifies the size of the margins in the order: bottom, left, top, and right (here c(4, 4, 4, 2)).This ensures there is enough space around the plot for labels and other elements.## Forest Plot for MEP_meta**Parameters in the forest() Function:**MIP_metaInput: This is the meta-analysis object generated using metacont(), which contains the meta-analytic data (effect sizes, confidence intervals, heterogeneity statistics, etc.).The plot will visualize the results stored in this meta-analysis object for MIP.**Left-Side Columns and Labels**leftcols: Specifies the data to be displayed in the left side of the forest plot. Here: studlab: Study labels (e.g., authors' names or study identifiers). n.e, mean.e, sd.e: Sample size (n), mean, and standard deviation (SD) for the experimental group. n.c, mean.c, sd.c: Sample size, mean, and SD for the control group. leftlabs: Customizes the labels for these columns in the plot: "Author" is the header for studlab. "n", "Mean", and "SD" are headers for the corresponding columns for both groups.**Right-Side Columns and Labels**rightcols: Specifies the data to be displayed on the right side of the forest plot.effect: The effect size (mean difference, MD, in this case). ci: The confidence interval of the effect size. rightlabs: Customizes the labels for these columns in the plot: "MD" is the label for the effect size. "95% CI" is the label for the confidence interval.**Results to Include**R comb.fixed = TRUE, comb.random = TRUE,comb.fixed and comb.random: Whether to show the pooled effect size (summary/combined effect size) for fixed- and random-effects models. Both are set to TRUE, so results for both models are displayed.R prediction = TRUE,prediction: Displays the prediction interval (the range in which the effect size of a new study is expected to lie). This is useful for assessing the generalizability of the results.Heterogeneity StatisticsR print.tau2 = TRUE, print.I2 = TRUE, print.H = TRUE,print.tau2: Displays the τ² (tau-squared) statistic, which quantifies the variance of the true effect sizes across studies (i.e., heterogeneity). print.I2: Displays the I² statistic, which represents the proportion of total variation in observed effects due to heterogeneity. print.H: Displays the H² statistic, which is another measure of heterogeneity (a ratio comparing total variability to sampling variability).**Colors and Customization**R col.predict = "red",col.predict: Color of the prediction interval. It is set to red, making it easily distinguishable from the confidence intervals.R col.diamond = "blue",col.diamond: Color of the diamonds that represent the pooled effect sizes for the fixed- and random-effects models. It is set to blue.**Options for Overall Statistics**R hetstat = TRUE, overall = TRUE, overall.hetstat = TRUE,hetstat: If TRUE, shows the heterogeneity statistics (I², τ², H²) below the forest plot. overall: If TRUE, includes the overall summary effect in the forest plot. overall.hetstat: If TRUE, includes heterogeneity statistics for the overall summary effect.**Hypothesis Tests**R test.overall.common = TRUE, test.overall.random = TRUE,test.overall.common: If TRUE, displays the p-value testing the null hypothesis that the pooled effect size (fixed-effects model) is zero. test.overall.random: If TRUE, displays the p-value testing the null hypothesis that the pooled effect size (random-effects model) is zero.**Title and Font Size**R main = "Maximum Inspiratory Pressure (MIP) generation in wind instrumentalists vs. controls", fontsize = 10main: The main title of the forest plot. fontsize: Adjusts the font size of text in the forest plot.## Forest Plot for MEP_metaR forest(MEP_meta, ... main = "Maximum Expiratory Pressure (MEP) generation in wind instrumentalists vs. controls", ...)The plot for MEP_meta is identical to the MIP_meta plot, except for the title (main). The process and parameters are the same as explained above, tailored for the maximum expiratory pressure dataset.##SummaryThese forest plots give a comprehensive visualization of the meta-analysis results for MIP and MEP:Individual studies with their effect sizes and confidence intervals are displayed.Pooled results for fixed- and random-effects models are included. Heterogeneity statistics (τ², I², H²) and prediction intervals are shown.Colors and labels are customized for clarity.The plots allow easy comparison of results across studies while assessing overall patterns and heterogeneity.```{r}#### 7. Summary Statistics ----# Display summary statistics for MIPcat("\nSummary for Maximum Inspiratory Pressure (MIP):\n")cat("Random-effects model (with Hartung-Knapp adjustment):\n")print(MIP_rma)cat("\nFixed-effect model:\n")print(MIP_fe)cat("\nHeterogeneity statistics:\n")cat("I² =", formatC(MIP_rma$I2, digits=1, format="f"), "%\n")cat("H² =", formatC(MIP_rma$H2, digits=2, format="f"), "\n")cat("τ² =", formatC(MIP_rma$tau2, digits=4, format="f"), "\n\n")# Display summary statistics for MEPcat("\nSummary for Maximum Expiratory Pressure (MEP):\n")cat("Random-effects model (with Hartung-Knapp adjustment):\n")print(MEP_rma)cat("\nFixed-effect model:\n")print(MEP_fe)cat("\nHeterogeneity statistics:\n")cat("I² =", formatC(MEP_rma$I2, digits=1, format="f"), "%\n")cat("H² =", formatC(MEP_rma$H2, digits=2, format="f"), "\n")cat("τ² =", formatC(MEP_rma$tau2, digits=4, format="f"), "\n")```# Summary Statistics## Section for MIPR cat("\nSummary for Maximum Inspiratory Pressure (MIP):\n")cat(): Displays custom text in the console without quotes. \n adds a newline at the beginning for better readability in the output. This prints the header for the MIP summary.**Random-Effects Model Statistics for MIP**R cat("Random-effects model (with Hartung-Knapp adjustment):\n") print(MIP_rma)Hartung-Knapp adjustment: This is a modification of the random-effects model to provide more accurate confidence intervals, especially when the number of studies is small. The random-effects model accounts for both within-study and between-study variations.print(MIP_rma): Outputs the summary of the random-effects meta-analysis object (MIP_rma).This typically includes: Pooled effect size (e.g., standardized mean difference or mean difference). Confidence intervals for the pooled effect size. Test statistics (e.g., z-value or p-value for the null hypothesis that the effect size is 0). Measures of heterogeneity (e.g., I², H², τ²). These results are generated by the metafor package in R.**Fixed-Effects Model Statistics for MIP**R cat("\nFixed-effect model:\n") print(MIP_fe)Fixed-effects models assume that all studies estimate the same underlying effect size (i.e., no between-study variance), unlike random-effects models.print(MIP_fe): Shows the summary of the fixed-effects meta-analysis object (MIP_fe). This includes the pooled effect size, confidence intervals, and test statistics based on the fixed-effects assumption.**Heterogeneity Statistics for MIP**R cat("\nHeterogeneity statistics:\n") cat("I² =", formatC(MIP_rma$I2, digits=1, format="f"), "%\n")cat("H² =", formatC(MIP_rma$H2, digits=2, format="f"), "\n") cat("τ² =", formatC(MIP_rma\$tau2, digits=4, format="f"), "\n\n")Heterogeneity statistics provide information on how much the effects vary between studies: I²: Percentage of total variability in effect sizes due to between-study variation (heterogeneity). Values range from 0% (no heterogeneity) to 100% (high heterogeneity). MIP_rma$I2: Extracts the I² value from the MIP_rma object.formatC(MIP_rma$I2, digits=1, format="f"): Formats the I² value to one decimal place as a percentage. H²: The ratio of total variability to within-study variability. Larger values of H² indicate more heterogeneity. Extracted using MIP_rma\$H2 and formatted to 2 decimal places. τ²: The estimated between-study variance. Indicates the extent of variability in effect sizes across the studies.Extracted using MIP_rma\$tau2 and formatted to 4 decimal places. These statistics help assess whether a random-effects model is appropriate.## Section for MEPThe code here mirrors the structure for MIP, but the operations are applied to the Maximum Expiratory Pressure (MEP) data.R cat("\nSummary for Maximum Expiratory Pressure (MEP):\n")Prints the header for the MEP summary.**Random-Effects Model Statistics for MEP**R cat("Random-effects model (with Hartung-Knapp adjustment):\n") print(MEP_rma)Displays the results of the random-effects meta-analysis for MEP, including pooled effect size, confidence intervals, and heterogeneity measures.**Fixed-Effects Model Statistics for MEP**R cat("\nFixed-effect model:\n") print(MEP_fe)Displays the results of the fixed-effects meta-analysis for MEP.**Heterogeneity Statistics for MEP**R cat("\nHeterogeneity statistics:\n") cat("I² =", formatC(MEP_rma$I2, digits=1, format="f"), "%\n")cat("H² =", formatC(MEP_rma$H2, digits=2, format="f"), "\n") cat("τ² =", formatC(MEP_rma\$tau2, digits=4, format="f"), "\n")Outputs heterogeneity statistics for the MEP random-effects meta-analysis: I²: Percentage of heterogeneity among studies. H²: Variance ratio for total and within-study variability. τ²: Between-study variance (absolute measure of heterogeneity).The values are extracted from the MEP_rma object and formatted similarly to those for MIP.**Key Functions and Outputs**cat(): Used to print custom text and formatted results. It ensures good readability for the console output. print(): Displays meta-analysis object summaries (created via metafor or similar packages). These include effect sizes, confidence intervals, and heterogeneity. formatC(): Formats numerical values with specific formatting options (digits, format): digits: Number of decimal places. format: Specifies the format of the output (e.g., "f" for fixed-point notation). Objects (MIP_rma, MIP_fe, MEP_rma, MEP_fe): Contain the results of meta-analyses (random- and fixed-effects models) and heterogeneity statistics.
