Fine-scale habitat modeling of a top marine predator: do prey data improve predictive capacity?

Predators and prey assort themselves relative to each other, the availability of resources and refuges, and the temporal and spatial scale of their interaction. Predictive models of predator distributions often rely on these relationships by incorporating data on environmental variability and prey availability to determine predator habitat selection patterns. This approach to predictive modeling holds true in marine systems where observations of predators are logistically difficult, emphasizing the need for accurate models. In this paper, we ask whether including prey distribution data in fine-scale predictive models of bottlenose dolphin (Tursiops truncatus) habitat selection in Florida Bay, Florida, U.S.A., improves predictive capacity. Environmental characteristics are often used as predictor variables in habitat models of top marine predators with the assumption that they act as proxies of prey distribution. We examine the validity of this assumption by comparing the response of dolphin distribution and fish catch rates to the same environmental variables. Next, the predictive capacities of four models, with and without prey distribution data, are tested to determine whether dolphin habitat selection can be predicted without recourse to describing the distribution of their prey. The final analysis determines the accuracy of predictive maps of dolphin distribution produced by modeling areas of high fish catch based on significant environmental characteristics. We use spatial analysis and independent data sets to train and test the models. Our results indicate that, due to high habitat heterogeneity and the spatial variability of prey patches, fine-scale models of dolphin habitat selection in coastal habitats will be more successful if environmental variables are used as predictor variables of predator distributions rather than relying on prey data as explanatory variables. However, predictive modeling of prey distribution as the response variable based on environmental variability did produce high predictive performance of dolphin habitat selection, particularly foraging habitat.     

Bridging gaps in demographic analysis with phylogenetic imputation

Phylogenetically informed imputation methods have rarely been applied to estimate missing values in demographic data but may be a powerful tool for reconstructing vital rates of survival, maturation, and fecundity for species of conservation concern. Imputed vital rates could be used to parameterize demographic models to explore how populations respond when vital rates are perturbed. We used standardized vital rate estimates for 50 bird species to assess the use of phylogenetic imputation to fill gaps in demographic data. We calculated imputation accuracy for vital rates of focal species excluded from the data set either singly or in combination and with and without phylogeny, body mass, and life-history trait data. We used imputed vital rates to calculate demographic metrics, including generation time, to validate the use of imputation in demographic analyses. Covariance among vital rates and other trait data provided a strong basis to guide imputation of missing vital rates in birds, even in the absence of phylogenetic information. Mean NRMSE for null and phylogenetic models differed by <0.01 except when no vital rates were available or for vital rates with high phylogenetic signal (Pagel's λ > 0.8). In these cases, including body mass and life-history trait data compensated for lack of phylogenetic information: mean normalized root mean square error (NRMSE) for null and phylogenetic models differed by <0.01 for adult survival and <0.04 for maturation rate. Estimates of demographic metrics were sensitive to the accuracy of imputed vital rates. For example, mean error in generation time doubled in response to inaccurate estimates of maturation time. Accurate demographic data and metrics, such as generation time, are needed to inform conservation planning processes, for example through International Union for Conservation of Nature Red List assessments and population viability analysis. Imputed vital rates could be useful in this context but, as for any estimated model parameters, awareness of the sensitivities of demographic model outputs to the imputed vital rates is essential.     

Tropical deforestation in Madagascar: analysis using hierarchical,spatially explicit,Bayesian regression models

Establishing cause–effect relationships for deforestation at various scales has proven difficult even when rates of deforestation appear well documented. There is a need for better explanatory models, which also provide insight into the process of deforestation. We propose a novel hierarchical modeling specification incorporating spatial association. The hierarchical aspect allows us to accommodate misalignment between the land-use (response) data layer and explanatory data layers. Spatial structure seems appropriate due to the inherently spatial nature of land use and data layers explaining land use. Typically, there will be missing values or holes in the response data. To accommodate this we propose an imputation strategy. We apply our modeling approach to develop a novel deforestation model for the eastern wet forested zone of Madagascar, a global rain forest “hot spot”. Using five data layers created for this region, we fit a suitable spatial hierarchical model. Though fitting such models is computationally much more demanding than fitting more standard models, we show that the resulting interpretation is much richer. Also, we employ a model choice criterion to argue that our fully Bayesian model performs better than simpler ones. To the best of our knowledge, this is the first work that applies hierarchical Bayesian modeling techniques to study deforestation processes. We conclude with a discussion of our findings and an indication of the broader ecological applicability of our modeling style.     

Confidence intervals for the mean of the delta-lognormal distribution

Data that are skewed and contain a relatively high proportion of zeros can often be modelled using a delta-lognormal distribution. We consider three methods of calculating a 95% confidence interval for the mean of this distribution, and use simulation to compare the methods, across a range of realistic scenarios. The best method, in terms of coverage, is that based on the profile-likelihood. This gives error rates that are within 1% (lower limit) or 3% (upper limit) of the nominal level, unless the sample size is small and the level of skewness is moderate to high. Our results will also apply to the delta-lognormal linear model, when we wish to calculate a confidence interval for the expected value of the response variable, given the value of one or more explanatory variables. We illustrate the three methods using data on red cod densities, taken from a fisheries trawl survey in New Zealand.

Environmental filtering and land-use history drive patterns in biomass accumulation in a mediterranean-type landscape

Aboveground biomass (AGB) reflects multiple and often undetermined ecological and land-use processes, yet detailed landscape-level studies of AGB are uncommon due to the difficulty in making consistent measurements at ecologically relevant scales. Working in a protected mediterranean-type landscape (Jasper Ridge Biological Preserve, California, USA), we combined field measurements with remotely sensed data from the Carnegie Airborne Observatory's light detection and ranging (lidar) system to create a detailed AGB map. We then developed a predictive model using a maximum of 56 explanatory variables derived from geologic and historic-ownership maps, a digital elevation model, and geographic coordinates to evaluate possible controls over currently observed AGB patterns. We tested both ordinary least-squares regression (OLS) and autoregressive approaches. OLS explained 44% of the variation in AGB, and simultaneous autoregression with a 100-m neighborhood improved the fit to an r2 = 0.72, while reducing the number of significant predictor variables from 27 variables in the OLS model to 11 variables in the autoregressive model. We also compared the results from these approaches to a more typical field-derived data set; we randomly sampled 5% of the data 1000 times and used the same OLS approach each time. Environmental filters including incident solar radiation, substrate type, and topographic position were significant predictors of AGB in all models. Past ownership was a minor but significant predictor, despite the long history of conservation at the site. The weak predictive power of these environmental variables, and the significant improvement when spatial autocorrelation was incorporated, highlight the importance of land-use history, disturbance regime, and population dynamics as controllers of AGB.     

Using habitat suitability index and particle dispersion models for early detection of marine invaders.

Eradication and control of invasive species are often possible only if populations are detected when they are small and localized. To be efficient, detection surveys should be targeted at locations where there is the greatest risk of incursions. We examine the utility of habitat suitability index (HSI) and particle dispersion models for targeting sampling for marine pests. Habitat suitability index models are a simple way to identify suitable habitat when species distribution data are lacking. We compared the performance of HSI models with statistical models derived from independent data from New Zealand on the distribution of two nonindigenous bivalves: Theora lubrica and Musculista senhousia. Logistic regression models developed using the HSI scores as predictors of the presence/absence of Theora and Musculista explained 26.7% and 6.2% of the deviance in the data, respectively. Odds ratios for the HSI scores were greater than unity, indicating that they were genuine predictors of the presence/ absence of each species. The fit and predictive accuracy of each logistic model were improved when simulated patterns of dispersion from the nearest port were added as a predictor variable. Nevertheless, the combined model explained, at best, 46.5% of the deviance in the distribution of Theora and correctly predicted 56% of true presences and 50% of all cases. Omission errors were between 6% and 16%. Although statistical distribution models built directly from environmental predictors always outperformed the equivalent HSI models, the gain in model fit and accuracy was modest. High residual deviance in both types of model suggests that the distributions realized by Theora and Musculista in the field data were influenced by factors not explicitly modeled as explanatory variables and by error in the environmental data used to project suitable habitat for the species. Our results highlight the difficulty of accurately predicting the distribution of invasive marine species that exhibit low habitat occupancy and patchy distributions in time and space. Although the HSI and statistical models had utility as predictors of the likely distribution of nonindigenous marine species, the level of spatial accuracy achieved with them may be well below expectations for sensitive surveillance programs.     

Interpolation performance of a spatio-temporal model with spatially varying coefficients: application to PM<Subscript>10</Subscript> concentrations in Rio de Janeiro

In this work we present a Bayesian analysis in linear regression models with spatially varying coefficients for modeling and inference in spatio-temporal processes. This kind of model is particularly appealing in situations where the effect of one or more explanatory processes on the response present substantial spatial heterogeneity. We describe for this model how to make inference about the regression coefficients and response processes under two scenarios: when the explanatory processes are known throughout the study region, and when they are known only at the sampling locations. Using a simulation experiment we investigate how parameter inference and interpolation performance are affected by some features of the data and prior distribution that is used. The proposed methodology is used to model the dataset on PM₁₀ levels in the metropolitan region of Rio de Janeiro presented in Paez and Gamerman (2003).     

Time-varying models for extreme values

We propose a new approach for modeling extreme values that are measured in time and space. First we assume that the observations follow a Generalized Extreme Value (GEV) distribution for which the location, scale or shape parameters define the space–time structure. The temporal component is defined through a Dynamic Linear Model (DLM) or state space representation that allows to estimate the trend or seasonality of the data in time. The spatial element is imposed through the evolution matrix of the DLM where we adopt a process convolution form. We show how to produce temporal and spatial estimates of our model via customized Markov Chain Monte Carlo (MCMC) simulation. We illustrate our methodology with extreme values of ozone levels produced daily in the metropolitan area of Mexico City and with rainfall extremes measured at the Caribbean coast of Venezuela.     

Landscape history, fragmentation, and patch occupancy: models for a forest bird with limited dispersal.

We developed and tested patch occupancy models for an endemic understory bird with limited dispersal ability, the Chucao Tapaculo (Scelorchilus rubecula), in two South American temperate rain forest landscapes that differed in levels and duration of forest loss. We assessed cover changes since 1961 in each landscape and surveyed patches for Chucao Tapaculo occupancy. We then developed incidence-based predictive models independently for each landscape and tested each model reciprocally in the alternative study area. We thereby assessed the domain of model applicability and identified those predictor variables with general effects and those that varied between the two landscapes. The two models were consistent regarding variable selection, and predictive accuracy of each model was high in the landscape where training data were collected. However, the models differed substantially in the magnitudes of effects related to patch size, with larger unoccupied patches observed in the landscape with the more advanced stage of fragmentation. Due to this discrepancy, each model performed poorly when applied to the alternative landscape, potentially reflecting the contrasting stages of habitat loss. Although it was impossible to dissociate effects of level and duration of forest loss, we viewed the landscapes as representing two extremes along a continuum of fragmentation, providing insights into potential trajectories for portions of the biome where deforestation is occurring. Further, our data suggest that static equilibrium models developed from distribution patterns in recently fragmented landscapes may overestimate persistence when used as a forecasting tool, or when extrapolated to alternative landscapes where fragmentation is more advanced. Thus, we suggest that landscapes used as standards for model building should be selected with caution. We recommend that distribution patterns be obtained from landscapes where fragmentation is advanced, preferably with histories of fragmentation long enough that time-delayed extinctions already would have occurred.     

Using data augmentation via the Gibbs Sampler to incorporate missing covariate structure in linear models for ecological assessments

Missing covariate values in linear regression models can be an important problem facing environmental researchers. Existing missing value treatment methods such as Multiple Imputation (MI), the EM algorithm and Data Augmentation (DA) have the assumption that both observed and unobserved data come from the same distribution, most commonly a multivariate normal or a conditionally multivariate normal family. These methods do try to incorporate the missing data mechanism and rely on the assumption of Missing At Random (MAR). We present a DA method which does not rely on the MAR assumption and can model missing data mechanisms and covariate structure. This method utilizes the Gibbs Sampler as a tool for incorporating these structures and mechanisms. We apply this method to an ecological data set that relates fish condition to environmental variables. Notice that the presented DA method detects relationships that are not detected when other missing data methods are employed.

A stochastic weather generation model for Alaska

A specific problem encountered in ecosystem-level simulation of Arctic ecosystems is the depth and extent of the driving variable record. Often, climate records are of short duration, gathered at locations different from the area to be simulated, or do not contain all the variables required by a given model. This paper addresses this problem for ecosystem simulation in Alaska with the development of a weather generator. The generator, called WGENAL, is based on the WGEN climate generator developed and validated in the 48 conterminous states. Because of the extreme variability of weather in Alaska that is not accommodated by the statistical metrics in the earlier model, a new climate generator was developed. WGENAL generates daily values of precipitation, maximum temperature, minimum temperature, solar radiation, and wind run. Precipitation is generated using a Markov chain-gamma model. A two-parameter gamma distribution is used to generate wind run. Temperatures and solar radiation are generated using procedures developed in the earlier study. Validation of the generator shows it provides adequate diurnal and seasonal weather records for Fairbanks. Other comparisons of synthetic weather with observed weather for sites north of the Brooks Range in Alaska are also within the error of the original data.     

The roles of climate, phylogenetic relatedness, introduction effort, and reproductive traits in the establishment of non-native reptiles and amphibians

We developed a method to predict the potential of non-native reptiles and amphibians (herpetofauna) to establish populations. This method may inform efforts to prevent the introduction of invasive non-native species. We used boosted regression trees to determine whether nine variables influence establishment success of introduced herpetofauna in California and Florida. We used an independent data set to assess model performance. Propagule pressure was the variable most strongly associated with establishment success. Species with short juvenile periods and species with phylogenetically more distant relatives in regional biotas were more likely to establish than species that start breeding later and those that have close relatives. Average climate match (the similarity of climate between native and non-native range) and life form were also important. Frogs and lizards were the taxonomic groups most likely to establish, whereas a much lower proportion of snakes and turtles established. We used results from our best model to compile a spreadsheet-based model for easy use and interpretation. Probability scores obtained from the spreadsheet model were strongly correlated with establishment success as were probabilities predicted for independent data by the boosted regression tree model. However, the error rate for predictions made with independent data was much higher than with cross validation using training data. This difference in predictive power does not preclude use of the model to assess the probability of establishment of herpetofauna because (1) the independent data had no information for two variables (meaning the full predictive capacity of the model could not be realized) and (2) the model structure is consistent with the recent literature on the primary determinants of establishment success for herpetofauna. It may still be difficult to predict the establishment probability of poorly studied taxa, but it is clear that non-native species (especially lizards and frogs) that mature early and come from environments similar to that of the introduction region have the highest probability of establishment.     

Predicting fish diet composition using a bagged classification tree approach: a case study using yellowfin tuna (<Emphasis Type="Italic">Thunnus albacares</Emphasis>)

We provided a classification tree modeling framework for investigating complex feeding relationships and illustrated the method using stomach contents data for yellowfin tuna (Thunnus albacares) collected by longline fishing gear deployed off eastern Australia between 1992 and 2006. The non-parametric method is both exploratory and predictive, can be applied to varying size datasets and therefore is not restricted to a minimum sample size. The method uses a bootstrap approach to provide standard errors of predicted prey proportions, variable importance measures to highlight important variables and partial dependence plots to explore the relationships between explanatory variables and predicted prey composition. Our results supported previous studies of yellowfin tuna feeding ecology in the region. However, the method provided a number of novel insights. For example, significant differences were noted in the prey of yellowfin tuna sampled north of 20°S in summer where oligotrophic waters dominate. The analysis also identified that sea-surface temperature, latitude and yellowfin size were the most important variables associated with dietary differences. The methodology is appropriate for delineating ecosystem-level trophic dynamics, as it can easily incorporate large datasets comprising multiple predators to explore trophic interactions among members of a community. Broad-scale relationships among explanatory variables (environmental, biological, temporal and spatial) and prey composition elucidated by this method then serve to focus and lend validity to subsequent fine-scale analyses of important parameters using standard diet methods and chemical tracers such as stable isotopes.     

Prediction of non-methane hydrocarbons in Kuwait using regression and Bayesian kriged Kalman model

This article describes the hierarchical Bayesian approach for predicting average hourly concentrations of ambient non-methane hydrocarbons (NMHC) in Kuwait where records of six monitor stations located in different sites are observed at successive time points. Our objective is to predict the concentration level of NMHC in unmonitored areas. Here an attempt is made for the prediction of unmeasured concentration of NMHC at two additional locations in Kuwait. We will implement a kriged Kalman filter (KKF) hierarchical Bayesian approach assuming a Gaussian random field, a technique that allows the pooling of data from different sites in order to predict the exposure of the NMHC in different regions of Kuwait. In order to increase the accuracy of the KKF we will use other statistical models such as imputation, regression, principal components, and time series analysis in our approach. We considered four different types of imputation techniques to address the missing data. At the primary level, the logarithmic field is modeled as a trend plus Gaussian stochastic residual model. The trend model depends on hourly meteorological predictors which are common to all sites. The residuals are then modeled using KKF, and the prediction equation is derived conditioned on adjoining hours. On this basis we developed a spatial predictive distribution for these residuals at unmonitored sites. By transforming the predicted residuals back to the original data scales, we can impute Kuwait’s hourly non-methane hydrocarbons field.     

The role of local and landscape level measures of greenness in modelling boreal plant species richness

Measurements of primary productivity and its heterogeneity based on satellite images can provide useful estimates of species richness and distribution patterns. However, species richness at a given site may depend not only on local habitat quality and productivity but also on the characteristics of the surrounding landscape. In this study we investigated whether the predictions of species richness of plant families in northern boreal landscape in Finland can be improved by incorporating greenness information from the surrounding landscape, as derived from remotely sensed data (mean, maximum, standard deviation and range values of NDVI derived from Landsat ETM), into local greenness models. Using plant species richness data of 28 plant families from 440 grid cells of 25 ha in size, generalized additive models (GAMs) were fitted into three different sets of explanatory variables: (1) local greenness only, (2) landscape greenness only, and (3) combined local and landscape greenness. The derived richness–greenness relationships were mainly unimodal or positively increasing but varied between different plant families, and depended also on whether greenness was measured as mean or maximum greenness. Incorporation of landscape level greenness variables improved significantly both the explanatory power and cross-validation statistics of the models including only local greenness variables. Landscape greenness information derived from remote sensing data integrated with local information has thus the potentiality to improve predictive assessments of species richness over extensive and inaccessible areas, especially in high-latitude landscapes. Overall, the significant relationship between plants and surrounding landscape quality detected here suggests that landscape factors should be considered in preserving species richness of boreal environments, as well as in conservation planning for biodiversity in other environments.     

Predictive vegetation modeling for conservation: impact of error propagation from digital elevation data.

The effect of digital elevation model (DEM) error on environmental variables, and subsequently on predictive habitat models, has not been explored. Based on an error analysis of a DEM, multiple error realizations of the DEM were created and used to develop both direct and indirect environmental variables for input to predictive habitat models. The study explores the effects of DEM error and the resultant uncertainty of results on typical steps in the modeling procedure for prediction of vegetation species presence/absence. Results indicate that all of these steps and results, including the statistical significance of environmental variables, shapes of species response curves in generalized additive models (GAMs), stepwise model selection, coefficients and standard errors for generalized linear models (GLMs), prediction accuracy (Cohen's kappa and AUC), and spatial extent of predictions, were greatly affected by this type of error. Error in the DEM can affect the reliability of interpretations of model results and level of accuracy in predictions, as well as the spatial extent of the predictions. We suggest that the sensitivity of DEM-derived environmental variables to error in the DEM should be considered before including them in the modeling processes.     

Quantifying parameter uncertainty in a coral reef model using Metropolis-Coupled Markov Chain Monte Carlo

Coral reefs are threatened ecosystems, so it is important to have predictive models of their dynamics. Most current models of coral reefs fall into two categories. The first is simple heuristic models which provide an abstract understanding of the possible behaviour of reefs in general, but do not describe real reefs. The second is complex simulations whose parameters are obtained from a range of sources such as literature estimates. We cannot estimate the parameters of these models from a single data set, and we have little idea of the uncertainty in their predictions.We have developed a compromise between these two extremes, which is complex enough to describe real reef data, but simple enough that we can estimate parameters for a specific reef from a time series. In previous work, we fitted this model to a long-term data set from Heron Island, Australia, using maximum likelihood methods. To evaluate predictions from this model, we need estimates of the uncertainty in our parameters. Here, we obtain such estimates using Bayesian Metropolis-Coupled Markov Chain Monte Carlo. We do this for versions of the model in which corals are aggregated into a single state variable (the three-state model), and in which corals are separated into four state variables (the six-state model), in order to determine the appropriate level of aggregation. We also estimate the posterior distribution of predicted trajectories in each case.In both cases, the fitted trajectories were close to the observed data, but we had doubts about the biological plausibility of some parameter estimates. We suggest that informative prior distributions incorporating expert knowledge may resolve this problem. In the six-state model, the posterior distribution of state frequencies after 40 years contained two divergent community types, one dominated by free space and soft corals, and one dominated by acroporid, pocilloporid, and massive corals. The three-state model predicts only a single community type. We conclude that the three-state model hides too much biological heterogeneity, but we need more data if we are to obtain reliable predictions from the six-state model. It is likely that there will be similarly large, but currently unevaluated, uncertainty in the predictions of other coral reef models, many of which are much more complex and harder to fit to real data.     

Invasive species distribution modeling (iSDM): Are absence data and dispersal constraints needed to predict actual distributions?

Species distribution models (SDMs) based on statistical relationships between occurrence data and underlying environmental conditions are increasingly used to predict spatial patterns of biological invasions and prioritize locations for early detection and control of invasion outbreaks. However, invasive species distribution models (iSDMs) face special challenges because (i) they typically violate SDM's assumption that the organism is in equilibrium with its environment, and (ii) species absence data are often unavailable or believed to be too difficult to interpret. This often leads researchers to generate pseudo-absences for model training or utilize presence-only methods, and to confuse the distinction between predictions of potential vs. actual distribution. We examined the hypothesis that true-absence data, when accompanied by dispersal constraints, improve prediction accuracy and ecological understanding of iSDMs that aim to predict the actual distribution of biological invasions. We evaluated the impact of presence-only, true-absence and pseudo-absence data on model accuracy using an extensive dataset on the distribution of the invasive forest pathogen Phytophthora ramorum in California. Two traditional presence/absence models (generalized linear model and classification trees) and two alternative presence-only models (ecological niche factor analysis and maximum entropy) were developed based on 890 field plots of pathogen occurrence and several climatic, topographic, host vegetation and dispersal variables. The effects of all three possible types of occurrence data on model performance were evaluated with receiver operating characteristic (ROC) and omission/commission error rates. Results show that prediction of actual distribution was less accurate when we ignored true-absences and dispersal constraints. Presence-only models and models without dispersal information tended to over-predict the actual range of invasions. Models based on pseudo-absence data exhibited similar accuracies as presence-only models but produced spatially less feasible predictions. We suggest that true-absence data are a critical ingredient not only for accurate calibration but also for ecologically meaningful assessment of iSDMs that focus on predictions of actual distributions.     

Sensitivity of species-distribution models to error,bias, and model design: An application to resource selection functions for woodland caribou

Models that predict distribution are now widely used to understand the patterns and processes of plant and animal occurrence as well as to guide conservation and management of rare or threatened species. Application of these methods has led to corresponding studies evaluating the sensitivity of model performance to requisite data and other factors that may lead to imprecise or false inferences. We expand upon these works by providing a relative measure of the sensitivity of model parameters and prediction to common sources of error, bias, and variability. We used a one-at-a-time sample design and GPS location data for woodland caribou (Rangifer tarandus caribou) to assess one common species-distribution model: a resource selection function. Our measures of sensitivity included change in coefficient values, prediction success, and the area of mapped habitats following the systematic introduction of geographic error and bias in occurrence data, thematic misclassification of resource maps, and variation in model design. Results suggested that error, bias and model variation have a large impact on the direct interpretation of coefficients. Prediction success and definition of important habitats were less responsive to the perturbations we introduced to the baseline model. Model coefficients, prediction success, and area of ranked habitats were most sensitive to positional error in species locations followed by sampling bias, misclassification of resources, and variation in model design. We recommend that researchers report, and practitioners consider, levels of error and bias introduced to predictive species-distribution models. Formal sensitivity and uncertainty analyses are the most effective means for evaluating and focusing improvements on input data and considering the range of values possible from imperfect models.     

Fourth-corner correlation is a score test statistic in a log-linear trait–environment model that is useful in permutation testing

Ecologists wish to understand the role of traits of species in determining where each species occurs in the environment. For this, they wish to detect associations between species traits and environmental variables from three data tables, species count data from sites with associated environmental data and species trait data from data bases. These three tables leave a missing part, the fourth-corner. The fourth-corner correlations between quantitative traits and environmental variables, heuristically proposed 20 years ago, fill this corner. Generalized linear (mixed) models have been proposed more recently as a model-based alternative. This paper shows that the squared fourth-corner correlation times the total count is precisely the score test statistic for testing the linear-by-linear interaction in a Poisson log-linear model that also contains species and sites as main effects. For multiple traits and environmental variables, the score test statistic is proportional to the total inertia of a doubly constrained correspondence analysis. When the count data are over-dispersed compared to the Poisson or when there are other deviations from the model such as unobserved traits or environmental variables that interact with the observed ones, the score test statistic does not have the usual chi-square distribution. For these types of deviations, row- and column-based permutation methods (and their sequential combination) are proposed to control the type I error without undue loss of power (unless no deviation is present), as illustrated in a small simulation study. The issues for valid statistical testing are illustrated using the well-known Dutch Dune Meadow data set.