Statistical Analyses
Percent-cover data were arcsine-square root transformed prior to analysis. In all tests, we examined plots of residuals for assumptions of normality, homogeneity of variance, and independence of error terms. In our experimental subordinate predator removals, Leptasterias spp. predators were very low in all treatment plots (zero in most plots). We used ANOVAs and Tukey’s HSD post-hoc tests to see if a) removals of only Leptasterias spp. (-L+W) were different from the no-removal treatment (+L+W), and b) removals of only whelks (-W+L) were different from removals of both predators (-W-L). In both cases there were no differences, so we combined treatments (at the replicate level) to make a +Whelks (i.e. +W+L and +W-L) and a -Whelks (-W+L and -W-L) comparison.

To test the effects of whelks on prey species, we used linear mixed effects models to analyze percent cover of prey species separately using the ‘nlme’ package in R (Pinheiro et al. 2014). Models included a random effect for plot to account for repeated measurements over time. The final model for each prey species was determined by removing terms from a full model based on the Akaike Information Criterion corrected for small sample size (AICc) from in the package ‘AICcmodavg’ (Mazerolle 2015). We used AIC to determine if adding a weighted variance structure was necessary in cases where plots of residuals revealed heterogeneity of variance. Time was treated as a categorical variable because of non-linear relationships with percent cover. Because plots started with different prey communities, we analyzed both percent cover and change in percent cover from initial cover with qualitatively similar results. We made post-hoc comparisons using the ‘phia’ package in R (De Rosario-Martinez 2015) with Holm’s-corrected p-values.
To determine whether population abundances of subordinate predators had changed from 2014 to 2015, we used negative binomial regression for survey count data, which were overdispersed, in the ‘MASS’ package (Venables & Ripley 2002) in R. In this case, a full model included a site x year interaction, which was dropped if non-significant. Because 2014 had many survey time points and 2015 had fewer surveys, we tried analyzing data from all time points and also from only the July 2014 and 2015 time points. We used all survey data as results were similar in both cases. To assess spatial changes in subordinate predator distributions, we ordered quadrats by measured tide height and analyzed differences in the cumulative distributions of whelk counts between the two years using two-sample Kolmogorov-Smirnov Tests with p-values bootstrapped using the package ‘Matching’ (Sekhon 2011). For these tests, we used only matched July time points in order to minimize distributional changes due to season.

BCO-DMO Processing:
- added conventional header with dataset name, PI name, version date
- renamed parameters to BCO-DMO standard
- reformatted date from m/d/yyyy to yyyymmdd
- sorted data by site, date, treatment