Can interaction coefficients be determined from cencus data?

Oecologia 9 Springer-Verlag1985

Oecologia (Berlin) (1985) 66:194-198

Can interaction coefficients be determined from census data? Michael L. Rosenzweig ~, Zvika Abramsky 2, Burt Kotler ~*, and William Mitchell ~ 1 Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, AZ 85721, USA 2 Department of Biology, Ben-Gurion University of the Negev, Beersheva, Israel Summary. The method of estimating interactions proposed independently by Pimm and Schoener is studied using field data from the community of rodents which lives in the arid, rocky habitats of Israel. One important problem the method addresses is how to remove the effects of habitat heterogeneity on the estimate. We tried six different variations of the analysis scheme outlined by Crowell and Pimm, and found their results qualitatively inconsistent. This was especially true when we compared the results produced from separate habitat variables with those produced from the principal components of the habitat variation. Another problem, this one not previously addressed, is great variation in the average abundance of the different species. We discovered that the ratio of the average abundances of two species is the best predictor of the value of their coefficients of interaction. Common species appear to have weak influence on rare ones; rare ones appear to have strong influence on common ones. The statistical mechanism which produces this relationship is clear, indicating that the relationship is an artifact.

A practical method of estimating coefficients in the field has been sought by ecologists since 1967 (MacArthur and Levins 1967). The hope has been to bypass perturbation analysis which is labor-intensive, time-consuming and expensive. A method which has recently been used is the PimmSchoener technique of repeated joint censuses (Schoener 1974; Crowell and Pimm 1976). The technique involves collecting independent census estimates of species at various times and places. Then the censuses are regressed against each other. The regression coefficients are supposed to be estimates of the coefficients of interaction near equilibrium. Hallett and Pimm (1979) did some simulation work on this technique and predicted it should give reasonable results for Lotka-Volterra systems in homogeneous environments with competitors that are not highly asymmetrical in their effects on each other and are also similar in carrying capacities. Its precision diminished as asymmetry increased. But Hallett et al. (1979) did not simulate non-linearity, heterogeneous environments or markedly unequal carrying capacities. * P e r m a n e n t a d d r e s s : Desert Research Institute, Ben-Gurion University, S'deh Boker, Israel

Previously, both Schoener (1975) and Crowell et al. (1976) had suggested methods to deal with habitat heterogeneity. In Schoener's, carrying capacities are estimated by finding sets of places where each species occurs in the absence of all others whose densities are being regressed against it. This may not be practical except in archipelagoes of many small islands with few species. The technique of Pimm (Crowell et al. 1976) promised to be more widely applicable because data entirely from areas of sympatry may be analyzed. Pimm's methods involves removing the effects of habitat heterogeneity by first regressing the censuses against a variety of habitat variables. Unfortunately there is more than one way to execute this method and the manner of its execution affects the results its yields (Rosenzweig et al. 1984). In this paper, we examine the results of censuses on 6 species of lithophilic rodents of the Negev Desert, Israel. In addition to our censuses, we collected a considerable amount of habitat information allowing us further to explore the consequences of habitat heterogeneity on the Pimm-Schoener coefficient estimates. Because of the rodent's great disparity in commonness, this set of species also allowed us to explore the effects of highly unequal densities on the coefficient estimates.

Sampling methods Seven sites were chosen in the rocky, mountainous sections of the Negev. They varied from the edge of the SyrianAfrican rift valley (Ein Gedi) to the extreme southern Judaean Hills (Yattir) to the central Negev's deepest canyon (Makhtesh Ramon). They varied in altitude, plant cover and composition, geology and productivity. See Abramsky et al. (in press) for a list of the sites and more description. At each site, four sampling grids were set up. Each consisted of 20 trapping stations at 20 m intervals, with two collapsible Sherman live traps per station. Grids could not be regular because of the extremely broken landscape. Thus they were each set up as a sinuous line following contours. We attempted at each site to distribute the grids in different macrohabitats. Usually this meant that two of the grids were along steep wadi walls or cliffs and two on much more level pebble-covered substrates. But the sites varied so much that it was not possible to find all macrohabitats at all sites. At each site during the academic year 1981-1982, the rodents were censused first during the cool, wet season and

195 then during the warm, dry season. Each census lasted 5 days. Captured mice were toe-clipped for identification and released. The census variable for a species is just the total number of individuals captured during the five days. (A more detailed description of census techniques is found in Abramsky et al. in press). The two censuses were averaged at each grid to yield an average for the year. A number of habitat variables were measured at each station and averaged to determine a grid's value. These are described fully elsewhere (Abramsky et al. in press). They included steepness, substrate rugosity, percent cover, compass heading (e.g. north-facing), number of plants in a sample area and distance between plants. The coefficient of variation of some of these was also computed for each grid, namely the CV of plant distance, compass heading, rugosity, and steepness. Hence, each grid was characterized by 10 habitat variables. There was essentially no annual cover in any grid and we did not use either the altitude or the estimated annual precipitation at a site.

are species (with each species taking its turn as the dependent variable). c) Free regression. Both the previous methods are biased in favor of the habitat variables. To see if this mattered, we also loaded the habitat and census variables all at once and allowed the regression to proceed without regard to whether a habitat or a census variable was chosen at a particular step. (We did, however, maintain the protocol about squared P C A terms and their crossproducts.) Subset 2: N o principal components. This subset is like the first except that the raw habitat variables were used. There were no squared terms or crossproducts. The percent cover was arcsine transformed. The compass heading, 0 was transformed to sin 0/2. Otherwise a reading of 1~ and another of 359 ~ would have appeared quite dissimilar instead of almost the same. Then the three types of regression analysis were performed as for subset 1.

Statistical methods

Results

We employed six different techniques to estimate the interaction coefficients. These differ only in the manner in which they deal with and attempt to account for habitat heterogeneity. The six are arranged into two subsets.

Table 1 shows the values of the regression coefficients when they met our criterion of significance. M a n y of the rows show considerable consistency, indicating that the particular method played only a small part in determining the coefficient. F o r example, all methods conclude that Acomys russatus does not affect Gerbillus dasyurus. A n d all agree that Sekeetamys calurus significantly depresses G. dasyurus with a coefficient of about - 1.5. Altogether, 14 of the 30 rows show complete consistency. Of the remaining 16, 6 show substantial consistency in that only one of the 6 columns differs from the others. For example, the effect of A. cahirinus on S. calurus is not significantly different from zero in 5 of 6 columns, but is - 0 . 0 6 for one. In sum, 20 of the 30 rows are consistent or nearly so. The remaining 10 rows fall into two categories. In one category, the output appears very variable and impossible to interpret. There are 4 such rows. All deal with effects on the 3 species which were least c o m m o n in these censuses (S. calurus, Meriones crassus and Eliomys melanurus). Quite probably larger sample sizes would have produced greater consistency. The remaining 6 rows show another pattern: regression with the raw habitat variables gave one result, but another result was obtained with principal components. This pattern appeared only in effects on the two spiny mice, A. russatus and A. cahirinus, but it was the predominant pattern for both of them.

Subset 1 : Reduction of dimensionality of the habitat variables by principal components analysis. The set of ten habitat variables was reduced to 3 principal components. These accounted for 88.4% of the habitat variation. In order to map nonlinearities on the principal component space, we generated the squares and crossproducts of the principal components and used these values as additional independent variables in the regression analyses. However, there was a protocol: no squared term was allowed to enter a regression equation until after its linear term had proved significant; and no crossproduct was allowed until both its linear components had proved significant. Significance was defined as F-ratios in excess of 2.0 both for entry into the equation and maintenance therein. This definition was maintained in regard to all variable types and for all the studies throughout the paper. It generally corresponds to 0 . 1 5 < P < 0 . 2 0 for this set of studies and so is a bit less conservative than usual. But it is only a bit less conservative than Crowell and Pimm's use of P = 0.10 for habitat variables and is more conservative than omitting a significance test which is what Pimm (pers. comm.) has recommended for interaction variables (i.e. censuses) themselves. a) Original Crowell and Pimm method. Stepwise multiple linear regression is performed first on the habitat variables (in this case the principal components) as independent variables and the censuses of one species as the dependent variable. Then the censuses of the other species are allowed to enter the equation. b) Residual analysis (Rosenzweig et al. 1984). After stepwise regression of habitat variables on all mouse species (separately) is performed as in a, the residuals of each species are calculated from their habitat regressions. Then all these residuals are entered into a data matrix and stepwise regression is performed on the matrix as many times as there

Discussion

The overall consistency of the results is, at first glance, reassuring. In particular, the results seem to point squarely at an asymmetrical competitive interaction between G. dasyurus and S. calurus, a similarly asymmetrical competition between M. crassus and A. cahirinus, and an unexpected mutualism between A. russatus and S. calurus. But appearances can be deceiving. In the first place, the consistency of the G. dasyurus results is not so reassuring after all because habitat does not account for this species' census. Using the free method, we did see some habitat effects, but only after the density of S. calurus was accounted for. In other words, in every

196 Table 1. Significant coefficients of interaction for lithophilic rodents in Israel according to 6 methods for removing habitat effects Species

Method

Receiver PCA Emiter Pimm

Non-PCA Resid- Free uals

G E R B I L L U S D A S Y U R US A. cahirinus -0.37 A. russatus S. calurus -1.32 -1.32 m . crassus E. melanurus ACOMYS CAHIRINUS G. dasyurus -0.28 -0.32 A. russatus --0.47 --0.49 S. calurus M. crassus --0.86 --0.98 E. melanurus A C O M YS R USSA TUS G. dasyurus A. cahirinus -0.26 -0.32 S. calurus +1.17 +0.98 M. crassus -0.63 -0.64 E. melanurus SEKEETA M Y S CAL UR US G. dasyurus -0.16 -0.11 A. cahirinus A. russatus +0.14 +0.17 M. crassus + 0.24 E. melanurus MERIONES CRASSUS G. dasyurus A. cahirinus -0.19 -0.21 A. russatus - 0.20 S. calurus +0.44 E. melanurus --1.28 ELIOMYS MELANURUS G. dasyurus A. cahirinus -- 0.04 A. russatus S. calurus M. crassus --0.15 -0.10

Pimm

Resid- Free uals

-1.32

-1.55

-1.66

--1.12 -- 3.09

-1.37 - 2.96

--1.19 -- 3.85

+0.77

+0.94

+1.08

+ 1.01

+ 1.46

+ 1.36

-0.1~l

-0.12

-0.10

+0.23

+0.16

+0.20

-0.11 --0.06 +0.19

-0.76

--0.59

-1.67

-0.25 -0.55 -1.16

-0.26 +1.16 -0.63

-0.22

-0.18

-0.17

-0.17

-1.56

+0.30 --1.33

-1.43

--1.84

+ 0.03 - 0.04

+ 0.03 - 0.05 + 0.06

-0.16

--0.14

-0.18

--0.15

case, the first independent variable chosen was the density o f S. calurus. The habitat matrix was not used, so it did not matter how we instructed the computer to treat it. The only reason for the m i n o r inconsistency (in effect o f A. cahirinus on G. dasyurus) is that habitat variables did affect the residuals of the other species. Otherwise an inconsistency would be literally impossible. W h a t is one to m a k e o f the results for the spiny mice? The principal components say they compete with each other. The raw habitat variables say they do not. They cannot both be correct. A n d that is only one o f the inconsistencies in the results for the Acornys. Altogether 18 o f the 60 estimates o f these ten coefficients are inconsistent. That represents an error rate o f 30% for the A c o m y s part o f Table 1. It means consistent estimates were p r o d u c e d for only 4 of the ten (and we do not yet know whether these four are right; consistency is only a first test). W e need

a m e t h o d which will give us an answer that does not depend on how we factor out the habitat relationships. The picture we are left with is rather discouraging. One species has a consistent set o f results because it is unaffected by habitat, two are inconsistent depending on whether raw variables or principal components are used, and three were not so consistent perhaps for lack o f enough data. We must conclude that properly accounting for habitat heterogeneity remains a m a j o r obstacle to using the Pimm-Schoener technique with confidence. Hallett (1982) has p r o p o s e d that his analyses using Pimm-Schoener reveal a general asymmetry in the structure o f communities. Those species which are habitat specialists exert strong competitive effects on those which are not. Those which are spread fairly evenly t h r o u g h o u t m a n y habitats exert only weak interactive effects on the specialists. A n o t h e r way o f putting the same thesis is this: the more a species census is determined by habitat, the less the magnitude of the competitive coefficients it suffers from other species. The index o f habitat specialization is the p r o p o r t i o n of a species census variance explained by habitat variables. This is readily calculated for each species. But the index of interactive strength is not so sharply defined. Are we to average all interaction coefficients including zeros ? W h a t a b o u t positive coefficients? There are three meaningful averages: with the zeros and positives; without positives; without positives or zeros. We tried all three. Moreover, we can look at species as receivers o f interaction (R) or as emitters o f interaction (E). Two vantage points times three coefficient definitions equals six ways o f looking for Hallett's pattern. W e tried all six, but found none to yield a significant result. A l t h o u g h there was no evidence o f Hallett's pattern in our data, it did seem that, when an interaction was significant, some species were emitting large coefficients consistently and others were emitting small ones just as consistently. In fact, there is a pattern in these results which has a very simple explanation. Let y be the average o f the absolute values o f the significant coefficients of species E's effect on species R. Let x be the ratio o f the variance of R ' s census to the variance o f E's. These variances are not residuals; they are the preregression variances. Figure 1 shows in y plotted against In x. Clearly In x is a fine predictor of in y (r = 0.976; P < 10 - 3 that regression coefficient is zero). Yet this is nothing more than a statistical artifact. A regression coefficient is the p r o d u c t of two numbers: the correlation coefficient and the ratio of standard deviations o f the dependent (numerator) and independent (denominator) variabes. This ratio is, of course, very nearly the same as x in Fig. 1. Thus, if all the correlation coefficients are the same in a set o f analyses such as we have done, any trend in the regression coefficients is determined solely by the ratio. N o w all correlation coefficients are not the same in our d a t a set, but they are severely constrained. To be significant, they h a d to have absolute values of at least 0.25 or 0.30 (depending on the size o f the d a t a set). Moreover, given the usual noise in natural systems, we are unlikely to encounter r's whose absolute values are greater than 0.7 or 0.8 (and we did not). So, in practice, the regression coeffid e n t will be the multiple of a number between 0.25 and

197

5 -2

-%

I

I -2

I

I 0

[

I 2

I

I 4

Ln [VAR#VARE] Fig. l. The magnitude of significant coefficients of interaction is determined by the ratio of variances in population size of the species. The x-axis is the natural logarithm of the ratio of the variances of the receiver's population to the emitter's. The y-axis is the natural logarithm of the average absolute value of the significant coefficients of interaction of E on R. Significance is determined by F = 2.0 (see text). The regression equation is: y=0.54x--l.07, for which R2=0.95, d.f.=21 and P <