--- title: "Diagnostics: when to trust an estimate" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Diagnostics: when to trust an estimate} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5) ``` ```{r setup} library(isoniche) ``` Isodars are estimated from the spread of a species' densities across replicates, so a species needs enough replicates, enough individuals and a consistent habitat preference for its niche breadth to be estimated well. Many species in a typical community data set will not meet these requirements. `niche_breadth()` never stops on such a species: it returns an estimate where it can, and records each potential problem as a logical `flag_*` column. This vignette explains each flag and what to do about it. ## All species of the Moorea data ```{r all} nb <- niche_breadth(moorea_fish, species = "species", site = "site", habitat = "habitat", abundance = "abundance", n_boot = 0, seed = 1) summary(nb) ``` In a data set of 13 sites, most rare species cannot be analysed. This is expected: the same species would also give unreliable classic niche breadth values, but without any warning. ## The flags Thresholds are arguments of `niche_breadth()`, set through `diagnostics` (e.g. `diagnostics = list(min_sites = 8)`); the defaults are given below. ### Data * **`few_sites`**: fewer than `min_sites` (10) replicates with the species. Isodars from few points are imprecise. Group nearby sites (`group_sites()`; see `vignette("data-preparation")`) or treat the species as unreliable. * **`few_pair_sites`**: some pair of used habitats was sampled together at fewer than `min_pair_sites` (5) sites, so its pairwise isodar is poorly estimated. Use the n-dimensional isodar with `ndim_method = "latent"`, or group sites. * **`low_abundance`**: the median site total is below `min_median_count` (5) individuals, or more than `max_zero_share` (half) of the counts in used habitats are zeros. Densities are then dominated by sampling noise. * **`undersampled`**: a habitat without records was sampled at too few sites to call it absent. Sample more sites, or decide on ecological grounds and name absent habitats with `absent`. ### Fit * **`fit_failed`**: the isodars could not be fitted; the reason is in the `note` column. Most often the species has too few sites. * **`disconnected`**: some habitats never occur together at a site, so their relationship cannot be estimated. * **`poor_fit`**: for the n-dimensional isodar, the line explains less than `min_var_explained` (50%) of the variance; for pairwise isodars, the isodars of three habitats contradict each other (the product of two slopes differs from the third by more than a factor `max_slope_ratio`, 2, or in sign). Habitat selection may not follow a single isodar, or the data are too noisy. Look at the isodars with `fit_isodars()`, and try the other method (`dim = "ndim"` or `"pairwise"`). * **`no_di`**: no total abundance up to `max_search` puts an individual in every used habitat: the least-favoured habitat is predicted never to fill. * **`intercept_noise`**: the isodar intercepts do not differ from those expected by sampling alone (see below). The data then give no evidence that niche breadth changes with abundance: with zero intercepts, the isodars pass through the origin and DI and DD niche breadth are equal. Report a single value, DD niche breadth, which rests on the slopes only and is the better estimated. With few sites the test has little power, so this flag does not show that the intercepts are zero. ### Extrapolation * **`di_extrapolation`**: the spillover abundance lies outside the observed range of site totals by more than a factor `extrapolation_ratio` (2). * **`dd_extrapolation`**: the largest site total is below twice the spillover abundance, so the data never reach the halfway point of the curve towards DD niche breadth. Both are extrapolations of the fitted isodars; the estimates may be right, but rest on the assumption that the isodars hold beyond the data. ### Uncertainty * **`unstable_bootstrap`** (only with `n_boot > 0`): more than `max_failed_share` (10%) of bootstrap replicates failed, or a confidence interval is wider than `max_ci_width` (0.5). Increase `n_boot`; if the interval stays wide, the data cannot resolve the species' niche breadth. ### Informational * **`reversal`**: two habitats swap rank as abundance increases; the abundances are in the `reversals` column. Niche breadth can then change non-monotonically with abundance, so compare species at a stated abundance. * **`abandoned`**: a habitat used at low abundance is predicted to be empty at high abundance; it counts for DI but not DD niche breadth. ## Is the intercept more than noise? When a species is rare, counts in its habitats are small and their ratios noisy, which alone produces isodar intercepts. `niche_breadth()` tests whether the observed intercepts exceed this: it simulates `n_null` (199) data sets in which each site's individuals are redistributed among its habitats at random, in proportion to the species' overall density in each habitat (no density dependence), refits the isodars, and compares the observed intercepts with the simulated ones. `intercept_p` is the smallest Holm-corrected p-value; the species is flagged when it exceeds `alpha` (0.05). The test needs integer counts. ```{r noise} noisy <- nb[nb$flag_intercept_noise & !nb$flag_few_sites, ] noisy[c("species", "nb_di", "nb_dd", "intercept_p")] ``` ## Looking at a flagged species Diagnostics are a starting point. For a flagged species, look at its isodars and curve: ```{r look} sp <- "Acanthurus nigrofuscus" nb[nb$species == sp, c("nb_di", "nb_dd", "var_explained", "flags")] fit <- fit_isodars(moorea_fish[moorea_fish$species == sp, ], site = "site", habitat = "habitat", abundance = "abundance") fit$fit[[1]] ``` ```{r look-curve, fig.height = 5} plot(iai_curve(fit)) nb$reversals[[which(nb$species == sp)]] ``` When rare, this species lives on the fringing and barrier reefs: both isodars against the outer slope have large intercepts. But the outer slope fills fastest as abundance rises (the reefs gain only 0.27 and 0.08 individuals per individual on the outer slope), so it overtakes both reefs at a few hundred individuals, the ranking of habitats changes and niche breadth first rises, then falls. Here the pairwise isodars are also somewhat inconsistent with each other (`poor_fit`). Neither DI nor DD niche breadth alone describes such a species well; the curve does. ## A workflow 1. Run `check_coverage()` and drop species or habitats that are clearly undersampled. 2. Run `niche_breadth()` with `n_boot = 0` on all species; look at `summary()`. 3. Decide how to handle each warning flag (for example: drop species with `fit_failed`, `few_sites` or `low_abundance`; report DI and DD niche breadth separately only where `intercept_noise` is not raised). 4. Rerun on the retained species with `n_boot = 999` or more, and check `unstable_bootstrap`. 5. Report the thresholds and how many species each rule removed.