--- title: "Preparing your data" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Preparing your data} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5) has_glmmTMB <- requireNamespace("glmmTMB", quietly = TRUE) ``` ```{r setup} library(isoniche) ``` This vignette covers the decisions to make before estimating niche breadth: the data format, sampling effort, what counts as a replicate, which habitats to include, and what to do with incomplete or small sites. ## Long format All the main functions (`check_coverage()`, `fit_isodars()`, `niche_breadth()`, `group_sites()`) take data in long format, with one row per species x site x habitat and columns named by arguments: | argument | column | required | |---|---|---| | `species` | species name | no (without it the data are one species) | | `site` | replicate (site, or survey) | yes | | `habitat` | habitat | yes | | `abundance` | number of individuals | yes | | `effort` | sampling effort | no (default: 1 for every row) | Rules: * A habitat not sampled at a site has no row. A habitat sampled where the species was not found has a row with abundance 0. If any species has a row for a site and habitat, that site and habitat count as sampled, and missing rows of other species there are treated as 0. * There is at most one row per species, site and habitat; pool repeated samples first (see below). * Sites where a species has no individuals at all are dropped for that species: they carry no information on how it distributes among habitats. * If `habitat` is a factor, its levels set the order of habitats. ## From a community matrix Community data are often stored as a matrix of samples x species, with a separate table of sample attributes. `community_to_long()` converts them, pooling samples from the same site and habitat and summing their effort: ```{r community} comm <- matrix(c(12, 0, 3, 8, 2, 0, 5, 1, 1, 0, 5, 9, 2, 7, 4), nrow = 5, byrow = TRUE, dimnames = list(NULL, c("sp_a", "sp_b", "sp_c"))) samples <- data.frame(site = c("north", "north", "north", "south", "south"), habitat = c("reef", "reef", "lagoon", "reef", "lagoon"), area_m2 = c(100, 50, 100, 100, 100), lon = c(-149.80, -149.80, -149.80, -149.85, -149.85), lat = c(-17.48, -17.48, -17.48, -17.55, -17.55)) long <- community_to_long(comm, samples, site = "site", habitat = "habitat", effort = "area_m2", keep = c("lon", "lat")) long ``` The result goes straight into the other functions, e.g. `niche_breadth(long, species = "species", site = "site", habitat = "habitat", abundance = "abundance", effort = "area_m2")`. Without an `effort` column, each sample has effort 1, so pooled cells have effort equal to their number of samples. ## Effort and the unit of abundance Isodars relate *densities*, so counts are divided by effort: density = abundance / effort. Effort can be area, transect length (with constant width), trap nights or any relative measure, as long as it is the same kind of measure in every row. Densities are then expressed per `reference_effort`, by default the median effort of the data. This is the sampling unit in which "one individual" and all reported abundances (spillover abundance, the x-axis of niche breadth curves) are measured. DI niche breadth is evaluated where the least favoured habitat holds one individual per reference effort, so **DI niche breadth and spillover abundance depend on this unit, while isodar slopes and DD niche breadth do not**. Use the same `reference_effort` for all species you compare, and report it. ```{r effort} fish <- moorea_fish[moorea_fish$species == "Chaetodon ornatissimus", ] fish$effort <- 1 unit1 <- niche_breadth(fish, site = "site", habitat = "habitat", abundance = "abundance", effort = "effort", n_boot = 0, diagnostics = list(n_null = 0)) unit2 <- niche_breadth(fish, site = "site", habitat = "habitat", abundance = "abundance", effort = "effort", reference_effort = 2, n_boot = 0, diagnostics = list(n_null = 0)) data.frame(reference_effort = c(1, 2), rbind(unit1[c("nb_di", "spillover_abundance", "nb_dd")], unit2[c("nb_di", "spillover_abundance", "nb_dd")])) ``` With `reference_effort = 2`, abundances are counted per two sampling units, so "one individual" in the least-favoured habitat is a different point along the isodars: the spillover abundance and DI niche breadth change (here the spillover abundance almost doubles, but the direction and size of the change depend on the isodars), while DD niche breadth is unchanged. The default `reference_effort`, the median effort of the pooled site x habitat cells, depends on how many samples happened to be pooled. Set it to a natural unit instead, for example the area of one standard transect, and use the same value for every species you compare. The unit is printed with the results of `niche_breadth()`. ## Replicates: sites or surveys Isodars are estimated across replicates that differ in abundance. A replicate can be a site, or one survey (a year, a season) of a site; put its identifier in the `site` column. With several sites surveyed repeatedly, each site x survey can be a replicate (combine the two into one identifier), but repeated surveys of a site are not independent, and the bootstrap, which resamples replicates, will then give intervals that are too narrow. Pooling surveys within a site, or analysing years at a single site, avoids this. ## Which habitats Include the habitats *available* to the species: those it could reach. A habitat the species never uses lowers its niche breadth, so species compared with each other should be assessed over the same set of habitats. Rare or poorly sampled habitats can be merged with similar ones before the analysis. ## Absent and undersampled habitats A habitat in which a species was never recorded is treated in one of two ways: * **absent**, when the habitat was sampled at enough sites where the species occurred (`absent_min_sites`, default 5) and the species was found there no more than `absent_max_count` times (default 0). It is fixed at 0, excluded from the isodar fit and counted as 0 in niche breadth. * **undersampled**, when it was sampled at too few such sites. The zeros may then be chance, so the habitat is kept and the species is flagged. ```{r absent} cov <- check_coverage(moorea_fish, species = "species", site = "site", habitat = "habitat", abundance = "abundance") cov[cov$species == "Chromis viridis", ] ``` The rule can be changed (`absent_max_count`, `absent_min_sites`), switched off (`absent = "none"`) or replaced by naming the absent habitats (`absent = "outer_slope"`), in `niche_breadth()` and `fit_isodars()`. With the default `absent_max_count = 0`, a single individual recorded in a habitat makes it used. That one record can move the spillover abundance a long way (the species must now reach one individual in that habitat), so if a few stray individuals should not count, raise `absent_max_count`. A habitat used at low abundance but predicted to be empty at high abundance is *not* absent: it counts for DI niche breadth and is listed in the `abandoned` column of `niche_breadth()`. ## Incomplete sites Not every habitat needs to be sampled at every site. Pairwise isodars (the default, `dim = "pairwise"`) use every site where both habitats of a pair were sampled. The n-dimensional isodar (`dim = "ndim"`) can use sites with at least two habitats through a one-factor model (`ndim_method = "latent"`, needs the **glmmTMB** package); with `ndim_method = "auto"` it is chosen whenever some sites are incomplete. Neither is preferred in general; when many sites are incomplete, fitting both and comparing them is a useful check. `n_sites` in the results counts the sites where the species was recorded; the latent fit uses those with at least two sampled habitats. ```{r latent, eval = has_glmmTMB} set.seed(1) inc <- simulate_isodars(slopes = c(1.5, 1, 0.8), intercepts = c(5, 0, 2), n = 60, p_missing = 0.3, format = "long") head(inc) niche_breadth(inc, site = "site", habitat = "habitat", abundance = "abundance", dim = "ndim", n_boot = 0, diagnostics = list(n_null = 0))[c("method", "n_sites", "nb_di", "nb_dd", "flags")] ``` ## Grouping nearby sites When single sites lack habitats or hold few individuals, nearby sites can be merged into larger sampling units with `group_sites()`. It uses complete linkage on great-circle distance, so every pair of sites in a group lies within the radius; counts and efforts are summed. ```{r group} set.seed(2) d <- simulate_isodars(1.5, 1.5, 5, 5, n = 30, p_missing = 0.3, format = "long") xy <- data.frame(site = unique(d$site), lon = stats::runif(30, 0, 0.1), lat = stats::runif(30, 0, 0.1)) d <- merge(d, xy) g <- group_sites(d, site = "site", lon = "lon", lat = "lat", radius_km = 2, habitat = "habitat", abundance = "abundance") head(g) ``` The radius is not chosen automatically: choosing it to improve the results would be circular. Instead, `radius_sensitivity()` shows how the number of units, their completeness and the estimates change with the radius; choose a radius on ecological grounds (e.g. the scale of movement of the species) and check that the conclusions do not hinge on it: ```{r radius} radius_sensitivity(d, site = "site", lon = "lon", lat = "lat", habitat = "habitat", abundance = "abundance", radii_km = c(0, 2, 5), diagnostics = list(n_null = 0)) ``` `radius_km` can also be a named vector with one radius per species. ## One species in wide format For a single species, data can also be given in wide format, one row per site and one column per habitat (`NA` for habitats not sampled), with the habitat columns named in `hab_cols`: ```{r wide} set.seed(1) sim <- simulate_isodars(1.5, 1.5, 5, 5, noise = 1, n = 30) head(sim) niche_breadth(sim, hab_cols = c("hab1", "hab2", "hab3"), n_boot = 0, diagnostics = list(n_null = 0))[c("nb_di", "nb_dd")] ``` The lower-level functions `iai()` and `isoderr()` work on fitted isodars of one species in this format; see their help pages.