######################################################## require(mvtnorm) # package for density and simulation of multivariate normals source("ecol-functions3.R") # functions to assist with ellipse and overlap calculations ######################################################## # figures and tables in manuscript ######################################################## # figure 1: area vs. probability overlap niche.sd <- 2 # niche size, ie number of sd's on each side of mean # Species A is N(0,1). # Species B is N(mu,sd) sd <- .4 eps <- -.1 mu <- (sd-1)*niche.sd - eps # this makes the left endpoints of each niche differ by eps #mu <- niche.sd*sd + eps # makes the mean of species A fall on the endpoint of species B #mu <- -1 # niches regions nicheA <- niche.sd * c(-1,1) nicheB <- mu + niche.sd * sd * c(-1,1) # plot the normals with their niches xlim <- range(c(-4,4,mu-4*sd,mu+4*sd)) # range of plot: 4 sds on each side of mean x <- seq(xlim[1], xlim[2], len = 200) yA <- dnorm(x) # niche A yB <- dnorm(x, mu, sd) # niche B ylim <- range(yA, yB) par(mai = c(.8, .82, .1, .8)) plot(x, dnorm(x), type = "l", xlab = "Isotope x", ylab = "p(x)", xlim = xlim, ylim = ylim) lines(x, dnorm(x, mu, sd), col = "red") abline(h = 0) xL <- c(-niche.sd, mu-niche.sd*sd) # lower niche endpoints xU <- c(niche.sd, mu+niche.sd*sd) # upper niche endpoints xM <- c(0, mu) # niche means segments(x0 = c(xL, xU, xM), y0 = 0, x1 = c(xL, xU, xM), y1 = dnorm(c(xL, xU, xM), mean = rep(c(0, mu), 3), sd = rep(c(1, sd), 3)), col = rep(c("black", "red"), 3), lty = rep(c(1,1,2), each = 2)) # legend legend(x = "topright", legend = c("Species A ", "Species B "), fill = c("black", "red")) # overlap calculations over <- matrix(NA, 2,2) dimnames(over) <- list(c("A on B", "B on A"), c("prob", "frac")) # area method over[,2] <- (min(nicheA[2],nicheB[2]) - max(nicheA[1],nicheB[1])) over[,2] <- over[,2]/c(diff(nicheB), diff(nicheA)) # probability method over[,1] <- c(diff(pnorm(nicheB, mean = 0, sd = 1)), diff(pnorm(nicheA, mean = mu, sd = sd))) signif(over,2) ################################################################################ # figure 2: ellipsoid projections require(rgl) level <- .95 d <- 3 Mu <- rep(0, d) V <- matrix(0, d, d) diag(V) <- c(3,1,2) # ellipsoid tmp <- ellipse3d(V, centre = Mu, level = level) plot3d(tmp, col = "yellow", type = "wire", axes = FALSE, aspect = "iso", xlab = "", ylab = "", zlab = "") # axes bbox <- par3d("bbox") segments3d(x = rep(bbox[2],2), y = bbox[4:3], z = rep(bbox[6],2)) text3d(x = bbox[2]+.5, y = bbox[3], z = bbox[6], texts = "U") segments3d(x = bbox[2:1], y = rep(bbox[4],2), z = rep(bbox[6],2)) text3d(x = bbox[1]+.5, y = bbox[4]+.5, z = bbox[6], texts = "V") segments3d(x = rep(bbox[2],2), y = rep(bbox[4],2), z = bbox[6:5], col = "blue", lwd = 2) text3d(x = bbox[2], y = bbox[4]+.5, z = .7*bbox[5] + .3*bbox[6], texts = "W") # geometric projection onto UV plane tmp2 <- tmp$vb[,which(tmp$vb[3,] == 0)] # order by angle tmp2 <- tmp2[,order(atan2(tmp2[2,], tmp2[1,]))] tmp2[3,] <- bbox[6] tmp2a <- polygon3d(tmp2[1,], tmp2[2,], tmp2[3,], plot = FALSE, fill = FALSE) n2a <- nrow(tmp2a)-1 # outline sapply(1:n2a, function(i) { invisible(segments3d(x = tmp2a[i+0:1,1], y = tmp2a[i+0:1,2], z = tmp2a[i+0:1,3], col = "black", lwd = 1)) }) s2a <- 18 # shade region apply(cbind((s2a-(1:(n2a/2))) %% n2a + 1, (s2a-1+(1:(n2a/2))) %% n2a + 1), 1, function(i) { invisible(segments3d(x = tmp2a[i,1], y = tmp2a[i,2], z = tmp2a[i,3], col = "black")) }) # probabilistic niche in UV plane tmp3 <- ellipse3d(V, centre = Mu, level = pchisq(qchisq(p = level, df = 3), df = 2)) tmp3 <- tmp3$vb[,which(tmp3$vb[3,] == 0)] tmp3 <- tmp3[,order(atan2(tmp3[2,], tmp3[1,]))] tmp3[3,] <- bbox[6] tmp3a <- polygon3d(tmp3[1,], tmp3[2,], tmp3[3,], plot = FALSE, fill = FALSE) sapply(1:(nrow(tmp3a)-1), function(i) { invisible(segments3d(x = tmp3a[i+0:1,1], y = tmp3a[i+0:1,2], z = tmp3a[i+0:1,3], col = "red", lwd = 2)) }) ######################################################################### # figure 3: ellipse plots # generate 10 parameter draws from the posteriors of each fish 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 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 plot.niche(niche.par = fish.par, niche.data = fish.data, col = clrs, xlab = expression("Isotope Ratio (\u2030)")) ######################################################################### # figure 4: overlap plots # generate 1e4 posterior parameter draws for each fish nsamples <- 1e4 system.time({ fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) }) # overlap calculation alpha <- .95 system.time({ over <- overlap(fish.par, nreps = 1e4, nprob = 1e4, alpha.level = alpha) }) # overlap plot clrs <- c("black", "red", "blue", "orange") # color for each species ii <- 3 # use 95% region plot.overlap(over, col = clrs, mean.cred.col = "turquoise", xlab = "Overlap Probability (%) -- Niche Size: 95\%") ############################################################## # figures and tables for supplementary material ############################################################## # figure S1: probabilistic niche regions alpha <- .95 # 95% coverage level alen <- .05 ttl <- c("Symmetric", "Skewed", "Bimodal") par(mfrow = c(1,3)) # normal distribution x <- seq(-5, 5, len = 100) f <- function(x) dnorm(x) # pdf y <- f(x) par(mai = c(.3, .3, .3, .05)) plot(x, y, type = "l", axes = TRUE, xlab = NA, ylab = NA, main = ttl[1]) xL <- qnorm((1-alpha)/2) xU <- -xL xseg <- c(0, xL, xU) yseg <- dnorm(xseg) segments(x0 = xseg, y0 = 0, x1 = xseg, y1 = yseg, lty = 2, col = c("black", "red", "red")) abline(h = 0) arrows(x0 = xL, x1 = xU, y0 = 0, col = "red", lwd = 2, angle = 90, code = 3, length = alen) points(x = c(xL, xU), y = c(0, 0), pch = 4, lwd = 2) # gamma distribution a <- 3 b <- 1 md <- (a-1)*b x <- seq(0, 15, len = 100) f <- function(x) dgamma(x, shape = a, scale = b) # pdf F <- function(x) pgamma(x, shape = a, scale = b) # cdf y <- f(x) xL <- .303 dL <- f(xL) xU <- uniroot(function(xu) f(xu)-dL, lower = md, upper = max(x))$root c(f(xL), f(xU)) # these should be identical diff(F(c(xL, xU))) # this should be 95% xseg <- c(md, xL, xU) yseg <- f(xseg) plot(x, y, type = "l", xlab = NA, ylab = NA, main = ttl[2]) segments(x0 = xseg, y0 = 0, x1 = xseg, y1 = yseg, lty = 2, col = c("black", "red", "red")) abline(h = 0) arrows(x0 = xL, x1 = xU, y0 = 0, col = "red", lwd = 2, angle = 90, code = 3, length = alen) points(x = qgamma(c(.025, .975), shape = a, scale = b), y = c(0, 0), pch = 4, lwd = 2) # bimodal distribution x <- seq(-8, 8, len = 100) f <- function(x) .5*(dnorm(x, -3) + dnorm(x, 3)) # pdf F <- function(x) .5*(pnorm(x, -3) + pnorm(x, 3)) # cdf y <- f(x) x1 <- -4.959 d1 <- f(x1) x2 <- uniroot(function(x) f(x) - d1, lower = x1+.1, upper = 0)$root xL <- -4.645 xU <- -xL c(f(x1), f(x2), f(-x2), f(-x1)) # these should be identical diff(F(c(x1, x2))) + diff(F(c(-x2, -x1))) # this should be 95% md <- c(-3, 3) xseg <- c(md, x1, x2, -x2, -x1) yseg <- f(xseg) plot(x, y, type = "l", xlab = NA, ylab = NA, main = ttl[3]) segments(x0 = xseg, y0 = 0, x1 = xseg, y1 = yseg, lty = 2, col = c("black", "black", "red", "red", "red", "red")) abline(h = 0) arrows(x0 = c(x1, -x2), x1 = c(x2, -x1), y0 = 0, col = "red", lwd = 2, angle = 90, code = 3, length = alen) points(x = c(xL, xU), y = c(0, 0), pch = 4, lwd = 2) ######################################################## # figure S2: niche overlap calculation ns <- .95 # niche size d <- 2 prob <- c(.1, .25, .5, .75, .95) MuA <- c(2.1, -.27) MuB <- c(.22, 1.3) VA <- matrix(c(2.8, -.80, -.80, .76), d, d) VB <- matrix(c(1.0, -.55, -.55, .70), d, d) sd <- 3 xA <- sapply(1:d, function(i) MuA[i] + sd*c(-1,1)*sqrt(diag(VA)[i])) xB <- sapply(1:d, function(i) MuB[i] + sd*c(-1,1)*sqrt(diag(VB)[i])) n <- 100 x <- seq(min(xA[1,1], xB[1,1]), max(xA[2,1], xB[2,1]), len = n) y <- seq(min(xA[1,2], xB[1,2]), max(xA[2,2], xB[2,2]), len = n) xy <- as.matrix(expand.grid(x=x,y=y)) zA <- matrix(dmvnorm(xy, mean = MuA, sigma = VA), n, n) zB <- matrix(dmvnorm(xy, mean = MuB, sigma = VB), n, n) # contour plots par(mai = c(.8, .82, .1, .8)) contour(x, y, zA, levels = norm.level(MuA, VA, prob = prob), labels = prob, xlab = "X", ylab = "Y", "Overlap of A on B") contour(x, y, zB, col = "red", add = TRUE, levels = norm.level(MuB, VB, prob = ns), drawlabels = FALSE) # add sample points nsim <- 1e2 xy.sim <- rmvnorm(nsim, mean = MuA, sigma = VA) over <- dmvnorm(xy.sim, mean = MuB, sigma = VB) > norm.level(MuB, VB, ns) points(xy.sim, pch = ifelse(over, 20, 1), cex = .7) # legend legend(x = "topright", legend = c("Species A ", "Species B "), fill = c("black", "red")) ################################################################## # figure S3: posterior densities with default prior, simulated data # fish data fish <- read.csv("ecol-fish2.csv") # in the package this will be loaded differently # estimate the fish data parameters species.names <- c("ARCS", "BDWF", "LKWF", "LSCS") ii <- 1 # use arctic cisco X0 <- fish[fish$species == species.names[ii], 2:4] N <- nrow(X0) mu <- colMeans(X0) Sigma <- var(X0) # simulate data X <- rmvnorm(N, mu, Sigma) # generate posterior draws from the default posterior nsamples <- 1e4 system.time({ Theta <- niw.post(nsamples = nsamples, X = X) }) # plot nbreaks <- 40 par(mfrow = c(4,3), mar = c(4.2, 4.2, .5, 1)+.1) # mu for(ii in 1:3) { hist(Theta$mu[,ii], breaks = nbreaks, freq = FALSE, main = "", xlab = parse(text = paste0("mu[", ii, "]")), ylab = parse(text = paste0("p(mu[", ii, "]*\" | \"*X)"))) abline(v = mu[ii], col = "red") } # Sigma for(ii in 1:3) { for(jj in 1:3) { hist(Theta$Sigma[ii,jj,], breaks = nbreaks, freq = FALSE, main = "", xlab = parse(text = paste0("Sigma[", ii, "*", jj, "]")), ylab = parse(text = paste0("p(Sigma[", ii, "*", jj, "]*\" | \"*X)"))) abline(v = Sigma[ii,jj], col = "red") } } ################################################################## # table S1: true coverage of 95% credible intervals, simulated data # obtain 1000 datasets and posteriors nreps <- 1e3 sim.res <- replicate(nreps, { # simulate data X <- rmvnorm(N, mu, Sigma) # generate posterior draws from the default posterior Theta <- niw.post(nsamples = 1e3, X = X) list(X = X, mu = Theta$mu, Sigma = Theta$Sigma) }) # coverage mu.cov <- matrix(NA, nreps, 3) colnames(mu.cov) <- names(mu) Sigma.cov <- array(NA, dim = c(3, 3, nreps)) dimnames(Sigma.cov) <- list(names(mu), names(mu), NULL) for(ii in 1:nreps) { mu.LU <- apply(sim.res[[2,ii]], 2, quantile, prob = c(.025, .975)) Sigma.LU <- apply(sim.res[[3,ii]], 1:2, quantile, prob = c(.025, .975)) mu.cov[ii,] <- (mu.LU[1,] <= mu) & (mu <= mu.LU[2,]) Sigma.cov[,,ii] <- (Sigma.LU[1,,] <= Sigma) & (Sigma <= Sigma.LU[2,,]) } # display in table f <- function() { mu.prob <- signif(colMeans(mu.cov)*100, 3) Sigma.prob <- signif(apply(Sigma.cov, 1:2, mean)*100, 3) message("\\begin{tabular}{ccclcccc}") message("& & & & \\multicolumn{4}{c}{$\\Sigma$} \\\\") message("\\hhline{~~~~====}") message("& & & & & $\\diso{15}{N}$ & $\\diso{13}{C}$ & $\\diso{34}{S}$ \\\\") message("\\cline{5-8}") message(paste(c("\\multicolumn{3}{c}{$\\mu$}", "", "$\\diso{15}{N}$", Sigma.prob[1,]), collapse = " & "), " \\\\") message("\\hhline{===~~~~~}") message(paste(c("$\\diso{15}{N}$", "$\\diso{13}{C}$", "$\\diso{34}{S}$", "", "$\\diso{13}{C}$", Sigma.prob[2,]), collapse = " & "), " \\\\") message("\\cline{1-3}") message(paste(c(mu.prob, "", "$\\diso{34}{S}$", Sigma.prob[3,]), collapse = " & "), " \\\\") message("\\cline{1-3}\\cline{5-8}") message("\\end{tabular}") } f() ################################################################## # table S2: summary of fish isotope data fish <- read.csv(file = "ecol-fish2.csv") # isotope data # summary table f <- function() { sp.names <- c("Arctic Cisco", "Broad Whitefish", "Lake Whitefish", "Least Cisco") tab <- tapply(1:length(fish$SPP), fish$SPP, function(i) { tmp <- signif(rbind(mean = colMeans(fish[i,2:4]), sd = apply(fish[i,2:4], 2, sd)), 2) paste(c(sp.names[fish[i[1],1]], length(i), apply(tmp, 2, function(x) paste0(x[1], "(", x[2], ")"))), collapse = " & ") }) message("\\begin{tabular}{ccccc}") message("\\hline \\hline") message(" & & \\multicolumn{3}{c}{Isotope -- mean(sd)} \\\\") message("Species & $N_{\\textrm{samples}}$ & $\\delta^{15}\\textrm{N}$ & $\\delta^{13}\\textrm{C}$ & $\\delta^{34}\\textrm{S}$ \\\\") message("\\hline") sapply(tab, function(x) message(x, " \\\\")) message("\\hline") message("\\end{tabular}") } f() ################################################################## # table S3: summary of overlap calculations # generate 1e4 posterior parameter draws for each fish nsamples <- 1e4 system.time({ fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) }) # overlap calculation alpha <- c(.8, .9, .95, .99) system.time({ over <- overlap(fish.par, nreps = 1e4, nprob = 1e4, alpha.level = alpha) }) # table. f <- function(over, alpha) { nspec <- dim(over)[1] spec.names <- dimnames(over)[[1]] over.mean <- apply(over, 1:2, mean)*100 over.cred <- apply(over, 1:2, quantile, prob = c(.025, .975), na.rm = TRUE)*100 over.cred <- array(over.cred, dim = c(2, nspec, nspec)) ti <- combn(nspec, 2) tmp2 <- t(ti) tmp <- rbind(apply(tmp2, 1, function(ii) paste(spec.names[ii[2:1]], collapse = " & ")), signif(over.mean[tmp2], 2), signif(over.cred[cbind(1,tmp2)], 2), signif(over.cred[cbind(2,tmp2)], 2)) tmp2 <- t(ti[2:1,]) tmp <- rbind(tmp, apply(tmp2, 1, function(ii) paste(spec.names[ii[2:1]], collapse = " & ")), signif(over.mean[tmp2], 2), signif(over.cred[cbind(1,tmp2)], 2), signif(over.cred[cbind(2,tmp2)], 2)) tmp <- t(tmp) message("\\begin{tabular}{cccclcccc}") message("\\multicolumn{9}{c}{Niche Size: $\\alpha = ", alpha, "$}\\\\") message("\\hhline{====~====}") message("\\multicolumn{2}{c}{$A$ onto $B$} & Mean & 95\\% C.I. & & \\multicolumn{2}{c}{$B$ onto $A$} & Mean & 95\\% C.I. \\\\") message("\\cline{1-4}\\cline{6-9}") message(paste(apply(tmp, 1, function(x) { paste0(x[1], " & ", x[2], " & (", x[3], ", ", x[4], ") & & ", x[5], " & ", x[6], " & (", x[7], ", ", x[8], ") \\\\") }), collapse = "\n")) message("\\cline{1-4}\\cline{6-9}") message("\\end{tabular}") } for(ii in 1:4) { f(over[,,,ii], alpha = alpha[ii]) } ######################################################################### # tables S4-7: raw data f <- function(lake.ind, ncolumns = 2) { sp.names <- c("Arctic Cisco", "Broad Whitefish", "Lake Whitefish", "Least Cisco") ind <- which(fish$SPP == levels(fish$SPP)[lake.ind]) br <- cut(1:length(ind), breaks = ncolumns) nr <- table(br) message("\\begin{table}[!ht]\n\\centering") message("\\caption{Stable isotope samples of ", sp.names[lake.ind], ".}\\label{tab:", levels(fish$SPP)[lake.ind], "}") for(ii in 1:ncolumns) { message("\\begin{tabular}{ccc}") message("$\\diso{15}N$ & $\\diso{13}C$ & $\\diso{34}S$ \\\\ \\hline") tmp <- apply(fish[ind[br==levels(br)[ii]],2:4], 1, function(x) { paste(round(x,2), collapse = " & ") }) message(paste(tmp, collapse = " \\\\ \n")) if(nr[ii] < max(nr)) { message("\\\\") message(paste(rep("&&", max(nr)-nr[ii]), collapse = "\\\\ \n")) } message("\\end{tabular}") if(ii < ncolumns) message("\\quad") } message("\\end{table}") } invisible(sapply(1:4, f, ncolumns = 3)) ################################################## # # R functions for ecological Niche Region and Niche Overlap metrics # this is the content of "ecol-functions3.R" is part of the package nicheROVER. # ################################################# # cutoff level of a multivariate Normal distribution for highest probability regions. # returns a value L such that if X ~ N(mu, V) has pdf f(x), then P(f(X) > L) = prob. # mu: mean of Normal # V: variance of Normal # prob: probability level norm.level <- function(mu, V, prob = .95) { d <- length(mu) q <- qchisq(prob, df = d) exp(-.5 * (q + determinant(V)$modulus[1] + d * log(2*pi))) } # coordinates of points on a 2-d ellipse # that is, X = c(x,y) satisfies (X-mu)' %*% V^{-1} %*% (X-mu) = C # mu: center of ellipse # V: scale of ellipse # C: value of contour. # n: number of points to return # returns a 2-column matrix of points on the ellipse in the xy plane. ellipse <- function(mu, V, C, n = 100) { tmp <- eigen(V) hlen <- sqrt(C*tmp$val) theta <- atan2(tmp$vec[2,1], tmp$vec[1,1]) t <- seq(0, 2*pi, len = n+1) x <- hlen[1] * cos(t) y <- hlen[2] * sin(t) alpha <- atan2(y, x) rad <- sqrt(x^2 + y^2) cbind(x = rad * cos(alpha + theta) + mu[1], y = rad * sin(alpha + theta) + mu[2]) } # random sample from a Wishart distribution W(Psi, nu). # setting inv = TRUE replaces Psi by Psi^{-1} and inverts the output random matrices, # such that they are being generated from an inverse Wishart iW(Psi, nu) distribution. # n: number of samples to return. # ***don't bother documenting this because it's going to end up being an internal # function, i.e. won't be directly callable by the user.*** # returns an array of wishart (or inverse wishart) draws of size c(dim(Psi), n) rwish <- function(n, Psi, nu, inv = FALSE, chol = TRUE) { if(inv) Psi <- solve(Psi) U <- chol(Psi) d <- nrow(Psi) ans <- array(0, dim = c(d, d, n)) if(!is.null(dimnames(Psi))) dimnames(ans) <- c(dimnames(Psi), list(NULL)) ans[rep(upper.tri(Psi), n)] <- rnorm(n*d*(d-1)/2) ans[rep(!lower.tri(Psi, diag = FALSE) & !upper.tri(Psi, diag = FALSE), n)] <- sqrt(rchisq(n*d, df = nu-1:d+1)) for(ii in 1:n) { tmp <- ans[,,ii] %*% U if(inv) tmp <- backsolve(tmp, diag(d), transpose = TRUE) ans[,,ii] <- crossprod(tmp) } ans } # random draws from a Normal-inverse-Wishart distribution. # returns a list with elements mu and Sigma, where # Sigma ~ iW(Psi, nu) # mu|Sigma ~ N(lambda, Sigma/kappa) # n: number of random samples # returns a list with components mu and Sigma of sizes c(n, length(lambda)) # and c(dim(Psi), n) rniw <- function(n, lambda, kappa, Psi, nu) { d <- length(lambda) Sigma <- rwish(n, Psi, nu, inv = TRUE) mu <- matrix(NA, n, d) colnames(mu) <- names(lambda) for(ii in 1:n) { mu[ii,] <- rmvnorm(1, mean = lambda, sigma = Sigma[,,ii]/kappa) } list(mu = mu, Sigma = Sigma) } # given X = rbind(X1, ..., XN) with Xi ~iid N(mu, Sigma), returns the coefficients of the # Normal-inverse-Wishart (NIW) posterior distribution p(mu, Sigma | X) for the NIW prior: # Sigma ~ iW(Psi, nu) # mu|Sigma ~ N(lambda, Sigma/kappa). # default values: # kappa = 0 uses the lebesgue prior on mu: p(mu) \propto 1. # all(Psi == 0) uses the scale invariant prior on Sigma: # p(Sigma) \propto |Sigma|^{-(nu+d+1)/2}. # nu = ncol(X)+1, with kappa =0 and all(Psi == 0) # makes E[mu | X] = colMeans(X) and E[Sigma | X] = var(X). # returns a list with values lambda, kappa, Psi, nu corresponding # to the NIW posterior p(mu, Sigma | X). # a list of posterior coefficients named lambda, kappa, Psi, and nu. niw.coeffs <- function(X, lambda, kappa, Psi, nu) { d <- ncol(X) N <- nrow(X) if(missing(kappa)) kappa <- 0 if(missing(Psi)) Psi <- 0 if(missing(nu)) nu <- ncol(X)+1 # sufficient statistics Xbar <- colMeans(X) S <- t(X)-Xbar S <- S %*% t(S) # posterior parameter values Psi2 <- Psi + S lambda2 <- N*Xbar if(kappa != 0) { Psi2 <- Psi2 + (N*kappa)/(N+kappa) * (Xbar-lambda) %*% t(Xbar-lambda) lambda2 <- lambda2 + kappa*lambda } lambda2 <- lambda2/(N+kappa) nu2 <- N+nu-(kappa==0) kappa2 <- N+kappa list(lambda = lambda2, kappa = kappa2, Psi = Psi2, nu = nu2) } # Monte Carlo (iid) draw from the posterior distribution of mu and Sigma for # X = rbind(X1, ..., XN) with Xi ~iid N(mu, Sigma) and the prior: # Sigma ~ iW(Psi, nu) # mu|Sigma ~ N(lambda, Sigma/kappa). # nsamples: the number of posterior samples # Uses the same default parameter values at niw.coeffs. # returns a list with components mu and Sigma (see rniw) niw.post <- function(nsamples, X, lambda, kappa, Psi, nu) { par <- niw.coeffs(X, lambda, kappa, Psi, nu) rniw(nsamples, par$lambda, par$kappa, par$Psi, par$nu) } # Gibbs sampling from the posterior distribution of mu and Sigma for # X = rbind(X1, ..., XN) with Xi ~iid N(mu, Sigma) and the prior: # Sigma ~ iW(Psi, nu) # mu|Sigma ~ N(lambda, V). # nsamples: the number of posterior samples. # burn: the burn-in, either a +ve integer or a fraction of nsamples. defaults to 1/10. # some default values: # all(V == 0) uses the lebesgue prior p(mu) \propto 1, though this choice of prior # reduces to an NIW for which iid draws can be obtained using niw.coeffs and rniw. # all(Psi == 0) uses the scale-invariant prior p(Sigma) \propto |Sigma|^{-(nu+d+1)/2}. # mu0: initial value of mu used to start the sampler. default is mu0 = lambda. # returns a list with components mu and Sigma (see rniw) niiw.post <- function(nsamples, X, lambda, V, Psi, nu, mu0 = lambda, burn) { # sufficient statistics d <- ncol(X) N <- nrow(X) Xbar <- colMeans(X) S <- t(X)-Xbar S <- S %*% t(S) # local variables mu <- rep(0, d) Sigma <- matrix(0, d, d) if(missing(burn)) burn <- .1 if(burn < 1) burn <- floor(nsamples*burn) mu <- mu0 # output mu.out <- matrix(NA, nsamples, d) Sigma.out <- array(NA, dim = c(d, d, nsamples)) # main for loop for(ii in (-burn+1):nsamples) { # sample Sigma Psi2 <- N * ((mu-Xbar) %o% (mu-Xbar)) + S + Psi nu2 <- N+nu Sigma <- matrix(rwish(1, Psi2, nu2, inv = TRUE), d, d) # sample mu Sigma2 <- Sigma/N if(!all(V == 0)) { B <- Sigma2 %*% solve(V + Sigma2) } else B <- matrix(0, d, d) IB <- diag(d)-B lambda2 <- c(IB %*% Xbar) if(!all(V == 0)) lambda2 <- lambda2 + c(B %*% lambda) V2 <- IB %*% Sigma2 mu <- c(rmvnorm(1, lambda2, V2)) # store if(ii > 0) { mu.out[ii,] <- mu Sigma.out[,,ii] <- Sigma } } list(mu = mu.out, Sigma = Sigma.out) } # default plot for 2-d niche regions # niche.par: a list of length nspecies, each element of which in turn is a list # of means and variance matrices. Each of these will correspond to an ellipse being # drawn for that species in the corresponding 2-dimensional plane. # niche.data: the raw data. this should be passed as a list of length # nspecies, with each element corresponding to a matrix with observations along the rows # and isotopes along the columns. # alpha.level: size of niche region to plot. default is 95% region. # species.names: names of the species. defaults to names(niche.par) # iso.names: names of the isotopes to plot. defaults to colnames(niche.par[[1]]$mu) # ndens: number of points at which to evaluate kernel density estimates # pfrac: fraction of the plot occupied by the 1-d raw data. pfrac = 0 means don't display # the 1-d raw data. # col: colors in which each species will be drawn # xlab: optional title of plot, located on the bottom # returns a plot plot.niche <- function(niche.par, niche.data, alpha.level = .95, species.names, iso.names, col, ndens = 512, pfrac = 0, xlab) { niso <- ncol(niche.par[[1]]$mu) nspec <- length(niche.par) npts <- 100 # number of points for each ellipse nell <- sapply(niche.par, function(x) nrow(x$mu)) # number of ellipses per species if(missing(species.names)) species.names <- names(niche.par) if(missing(iso.names)) iso.names <- colnames(niche.par[[1]]$mu) # create all the ellipses. need to do this straight up in order to get # the plot limits right. ell <- vector("list", nspec) names(ell) <- names(species.names) D <- combn(niso, 2) for(ii in 1:nspec) { ell.tmp <- array(NA, c(nell[ii], ncol(D), npts+1, 2)) for(jj in 1:nell[ii]) { for(kk in 1:ncol(D)) { ell.coord <- ellipse(niche.par[[ii]]$mu[jj, D[,kk]], V = niche.par[[ii]]$Sigma[D[,kk], D[,kk], jj], c = qchisq(alpha.level, df = 2), n = npts) ell.tmp[jj,kk,,] <- ell.coord } } ell[[ii]] <- ell.tmp } # plot limits. lims <- array(sapply(niche.data, function(x) apply(x, 2, range)), dim = c(2, niso, nspec)) lims <- apply(lims, 2, range) # plot par(mfcol = c(niso,niso), mar = rep(.5, 4), oma = rep(4,4)) # plots for(ci in 1:niso) { for(ri in 1:niso) { # initialize plot plot.new() plot.window(lims[,ci], lims[,ri]) if (ci == ri) { # diagonals: density plots xdens <- matrix(NA, ndens, nspec) ydens <- xdens for(ii in 1:nspec) { den <- density(niche.data[[ii]][,ci], n = ndens) xdens[,ii] <- den$x ydens[,ii] <- den$y } for(ii in 1:nspec) { ly <- par("usr")[1:2] ly[2] <- ly[1] + pfrac*(ly[2]-ly[1]) ly[3] <- (ly[2]-ly[1])/nspec segments(x0 = niche.data[[ii]][,ci], y0 = ly[1]+(ii-1)*ly[3], y1 = ly[1]+ii*ly[3], col = col[ii]) ly <- ly[2] + ydens[,ii]/max(ydens)*(lims[2,ci]-ly[2]) lines(xdens[,ii], ly, col = col[ii]) } } if (ri > ci) { # lower triangle: point plots for(ii in 1:nspec) { points(niche.data[[ii]][,c(ci,ri)], col = col[ii], pch = 16) } } if (ri < ci) { # upper triangle: ellipses for(ii in 1:nspec) { for(jj in 1:nell[ii]) { lines(ell[[ii]][jj,which(D[1,] == ri & D[2,] == ci),,2:1], col = col[ii]) } } } box() if(ci == niso) axis(side = 4) else axis(side = 4, labels = FALSE) if(ri == niso) axis(side = 1) else axis(side = 1, labels = FALSE) if(ri == 1) mtext(text = iso.names[ci], side = 3, line = 1) if(ci == 1) mtext(text = iso.names[ri], side = 2, line = 1) } } if(!missing(xlab)) { mtext(text = xlab, side = 1, outer = TRUE, line = 2.2, cex = .9) } legend(x = "topleft", legend = species.names, fill = col, bty = "n", cex = 1.25) } # bayesian Monte Carlo calculation of the overlap metric. # niche.par: a list of length nspecies, each element of which in turn is a list # of means and variance matrices. # nreps: the number of overlap calculations for each pair of species # nprob: the number of normal draws for each overlap calculation # alpha.level: scalar or vector of probability levels for the overlap calculation. # defaults to 0.95. # species.names: names of the species. defaults to names(niche.par) # norm.redraw: T/F. if False, uses the same nprob*nspecies iid N(0,1) draws for each # calculation of the overlap metric. increases the Monte Carlo error but # about 50% faster. # verbose: T/F. calculations can be long, so this produces output which can give # a sense of the timing. # returns an array of size c(nspecies, nspecies, nreps, nlevels), where nspecies is the # number of species and nlevels is the number of alpha levels to calculate. # if nlevels == 1 then the output has dimension c(nspecies, nspecies, nreps). overlap <- function(niche.par, nreps, nprob, alpha = 0.95, species.names, norm.redraw = TRUE) { niso <- ncol(niche.par[[1]]$mu) # number of isotopes nspec <- length(niche.par) # number of species nlevels <- length(alpha) # number of levels if(missing(species.names)) species.names <- names(niche.par) # temporary variables mu <- matrix(NA, niso, nspec) x <- array(NA, dim = c(niso, nprob, nspec)) C <- array(NA, dim = c(niso, niso, nspec)) Zsq <- array(NA, dim = c(nprob, nspec-1, nlevels)) # constants qlevels <- qchisq(alpha, df = niso) qlevels <- array(rep(qlevels, each = nprob*(nspec-1)), dim = c(nprob, nspec-1, nlevels)) if(!norm.redraw) z <- matrix(rnorm(niso*nprob), niso, nprob) # output over <- array(NA, dim = c(nspec, nspec, nreps, nlevels)) dimnames(over) <- list(species.names, species.names, NULL, paste0(alpha*100, "%")) names(dimnames(over)) <- c("Species A", "Species B", "", "alpha") # subsample the parameters posteriors of each niche nreps times ind <- sapply(niche.par, function(mv) { sample(nrow(mv$mu), nreps, replace = TRUE) }) # calculate each overlap for(jj in 1:nreps) { # get means, t(chol(variance)), and draws for each species for(ii in 1:nspec) { mu[,ii] <- niche.par[[ii]]$mu[ind[jj,ii],] C[,,ii] <- t(chol(niche.par[[ii]]$Sigma[,,ind[jj,ii]])) if(norm.redraw) z <- matrix(rnorm(niso*nprob), niso, nprob) x[,,ii] <- C[,,ii] %*% z + mu[,ii] } # check whether the draw of each species is within the ellipse of every other for(ii in 1:nspec) { Z <- backsolve(C[,,ii], matrix(x[,,-ii], nrow = niso)-mu[,ii], upper.tri = FALSE) Zsq[] <- colSums(Z^2) over[-ii,ii,jj,] <- colMeans(Zsq < qlevels) } } if(dim(over)[4] == 1) over <- over[,,,1] over } # plot the overlap metric. # over.stat: an array of size c(nspecies, nspecies, nreps) containing nreps calculations # of the overlap metric for each pair of species (see overlap). # nbreaks: number of histogram breaks. # equal.axis: T/F, if True puts all histograms on the same range of 0-100%. # species.names: names of the species. defaults to dimnames(over.stat)[[1]] # col: colors in which each species will be drawn # mean.cred: T/F include vertical lines for mean and 95% credible intervals # for the histogram of each overlap metric # mean.cred.col: color of the mean and credible interval lines # xlab: optional title of plot, located on the bottom # returns a plot overlap.plot <- function(over.stat, nbreaks = 50, equal.axis = TRUE, species.names, col, mean.cred = TRUE, mean.cred.col = "green", xlab) { if(length(dim(over.stat)) != 3 || dim(over.stat)[1] != dim(over.stat)[2]) stop("Incorrect specification of over.stat.") nspec <- dim(over.stat)[1] if(missing(species.names)) species.names <- dimnames(over.stat)[[1]] # histograms over.hist <- apply(over.stat, 1:2, function(x) { if(any(is.na(x))) return(NULL) else { tmp <- hist(x*100, breaks = nbreaks, plot = FALSE) tmp$density <- tmp$density/max(tmp$density) tmp$counts <- tmp$counts/max(tmp$counts) } tmp }) # x-axis limits xlim <- lapply(over.hist, function(h) { if(is.null(h)) tmp <- matrix(NA, 2, 2) else { tmp <- cbind(range(h$breaks), range(h$density)) } tmp }) xlim <- matrix(xlim, nspec, nspec) if(equal.axis) { for(ii in 1:nspec) { xlim[,ii] <- rep(list(cbind(range(sapply(xlim[,ii], function(ll) ll[,1]), na.rm = TRUE), range(sapply(xlim[,ii], function(ll) ll[,2]), na.rm = TRUE))), nspec) } } # mean and credible intervals if(mean.cred) { over.mean <- apply(over.stat, 1:2, mean)*100 over.cred <- apply(over.stat, 1:2, quantile, prob = c(.025, .975), na.rm = TRUE)*100 over.cred <- array(over.cred, dim = c(2, nspec, nspec)) } # plot par(mfcol = c(nspec,nspec), mar = c(1.5,rep(.5, 3)), oma = rep(4,4)) for(ci in 1:nspec) { for(ri in 1:nspec) { # initialize plot plot.new() if (ri != ci) { # off diagonals: overlap histograms plot.window(xlim[ri,ci][[1]][,1], xlim[ri,ci][[1]][,2]) if(equal.axis) { # recalculate histograms with new limits tmp <- hist(over.stat[ri,ci,]*100, breaks = seq(xlim[ri,ci][[1]][1,1], xlim[ri,ci][[1]][2,1], len = nbreaks+1), plot = FALSE) tmp$density <- tmp$density/max(tmp$density) over.hist[[ri,ci]] <- tmp } plot(over.hist[[ri,ci]], add = TRUE, freq = FALSE, col = col[ci], border = "white") if(mean.cred) { # mean and 95% credible intervals abline(v = c(over.mean[ri,ci], over.cred[,ri,ci]), col = mean.cred.col, lty = c(1,2,2), lwd = 2) } } else { # diagonals: empty plots plot.window(xlim = c(0,1), ylim = c(0,1)) } if(ri == 1 && ci == 1) { text(x = .5, y = .5, labels = expression("Niche Overlap: "*bgroup("(", atop("Row", "Column"), ")")), adj = c(.5, .5), cex = 1) } box() if(ci != ri) axis(side = 1) if(ci > 1) axis(side = 2, labels = FALSE) if(ci < nspec) axis(side = 4, labels = FALSE) if(ri == 1) mtext(text = species.names[ci], side = 3, line = 1, col = col[ci]) if(ci == 1) mtext(text = species.names[ri], side = 2, line = 1) if(mean.cred && ri == nspec && ci == nspec) { legend(x = "center", legend = c("Mean", "95% Credible Interval"), lty = c(1, 2), lwd = 2, col = mean.cred.col) } } } if(!missing(xlab)) { mtext(text = xlab, side = 1, line = 1.5, cex = 1, outer = TRUE) } }