Changes in streamflow extremes can affect the economy and are likely to impact the most vulnerable in society. Estimating these changes is crucial to develop rational adaptation strategies and to protect society. Streamflow data from 1106 gauges were used to provide a comprehensive analysis of change in eight different extreme indices. The modified trend-free prewhitening and the false discovery rate were used to account for serial correlation and multiplicity in regional analysis, issues shown here to distort the results if not properly addressed. The estimated proportion of gauges with significant trends in low and high flows was about 23% and 15%, respectively. Half of these significant gauges had more than 60 years of data and were associated with changes greater than 5% per decade. A clear spatial pattern was identified, where most increasing trends in both low and high flows were observed in Southern Brazil, and decreasing trends in the remaining regions, except for the Amazon, where a pattern is not clear, and the proportion of significant gauges is low. Results based only on gauges unaffected by reservoirs suggest that reservoirs alone do not explain the increasing trends of low flows in the southern regions nor the decreasing trends in high flows in the remaining hydrographic regions.

This paper develops the quasi‐analytic Bayesian analysis of the generalized least squares (GLS) (B‐GLS) model introduced by Reis et al. (2005, https://doi.org/10.1029/2004WR003445) into an operational and statistically comprehensive GLS regional hydrologic regression methodology to estimate flood quantiles, regional shape parameters, low flows, and other statistics with spatially correlated flow records. New GLS regression diagnostic statistics include a Bayesian plausibility value, pseudo adjusted R², pseudo analysis of variance table, and two diagnostic error variance ratios. Traditional leverage and influence are extended to identify rogue observations, address lack of fit, and support gauge network design and region‐of‐influence regression. Formulas are derived for the Bayesian computation of estimators, standard errors, and diagnostic statistics. The use of B‐GLS and the new diagnostic statistics are illustrated with a regional log‐space skew regression analysis for the Piedmont region in the Southeastern U.S. A comparison of ordinary, weighted, and GLS analyses documents the advantages of the Bayesian estimator over the method‐of‐moment estimator of model error variance introduced by Stedinger and Tasker (1985, https://doi.org/10.1029/WR021i009p01421). Of the three types of analyses, only GLS considers the covariance among concurrent flows. The example demonstrates that GLS regional skewness models can be highly accurate when correctly analyzed: The B‐GLS average variance of prediction is 0.090 for South Carolina (92 stations), whereas a traditional ordinary least squares analysis published by the U.S. Geological Survey yielded 0.193 (Feaster & Tasker, 2002, https://doi.org/10.3133/wri024140). B‐GLS provides a statistical valid framework for the rigorous analysis of spatially correlated hydrologic data, allowing for the estimation of parameters and their actual precision and computation of several diagnostic statistics, as well as correctly attributing variability to the three key sources: time sampling error, model error, and signal.