| Title: | Niche Region and Niche Overlap Metrics for Multidimensional Ecological Niches |
|---|---|
| Description: | Implementation of a probabilistic method to calculate 'nicheROVER' (_niche_ _r_egion and niche _over_lap) metrics using multidimensional niche indicator data (e.g., stable isotopes, environmental variables, etc.). The niche region is defined as the joint probability density function of the multidimensional niche indicators at a user-defined probability alpha (e.g., 95%). Uncertainty is accounted for in a Bayesian framework, and the method can be extended to three or more indicator dimensions. It provides directional estimates of niche overlap, accounts for species-specific distributions in multivariate niche space, and produces unique and consistent bivariate projections of the multivariate niche region. The article by Swanson et al. (2015) <doi:10.1890/14-0235.1> provides a detailed description of the methodology. See the package vignette for a worked example using fish stable isotope data. |
| Authors: | Martin Lysy [aut, cre], Ashley D. Stasko [aut], Heidi K. Swanson [aut] |
| Maintainer: | Martin Lysy <[email protected]> |
| License: | GPL-3 |
| Version: | 1.1.2 |
| Built: | 2024-11-19 05:40:44 UTC |
| Source: | https://github.com/mlysy/nicherover |
This package uses a probabilistic method to calculate niche regions and pairwise niche overlap using multidimensional niche indicator data (e.g., stable isotopes, environmental variables, etc.). The niche region is defined as the joint probability density function of the multidimensional niche indicators at a user-defined probability alpha (e.g., 95%). Uncertainty is accounted for in a Bayesian framework, and the method can be extended to three or more indicator dimensions. It provides directional estimates of niche overlap, accounts for species-specific distributions in multivariate niche space, and produces unique and consistent bivariate projections of the multivariate niche region. See Swanson et al. (2014) for a detailed description and worked example below using fish stable isotope data.
Maintainer: Martin Lysy [email protected]
Authors:
Ashley D. Stasko
Heidi K. Swanson
Swanson, H.K., Lysy, M., Stasko, A.D., Power, M., Johnson, J.D., and Reist, J.D. "A new probabilistic method for quantifying n-dimensional ecological niches and niche overlap." Ecology: Statistical Reports 96:2 (2015): 318-324. doi:10.1890/14-0235.1.
Useful links:
# analysis for fish data data(fish) # 4 fish, 3 isotopes aggregate(fish[2:4], fish[1], mean) # isotope means per fish # random draws from posterior distribution with default prior nsamples <- 500 fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) # display p(mu | X) and p(Sigma | X) for each fish clrs <- c("black", "red", "blue", "orange") # colors for each species par(mar = c(4.2, 4.2, 2, 1)+.1) niche.par.plot(fish.par, col = clrs) legend(x = "topright", legend = names(fish.par), fill = clrs) # 2-d projections of 10 niche regions nsamples <- 10 fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) # format data for plotting function fish.data <- tapply(1:nrow(fish), fish$species, function(ii) X = fish[ii,2:4]) niche.plot(niche.par = fish.par, niche.data = fish.data, pfrac = .05, iso.names = expression(delta^{15}*N, delta^{13}*C, delta^{34}*S), col = clrs, xlab = expression("Isotope Ratio (per mille)")) # niche overlap plots for 95% niche region sizes # overlap calculation. use nsamples = nprob = 1e4 for higher accuracy. nsamples <- 500 over.stat <- overlap(fish.par, nreps = nsamples, nprob = nsamples, alpha = .95) # overlap plot overlap.plot(over.stat, col = clrs, mean.cred.col = "turquoise", equal.axis = TRUE, xlab = "Overlap Probability (%) -- Niche Region Size: 95%")# analysis for fish data data(fish) # 4 fish, 3 isotopes aggregate(fish[2:4], fish[1], mean) # isotope means per fish # random draws from posterior distribution with default prior nsamples <- 500 fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) # display p(mu | X) and p(Sigma | X) for each fish clrs <- c("black", "red", "blue", "orange") # colors for each species par(mar = c(4.2, 4.2, 2, 1)+.1) niche.par.plot(fish.par, col = clrs) legend(x = "topright", legend = names(fish.par), fill = clrs) # 2-d projections of 10 niche regions nsamples <- 10 fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) # format data for plotting function fish.data <- tapply(1:nrow(fish), fish$species, function(ii) X = fish[ii,2:4]) niche.plot(niche.par = fish.par, niche.data = fish.data, pfrac = .05, iso.names = expression(delta^{15}*N, delta^{13}*C, delta^{34}*S), col = clrs, xlab = expression("Isotope Ratio (per mille)")) # niche overlap plots for 95% niche region sizes # overlap calculation. use nsamples = nprob = 1e4 for higher accuracy. nsamples <- 500 over.stat <- overlap(fish.par, nreps = nsamples, nprob = nsamples, alpha = .95) # overlap plot overlap.plot(over.stat, col = clrs, mean.cred.col = "turquoise", equal.axis = TRUE, xlab = "Overlap Probability (%) -- Niche Region Size: 95%")
Calculates coordinates of points for plotting a 2-dimensional ellipse based on user-defined parameters. Can be used for exploratory data analysis to produce ellipses at a given niche region size (e.g., ).
ellipse(mu, V, alpha = 0.95, n = 100)ellipse(mu, V, alpha = 0.95, n = 100)
mu |
Centre of ellipse. A vector of length 2. |
V |
Scale of ellipse. A 2x2 matrix. See 'Details'. |
alpha |
Niche region size. See 'Details'. |
n |
Number of points to return for plotting. |
This function provides the coordinates needed to plot a 2-dimensional ellipse based on user-defined parameters, such that X = c(x,y) satisfies the equation
where .
Returns a matrix of coordinates cbind(x,y) to plot a 2-dimensional ellipse.
niche.plot() for plotting.
mu <- rnorm(2) V <- crossprod(matrix(rnorm(4), 2, 2)) ell.pts <- ellipse(mu = mu, V = V, alpha = .9, n = 100) plot(ell.pts, col = rainbow(110)[1:100], type = "o") points(mu[1], mu[2], pch = "+")mu <- rnorm(2) V <- crossprod(matrix(rnorm(4), 2, 2)) ell.pts <- ellipse(mu = mu, V = V, alpha = .9, n = 100) plot(ell.pts, col = rainbow(110)[1:100], type = "o") points(mu[1], mu[2], pch = "+")
A dataset containing values for three stable isotopes measured in the muscle tissue of four species of arctic fish. For use in examples.
fishfish
A data frame with 278 rows (observations) and 4 columns (species, , , and ).
This dataset contains , , and values for the following fish species:
ARCS - Arctic Cisco (Coregonus autumnalis), .
BDWF - Broad Whitefish (Coregonus nasus), .
LKWF - Lake Whitefish (Coregonus clupeaformis), .
LSCS - Least Cisco (Coregonus sardinella),
Fish were collected between 2007 and 2008 from an estuarine area of the Beaufort Sea, North and West of the Mackenzie Delta at Phillips Bay, Yukon Territory, Canada (69.28 N, 138.49 W).
head(fish) aggregate(fish[2:4], fish[1], mean)head(fish) aggregate(fish[2:4], fish[1], mean)
For one or more species, plots some or all of the niche parameters and .
niche.par.plot( niche.par, plot.mu = TRUE, plot.Sigma = TRUE, plot.index, col, ndens = 512, ylab )niche.par.plot( niche.par, plot.mu = TRUE, plot.Sigma = TRUE, plot.index, col, ndens = 512, ylab )
niche.par |
List with |
plot.mu |
Logical. If |
plot.Sigma |
Logical. If |
plot.index |
Either a scalar of a numeric vector of length 2. If |
col |
Vector of colors in which to plot each species. |
ndens |
Number of points at which to evaluate density estimates. |
ylab |
Optional label for |
Each element of the list niche.par is a distribution of niche parameters. That is, names(niche.par[[1]]) = c("mu", "Sigma"), and if niso is the number of niche indicators (e.g., stable isotopes), then dim(niche.par[[1]]$mu) = c(nsamples, niso) and dim(niche.par[[1]]$Sigma) = c(niso, niso, nsamples).
Returns a plot of the distribution of some or all niche parameters.
niw.post(), niiw.post() for niche parameter output, stats::density() for density estimation from sample data.
# fish data data(fish) # generate parameter draws from the "default" posteriors of each fish nsamples <- 1e3 system.time({ fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) }) # various parameter plots clrs <- c("black", "red", "blue", "orange") # colors for each species # mu1, mu2, and Sigma12 par(mar = c(4, 4, .5, .1)+.1, mfrow = c(1,3)) niche.par.plot(fish.par, col = clrs, plot.index = 1) niche.par.plot(fish.par, col = clrs, plot.index = 2) niche.par.plot(fish.par, col = clrs, plot.index = 1:2) legend("topright", legend = names(fish.par), fill = clrs) # all mu niche.par.plot(fish.par, col = clrs, plot.mu = TRUE, plot.Sigma = FALSE) legend("topright", legend = names(fish.par), fill = clrs) # all mu and Sigma par(mar = c(4.2, 4.2, 2, 1)+.1) niche.par.plot(fish.par, col = clrs, plot.mu = TRUE, plot.Sigma = TRUE) legend("topright", legend = names(fish.par), fill = clrs)# fish data data(fish) # generate parameter draws from the "default" posteriors of each fish nsamples <- 1e3 system.time({ fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) }) # various parameter plots clrs <- c("black", "red", "blue", "orange") # colors for each species # mu1, mu2, and Sigma12 par(mar = c(4, 4, .5, .1)+.1, mfrow = c(1,3)) niche.par.plot(fish.par, col = clrs, plot.index = 1) niche.par.plot(fish.par, col = clrs, plot.index = 2) niche.par.plot(fish.par, col = clrs, plot.index = 1:2) legend("topright", legend = names(fish.par), fill = clrs) # all mu niche.par.plot(fish.par, col = clrs, plot.mu = TRUE, plot.Sigma = FALSE) legend("topright", legend = names(fish.par), fill = clrs) # all mu and Sigma par(mar = c(4.2, 4.2, 2, 1)+.1) niche.par.plot(fish.par, col = clrs, plot.mu = TRUE, plot.Sigma = TRUE) legend("topright", legend = names(fish.par), fill = clrs)
For one or more species, creates a series of plots: (i) the raw niche indicators (e.g., stable isotope) data, (ii) their density estimates, and (iii) 2-dimensional projections of probabilistic niche regions based on -dimensionsional data.
niche.plot( niche.par, niche.data, alpha = 0.95, species.names, iso.names, lims, col, ndens = 512, pfrac = 0, xlab, legend = TRUE )niche.plot( niche.par, niche.data, alpha = 0.95, species.names, iso.names, lims, col, ndens = 512, pfrac = 0, xlab, legend = TRUE )
niche.par |
A list of length |
niche.data |
A list of length |
alpha |
Size of the niche region to plot. Defaults to 0.95. |
species.names |
Names of the species. Defaults to |
iso.names |
Names of the niche indicators (or isotopes) to plot. Defaults to |
lims |
Two-row matrix of range limits for each niche indicator. Defaults to include all ellipses. |
col |
Vector of colours in which each species will be drawn. |
ndens |
Number of points at which to evaluate kernel density estimates. |
pfrac |
Fraction of the plot on which to display 1-dimensional raw niche indicator data. |
xlab |
Title of plot, located on the bottom. Defaults to no title. |
legend |
Whether or not to add a legend. Defaults to |
A set of plots is created for each pairwise combination of niche indicators. Below the diagonal are scatterplots for each species, above the diagonal are ellipses corresponding to 2-d projections of the probabilistic niche regions. The diagonal displays density estimates for each indicator, and optionally the raw 1-d data. See Swanson et al. (2015) for detailed description of the probabilistic niche region.
Returns a series of plots displaying niche indicator data and their probabilistic niche projections.
Swanson, H.K., Lysy, M., Stasko, A.D., Power, M., Johnson, J.D., and Reist, J.D. "A new probabilistic method for quantifying n-dimensional ecological niches and niche overlap." Ecology: Statistical Reports 96:2 (2015): 318-324. doi:10.1890/14-0235.1.
overlap.plot(), niw.post(), niiw.post().
data(fish) # 4 fish, 3 isotopes # generate 10 parameter draws from the posteriors # of each fish with default prior nsamples <- 10 fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) # format data for plotting function fish.data <- tapply(1:nrow(fish), fish$species, function(ii) X = fish[ii,2:4]) clrs <- c("black", "red", "blue", "orange") # colors for each species niche.plot(niche.par = fish.par, niche.data = fish.data, pfrac = .1, iso.names = expression(delta^{15}*N, delta^{13}*C, delta^{34}*S), col = clrs, xlab = expression("Isotope Ratio (per mille)"))data(fish) # 4 fish, 3 isotopes # generate 10 parameter draws from the posteriors # of each fish with default prior nsamples <- 10 fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) # format data for plotting function fish.data <- tapply(1:nrow(fish), fish$species, function(ii) X = fish[ii,2:4]) clrs <- c("black", "red", "blue", "orange") # colors for each species niche.plot(niche.par = fish.par, niche.data = fish.data, pfrac = .1, iso.names = expression(delta^{15}*N, delta^{13}*C, delta^{34}*S), col = clrs, xlab = expression("Isotope Ratio (per mille)"))
Uniform sampling from an elliptical niche region.
niche.runif(n, mu, Sigma, alpha = 0.95)niche.runif(n, mu, Sigma, alpha = 0.95)
n |
Number of random draws. |
mu |
Mean vector. |
Sigma |
Variance matrix. |
alpha |
Probabilistic niche size |
IID draws from a uniform distribution on the elliptical niche region.
ellipse() and niche.size() for the definition of the elliptical niche region.
# 2d example d <- 2 # number of dimensions V <- crossprod(matrix(rnorm(4),d,d)) mu <- rnorm(d) plot(ellipse(mu, V), type = "l") points(niche.runif(1e4, mu, V), col = "brown", pch = ".")# 2d example d <- 2 # number of dimensions V <- crossprod(matrix(rnorm(4),d,d)) mu <- rnorm(d) plot(ellipse(mu, V), type = "l") points(niche.runif(1e4, mu, V), col = "brown", pch = ".")
Calculate the size of an elliptical niche region.
niche.size(Sigma, alpha = 0.95)niche.size(Sigma, alpha = 0.95)
Sigma |
Variance matrix for normally distributed niche axes. |
alpha |
Probabilistic niche size. |
For a given niche region , the niche size is defined as the hypervolume of this region: .
Hypervolume niche size.
# for each species, size of 95% niche region using sample variance tapply(1:nrow(fish), fish$species, function(ind) { X <- fish[ind,2:4] # all measurements for given species Sighat <- var(X) # sample variance niche.size(Sigma = Sighat, alpha = .95) })# for each species, size of 95% niche region using sample variance tapply(1:nrow(fish), fish$species, function(ind) { X <- fish[ind,2:4] # all measurements for given species Sighat <- var(X) # sample variance niche.size(Sigma = Sighat, alpha = .95) })
Given iid -dimensional niche indicators with , this function generates random draws from for the Normal-Independent-Inverse-Wishart (NIIW) prior.
niiw.post(nsamples, X, lambda, Omega, Psi, nu, mu0 = lambda, burn)niiw.post(nsamples, X, lambda, Omega, Psi, nu, mu0 = lambda, burn)
nsamples |
The number of posterior draws. |
X |
A data matrix with observations along the rows. |
lambda |
Mean of |
Omega |
Variance of |
Psi |
Scale matrix of |
nu |
Degrees of freedom of |
mu0 |
Initial value of |
burn |
Burn-in for the MCMC sampling algorithm. Either an integer giving the number of initial samples to discard, or a fraction with |
The NIIW distribution is defined as
The default value Omega = 0 uses the Lebesque prior on : . In this case the NIW and NIIW priors produce identical resuls, but niw.post() is faster.
The default value Psi = 0 uses the scale-invariant prior on : .
The default value nu = ncol(X)+1 for Omega = 0 and Psi = 0 makes and .
Random draws are obtained by a Markov chain Monte Carlo (MCMC) algorithm; specifically, a Gibbs sampler alternates between draws from and , which are Normal and Inverse-Wishart distributions respectively.
Returns a list with elements mu and Sigma of sizes c(nsamples, length(lambda)) and c(dim(Psi), nsamples).
# simulate normal data with mean and variance (mu0, Sigma0) d <- 4 mu0 <- rnorm(d) Sigma0 <- matrix(rnorm(d^2), d, d) Sigma0 <- Sigma0 %*% t(Sigma0) N <- 1e2 X <- matrix(rnorm(N*d), N, d) # iid N(0,1) X <- t(t(X %*% chol(Sigma0)) + mu0) # each row is N(mu0, Sigma) # prior parameters # flat prior on mu lambda <- 0 Omega <- 0 # informative prior on Sigma Psi <- crossprod(matrix(rnorm(d^2), d, d)) nu <- 5 # sample from NIIW posterior nsamples <- 2e3 system.time({ siiw <- niiw.post(nsamples, X, lambda, Omega, Psi, nu, burn = 100) }) # sample from NIW posterior kappa <- 0 system.time({ siw <- niw.post(nsamples, X, lambda, kappa, Psi, nu) }) # check that posteriors are the same # p(mu | X) clrs <- c("black", "red") par(mar = c(4.2, 4.2, 2, 1)+.1) niche.par.plot(list(siiw, siw), col = clrs, plot.mu = TRUE, plot.Sigma = FALSE) legend(x = "topright", legend = c("NIIW Prior", "NIW Prior"), fill = clrs) # p(Sigma | X) par(mar = c(4.2, 4.2, 2, 1)+.1) niche.par.plot(list(siiw, siw), col = clrs, plot.mu = FALSE, plot.Sigma = TRUE) legend(x = "topright", legend = c("NIIW Prior", "NIW Prior"), fill = clrs)# simulate normal data with mean and variance (mu0, Sigma0) d <- 4 mu0 <- rnorm(d) Sigma0 <- matrix(rnorm(d^2), d, d) Sigma0 <- Sigma0 %*% t(Sigma0) N <- 1e2 X <- matrix(rnorm(N*d), N, d) # iid N(0,1) X <- t(t(X %*% chol(Sigma0)) + mu0) # each row is N(mu0, Sigma) # prior parameters # flat prior on mu lambda <- 0 Omega <- 0 # informative prior on Sigma Psi <- crossprod(matrix(rnorm(d^2), d, d)) nu <- 5 # sample from NIIW posterior nsamples <- 2e3 system.time({ siiw <- niiw.post(nsamples, X, lambda, Omega, Psi, nu, burn = 100) }) # sample from NIW posterior kappa <- 0 system.time({ siw <- niw.post(nsamples, X, lambda, kappa, Psi, nu) }) # check that posteriors are the same # p(mu | X) clrs <- c("black", "red") par(mar = c(4.2, 4.2, 2, 1)+.1) niche.par.plot(list(siiw, siw), col = clrs, plot.mu = TRUE, plot.Sigma = FALSE) legend(x = "topright", legend = c("NIIW Prior", "NIW Prior"), fill = clrs) # p(Sigma | X) par(mar = c(4.2, 4.2, 2, 1)+.1) niche.par.plot(list(siiw, siw), col = clrs, plot.mu = FALSE, plot.Sigma = TRUE) legend(x = "topright", legend = c("NIIW Prior", "NIW Prior"), fill = clrs)
Given iid -dimensional niche indicators with , this function calculates the coefficients of the Normal-Inverse-Wishart (NIW) posterior for a conjugate NIW prior. Together with niw.mom(), this can be used to rapidly compute the point estimates and .
niw.coeffs(X, lambda, kappa, Psi, nu)niw.coeffs(X, lambda, kappa, Psi, nu)
X |
A data matrix with observations along the rows. |
lambda |
Location parameter. See 'Details'. |
kappa |
Scale parameter. Defaults to |
Psi |
Scale matrix. Defaults to |
nu |
Degrees of freedom. Defaults to |
The NIW distribution is defined as
The default value kappa = 0 uses the Lebesque prior on : .
The default value Psi = 0 uses the scale-invariant prior on : .
The default value nu = ncol(X)+1 for kappa = 0 and Psi = 0 makes and .
Returns a list with elements lambda, kappa, Psi, nu corresponding to the coefficients of the NIW posterior distribution .
rniw(), niw.mom(), niw.post().
# NIW prior coefficients d <- 3 lambda <- rnorm(d) kappa <- 5 Psi <- crossprod(matrix(rnorm(d^2), d, d)) nu <- 10 # data data(fish) X <- fish[fish$species == "ARCS",2:4] # NIW posterior coefficients post.coef <- niw.coeffs(X, lambda, kappa, Psi, nu) # compare mu.mean <- niw.mom(post.coef$lambda, post.coef$kappa, post.coef$Psi, post.coef$nu)$mu$mean mu.est <- rbind(prior = niw.mom(lambda, kappa, Psi, nu)$mu$mean, data = colMeans(X), post = mu.mean) round(mu.est, 2)# NIW prior coefficients d <- 3 lambda <- rnorm(d) kappa <- 5 Psi <- crossprod(matrix(rnorm(d^2), d, d)) nu <- 10 # data data(fish) X <- fish[fish$species == "ARCS",2:4] # NIW posterior coefficients post.coef <- niw.coeffs(X, lambda, kappa, Psi, nu) # compare mu.mean <- niw.mom(post.coef$lambda, post.coef$kappa, post.coef$Psi, post.coef$nu)$mu$mean mu.est <- rbind(prior = niw.mom(lambda, kappa, Psi, nu)$mu$mean, data = colMeans(X), post = mu.mean) round(mu.est, 2)
This function computes the mean and variance of the Normal-Inverse-Wishart (NIW) distribution. Can be used to very quickly compute Bayesian point estimates for the conjugate NIW prior.
niw.mom(lambda, kappa, Psi, nu)niw.mom(lambda, kappa, Psi, nu)
lambda |
Location parameter. See 'Details'. |
kappa |
Scale parameter. See 'Details'. |
Psi |
Scale matrix. See 'Details'. |
nu |
Degrees of freedom. See 'Details'. |
The NIW distribution is defined as
Note that cov.
Returns a list with elements mu and Sigma, each containing lists with elements mean and var. For mu these elements are of size length(lambda) and c(length(lambda),length(lambda)). For Sigma they are of size dim(Psi) and c(dim(Psi), dim(Psi)), such that covSigma$var[i,j,k,l].
rniw(), niw.coeffs(), niw.post().
# NIW parameters d <- 3 # number of dimensions lambda <- rnorm(d) kappa <- 2 Psi <- crossprod(matrix(rnorm(d^2), d, d)) nu <- 10 # simulate data niw.sim <- rniw(n = 1e4, lambda, kappa, Psi, nu) # check moments niw.mV <- niw.mom(lambda, kappa, Psi, nu) # mean of mu ii <- 1 c(true = niw.mV$mu$mean[ii], sim = mean(niw.sim$mu[,ii])) # variance of Sigma II <- c(1,2) JJ <- c(2,3) c(true = niw.mV$var[II[1],II[2],JJ[1],JJ[2]], sim = cov(niw.sim$Sigma[II[1],II[2],], niw.sim$Sigma[JJ[1],JJ[2],]))# NIW parameters d <- 3 # number of dimensions lambda <- rnorm(d) kappa <- 2 Psi <- crossprod(matrix(rnorm(d^2), d, d)) nu <- 10 # simulate data niw.sim <- rniw(n = 1e4, lambda, kappa, Psi, nu) # check moments niw.mV <- niw.mom(lambda, kappa, Psi, nu) # mean of mu ii <- 1 c(true = niw.mV$mu$mean[ii], sim = mean(niw.sim$mu[,ii])) # variance of Sigma II <- c(1,2) JJ <- c(2,3) c(true = niw.mV$var[II[1],II[2],JJ[1],JJ[2]], sim = cov(niw.sim$Sigma[II[1],II[2],], niw.sim$Sigma[JJ[1],JJ[2],]))
Given iid -dimensional niche indicators with , this function generates random draws from for the Normal-Inverse-Wishart (NIW) prior.
niw.post(nsamples, X, lambda, kappa, Psi, nu)niw.post(nsamples, X, lambda, kappa, Psi, nu)
nsamples |
The number of posterior draws. |
X |
A data matrix with observations along the rows. |
lambda |
Location parameter. See 'Details'. |
kappa |
Scale parameter. Defaults to |
Psi |
Scale matrix. Defaults to |
nu |
Degrees of freedom. Defaults to |
The NIW distribution is defined as
The default value kappa = 0 uses the Lebesque prior on : .
The default value Psi = 0 uses the scale-invariant prior on : .
The default value nu = ncol(X)+1 for kappa = 0 and Psi = 0 makes and .
Returns a list with elements mu and Sigma of sizes c(nsamples, length(lambda)) and c(dim(Psi), nsamples).
# compare the default non-informative prior to an arbitrary informative prior # for simulated data # simulate normal data with mean and variance (mu0, Sigma0) d <- 4 mu0 <- rnorm(d) Sigma0 <- matrix(rnorm(d^2), d, d) Sigma0 <- Sigma0 %*% t(Sigma0) N <- 1e2 X <- matrix(rnorm(N*d), N, d) # iid N(0,1) X <- t(t(X %*% chol(Sigma0)) + mu0) # each row is N(mu0, Sigma) # informative prior parameters lambda <- rnorm(d) kappa <- 20 Psi <- crossprod(matrix(rnorm(d^2), d, d)) nu <- 5 # iid draws from informative prior pi(mu, Sigma) nsamples <- 2e3 siw0 <- rniw(nsamples, lambda, kappa, Psi, nu) # iid draws from posterior p(mu, Sigma | X) with informative prior siw1 <- niw.post(nsamples, X, lambda, kappa, Psi, nu) # iid draws from posterior p(mu, Sigma | X) with default noninformative prior siw2 <- niw.post(nsamples, X) # compare # prior and posterior densities of mu clrs <- c("orange", "red", "blue", "black") ii <- 1 par(mar = c(4.2, 4.2, 2, 1)+.1) niche.par.plot(list(siw0, siw1, siw2), col = clrs[1:3], plot.index = ii, ylab = "Density") abline(v = mu0[ii], col = clrs[4]) # true value of mu legend(x = "topright", legend = c(parse(text = paste0("pi(mu[", ii, "])")), parse(text = paste0("p(mu[", ii, "]*\" | \"*X)*\", Informative Prior\"")), parse(text = paste0("p(mu[", ii, "]*\" | \"*X)*\", Noninformative Prior\"")), parse(text = paste0("\"True value of \"*mu[", ii, "]"))), fill = clrs) # prior and posterior densities of Sigma ii <- 1 jj <- 2 par(mar = c(4.2, 4.2, 2, 1)+.1) niche.par.plot(list(siw0, siw1, siw2), col = clrs[1:3], plot.index = c(ii,jj), ylab = "Density") abline(v = Sigma0[ii,jj], col = clrs[4]) legend(x = "topright", legend = c(parse(text = paste0("pi(Sigma[", ii, "*", jj, "])")), parse(text = paste0("p(Sigma[", ii, "*", jj, "]*\" | \"*X)*\", Informative Prior\"")), parse(text = paste0("p(Sigma[", ii, "*", jj, "]*\" | \"*X)*\", Noninformative Prior\"")), parse(text = paste0("\"True value of \"*Sigma[", ii, "*", jj, "]"))), fill = clrs)# compare the default non-informative prior to an arbitrary informative prior # for simulated data # simulate normal data with mean and variance (mu0, Sigma0) d <- 4 mu0 <- rnorm(d) Sigma0 <- matrix(rnorm(d^2), d, d) Sigma0 <- Sigma0 %*% t(Sigma0) N <- 1e2 X <- matrix(rnorm(N*d), N, d) # iid N(0,1) X <- t(t(X %*% chol(Sigma0)) + mu0) # each row is N(mu0, Sigma) # informative prior parameters lambda <- rnorm(d) kappa <- 20 Psi <- crossprod(matrix(rnorm(d^2), d, d)) nu <- 5 # iid draws from informative prior pi(mu, Sigma) nsamples <- 2e3 siw0 <- rniw(nsamples, lambda, kappa, Psi, nu) # iid draws from posterior p(mu, Sigma | X) with informative prior siw1 <- niw.post(nsamples, X, lambda, kappa, Psi, nu) # iid draws from posterior p(mu, Sigma | X) with default noninformative prior siw2 <- niw.post(nsamples, X) # compare # prior and posterior densities of mu clrs <- c("orange", "red", "blue", "black") ii <- 1 par(mar = c(4.2, 4.2, 2, 1)+.1) niche.par.plot(list(siw0, siw1, siw2), col = clrs[1:3], plot.index = ii, ylab = "Density") abline(v = mu0[ii], col = clrs[4]) # true value of mu legend(x = "topright", legend = c(parse(text = paste0("pi(mu[", ii, "])")), parse(text = paste0("p(mu[", ii, "]*\" | \"*X)*\", Informative Prior\"")), parse(text = paste0("p(mu[", ii, "]*\" | \"*X)*\", Noninformative Prior\"")), parse(text = paste0("\"True value of \"*mu[", ii, "]"))), fill = clrs) # prior and posterior densities of Sigma ii <- 1 jj <- 2 par(mar = c(4.2, 4.2, 2, 1)+.1) niche.par.plot(list(siw0, siw1, siw2), col = clrs[1:3], plot.index = c(ii,jj), ylab = "Density") abline(v = Sigma0[ii,jj], col = clrs[4]) legend(x = "topright", legend = c(parse(text = paste0("pi(Sigma[", ii, "*", jj, "])")), parse(text = paste0("p(Sigma[", ii, "*", jj, "]*\" | \"*X)*\", Informative Prior\"")), parse(text = paste0("p(Sigma[", ii, "*", jj, "]*\" | \"*X)*\", Noninformative Prior\"")), parse(text = paste0("\"True value of \"*Sigma[", ii, "*", jj, "]"))), fill = clrs)
Calculates the distribution of a niche region overlap metric for each pairwise species combination and user-specified niche region sizes.
overlap( niche.par, nreps, nprob, alpha = 0.95, species.names, norm.redraw = TRUE )overlap( niche.par, nreps, nprob, alpha = 0.95, species.names, norm.redraw = TRUE )
niche.par |
A list with |
nreps |
The number of overlap metric calculations for each species. Defaults to the smallest number of parameter samples supplied by |
nprob |
The number of normal draws for each Monte Carlo overlap metric calculation. See 'Details'. |
alpha |
Scalar or vector of niche region sizes for calculating the niche overlap metric. Defaults to 0.95. |
species.names |
Names of the species. Defaults to |
norm.redraw |
Logical. If |
The overlap metric is the probability that a randomly drawn individual from species will be found within the niche region of species (for a given niche region size, e.g., alpha = .95). It is a single number which is a function of the parameters for each species, and . This number is difficult to calculate directly, but easy to approximate stochastically by generating nprob draws from the distribution of species and counting the fraction of them which fall in the niche region of species .
Typically the true values of and are unknown and must be estimated from the data. Thus, the overlap metric is calculated for nreps combinations of samples from and which are supplied in niche.par.
See Swanson et al. (2015) for a detailed description of niche overlap and its calculation.
Returns an array of size c(nspecies, nspecies, nreps, nlevels), where nlevels is the number of alpha levels at which to calculate the overlap metric. For each of the last two dimensions of the output array, the first two dimensions form an nspecies by nspecies matrix giving each pairwise calculation of overlap metric between two species for given , , and alpha. In each of these matrices, Species is along the rows of this matrix and Species is along the columns.
Swanson, H.K., Lysy, M., Stasko, A.D., Power, M., Johnson, J.D., and Reist, J.D. "A new probabilistic method for quantifying n-dimensional ecological niches and niche overlap." Ecology: Statistical Reports 96:2 (2015): 318-324. doi:10.1890/14-0235.1.
overlap.plot(), niw.post(), niiw.post().
# fish data data(fish) # generate parameter draws from the "default" posteriors of each fish nsamples <- 500 system.time({ fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) }) # overlap calculation. use nsamples = nprob = 1e4 for more accurate results. system.time({ over <- overlap(fish.par, nreps = nsamples, nprob = nsamples, alpha = c(.95, .99)) }) # posterior expectations of overlap metrics over.mean <- apply(over*100, c(1:2, 4), mean) round(over.mean) # posterior 95% credible intervals of overlap metrics over.cred <- apply(over*100, c(1:2, 4), quantile, prob = c(.025, .975), na.rm = TRUE) round(over.cred[,,,1]) # display alpha = .95 niche region# fish data data(fish) # generate parameter draws from the "default" posteriors of each fish nsamples <- 500 system.time({ fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) }) # overlap calculation. use nsamples = nprob = 1e4 for more accurate results. system.time({ over <- overlap(fish.par, nreps = nsamples, nprob = nsamples, alpha = c(.95, .99)) }) # posterior expectations of overlap metrics over.mean <- apply(over*100, c(1:2, 4), mean) round(over.mean) # posterior 95% credible intervals of overlap metrics over.cred <- apply(over*100, c(1:2, 4), quantile, prob = c(.025, .975), na.rm = TRUE) round(over.cred[,,,1]) # display alpha = .95 niche region
Plots the posterior distribution of the niche region overlap metric calculated for each pairwise combination of species.
overlap.plot( over.stat, nbreaks = 50, equal.axis = FALSE, species.names, col, mean.cred = TRUE, mean.cred.col = "green", xlab )overlap.plot( over.stat, nbreaks = 50, equal.axis = FALSE, species.names, col, mean.cred = TRUE, mean.cred.col = "green", xlab )
over.stat |
An array with |
nbreaks |
Number of breaks in the histogram. Defaults to 50. |
equal.axis |
Logical. If |
species.names |
A vector of species names. Defaults to |
col |
A vector of the colours in which each species will be drawn. |
mean.cred |
Logical. If |
mean.cred.col |
Colour of the mean and credible interval lines in the histogram. |
xlab |
Optional plot title, located on the bottom. Default is no title. |
This function uses the overlap metric information in over.stat calculated by overlap() to create 2-dimensional plots of interspecific niche region overlap.
Returns a series of histograms illustrating the probability of pairwise species overlap.
overlap(), niw.post(), niiw.post().
# fish data data(fish) # parameter draws from the "default" posteriors of each fish nsamples <- 500 system.time({ fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) }) # overlap calculation system.time({ over <- overlap(fish.par, nreps = nsamples, nprob = nsamples, alpha = c(.95, .99)) }) # overlap plot clrs <- c("black", "red", "blue", "orange") # color for each species ii <- 1 # which niche region alpha level to use overlap.plot(over[,,,ii], col = clrs, mean.cred.col = "turquoise", xlab = paste0("Overlap Probability (%) -- Niche Region Size: ", dimnames(over)[[4]][ii]))# fish data data(fish) # parameter draws from the "default" posteriors of each fish nsamples <- 500 system.time({ fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) }) # overlap calculation system.time({ over <- overlap(fish.par, nreps = nsamples, nprob = nsamples, alpha = c(.95, .99)) }) # overlap plot clrs <- c("black", "red", "blue", "orange") # color for each species ii <- 1 # which niche region alpha level to use overlap.plot(over[,,,ii], col = clrs, mean.cred.col = "turquoise", xlab = paste0("Overlap Probability (%) -- Niche Region Size: ", dimnames(over)[[4]][ii]))
Overlap calculation for uniform niche regions.
overlap.unif(muA, SigmaA, muB, SigmaB, alphaA = 0.95, alphaB = 0.95, nprob) overlap.sphere(muA, sigmaA, muB, sigmaB, alphaA = 0.95, alphaB = 0.95)overlap.unif(muA, SigmaA, muB, SigmaB, alphaA = 0.95, alphaB = 0.95, nprob) overlap.sphere(muA, sigmaA, muB, sigmaB, alphaA = 0.95, alphaB = 0.95)
muA, muB
|
Mean of niche regions. |
SigmaA, SigmaB
|
Variance matrix of elliptical niche regions. |
alphaA, alphaB
|
Probabilistic size of niche regions. |
nprob |
Number of uniform draws from niche region |
sigmaA, sigmaB
|
standard deviations (scalars) of spherical niche regions. |
The overlap between niche regions and is defined as , where the hypervolume of an -dimensional region is . For elliptical niche regions, there are simple formulas for and . Thus, we need only determine the volume of the intersection , as the volume of the union is given by the formula .
For spherical niche regions, has a closed-form expression (see 'References'). For elliptical regions, no such formula exists and a Monte Carlo method is used instead. That is, is calculated by sampling uniformly from , then multiplying by the fraction of sampled points which fall into .
While the uniform overlap metric is invariant to permutation of niche regions and , the accuracy of the Monte Carlo calculation of is not: higher accuracy is obtained when a higher fraction of sampled points are in the opposite niche region. overlap.unif() does not attempt to determine for which region this is the case, though the choice can be informed by plotting the niche regions, e.g., with niche.plot().
A Monte Carlo estimate of the niche overlap for overlap.unif(), and an analytic calculation for overlap.sphere().
Li, S. "Concise formulas for the area and volume of a hyperspherical cap." Asian Journal of Mathematics & Statistics 4.1 (2011): 66-70. doi:10.3923/ajms.2011.66.70.
# spherical case: compare Monte Carlo method to analytic formula d <- 2 # 2D example mA <- rnorm(d) mB <- rnorm(d) sigA <- rexp(1) SigA <- sigA^2 * diag(d) sigB <- rexp(1) SigB <- sigB^2 * diag(d) # plot circles ellA <- ellipse(mA, SigA) ellB <- ellipse(mB, SigB) plot(0, type = "n", xlim = range(ellA[,1], ellB[,1]), ylim = range(ellA[,2], ellB[,2]), xlab = "x", ylab = "y") lines(ellA, col = "red") lines(ellB, col = "blue") legend("topright", legend = c("niche A", "niche B"), fill = c("red", "blue"), bg = "white") # compare niche calculations overlap.sphere(mA, sigA, mB, sigB) overlap.unif(mA, SigA, mB, SigB, nprob = 1e5)# spherical case: compare Monte Carlo method to analytic formula d <- 2 # 2D example mA <- rnorm(d) mB <- rnorm(d) sigA <- rexp(1) SigA <- sigA^2 * diag(d) sigB <- rexp(1) SigB <- sigB^2 * diag(d) # plot circles ellA <- ellipse(mA, SigA) ellB <- ellipse(mB, SigB) plot(0, type = "n", xlim = range(ellA[,1], ellB[,1]), ylim = range(ellA[,2], ellB[,2]), xlab = "x", ylab = "y") lines(ellA, col = "red") lines(ellB, col = "blue") legend("topright", legend = c("niche A", "niche B"), fill = c("red", "blue"), bg = "white") # compare niche calculations overlap.sphere(mA, sigA, mB, sigB) overlap.unif(mA, SigA, mB, SigB, nprob = 1e5)
Generates random draws from a Normal-Inverse-Wishart (NIW) distribution. Can be used to compare prior to posterior parameter distributions.
rniw(n, lambda, kappa, Psi, nu)rniw(n, lambda, kappa, Psi, nu)
n |
Number of samples to draw. |
lambda |
Location parameter. See 'Details'. |
kappa |
Scale parameter. See 'Details'. |
Psi |
Scale matrix. See 'Details'. |
nu |
Degrees of freedom. See 'Details'. |
The NIW distribution is defined as
Returns a list with elements mu and Sigma of sizes c(n,length(lambda)) and c(nrow(Psi),ncol(Psi),n).
rwish(), niw.mom(), niw.coeffs().
d <- 4 # number of dimensions nu <- 7 # degrees of freedom Psi <- crossprod(matrix(rnorm(d^2), d, d)) # scale lambda <- rnorm(d) kappa <- 2 n <- 1e4 niw.sim <- rniw(n, lambda, kappa, Psi, nu) # diagonal elements of Sigma^{-1} are const * chi^2 S <- apply(niw.sim$Sigma, 3, function(M) diag(solve(M))) ii <- 2 const <- solve(Psi)[ii,ii] hist(S[ii,], breaks = 100, freq = FALSE, main = parse(text = paste0("\"Histogram of \"*(Sigma^{-1})[", ii,ii,"]")), xlab = parse(text = paste0("(Sigma^{-1})[", ii,ii,"]"))) curve(dchisq(x/const, df = nu)/const, from = min(S[ii,]), to = max(S[ii,]), col = "red", add = TRUE) # elements of mu have a t-distribution mu <- niw.sim$mu ii <- 4 const <- sqrt(Psi[ii,ii]/(kappa*(nu-d+1))) hist(mu[,ii], breaks = 100, freq = FALSE, main = parse(text = paste0("\"Histogram of \"*mu[", ii, "]")), xlab = parse(text = paste0("mu[", ii, "]"))) curve(dt((x-lambda[ii])/const, df = nu-d+1)/const, add = TRUE, col = "red")d <- 4 # number of dimensions nu <- 7 # degrees of freedom Psi <- crossprod(matrix(rnorm(d^2), d, d)) # scale lambda <- rnorm(d) kappa <- 2 n <- 1e4 niw.sim <- rniw(n, lambda, kappa, Psi, nu) # diagonal elements of Sigma^{-1} are const * chi^2 S <- apply(niw.sim$Sigma, 3, function(M) diag(solve(M))) ii <- 2 const <- solve(Psi)[ii,ii] hist(S[ii,], breaks = 100, freq = FALSE, main = parse(text = paste0("\"Histogram of \"*(Sigma^{-1})[", ii,ii,"]")), xlab = parse(text = paste0("(Sigma^{-1})[", ii,ii,"]"))) curve(dchisq(x/const, df = nu)/const, from = min(S[ii,]), to = max(S[ii,]), col = "red", add = TRUE) # elements of mu have a t-distribution mu <- niw.sim$mu ii <- 4 const <- sqrt(Psi[ii,ii]/(kappa*(nu-d+1))) hist(mu[,ii], breaks = 100, freq = FALSE, main = parse(text = paste0("\"Histogram of \"*mu[", ii, "]")), xlab = parse(text = paste0("mu[", ii, "]"))) curve(dt((x-lambda[ii])/const, df = nu-d+1)/const, add = TRUE, col = "red")
Generates a random samples from a Wishart distribution defined as , or an Inverse-Wishart distribution defined as .
rwish(n, Psi, nu, inv = FALSE)rwish(n, Psi, nu, inv = FALSE)
n |
Number of samples to draw. |
Psi |
Scale matrix. |
nu |
Degrees of freedom. |
inv |
Logical. Setting |
Setting inv = TRUE replaces by and inverts the output random matrices, such that they are being generated from an Inverse-Wishart distribution.
Returns an array of Wishart (or Inverse-Wishart) draws of size c(nrow(Psi),ncol(Psi),n).
d <- 4 # number of dimensions nu <- 7 # degrees of freedom Psi <- crossprod(matrix(rnorm(d^2), d, d)) # scale matrix n <- 1e4 Sigma <- rwish(n, Psi, nu) # for any vector a, X = (a' Sigma a) has a const * chi^2 distribution a <- rnorm(d) X <- apply(Sigma, 3, function(S) crossprod(a, S %*% a)) const <- c(a %*% Psi %*% a) hist(X, breaks = 100, freq = FALSE, main = parse(text = "\"Histogram of \"*X==a*minute*Sigma*a"), xlab = parse(text = "X==a*minute*Sigma*a")) curve(dchisq(x/const, df = nu)/const, from = min(X), to = max(X), col = "red", add = TRUE)d <- 4 # number of dimensions nu <- 7 # degrees of freedom Psi <- crossprod(matrix(rnorm(d^2), d, d)) # scale matrix n <- 1e4 Sigma <- rwish(n, Psi, nu) # for any vector a, X = (a' Sigma a) has a const * chi^2 distribution a <- rnorm(d) X <- apply(Sigma, 3, function(S) crossprod(a, S %*% a)) const <- c(a %*% Psi %*% a) hist(X, breaks = 100, freq = FALSE, main = parse(text = "\"Histogram of \"*X==a*minute*Sigma*a"), xlab = parse(text = "X==a*minute*Sigma*a")) curve(dchisq(x/const, df = nu)/const, from = min(X), to = max(X), col = "red", add = TRUE)