Climate and hydrologic variables such as temperature, precipitation, streamflow, and baseflow generally do not follow Gaussian distribution due to the presence of outliers and heavy tails. Therefore, they are usually analyzed using the nonparametric Wilcoxon rank sum test rather than parametric methods like classical t tests and analysis of variance. Furthermore, in addition to having a non-Gaussian distribution, these data exhibit monthly/seasonal variability, which leads to within month/season cluster correlation. In this study, a nonparametric procedure, called joint rank fit (JRFit), for analyzing cluster-correlated data was implemented and compared against traditional methods such as restricted maximum likelihood, least absolute deviations, and rank-based fit (a model-based extension of Wilcoxon rank sum) for studying the coupled effect of the phases of El Niño–Southern Oscillation and Atlantic Multidecadal Oscillation on baseflow levels. The results from a large Monte Carlo simulation experiment showed that JRFit was more efficient than the other three methods for data with (i) high variability, (ii) outliers due to contamination, or (iii) strong monthly/seasonal correlation. The efficiency gain of JRFit was up to 50% compared to restricted maximum likelihood for heavy-tailed and highly correlated data. Predictive performance evaluated using the mean absolute prediction error and mean prediction standard error from an out-of-sample cross-validation study showed JRFit to be optimal for providing predictions of baseflow on the basis of the phases of El Niño–Southern Oscillation and Atlantic Multidecadal Oscillation. Thus, it is recommended that JRFit be implemented in hydroclimatic studies to provide powerful inference when there is evidence of clustering in the data. ©2018. American Geophysical Union. All Rights Reserved.