# Scalability of Water Property Measurements in Space and Time on a Brackish Archipelago Coast
## 2 Materials and Methods
### Statistical Analyses
The more complex intercept-only models (Equation (2)) included nested random intercepts for the site and sub-site structures:  
\[\text{Response}\sim(1\,|\,\text{Site}/\text{Sub-site}) \tag{2}\]  
A simpler variant of the intercept-only multilevel model included a random intercept for site (Equation (3)). Thus, the intercept-only models allowed the overall level of the water property variable to vary between spatial locations. \[\text{Response}\sim(1\,|\,\text{Site})\text{ (or Y}_{\text{ij}}=\gamma_{00}+\text{u}_{ \text{ij}}+\text{e}_{\text{ij}}) \tag{3}\]  
where Yij = the ith observation in the jth group, \(\gamma_{00}\) = the overall mean, u0 = group-level residuals, and eij = observation-level residuals (residual error). We calculated intraclass correlation (\(\rho\); Equation (4)) for these intercept-only models, to quantify the proportion of the variance explained by the spatial grouping structure in the observations. Thus, intraclass correlation examined how much of the variation in one water property variable depended on location. \[\rho=\sigma^{2}\text{u}_{0}/(\sigma^{2}\text{u}_{0}+\sigma^{2}\text{e}) \tag{4}\]  
where \(\sigma^{2}\text{u}_{0}\) = variance of the group-level residuals u0, and \(\sigma^{2}\text{e}\) = variance of the observation-level residuals eij [33]. Since no variance was found at the sub-site level, the following procedures are described for the simpler site-level models only. Fourth, we examined the interaction of temporal and vertical gradients of the four water quality variables with spatial location. In other words, we examined how the temporal and vertical gradients vary between sites and sub-sites. We partitioned the variation in water property variables to temporal (week), vertical (depth) and horizontal (site) scales by fitting full multilevel models. The multilevel models were fitted using the restricted maximum likelihood method (REML). These models combined observation-level terms of explanatory variables with random intercepts (Equation (2) or (3)) and random slopes for \"week\" and \"depth\". The multilevel models were fitted sequentially (following [33]) to examine the contribution of each term to the BIC and residual variance. The sequential fitting started from random intercept models, including observation-level terms and nested random intercepts for \"site\" and \"sub-site\":  
\[\text{Response}\sim\text{Week}+\text{Depth}+\text{Week}^{2}+\text{Depth}^{2} \ldots+\text{Week}^{10}+\text{Depth}^{10}+\text{Week:Depth}+(1\,|\,\text{Site}/ \text{Sub-site}) \tag{5}\]  
Finally, random slope models were fitted, including site-level random terms for the first- and second-order terms of the explanatory variables:  
\[\text{Response}\sim\text{Week}+\text{Depth}+\text{Week}^{2}+\text{Depth}^{2} \ldots+\text{Week}^{10}+\text{Depth}^{10}+\text{Week:Depth}+(1+\text{Week}+ \text{Week}^{2}+\text{Depth}+\text{Depth}^{2}|\,\text{Site}) \tag{6}\]The random structure of the multilevel model thus allowed both the intercept and slope coefficients of all fixed terms to vary between sites. Stepwise regression was again applied to identify an \"optimal\" model structure, with the lowest BIC (Appendix A). All statistical analyses and modelling were performed with the R statistical software (version 4.0.3; [35]). Linear mixed models were fitted with the package lme4 (1.1-26; [36]) and examined with package lmerTest (3.1-3; [37]). The effects of the explanatory variables were calculated with package effects (4.2-0; [34]). The main and interactive effects of \"week\" and \"depth\" were calculated with the effect function by absorbing the lower-order terms marginal to the higher-order term in question and averaging over other model terms [38]. As a result, a table of predicted values at different levels of the two interacting factors was calculated. The package ggeffects (1.1.0; [39]) was utilised to calculate random effects for visualisation. Maps were created in QGIS (3.16; [40]) and all visualisations were finalised in R with packages ggplot2 (3.3.3; [41]) and gridExtra (2.3; [42]).