Refining predictions of metacommunity dynamics by modeling species non‐independence

Predicting the dynamics of biotic communities is difficult because species’ environmental responses are not independent, but covary due to shared or contrasting ecological strategies and the influence of species interactions. We used latent-variable joint species distribution models to analyze paired historical and contemporary inventories of 585 vascular plant species on 471 islands in the southwest Finnish archipelago. Larger, more heterogeneous islands were characterized by higher colonization rates and lower extinction rates. Ecological and taxonomical species groups explained small but detectable proportions of variance in species’ environmental responses. To assess the potential influence of species interactions on community dynamics, we estimated species associations as species-to-species residual correlations for historical occurrences, for colorizations, and for extinctions. Historical species associations could to some extent predict joint colonization patterns, but the overall estimated influence of species associations on community dynamics was weak. These results illustrate the benefits of considering metacommunity dynamics within a joint framework, but also suggest that any influence of species interactions on community dynamics may be hard to detect from observational data.


INTRODUCTION
Predicting the dynamics of biotic communities over time is a fundamental challenge of ecology, and of undeniable importance in a time of rapid environmental change. Demographic studies of single species allow projection of population dynamics through time (e.g., Caswell 2001, Merow et al. 2014), but are not directly informative about community dynamics. An important reason for this is that species are not independent units responding to external environmental drivers, but are dependent on the dynamics of the community to which they belong. Consequently, explicit consideration of community-level patterns and processes may improve our ability to predict how communities change.
As a simple example, ecologically similar species may respond similarly to environmental drivers, leading to synchronized (meta-)population dynamics. This view is central to the research program centered on species' functional traits (McGill et al. 2006), where "response traits" are those attributes of species that are predictive of environmental responses Garnier 2002, Suding et al. 2008). Ecologically similar species may also tend to be closely related, so that functional traits exhibit a strong phylogenetic signal (e.g., Webb et al. 2002, Cavender-Bares et al. 2009; but see Losos 2008). An advantage of leveraging phylogenetic information as a complement to trait data is that species relatedness may convey information about the net phenotypic similarities among species (aspects of ecological similarity not captured by individual traits), or traits that are ecologically relevant but difficult to measure.
Changes in species distributions and abundances may depend on or modify species interactions, with unknown consequences for community dynamics (Holt 1997, Blois et al. 2013, Wisz et al. 2013, Alexander et al. 2015, Magurran 2016, Schleuning et al. 2016, Urban et al. 2016. Interactions can occur both among trophic levels (e.g., pollination, herbivory), and within trophic levels (e.g., plant-plant interactions; Table 1). As an example of the latter, the colonization of a site by one species may lead to the extinction of other species through competition, as when woody shrub expansion suppresses populations of understory plants (Pajunen et al. 2011, Mod andLuoto 2016). Similarly, colonization by foundation species may facilitate colonization by other species, as observed for cushion plants in harsh alpine environments (K€ orner 2003). Therefore, joint analysis of the spatial distribution and temporal dynamics of entire communities may lead to a more predictive understanding of biotic community responses to environmental change.
We can think about the processes discussed above as aspects of species non-independence, either in the statistical sense of species falling into groups that respond similarly to environmental variation, or in the ecological sense of species directly interacting with each other. Our overall hypothesis is that explicit consideration of such patterns and processes in community-wide analyses can refine predictions of community dynamics. Recent developments in the field of joint species distribution modeling now allows estimating shared responses to environmental drivers among predefined ecological, phenotypic, and phylogenetic species groups . Furthermore, these models yield estimates of residual correlations among species after accounting for shared environmental responses (Pollock et al. 2014, Warton et al. 2015. Species associations quantified by residual correlations may represent shared responses to unmeasured aspects of the environment, but also the potential influence of species interactions. By considering entire communities within a single model, we can therefore start to assess quantitatively the potential influence of species interactions on community structure and dynamics. Here, we analyze changes in the composition of the vascular plant communities of 471 islands in the southwest Finnish archipelago over 80 yr. Because islands represent distinct sampling units, insular metacommunities provide unique systems for testing biological hypotheses related to community assembly, diversity, and turnover (e.g., MacArthur and Wilson 1967, Simberloff 1974, Diamond 1975, Brown and Kodric-Brown 1977, Kohn and Walsh 1994, Cody 2006, Hanski 2016. In the study area, previous analyses of individual species have identified a set of environmental drivers of species occurrence, extinction, and colonization (e.g., von Numers and Korvenpaa 2007, Hannus and von Numers 2010, von Numers 2011, 2015. Here, we extend the scope of analysis by jointly modeling extinction and colonization rates of all 585 vascular plant species that occur in the study area, allowing us to assess patterns of shared dynamics among species. To assess whether quantifiable axes of species similarity are informative about changes in communities over time, we asked whether and to what extent changes in the metacommunity are structured by (1) species ecological groups (shoreline specialists vs. ordinary species) and (2) taxonomic relatedness. To further assess the potential role of species interactions in structuring community dynamics, we also asked (3) whether statistical species associations inferred from historical snapshot data are informative about species' shared dynamics over time, and (4) whether these associations are taxonomically structured.

Study area and plant community data
The southwest Finnish archipelago comprises at least 22,000 islands ranging from small skerries to large,  Ecology,Vol. 101,No. 8 forested islands (Gran€ o et al. 1999). The climate is maritime, but the early part of the growing season often has continental characteristics, with low precipitation and higher insolation compared to the mainland. Further description of the study area is provided in von Numers and van der Maarel (1998) andvon Numers (2015). The species occurrence data consist of paired historical and recent inventories of all vascular plants on 471 islands (for details, see von Numers 2011Numers , 2017. The historical inventories were conducted between 1925 and 1946 (mainly in the 1930s) by Eklund (1958), and in the 1940s by Skult (1960), and the recent inventories were conducted between 1996 and 2017 by M. von Numers. In both inventories, species lists were compiled for each island, but no data on abundance were collected. In the recent inventory, the islands were surveyed during the same part of the growing season as in the historical inventories, usually between 10 June and 30 July. The islands were surveyed systematically using a species checklist. Some islands were surveyed twice, and some of the largest islands demanded several days of work. Both Eklund and Skult similarly compiled the species lists in a systematic way using checklists. Eklund was known to work systematically and scrupulously in the field (H. Skult, personal communication), devoting his life to surveying the flora of the archipelago of southwest Finland. Skult worked in a similar systematic way. While both the historical and recent inventories were conducted by highly experienced field botanists, it is possible that a few species remained undetected (see von Numers and Korvenpaa 2007 for examples and discussion). Importantly, we except the detectability of individual species to have been similar between the two inventories.
A total of 585 species were recorded on the study islands, with a median species richness of 113 species per island (range = 22-279 species). Species´life histories range from short-lived annuals (e.g., Atriplex prostrata, Moehringia trinervia, Saxifraga tridactylites, Linum catharticum) and biennials (e.g., Isatis tinctoria, Peucedanum palustre, Cirsium palustre) to long-lived perennials (e.g., Vincetoxicum hirundinaria, Juniperus communis, Pinus sylvestris). Thus, while most species will have gone through several to many generations between the inventories, a few individuals of the more long-lived species may have persisted on a given island for the entire duration of the study. We defined shoreline species, i.e., species of the geolittoral zone (shore between the mean water level and the high water level), the hydrolittoral zone (between mean water level and low water level), and the emergent species of the sublittoral zone (always under water), as in von Numers (2011). We expected these species to exhibit shared patterns of change because they inhabit a shared and distinct habitat.
We assigned colonization events when a species was absent on an island in the historical inventory but present in the recent inventory, and extinction events when a species was present in the historical inventory but absent in the recent inventory. Thus, "colonizations" and "extinctions" can include cases of repeated change in presence/absence status, if a species, say, colonised, went extinct, and then recolonised during the approximately 80 yr between the two inventories. Using these definitions, 15,008 colonization events and 9,384 extinction events were recorded across all islands between the two inventories.

Environmental covariates
Island diversity patterns are known to depend on island area, isolation, and environmental heterogeneity (MacArthur and Wilson 1967, Brown and Kodric-Brown 1977, Kohn and Walsh 1994. To account for these underlying gradients in diversity, we used ArcGIS 10.1 to measure a set of geographical covariates describing the position and structure of each island, namely island area (natural logarithm of island area in m 2 ), topographic complexity (standard deviation of elevation points of 10 9 10 m grid cells of a digital elevation model covering each island), and neighborhood size (natural logarithm of total area of land within 5 km of the shoreline of each island in m 2 ). We assumed that the topography variable would capture variation in microclimate due, for example, to differences in solar radiation and drainage of moisture (Opedal et al. 2015). The neighborhood variable captures variation in the probability of colonization due to distance to source populations, as well as variation in exposure to wind, waves, and sea currents.

Joint species distribution models
To assess community-wide patterns of vegetation change between the two inventories, and to answer our specific research questions derived from the overall hypothesis of species non-independence in community dynamics, we fitted latent-variable joint species distribution models using the Hierarchical Modeling of Species Communities (HMSC) framework of Ovaskainen et al. (2017) implemented in the Bayesian framework in the Hmsc 3.0 R package (Tikhonov et al. 2020).
We fitted a series of models aiming to answer distinct research questions. In the first analysis (repeated-measures model), we treated the historical and recent inventories as repeated observations, and used probit regression to model the presence/absence of species (Y matrix in HMSC) as a function of time (historical vs. recent) and the geographic covariates island area, topographic complexity, and neighborhood size (X matrix). This approach allowed us to estimate the net change in occurrence of each species between the two inventories, and to estimate and compare temporal and spatial sources of variation in species occurrences and distributions. To assess whether ecologically similar species exhibit similar patterns of change in occurrence over time (research question i), we included in the model trait matrix (T) an indicator of species that are shoreline specialists, as a simple proxy of ecological similarity.
In the second analysis (colonization-extinction model), we explicitly modeled the colonization and extinction probabilities of each species. By modeling historical occurrences, colonizations and extinctions of all species jointly (in a single model), this approach allowed us to assess correlated probabilities of colonization and extinction across species. Each species occurred three times in the response matrix (Y): first as presence/absence in the historical data, and then as presence/absence for colonization and extinction, conditional on absence (for colonization) or presence (for extinction) in the historical data. Covariates were the same as in the repeated-measures model, except there was no time variable. We included response variable type (historical occurrence vs. colonization vs. extinction), species ecological type (shoreline specialist vs. ordinary species), and their interaction in the trait matrix (T). To assess whether species that were positively or negatively associated in the historical data also tended to colonize or go extinct together (research question iii), we correlated the estimated residual correlations for historical occurrences, colonization probabilities, and extinction probabilities for each species pair. We repeated this analysis for all species and for the shoreline species only.
In the third analysis (taxonomic-signal models), we asked whether related species responded similarly to covariates (research question ii), and whether positive and negative species associations were taxonomically structured (research question iv). We fitted separate models for the historical inventory, for colonization conditional on historical absence, and for extinction conditional on historical presence, excluding those species missing in each of these subsets of the data. Covariates and traits were the same as in the repeated-measures model, except there was no time variable. We modeled the matrix b describing species' responses to the environmental covariates as b $ N l; V qC þ 1 À q ð ÞI ½ ð Þ , where is the Kronecker product, C is a phylogenetic correlation matrix, and I is an identity matrix, and the mean l is modeled as a linear regression on species' traits . The strength of the phylogenetic signal is measured by the parameter q, with q = 0 indicating phylogenetic independence. In the absence of a dated phylogeny for the diverse plant communities in the study area, we assembled an undated phylogeny from Phylomatic (Webb and Donoghue 2005), based on the default Phylomatic R20120829 megatree for plants (program available online). 2 We set all branch lengths to one, and thus refer to the inferred signal as taxonomic rather than phylogenetic. To explore whether species associations were taxonomically structured, we correlated the taxonomic relatedness of each pair of species with their estimated residual cooccurrences.

Model fitting, evaluation, and cross-validation
We sampled the posterior distributions with two independent Monte Carlo Markov chains (MCMCs) for an increasing number of iterations until results stabilized across chains, and posterior trace plots revealed adequate chain mixing and convergence (the final number of MCMC iterations per chain ranged from 60,000 to 300,000; Appendix S1: Table S1). To evaluate chain mixing and convergence quantitatively, we computed effective sample sizes and potential scale-reduction factors for those model parameters used to make inferences (Gelman et al. 2014). We assessed explanatory power as the squared correlation coefficient (r 2 ) between the observed and predicted numbers of species, colonization events and extinction events per island. We also computed area-under-curve (AUC) statistics and coefficients of discrimination (Tjur r 2 ) for each species. The latter is defined as the difference in average model-predicted probability of occurrence, colonization, or extinction for islands where a species is present vs. absent (Tjur 2009). To assess predictive power, we cross-validated the repeated-measures and colonization-extinction models by refitting sequentially to one half of the data and obtaining predictions for those islands not used in model training (i.e., two-fold cross-validation).

Variance partitioning
We partitioned the variance explained by each model into contributions of each fixed and random effect based on the relative proportion of variances for each of the Gaussian components of the model's latent predictor . We further computed the trait r 2 , i.e., the percentage of the variance explained by fixed effects attributable to traits included in the model (see Abrego et al. [2017] for details).

Quantifying the influence of environment and species associations on predicted current occurrences
A species is present on an island in the recent inventory if it was historically present and did not go extinct, or if it colonized between the inventories. Hence, the probability of current occurrence of species i can be written as where p i is the probability of current occurrence, q i is the probability of historical occurrence, e i is the extinction probability, and c i is the colonization probability for species i. To assess how current occurrence probabilities (p i ) depend on environmental covariates vs. the occurrence and colonization-extinction dynamics of the complete metacommunity (Table 1), we performed a series of predictions using the colonization-extinction model where we sequentially excluded subsets of information, and computed current occurrence probabilities using the equation above. As a benchmark, we performed full-model conditional predictions including all available information, i.e., the covariates, the historical species occurrences, and the colonizations and extinctions of all non-focal species. Conditional predictions in HMSC use the model-inferred species associations to inform predictions, thus incorporating information about local community structure and dynamics. To assess the overall (potential) influence of historical occurrences, colonizations, and extinctions, we made predictions while holding q i , c i , and e i constant at their mean.
To assess the (realized) influence of each factor (covariates, historical occurrences, colonizations, and extinctions), we compared the full-model predictions to predictions obtained after excluding each focal factor. For example, we assessed the influence of historical occurrences of all non-focal species on the colonization probability of a focal species by comparing predictions from the full model vs. predictions made with all historical occurrences set to NA. We quantified the influence of each factor as the mean absolute deviation of the reduced-model predictions from the full-model predic- where p F ij is the fullmodel prediction for species i on island j, p R ij is the reduced-model prediction, and n is the number of islands.

Repeated-measures model: Patterns of net change in the metacommunity
The model treating the historical and recent inventories as repeated measures explained very well the number of species per island (r 2 = 1, cross-validation r 2 CV = 0.57), but less well the occurrences of individual species (mean coefficient of discrimination [Tjur r 2 ] = 0.30, range = 0-0.71, mean AUC = 0.94). Thus, the average predicted occurrence probabilities of species on islands where they were present were on average 0.30 greater than for islands where they were not present. Of the total variance explained by the model, time (historical vs. recent) explained on average 4.7%, the geographical covariates explained 31.7%, and the spatial latent factors representing unmeasured environmental variation as well as species interactions explained 63.6% ( Table 2). The limited species-specific explanatory power arising from the geographical covariates was also apparent as low species-specific predictive power as evaluated by unconditional cross-validation, which does not leverage the information captured by the latent variables ( Table 2).
The average predicted species richness per island increased by 12 species between the two inventories, from 107.8 species in the historical data to 119.8 species in the recent data. Proportional species gains and losses occurred rather randomly across the study area (Fig. 1).
The species ecological type (shoreline specialist vs. ordinary species) explained 9.5% of the variance due to fixed effects. The shoreline species increased more over time than did other species (

Colonization-extinction model: Model performance and effects of covariates
The explicit colonization-extinction model explained well the historical species richness (r 2 = 0.99, cross- Notes: Trait r 2 is the percentage of the variance explained by fixed effects explained by traits, i.e., the plant ecological type for the repeated-measures and taxonomic-signal models, and the plant ecological type 9 response variable type interaction for the colonization-extinction model. Parameters are r 2 CV , cross-validated r 2 ; AUC, area under the curve; CV, cross-validation; q, taxonomic signal (95% credible interval); SU, sampling unit.   validated r 2 CV = 0.45), number of colonization events (r 2 = 0.90, r 2 CV = 0.39), and number of extinction events per island (r 2 = 0.92, r 2 CV = 0.65). The ability to discriminate between presences and absences of individual species was reasonably good for historical occurrences (mean Tjur r 2 = 0.31, range = 0 -0.73, mean AUC = 0.94), colonizations (mean Tjur r 2 = 0.20, range = 0 -0.63, mean AUC = 0.92), and extinctions (mean Tjur r 2 = 0.24, range = 0-0.64, mean AUC = 0.92). The island-level latent factors explained more variance (68.9-77.3%) than did the geographical covariates, which led to low species-specific predictive power when cross-validating the model (Table 2).
Larger islands contained more species, were more likely to gain new species through colonization, and were less likely to suffer extinctions (Fig. 2). A similar yet weaker effect was observed for island topographic complexity, but only when allowing island area and neighborhood size to covary with topographic complexity (Fig. 2, upper panels). Neighborhood size had no detectable effects on species occurrences or extinction probabilities, but islands with larger neighborhood sizes had somewhat greater colonization probabilities (Fig. 2).
The response variable type (historical occurrence vs. colonization probability vs. extinction probability) and species ecological type (shoreline specialist vs. ordinary species) jointly explained 39.6% of the variance due to fixed effects (Table 2), largely reflecting the contrasting effects of the covariates on historical occurrences, colonizations, and extinctions (Fig. 2).

Colonization-extinction model: Patterns of cocolonization and co-extinction
To assess whether species tended to colonize or go extinct together, we first asked whether species that were positively or negatively associated in the historical inventory also exhibited shared patterns of colonization and extinction (Fig. 3). Overall, colonization patterns were more predictable from historical associations (historical-historical associations vs. colonization-colonization associations, r 2 = 22.3%) than were extinction patterns (historical-historical associations vs. extinction-extinction associations, r 2 = 5.6%). If only those associations that were either positive with posterior probability> 0.75 or negative with posterior probability> 0.75 were included, these values rose to 46.1% and 19.0%, respectively, suggesting that stronger associations were more predictive than were weaker associations. The relationships were also somewhat stronger when restricting the analysis to the shoreline species only (r 2 = 34.5% and 22.2%, respectively; Fig. 3).
To further explore patterns of residual species associations and their influence on model predictions, we first tabulated the direction and mean strength of species , topographic complexity (standard deviation of elevation points), and local neighborhood size (in log m 2 of land within 5km) on species occurrence probability (gray), conditional extinction rate (red), and conditional colonization rate (blue). Colonization and extinction rates are conditional on historical absence and presence, respectively. Upper panels show net effects of each covariate and lower panels show marginal effects with the non-focal covariates held constant at their mean value. Lines show the mean estimate, and shading shows the 95% credible intervals. associations within and across historical occurrence, colonization, and extinction probabilities across islands. Overall, positive and negative associations were about equally common and equally strong (Fig. 4a). Positive historical associations, colonization-colonization associations and historical-colonization associations were somewhat more common than were the corresponding negative association types, and these association types were slightly stronger than were associations involving extinction probabilities.
We next performed a series of predictions of the current occurrence of each species on each island, where we sequentially excluded different components of the complete model. Predicted current occurrence probabilities were less sensitive to variation in extinction probabilities (mean absolute deviation = 0.022) than to variation in historical occurrence probabilities (0.033) and colonization probabilities (0.045). In turn, holding environmental covariates constant at their mean value had detectable effects on model predictions, as indicated by non-zero deviations of the reduced-model predictions from the full-model predictions (Fig. 4b) and reduced ability of the model to discriminate between true presences and absences when environmental covariates were  ) and their influence on model predictions. (a) Positive (red) and negative (blue) residual correlations between species' historical occurrence (Hist), colonization (Col), and extinction (Ext) probabilities across islands. Solid lines indicate medians, boxes range from the first to third quartile, and whiskers range over 1.59 the interquartile range. Box widths are proportional to the number of coefficients in each category. (b) Influence of environmental covariates and species associations on predicted current occurrence of species across islands, as measured by the mean absolute deviation of reduced-model predictions from full-model predictions. Each plotting symbol gives the mean across species for the influence of the parameter on the x-axis on the focal parameter (historical occurrence, colonization, or extinction), and error bars span values within one standard deviation computed over species. (c) Influence of environmental covariates and species associations on model-predicted current occurrence of species across islands, as measured by the coefficient of discrimination (Tjur r 2 ). Error bars indicate standard errors. held constant (Fig. 4c). In contrast, excluding information about historical species occurrences, colonizations and extinctions influenced model predictions only weakly (Fig. 4b, c).

Taxonomic-signal models
Species' responses to the covariates exhibited a moderately strong yet well-supported taxonomic signal, which was broadly consistent for historical occurrences, colonization probabilities, and extinction probabilities (q = 0.21-0.31, Table 2). Hence, related species tended to respond in a similar manner to island area, topographic complexity and neighborhood size. To explore whether species associations were taxonomically structured, we correlated the taxonomic relatedness of each pair of species with their estimated residual co-occurrences. These correlations were very weak (r < 0.04), suggesting no detectable taxonomic structure to species associations.

DISCUSSION
The explanatory power of our community-wide models was higher for species richness, colonization rates, and extinction rates at the level of islands than for island-specific probabilities of occurrence, colonization, and extinction of individual species. At the species level, the models explained historical and current species distributions better than colonization and extinction probabilities. At face value, these observations suggest either dissimilar environmental drivers of species distributions vs. colonization/extinction dynamics, or that colonizations and extinctions are stochastic and "unpredictable" events, as would indeed be the implicit assumption of "neutral" models of island biogeography (MacArthur and Wilson 1967). Our species occurrence data were collected as presence/absence at the level of islands, ignoring variation in abundances and microhabitats of species, thus adding noise to the species-specific results. Furthermore, we focused on broad-scale geographical covariates, and assumed that additional environmental variation among islands would be captured by the island-level latent factors included in our models (Ovaskainen et al. 2016b). The "success" of the latent factors in capturing unmeasured environmental variation was apparent from the high explanatory power at the level of islands. The latent factors also represent potential species interactions (Table 1), yet we failed to detect a strong signal of species associations on metacommunity dynamics in this system. In the following, we discuss the implications of our results for the predictability of community-level responses to environmental change.

Spatial and temporal diversity patterns
Before turning to patterns of vegetation change in the study area over the last 80 yr, some notes on spatial diversity patterns are warranted. The species richness of islands increased with area following a simple power relationship, a well-known biogeographic pattern. However, species diversity is not only a function of area per se, but also of the increased environmental heterogeneity normally associated with increasing area (MacArthur and Wilson 1967, Kohn and Walsh 1994, Lundholm 2009, Stein et al. 2014. Larger islands typically contain a greater number of distinct habitat types, and consequently more species (Kohn andWalsh 1994, Hannus andvon Numers 2008). Furthermore, recent work adopting a more fine-scale view of environmental heterogeneity has documented that topographic complexity creates mosaics of variation in microclimate at small scales (Scherrer and K€ orner 2010, Lenoir et al. 2013, Opedal et al. 2015. The size and environmental heterogeneity of islands may also affect community dynamics over time. Larger islands are likely to sustain larger average population sizes and consequently to exhibit reduced extinction risk, and are bigger "targets" for colonization (MacArthur andWilson 1967, Cody 2006). Our results are largely consistent with these expectations, as larger islands had both greater colonization rates and lower extinction rates (Fig. 2).
Topographic complexity may affect community dynamics by altering the probabilities of extinction and colonization . First, in the event of climate change, fine-scale environmental variation may allow species to track changing environmental conditions by migrating into local microrefugia within complex landscapes, thus reducing extinction risk. Second, the greater environmental heterogeneity of topographically more complex landscapes may increase the probability of incoming propagules arriving in a suitable site for establishment. These expectations were met in our results, but only when island area was allowed to covary with topographic complexity (Fig. 2, upper row). Thus, although we cannot conclusively separate the area effect from effects of unmeasured environmental heterogeneity (Kohn and Walsh 1994), we failed to detect a strong effect of microclimatic variation independently of island area.
Colonization rates may also depend on the local neighborhood size (i.e., the distance to and size of nearby islands), and hence the intensity of propagule influx. As expected, we detected a weak positive effect of neighborhood size on colonisation rates (Fig. 2). Neighborhood size has also been hypothesized to affect observed extinction rates through the rescue effect (Brown and Kodric-Brown 1977), but we detected no such effects in the archipelago plant metacommunity.

Ecological similarity, taxonomic relatedness, and species' environmental responses
If similarities in species' environmental responses are to some extent structured by their ecological strategies, phylogenetic histories, or both, this would improve our ability to predict how communities change. Instead of focusing on individual functional traits, we chose to focus on a group of species that are ecologically similar by using the same general habitat. A main finding was that the occurrence probabilities of shoreline species had increased more over the last 80 yr than had those of other species, and including the separation of these broad species groups into the analysis therefore proved meaningful. von Numers (2011) discusses possible causes of increased abundance of shoreline species in the study area, emphasizing the roles of changes in grazing pressure and increased nutrient availability caused by the eutrophicated sea.
The explanatory power arising from ecological similarity may rise further by adopting a more mechanistic approach focused on functional traits explicitly hypothesized to mediate species´responses to specific environmental variables (e.g., Blonder et al. 2018). For example, working in the Barkley sound archipelago in Canada (see Cody 2006), Burns and Neufeld (2009) reported reduced extinction probabilities of species with thicker leaves, consistent with an important role of physical stress tolerance in the rough environment characterizing small islands. Similarly, we may expect population dynamics to depend on generation times, an effect that could be easily assessed by including generation time or at least the separation of annual, biennial, and perennial species as a trait in joint models.
The shoreline species are an ecological group representing numerous families and genera, many of which also comprise species confined to different habitat types. Our taxonomic analysis therefore assessed a largely independent axis of species similarity, and revealed moderately strong similarities in the responses of related species to island size, topographic complexity, and neighborhood size. This finding is encouraging, especially because the strength of the signal was almost certainly reduced by the coarse resolution of our analysis, which was restricted to separating higher taxa, and more detailed phylogenetic data could reveal patterns hidden at the current level of analysis. Nevertheless, these results illustrate the utility of phylogenetic information in understanding the dynamics of communities.

The role of species associations in community structure and dynamics
There is growing consensus that species interactions will affect the dynamics of communities over time, yet there is substantial uncertainty in exactly how this will occur (e.g., Holt 1997, Blois et al. 2013, Magurran 2016. For example, an important yet largely unexplored question concerning the temporal dynamics of biotic communities is whether associated species tend also to exhibit shared dynamics over time. Our results suggest that the answer to this question depends on the specific component of community dynamics considered, because the direction and strength of species associations in the historical data explained patterns of joint colonizations better than joint extinctions (Fig. 3). The difficulty in identifying drivers of colonizations and particularly extinctions was further confirmed by our assessment of species associations within and among historical occurrences, colonizations, and extinctions (Fig. 4). Despite our focus on estimating the overall influence of species associations on model predictions, rather than inferring specific associations, we failed to detectably improve predictions of community dynamics by considering species associations. The difficulty of detecting net effects of extinctions on community dynamics is compounded by their relative infrequency, which reduces the potential influence of extinctions on current community structure.
We cannot at this time ascertain whether the weak influence of species associations on model predictions in our study resulted from limited importance of species interactions for community dynamics in the study area, or from the coarse resolution of our data precluding the detection of such effects. As part of his long-term study of the plant communities of Barkley sound islands, Cody (2006:Chapter 5) collected fine-scale data on the distributions and dynamics of species hypothesized to compete for a certain resource. While some closely related or ecologically similar species did indeed appear to exhibit disjunct distributions, these cases appeared to be exceptions, suggesting that overall effects of species interactions on community dynamics may be hard to detect even if they exist.
A related issue is whether joint species distribution models are indeed able to detect species interactions (Barner et al. 2018, Zurell et al. 2018). There are several ways to interpret residual species associations inferred from joint species distribution models, including joint responses to unmeasured environmental variables and direct interactions between species (Table 1, and see Ovaskainen et al. 2016a). Because the covariates included in our colonisation-extinction model described island geography (island size, topographic complexity, and neighborhood size) rather than island-specific aspects of the abiotic environments, the conservative interpretation is that species exhibiting positive associations respond to some extent to the same environmental variables, while those species exhibiting negative associations have contrasting responses. Furthermore, because the analysis of residual species associations in joint species distribution models is essentially a correlation analysis, we cannot with current methods distinguish the effect of species A going extinct on the probability of species B colonizing from the effect of species B colonizing on the probability of species A going extinct. While these scenarios may differ in their ecological interpretation (Table 1), they are indistinguishable from the symmetric association matrices produced by current joint species distribution models.
Finally, we detected essentially no taxonomic structure to residual species associations beyond those Article e03067; page 10 ØYSTEIN H. OPEDAL ET AL. Ecology,Vol. 101,No. 8 explained by shared environmental responses. Hence, there was no strong overall tendency for related species to be positively or negatively associated on the study islands, or to affect each other's probabilities of colonisation and extinction. Despite the long-standing notion that related species should be ecologically similar and thus tend to enter competition with each other (e.g., Darwin 1859), our results are consistent with several studies reporting weak effects of phylogenetic relatedness on competitive interactions (Cahill et al. 2008, Kunstler et al. 2012. We cannot, however, ascertain whether this pattern would hold also with more detailed phylogenetic data.

CONCLUSIONS: SPECIES NON-INDEPENDENCE AND THE PRE-DICTABILITY OF BIOTIC COMMUNITY DYNAMICS
We set out to assess the utility of several quantifiable components of species non-independence in refining predictions of community dynamics. Despite the coarse resolution of our data, both species' ecological groups and taxonomic relatedness explained detectable proportions of the variance in community structure and dynamics. Our results also suggest that species associations inferred from snapshot data can be informative about community dynamics, because historically associated species tended also to colonize the same islands. However, explicit modeling of joint colonisation and extinction dynamics had only minor influences on model predictions, suggesting limited influence of species interactions on community dynamics in this system. These results illustrate both the complexity of biotic community dynamics, and the usefulness of considering community dynamics within a joint framework, rather than treating species as independent entities.