The t-test, introduced by William Sealy Gosset in 1908 under the pseudonym “Student,” is a parametric statistical test designed to compare the means of two groups under specific assumptions [132]. The standard t-test, also called the independent t-test, evaluates the difference between two independent group means. The test statistic is calculated as: t = X ¯ 1 − X ¯ 2 / s p 2 ( n 1 − 1 + n 2 − 1 ) , where X ¯ 1 and X ¯ 2 are the sample means, n 1 and n 2 are the sample sizes, and s p 2 = ( n 1 − 1 ) s 1 2 + ( n 2 − 1 ) s 2 2 / n 1 + n 2 − 2 is the pooled variance. Here, s 1 2 and s 2 2 represent the variances of the two groups. The paired t-test, a variation of the standard t-test, is used for dependent or matched samples, analyzing the mean differences between paired observations using the formula: t = d ¯ / s d / n , where d ¯ is the mean of the differences, s d is the standard deviation of the differences, and n is the number of paired observations. Another important extension, Welch’s t-test, modifies the standard t-test to account for unequal variances and sample sizes between groups, with its test statistic computed as: t = X ¯ 1 − X ¯ 2 / ( s 1 2 / n 1 ) + ( s 2 2 / n 2 ) . These methods, which rely on the assumption of normally distributed data, have become foundational tools for statistical inference in clinical trials and experimental research [150].

T-tests are widely recognized for their simplicity, computational efficiency, and effectiveness in comparing means, particularly in small to moderate-sized datasets (Cf. Table 1). The paired t-test enhances statistical power by focusing on within-subject differences, while Welch’s t-test provides a robust alternative for scenarios involving unequal variances or sample sizes. However, these methods rely on assumptions of normality and are sensitive to outliers, which can compromise their validity in non-normal or heavily skewed data. Additionally, t-tests are limited to mean-based comparisons, rendering them unsuitable for ordinal data or high-dimensional datasets that require more sophisticated statistical approaches, such as mixed-effects models or machine learning techniques.

Analysis of Variance (ANOVA), introduced by Fisher in 1918, is a parametric method for comparing the means of multiple groups. It partitions total variation into between-group variation ( S S between ) and within-group variation ( S S within ), with the F-statistic calculated as F = S S between ( N − k ) / S S within ( k − 1 ) where k is the number of groups and N is the total number of observations. Under the null hypothesis ( H 0 ), the F-statistic follows an F-distribution with k − 1 and N − k degrees of freedom. ANOVA assumes normality within groups, homoscedasticity, and independence of observations. Extensions like repeated measures and two-way ANOVA have been developed for longitudinal and factorial designs. While ANOVA controls Type I error better than multiple pairwise t-tests, violations of its assumptions can lead to misleading results. Nonparametric alternatives like the Kruskal-Wallis test address these limitations. ANOVA is also less powerful with unequal group sizes.

Analysis of Covariance (ANCOVA), introduced by Fisher in 1925, combines ANOVA and linear regression to compare group means while adjusting for one or more covariates. ANCOVA is particularly useful in clinical trials to control for baseline imbalances (Cf. Table 2). The model is Y i j = μ + α i + β X i j + ϵ i j where Y i j is the outcome, X i j is the covariate, and β is the regression coefficient for the covariate. ANCOVA adjusts the outcome by removing the effect of the covariate, enabling a comparison of group means for a common covariate value. It assumes linearity between the covariate and the outcome, homogeneity of regression slopes, and normality and homoscedasticity of residuals. ANCOVA has been extended to handle multiple covariates and more complex designs. However, violations of assumptions or failure to account for covariate-treatment interactions can lead to biased results.

LMMs, introduced by [61], extend linear regression by incorporating both fixed and random effects, making them particularly suitable for analyzing hierarchical or longitudinal data, such as repeated measures in clinical trials. Fixed effects represent population-level relationships, while random effects capture subject-specific variability, allowing LMMs to effectively model within-subject or within-cluster correlations. The general model is expressed as: Y = X β + Z b + ϵ, where X β denotes the fixed effects, Z b represents the random effects ( b ∼ N ( 0 , G )), and ϵ ∼ N ( 0 , R ) accounts for residual errors. One of the significant strengths of LMMs is their ability to handle missing data under the Missing at Random (MAR) assumption, where the probability of missingness depends only on observed data. Unlike traditional methods that exclude incomplete cases, LMMs use all available data through maximum likelihood or restricted maximum likelihood estimation, reducing bias and improving efficiency. Additionally, LMMs accommodate complex covariance structures, such as autoregressive or unstructured correlations, providing flexibility in modeling dependencies.

LMMs, while versatile, have limitations. They assume linear relationships, normality of random effects and residuals, and independence between random and fixed effects. While diagnostic checks, transformations, and robust estimation techniques can address some violations, these assumptions inherently restrict their use in datasets with non-linear relationships or non-normal distributions. Computational challenges with large datasets or complex random-effects structures can be mitigated using efficient algorithms like Expectation-Maximization or parallel computing. However, LMMs struggle with data missing not at random (MNAR), as the mechanism depends on unobserved variables, and biased estimates may result when the independence of random and fixed effects is violated. Additionally, interpreting random effects in complex hierarchical models can be challenging, limiting subject-specific inferences. These limitations underscore the need to carefully evaluate assumptions and choose appropriate methods for specific research contexts.

Generalized Estimating Equations (GEEs), introduced by Liang and Zeger in 1986, are an extension of generalized linear models (GLMs) designed to handle correlated or clustered data, such as repeated measures or longitudinal observations. Unlike LMMs, which explicitly include random effects to account for individual-level variability, GEEs focus on estimating population-averaged effects, making them ideal for studies where the primary interest lies in overall population trends rather than subject-specific inferences. The general form of a GEE is: g ( μ i j ) = X i j β, where g ( · ) is the link function (e.g., logit for binary outcomes, log for count data, identity for continuous data), μ i j represents the mean response for the j-th observation of the i-th cluster, X i j is the covariate matrix, and β denotes the regression coefficients. GEEs account for within-cluster correlation by specifying a working correlation structure, such as exchangeable, autoregressive, or unstructured. While the correct specification of this structure enhances efficiency, GEEs remain robust even if it is misspecified, providing consistent estimates of regression coefficients.

GEEs are highly flexible, accommodating various outcome types—binary, count, and continuous—via appropriate link functions. They are computationally efficient as they avoid explicit random-effects modeling and are robust to correlation structure misspecification, ensuring reliable population-level inferences. However, they assume data are Missing Completely at Random (MCAR), a stricter condition than the MAR assumption used in LMMs. GEEs focus on population-averaged effects, limiting their ability to provide individual-level inferences or subject-specific predictions. Additionally, their efficiency relies on correctly specifying the working correlation matrix, and they cannot naturally handle complex random-effect structures, restricting their use in hierarchical settings. Despite these limitations, GEEs are widely applied in clinical and epidemiological studies for analyzing correlated data and modeling average treatment effects.

Linear Mixed-Effects Models (LMMs) and Generalized Estimating Equations (GEEs) are widely used for analyzing longitudinal and hierarchical data. However, complex relationships often require extensions such as Nonlinear Mixed-Effects Models (NLMMs) and Generalized Additive Mixed Models (GAMMs). NLMMs extend LMMs by incorporating nonlinear functions to model complex biological processes like drug absorption and elimination in pharmacokinetics. These models take the form Y i j = f ( X i j , β , b i ) + ϵ i j , where f ( · ) is a nonlinear function, β represents fixed effects, b i accounts for random effects, and ϵ i j captures residual errors. Tools such as NONMEM, Monolix, and the R package nlme facilitate the application of NLMMs, although their computational demands are significant [104, 27]. GAMMs combine smooth, nonparametric functions with mixed effects to model unknown or varying relationships. The general form of a GAMM is: Y i j = X i j β + Z i j b i + ∑ k f k ( X i j k ) + ϵ i j , where ∑ k f k ( X i j k ) represents smooth functions of predictors, and Z i j b i captures the random effects. GAMMs are commonly used in environmental and disease progression studies but may face overfitting and computational challenges in large datasets [154].

NLMMs and GAMMs expand the flexibility of LMMs and GEEs by accommodating nonlinear and nonparametric relationships, making them invaluable for modeling complex processes in clinical and epidemiological research. However, careful consideration is needed to address computational demands and ensure model interpretability.

Accelerated Failure Time (AFT) models are commonly used to analyze time-to-event data by directly modeling the logarithm of survival times as a linear function of covariates. The general form of an AFT model is: log ( T i ) = X i β + ϵ i , where T i is the survival time for the i-th individual, X i is the covariate vector, β is the vector of regression coefficients and ϵ i is the error term. AFT models assume that covariates accelerate or decelerate the time to the event by a constant factor. For example, if β > 0, the covariate increases survival time by a factor e β ; if β < 0, it decreases survival time by the same factor. Common distributions for ϵ yield specific AFT models, such as the Weibull AFT or log-logistic AFT.

When longitudinal data is integrated with survival response, Joint Models (JMs) are used to address the dependency between longitudinal trajectories and time-to-event processes, offering a unified framework for analyzing such data. The general joint modeling framework consists of two submodels: the longitudinal submodel, given by Y i ( t ) = X i ⊤ ( t ) β + Z i ⊤ ( t ) b i + ϵ i ( t ), where Y i ( t ) is the longitudinal outcome for subject i at time t, X i ( t ) and Z i ( t ) are design matrices for fixed (β) and random ( b i ) effects, respectively, ϵ i ( t ) ∼ N ( 0 , σ 2 ) represents residual errors, and the survival submodel h i ( t | M i ( t ) , w i ) = h 0 ( t ) exp ( w i ⊤ γ + α M i ( t ) ), where h i ( t ) is the hazard function for subject i at time t, h 0 ( t ) is the baseline hazard, w i is a vector of baseline covariates with coefficients γ, M i ( t ) is a function of the subject-specific longitudinal trajectory, linking the two submodels via the association parameter α. JMs provide more accurate estimates of survival outcomes and dynamic predictions of survival probabilities by incorporating longitudinal biomarker information, thus reducing bias compared to separate analyses. However, they rely on strong parametric assumptions about the underlying processes, which, if misspecified, can lead to biased results. Fitting JMs is computationally intensive, requiring specialized algorithms like the Expectation-Maximization (EM) algorithm or Bayesian Markov Chain Monte Carlo (MCMC) methods. Recent advancements have addressed these challenges by developing efficient algorithms, such as the R package JSM [156], which provides a semiparametric joint modeling framework and offers a maximum likelihood approach to fit these models.

While AFT models and the proportional hazards framework remain the most widely used for survival data, several extensions have been proposed to address more complex event structures. When the proportional hazards assumption is violated, alternatives such as time-varying coefficient models, stratified Cox models, Random Survival Forest and Gradient Boosting type methods provide flexible inference by allowing hazard ratios to change over time [62, 137, 50]. In studies with complex censoring patterns, including interval-censored, left-truncated, or competing-risk data, specialized likelihood-based and nonparametric estimators have been developed, such as the Fine–Gray subdistribution model for competing risks [41]. Furthermore, when participants may experience multiple or repeated events, recurrent-event models such as the Andersen–Gill counting-process formulation, Prentice–Williams–Peterson (PWP) total- or gap-time models, and Wei–Lin–Weissfeld (WLW) marginal models offer distinct strategies for capturing dependence among event recurrences [100]. Together, these approaches expand survival modeling beyond proportional hazards assumptions and accommodate the complexities frequently encountered in modern clinical trials.

The Wilcoxon Signed-Rank Test [152] is a non-parametric alternative to the paired t-test, used to compare paired observations when parametric assumptions are violated. It tests whether the median of the differences between paired values is zero by ranking the absolute values of the differences and summing the signed ranks. The test statistic is calculated as W = ∑ i = 1 n R i · sgn ( D i ), where D i = X i − Y i and sgn ( D i ) is the sign of the difference. For small samples, critical values are from exact tables, while for large samples, the statistic approximates a normal distribution. This test is robust for small sample sizes and ordinal data but assumes symmetry in the differences, which if violated, can lead to biased results. It is less powerful than the paired t-test for normally distributed data.

The Mann-Whitney U Test [88] compares two independent groups, assessing if one group tends to have larger values than the other. It is particularly useful for ordinal data or skewed distributions. The U statistic is calculated based on ranks of the pooled data, with the formula: U = n 1 n 2 + n 1 ( n 1 + 1 ) 2 − ∑ i = 1 n 1 R ( X i ), where R ( X i ) is the rank of X i in the combined sample. The test assumes that the groups are independent and that the distributions have the same shape and scale. For small sample sizes, critical values are obtained from exact tables, while larger samples use a normal approximation. It is less efficient than the t-test for normally distributed data, trading power for robustness. Furthermore, it assumes that the groups are independent and the distributions of the two groups have the same shape and scale; otherwise, the results may reflect differences in spread rather than central tendency.

Nonparametric tests for efficacy in a two-group setup has been widely studied in the literature, including [5], who evaluated blood pressure changes in antihypertensive drug trials, while [53] compared pain scores and skin improvement in crossover and dermatology trials, respectively. Other studies assessed the accuracy of diagnostic tools [82], cholesterol level changes in nutritional research [65], genetic and biomarker disclosure process [92] and progression-free survival in oncology [53].

The Kruskal-Wallis (KW) test, developed by Kruskal and Wallis in 1952, is a non-parametric alternative to one-way Analysis of Variance (ANOVA). It is used to compare medians across two or more independent groups, particularly when the assumption of normality or homogeneity of variances is violated. Unlike ANOVA, the Kruskal-Wallis test ranks the data instead of analyzing raw values, making it robust against outliers and non-normal distributions. Consider k independent groups with sample sizes n 1 , n 2 , … , n k and combined total observations N = ∑ i = 1 k n i . Let R i j denote the rank of the j-th observation in group i, where all observations are ranked jointly across groups. The test statistic H is calculated as H = 12 ( N ( N + 1 ) ) − 1 ∑ i = 1 k ( T i 2 / n i ) − 3 ( N + 1 ), where T i = ∑ j = 1 n i R i j is the sum of ranks for the i-th group. Under the null hypothesis ( H 0 ), which assumes that all groups are sampled from the same distribution, H approximately follows a chi-squared distribution with k − 1 degrees of freedom for large N. When H 0 is rejected, it indicates significant differences among the group medians. While the test detects differences among groups, it does not identify which specific groups differ. Post-hoc pairwise comparisons must be conducted to pinpoint group differences, often using adjusted rank-based tests or Bonferroni corrections.

The Mantel-Haenszel (MH) test, introduced by Mantel and Haenszel in 1959, is a non-parametric statistical method used to evaluate associations between two categorical variables while controlling for a confounding variable. It is widely used in clinical trials and epidemiology to analyze stratified data, particularly in cases where the data are organized into contingency tables across strata. Given K strata, each represented by a 2 × 2 table, the MH test combines information across all strata to compute a pooled odds ratio and test for homogeneity or association. For each stratum k, the contingency table is represented as:

The Log-Rank test, first introduced by Mantel in 1966, is a nonparametric statistical method used to compare survival distributions between two or more groups. The test is based on the null hypothesis ( H 0 ) that there is no difference in the survival distributions between groups. Consider k groups, and let n i ( t ) and d i ( t ) denote the number of subjects at risk and the number of events observed at time t in group i, respectively. The observed ( O i ) and expected ( E i ) event counts for each group are calculated as follows: E i ( t ) = ( n i ( t ) / n ( t ) ) · d ( t ), where n ( t ) = ∑ i = 1 k n i ( t ) is the total number of subjects at risk at time t, and d ( t ) = ∑ i = 1 k d i ( t ) is the total number of events at time t. The observed and expected values are then summed across all time points. The test statistic is given by χ 2 = ( ∑ i = 1 k ( O i − E i ) ) 2 / ∑ i = 1 k V i , where V i is the variance of O i under the null hypothesis, calculated as V i = ( ∑ t n i ( t ) · ( n ( t ) − n i ( t ) ) · d ( t ) · ( n ( t ) − d ( t ) ) ) / ( n ( t ) 2 · ( n ( t ) − 1 ) ). Under H 0 , the test statistic χ 2 follows a chi-squared distribution with k − 1 degrees of freedom. For k = 2, the test reduces to a single degree of freedom test comparing two groups. The Log-Rank is a robust, distribution-free method for comparing survival curves and it retains good power when the hazard ratios between groups are proportional. However, violations of the proportional hazards assumption can lead to biased results. The test is also sensitive to censoring patterns; unbalanced censoring across groups can distort the results. Moreover, it does not adjust for covariates, requiring the use of stratified methods for more complex analyses.

Gray’s test, introduced in [55], is a non-parametric method for comparing cumulative incidence functions in the presence of competing risks. In clinical trials and epidemiological studies, competing risks occur when an individual is at risk of experiencing more than one mutually exclusive event, such as death from different causes. Unlike the Log-Rank test, which assumes all events are of the same type, Gray’s test accounts for the sub-distribution of competing risks. The cumulative incidence function (CIF), denoted as F j ( t ), represents the probability of experiencing event j by time t, considering the presence of other competing events. Gray’s test assesses whether the CIFs differ significantly between groups for a specific event type j. The test statistic for Gray’s test is derived from weighted differences in observed and expected event counts for each group. For k groups, let d i j ( t ) denote the number of events of type j at time t in group i, and n i ( t ) denote the number of individuals at risk in group i at time t. The weighted observed-minus-expected difference is computed as S ( t ) = ∑ i = 1 k w ( t ) · d i j ( t ) − E i j ( t ) , where E i j ( t ) is the expected number of events in group i under the null hypothesis, and w ( t ) is a weight function. The variance of S ( t ) is estimated to compute the test statistic χ 2 = S ( t ) 2 / Var [ S ( t ) ]. Under the null hypothesis, χ 2 follows a chi-squared distribution with k − 1 degrees of freedom. This test is found to be particularly valuable in clinical trials where competing risks are prominent, such as cancer studies with multiple causes of mortality. However, Gray’s test is sensitive to the choice of weights, and assumes independence between competing events and censoring, which may not always hold in practice. Additionally, the test does not adjust for covariates, requiring extensions such as Fine and Gray’s regression model for more complex analyses.

The Friedman Test, introduced by [45], is a nonparametric method designed for analyzing repeated measures or matched group data across multiple treatments. The test ranks observations within each subject and evaluates differences in ranks across treatments, making it robust to non-normal distributions. For n subjects and k treatments, let R i j denote the rank of the j-th treatment for the i-th subject. The test statistic is computed as: Q = 12 n k ( k + 1 ) ∑ j = 1 k R ¯ j − k + 1 2 2 , where R ¯ j = 1 n ∑ i = 1 n R i j is the average rank of treatment j. Under the null hypothesis that all treatments have identical distributions, Q approximately follows a chi-squared distribution with k − 1 degrees of freedom when n is large. A significant result suggests differences in treatment effects. The Friedman Test is particularly suitable for ordinal data or non-normally distributed outcomes, where traditional parametric methods like repeated measures ANOVA are inappropriate.

The Friedman Test offers simplicity and robustness as its primary advantages. It does not assume normality and is straightforward to implement, making it accessible for small datasets and early-stage studies. However, it comes with notable limitations. The test assumes balanced data, which means every subject must have observations across all time points or treatments—a condition rarely met in real-world longitudinal studies with missing data. Additionally, it treats measurements as independent ranks within subjects, ignoring temporal trends or within-subject variability, which are critical in longitudinal studies. Despite these limitations, the Friedman Test remains a foundational method for small-scale repeated measures experiments, especially in fields like nutrition and behavioral studies.

The Brunner-Munzel Test, introduced by [15], is a rank-based nonparametric method designed to compare two groups while allowing for unequal variances and non-normal distributions. Unlike the Wilcoxon Rank-Sum Test, the Brunner-Munzel Test does not assume homogeneity of variances, making it robust in situations where group variances differ. For independent samples X 1 , … , X n from group 1 and Y 1 , … , Y m from group 2, the test statistic evaluates the probability P ( X < Y ) + 0.5 P ( X = Y ), interpreted as the stochastic dominance of one group over the other. The test statistic T is calculated as: T = Δ ˆ σ ˆ 2 , where Δ ˆ is the difference in rank-based means between the two groups, and σ ˆ 2 estimates the variance of the rank means. Under the null hypothesis of no difference between the groups, T follows a standard normal distribution asymptotically. The Brunner-Munzel Test is particularly useful in comparing treatment effects where variability differs substantially between groups, such as in clinical trials with heterogeneous populations.

The Brunner-Munzel Test offers significant advantages over traditional rank-based methods. Its ability to handle unequal variances makes it particularly robust in practical scenarios, such as comparing treatment efficacy in diverse patient populations or analyzing data with extreme outliers. Additionally, it retains the benefits of nonparametric methods, including robustness to non-normality and suitability for ordinal outcomes. However, the test is limited to two-group comparisons, making it less versatile for studies involving more than two treatments or time points. Moreover, like many rank-based tests, it does not explicitly incorporate temporal dependencies or repeated measures, which restricts its applicability in longitudinal settings. Despite these limitations, the Brunner-Munzel Test remains a valuable tool in scenarios requiring robust, nonparametric two-sample comparisons.

The nparLD framework, introduced by [98] is a nonparametric rank-based approach specifically designed for the analysis of longitudinal and repeated measures data. It extends traditional rank-based methods by using pseudo-rank transformations, which preserve the ordinal structure of the data while accounting for repeated measures. For a dataset with n subjects, t time points, and k treatment groups, let Y i j represent the observation for subject i at time j. The pseudo-rank R i j ∗ of Y i j is calculated based on its rank relative to all observations across time points and groups. Using these pseudo-ranks, the framework evaluates treatment, time, and interaction effects through F-type test statistics. The F-type statistic for a given effect is expressed as: F = Between-group pseudo-rank variability/Residual pseudo-rank variability, where the numerator captures the variability explained by the effect (e.g., group, time, or interaction), and the denominator accounts for residual variability. For example, the test for a group effect compares the average pseudo-ranks across treatment groups, while controlling for time and subject variability. The null hypothesis of no effect is tested using permutation or large-sample asymptotic methods, with the F-statistic approximated by an F-distribution under the null.

The nparLD framework offers several advantages. Its ability to handle unbalanced data and missing observations makes it particularly suitable for real-world longitudinal studies, where dropout and irregular follow-up times are common. Additionally, it accommodates both ordinal and continuous outcomes, providing robust results even in the presence of non-normal data or outliers. The inclusion of interaction effects, such as group-by-time interactions, further enhances its utility for complex study designs. However, there are notable limitations. The need for multiplicity adjustments when testing multiple hypotheses can reduce statistical power, and the computational demands of pseudo-rank transformations and F-statistic calculations increase with larger datasets. While nparLD is robust for moderate-to-large sample sizes, its power may be limited in small-sample settings. Despite these challenges, nparLD has been widely adopted in clinical trials, with applications ranging from the evaluation of repeated measures of biomarkers to studies of vaccine efficacy and weight-loss interventions.

The Longitudinal Rank-Sum Test (LRST), introduced by [158], is a nonparametric method developed to evaluate treatment effects across multiple longitudinal endpoints in clinical trials. Unlike traditional rank-based methods, the LRST accounts for the complexity of repeated measures and multi-endpoint designs while maintaining robustness against non-normal distributions and outliers. Let x i t k and y j t k represent the changes from baseline for subjects in the control and treatment groups, respectively, where i and j index subjects, t denotes time, and k represents outcomes. Observations are ranked across all subjects, time points, and outcomes, with ranks R i t k and R j t k . The test statistic is calculated as T LRST = ( R ¯ y · · · − R ¯ x · · · ) / Var ˆ ( R ¯ y · · · − R ¯ x · · · ) , where R ¯ y · · · and R ¯ x · · · are the mean ranks for the treatment and control groups, respectively, aggregated across all time points and outcomes. Under the null hypothesis of no treatment effect, T LRST asymptotically follows a standard normal distribution. LRST has also been developed for multi-arm clinical trials by [49, 51].

The LRST offers several advantages, making it a valuable tool in modern clinical trials. Its rank-based approach is robust to outliers and non-normal data distributions, which are common in real-world datasets. By simultaneously evaluating multiple longitudinal endpoints, the LRST reduces the need for multiplicity adjustments, enhancing statistical power while maintaining control of Type I error rates. Additionally, it handles missing data and irregular follow-ups efficiently, ensuring flexibility in complex trial designs. However, the test relies on large-sample approximations for its validity, and its performance in small-sample studies may require further evaluation. Despite these limitations, the LRST has shown strong applicability in neurodegenerative disease and oncology trials, where multiple outcomes such as motor function and cognitive performance are monitored over time, providing a comprehensive evaluation of treatment efficacy.