VARIATION PARTITIONING: PCA, RDA AND PRDA

To tackle collinearity issues, many researchers have used principal components analysis (PCA) to reduce a number of correlated variables into a set of uncorrelated variables, which reserves its total variance and uncovers its hidden patterns (Varanka et al., 2012). In addition, in order to realize the relative importance of different explanatory variables and their interactive effects, direct gradient analysis such as redundancy analysis (RDA) and its successive partial constrained ordinations, i.e., partial RDA (pRDA) have been commonly proposed (ter Braak, 1998; Borcard et al., 2014). This allows researchers to explore the relationships between predictor variables and dependent variables by removing the intertwined effects among them (Liu, 1997). However, most previous studies utilized PCA or RDA methods focusing on the relationship between biological phenomena and environmental influence, and there are only a few studies on water quality (e.g., Nava-López et al., 2016). Therefore, exploring the likely collinear controlling factors and their internal relationships to riverine DIN transport based on PCA and RDA will be valuable.

In order to interpret the spatial and temporal patterns of riverine NO3− and NH4+

export and unravel the dependency among controlling factors, i.e., the human disturbance, climatic factors and landscape settings (Table 1), the whole analysis was carried out in three steps (Liu, 1997): (i) the PCA was applied to find out a set of uncorrelated variables, (ii) detailed relations between export and each one of the controlling factors were displayed using a scatterplot matrix, and (iii) the RDA and pRDA were conducted to disentangle the contribution of the major variables.

Table 1. Definition of different variables used in the three dimensions. Abbre. =

Dimension Variables Abbre. Definition

Climatic Factors

Rainfall (mm)

RDry Rainfall in dry season of the year RWet Rainfall in wet season of the year Streamflow

(mm)

SFDry Discharge rate in dry season of the year SFWet Discharge rate in wet season of the year Temperature

(°C) T The degree of hotness or coldness of environment

Landscape Settings

Channel

length (km) CL Total length of the stream channel Longest

channel length (km)

LCL The length of the longest stream channel in watershed

Relief Rel The difference between the highest and lowest elevations in watershed

Area (km

) A Drainage area of watershed

Slope (%)

SLP100 The average slope in the 100 m buffer zone SLP200 The average slope in the 200 m buffer zone SLP500 The average slope in the 500 m buffer zone SLP1000 The average slope in the 1000 m buffer zone SLP2000 The average slope in the 2000 m buffer zone

SLP The average slope in watershed Drainage

density (km

)

DD Total channel length over drainage area

L/G (m) The ratio of median flow path length to median flow path gradient

PD100 Population density in the 100 m buffer zone PD200 Population density in the 200 m buffer zone PD500 Population density in the 500 m buffer zone PD1000 Population density in the 1000 m buffer zone PD2000 Population density in the 2000 m buffer zone

PD Population density in watershed

Buildup (%)

BD100 The percentage of buildup area in the 100 m buffer zone BD200 The percentage of buildup area in the 200 m buffer zone BD500 The percentage of buildup area in the 500 m buffer zone BD1000 The percentage of buildup area in the 1000 m buffer zone BD2000 The percentage of buildup area in the 2000 m buffer zone

BD The percentage of the buildup area in watershed

Agriculture (%)

AGR100 The percentage of agriculture in the 100 m buffer zone AGR200 The percentage of agriculture in the 200 m buffer zone AGR500 The percentage of agriculture in the 500 m buffer zone AGR1000 The percentage of agriculture in the 1000 m buffer zone AGR2000 The percentage of agriculture in the 2000 m buffer zone

AGR The percentage of the agriculture in the watershed

PCA was applied to reduce redundant information and to transform the original correlated data into another set of uncorrelated variables. The PCA keeps only a few independent sets (patterns) of environmental data that are distinct from each other, which

will help to realize the effects of various characteristics of watersheds on NO3− and NH4+

exports in our study (Liu, 1997; Johnson et al., 2007). The varimax rotation was selected to better separate divergent groups of variables, as suggested (Jolliffe and Cadima, 2016).

The environmental variables were centered and standardized in order to approximate normally distributed random errors and then were derived from the PCs via a standardized linear projection which maximizes the variance in the projected space (Hotelling, 1933).

For a set of observed-dimensional data vectors, { the q principal axes { could be derived as the orthonormal axes onto which the retained variance under projection is maximal. It can be shown that the vectors are given by the q dominant eigenvectors (i.e., those with the largest associated eigenvalues ) of the sample covariance matrix. The outcomes of PCA help us to identify relationships between these variables and determine which variables require further investigation. The variables with loading higher than 0.1 in the first and second PCs were kept for the following RDA and pRDA analysis to constrain the ordination of environmental variables and to avoid the collinearity problem (Sutter and Kalivas, 1993; Johnson et al., 2007).

Moreover, we know water quality is regulated by riparian zones along the river and stream networks, but what needs to be clarified is spatially to what extent their individual effect is (Uriarte et al., 2011). Here, we delineated the buffer zones of 100, 200, 500, 1000 and 2000 m along the stream network using the buffer tool in ArcGIS v.10.7. (ESRI Inc., Redlands, CA, USA) The environmental variables within the entire watershed and five buffer zones were also retrieved as previous studies suggested (Nielsen et al., 2012;

Nava-López et al., 2016; Xiao et al., 2016). The land cover/land use data were acquired from the Ministry of the Interior of Taiwan (Figure 1b; https://www.moi.gov.tw), and the digital elevation model (DEM) data were derived from the open data platform in Taiwan

(https://data.gov.tw/), which were provided as input for calculations of landscape settings and human disturbance variables (Table 1).

The RDA and pRDA were further applied to quantify the individual effect and integrative contribution among human disturbance, climatic factors and landscape setting on riverine NO3− and NH4+

exports (ter Braak 1988). RDA extends the algorithm of PCA with a response matrix Y (with n objects and p variables) by an explanatory matrix X (with n objects and m variables). First, RDA produces a matrix of fitted values Yˆ through Equation (2),

  X X X Y X

Y ˆ  



and second, runs a PCA based on Yˆ (legendre and legendre, 2012).

For pRDA, the additional explanatory variables, called covariables, are assembled in matrix W; the linear effects of the explanatory variables in X on the response variables in Y are adjusted for the effects of the covariables in W (Legendre et al., 2011). In our study, the total variance of riverine NO3− and NH4+ exports could be explained by the variables derived from human disturbance, climatic factors and landscape setting, and their individual contribution of NO3− and NH4+ export can be finally figured out. We further partitioned the total variation of the riverine NO3− and NH4+ response variables using three steps (Table 2). First, canonical ordination with no covariables was used to estimate the total amount of variance explained (as sum of canonical eigenvalues) in the NO3− and NH4+

export attributed to all explanatory variables, human disturbances (H), climatic factors (C) and landscape setting (L), and the total unexplained variance (1HCL). Second, the combinations of various covariables were considered to calculate the separate effect of each variable (H, C or L), in which an individual predictor variable was run (e.g., H) with the remaining other two as covariables (e.g., C&L). Third, a series of partial canonical ordinations were used to calculate the unique and interactive effects for each set of

predictors (e.g., C&LH) by considering the interaction term of interest as explanatory (C&L) and excluding the effect of not interest (e.g., H). For more details of calculations, please refer to Borcard et al. (1992) and Liu (1997).

Table 2. Eigenvalues of partial RDA (pRDA) of NO

and separate climatic (C), landscape setting (L) and human disturbance (H) and interactive effects among C, L and H.

Environmental factor Covariable NO

NH

Wet Dry Wet Dry

Unexplained variable 0.27 0.14 0.31 0.21

CLH None 0.73 0.86 0.69 0.79

C L&H 0.31 0.27 0.03 0.02

L&H C 0.44 0.27 0.68 0.77

L C&H 0.00 0.00 0.05 0.00

C&H L 0.41 0.74 0.06 0.18

H C&L 0.07 0.09 0.02 0.13

C&L H 0.31 0.2 0.08 0.02