Primary outcome measure of the The Immediate Effects of Impact Massage at Different Frequencies on Lower Limb Explosive Power and Neuromuscular Performance in Volleyball Players During Exercise

Given the small sample size, correlation analysis was used solely as an exploratory indicator of individual response stability. Results are reported as the r value and its 95% confidence interval, without determining the existence of true individual differences based on individual P-values or fixed thresholds for correlation strength. The 95% confidence intervals for correlation coefficients were estimated using Fisher's z-transformation and validated with non-parametric bootstrap methods.

This analysis included only subjects with complete corresponding data across both repeated cycles. （2）Calculate the standard deviation of individual responses (SDIR) using the following formula, which is used to describe the additional variability in observed values for the 40 Hz condition relative to the 30 Hz condition: Among these, SDIR represents the standard deviation of the additional variability values relative to the 30 Hz condition under the 40 Hz condition; SD40 and SD30  denote the standard deviations of the variability values relative to the corresponding baseline periods under the 40 Hz and 30 Hz conditions, respectively.

The variances of the variability values for each condition are first calculated within repeated periods and then aggregated according to a pre-specified combination method. This estimate relies on the assumption that the two conditions have comparable random measurement errors, natural fluctuations, and within-period variations. Since the 30 Hz condition is an active comparison condition of low intensity rather than a no-treatment control, it may itself exhibit genuine individual responses.

Therefore, SDIR serves only as a supportive and exploratory indicator and cannot be directly interpreted as an unbiased estimate of individual true treatment effect differences. If the variance difference within the square root is less than or equal to 0, the original variance difference and its uncertainty are reported, and SDIR is considered to lie at the 0 boundary. This result indicates that the existing data are insufficient to support the presence of additional variability in the 40 Hz condition beyond that in the 30 Hz condition, and it should not be interpreted as evidence proving the absence of individual response heterogeneity. （3）A linear mixed-effects model was used to analyze the changes in outcome measures over time under different tapping frequency conditions.

The primary analysis employed the change from baseline values during corresponding periods within each cycle as the dependent variable, with the same-day baseline values included as a covariate in the model. Fixed effects included tapping frequency condition, post-treatment time point, condition × time interaction, and repeated cycles coded according to the actual randomization structure. Time points were treated as categorical variables to avoid forcing an ungrounded linear time trend assumption.

Body fat percentage was incorporated as a centralized continuous covariate into the adjustment model only when predefined and biologically justified. Subjects were specified as random intercepts to account for the correlation of repeated measurements within the same subject; for multiple time points within the same treatment cycle, an appropriate residual correlation structure was pre-specified based on measurement order, with sensitivity analyses conducted using alternative correlation structures.

An extended model incorporating subject × condition random slopes was also built to explore individual treatment effect heterogeneity; if the correlation parameter between random intercepts and random slopes could not be stably estimated due to small sample size, an uncorrelated random effects structure was adopted, with singular fits or boundary estimates reported. Due to the small sample size, specific model parameters were estimated using Restricted Maximum Likelihood (REML).

Standard errors of fixed effects, denominator degrees of freedom, and overall tests were corrected for small samples using the Kenward-Roger method. The condition × time interaction was tested using a corrected F-test, and estimated marginal means, inter-condition differences, 95% confidence intervals, and two-sided P-values were reported. （4）Using the method proposed by Senn[61], the pooled within-subject variance was estimated with all repeated crossover data.

The average treatment effect, standard error, and 95% confidence interval for repeated measurements were calculated based on the number of repetitions completed by each subject and the corresponding degrees of freedom. Each subject was then treated as an independent unit for crossover analysis, and the average treatment effects and their standard errors at the subject level were meta-analyzed using the metagen function with a random-effects model.

The overall average treatment effect was weighted by inverse variance, and the between-subject heterogeneity variance τ² was estimated using the REML method. τ, τ², and their 95% confidence intervals were reported. The 95% confidence interval for the overall average effect was adjusted for small samples using the Hartung-Knapp method, and the confidence interval for the heterogeneity variance was calculated using the Q-profile method.

Additionally, a 95% prediction interval was reported to describe the range in which the true average treatment effect of a new, similar subject would lie, assuming the model assumptions hold and the target population is exchangeable with the current sample. Given the small number of subjects, τ, τ², and the prediction interval were interpreted as exploratory results, and individual responses were not determined solely based on the P-value from heterogeneity tests.

All statistical analyses were performed in R 4.6.0. The statistical analysis of this study is divided into two levels: overall average effect analysis and exploration of individual response heterogeneity. First, a linear mixed-effects model is used as the primary analysis to examine the overall effects of tapping frequency, time, and their interaction on each outcome measure.

A mixed model incorporating random slopes for participants × conditions is employed as the primary model for exploring individual heterogeneity, with supportive interpretation provided by participant-level random effects summary results. Pearson correlation coefficients, Bland-Altman agreement analysis, and SDIR are used solely as supplementary, hypothesis-generating analyses. To investigate whether outcome measures exhibit systematic changes with measurement order or repetition count, a hierarchical linear mixed-effects model was employed for analysis.

The model treated outcome measures as the dependent variable, with measurement order or repetition number as continuous independent variables, and included the necessary interaction terms between tapping frequency condition, repetition cycle, and time trend. Participants were specified as random intercepts, and participant-specific time random slopes were added when the model could stably estimate them, allowing for different change trends across participants.

Regression coefficients of the time term along with their 95% confidence intervals were used to evaluate overall learning, familiarity, or fatigue trends. Independent trend regressions constructed for individual participants were not used as primary statistical inferences; if retained as supplementary analyses, they were only used for descriptive purposes and hypothesis generation, with P-values for the time term controlled for false discovery rate (FDR) using the Benjamini-Hochberg method.

Participants were not simply categorized as “trend-positive” or “trend-negative” based on this. Furthermore, utilize the difference between 40 Hz and 30 Hz from individual observations in the first crossover period as the X-axis value, and place the corresponding difference from the second crossover period on the Y-axis to plot a scatter plot for evaluating the consistency of intervention effects across repeated test subjects[63-65].

The data from this study are described as mean ± standard deviation. Primary results are reported with effect estimates, 95% confidence intervals, and exact two-sided P values. The statistical significance criterion for primary analyses is a two-sided P < 0.05; however, P values from exploratory analyses are provided for descriptive purposes only.

Linear mixed-effects models utilize all available observations under the assumption of missing at random, without single-value imputation for missing data; correlation analyses and inter-period consistency analyses are conducted using only participants with complete data for the respective periods.