## ----echo = FALSE------------------------------------------------------------- knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 5) ## ----------------------------------------------------------------------------- # remove previously loaded items from the current environment and remove previous graphics. rm(list=ls()) graphics.off() # Here, I set the seed each time so that the results are comparable. # This is useful as it means that anyone that runs your code, *should* # get the same results as you, although random number generators change # from time to time. set.seed(1) # load SIBER library(SIBER) library(viridis) # set a new three-colour palette from the viridis package palette(viridis::viridis(3)) # load in the included demonstration dataset data("demo.siber.data") # # create the siber object siber.example <- createSiberObject(demo.siber.data) # Or if working with your own data read in from a *.csv file, you would use # This *.csv file is included with this package. To find its location # type # fname <- system.file("extdata", "demo.siber.data.csv", package = "SIBER") # in your command window. You could load it directly by using the # returned path, or perhaps better, you could navigate to this folder # and copy this file to a folder of your own choice, and create a # script from this vingette to analyse it. This *.csv file provides # a template for how your own files should be formatted. # mydata <- read.csv(fname, header=T) # siber.example <- createSiberObject(mydata) # Create lists of plotting arguments to be passed onwards to the # plotting functions. With p.interval = NULL, these are SEA. NB not SEAc though # which is what we will base our overlap calculations on. This implementation # needs to be added in a future update. For now, the best way to plot SEAc is to # add the ellipses manually following the vignette on this topic. group.ellipses.args <- list(n = 100, p.interval = NULL, lty = 1, lwd = 2) par(mfrow=c(1,1)) plotSiberObject(siber.example, ax.pad = 2, hulls = F, community.hulls.args, ellipses = T, group.ellipses.args, group.hulls = F, group.hull.args, bty = "L", iso.order = c(1,2), xlab = expression({delta}^13*C~'permille'), ylab = expression({delta}^15*N~'permille') ) ## ----MLoverlap---------------------------------------------------------------- # In this example, I will calculate the overlap between ellipses for groups 2 # and 3 in community 1 (i.e. the green and yellow open circles of data). # The first ellipse is referenced using a character string representation where # in "x.y", "x" is the community, and "y" is the group within that community. # So in this example: community 1, group 2 ellipse1 <- "1.2" # Ellipse two is similarly defined: community 1, group3 ellipse2 <- "1.3" # The overlap of the maximum likelihood fitted standard ellipses are # estimated using sea.overlap <- maxLikOverlap(ellipse1, ellipse2, siber.example, p.interval = NULL, n = 100) # the overlap betweeen the corresponding 95% prediction ellipses is given by: ellipse95.overlap <- maxLikOverlap(ellipse1, ellipse2, siber.example, p.interval = 0.95, n = 100) # so in this case, the overlap as a proportion of the non-overlapping area of # the two ellipses, would be prop.95.over <- ellipse95.overlap[3] / (ellipse95.overlap[2] + ellipse95.overlap[1] - ellipse95.overlap[3]) ## ----bayesOverlap------------------------------------------------------------- # options for running jags parms <- list() parms$n.iter <- 2 * 10^4 # number of iterations to run the model for parms$n.burnin <- 1 * 10^3 # discard the first set of values parms$n.thin <- 10 # thin the posterior by this many parms$n.chains <- 2 # run this many chains # define the priors priors <- list() priors$R <- 1 * diag(2) priors$k <- 2 priors$tau.mu <- 1.0E-3 # fit the ellipses which uses an Inverse Wishart prior # on the covariance matrix Sigma, and a vague normal prior on the # means. Fitting is via the JAGS method. ellipses.posterior <- siberMVN(siber.example, parms, priors) # and teh corresponding Bayesian estimates for the overlap between the # 95% ellipses is given by: bayes95.overlap <- bayesianOverlap(ellipse1, ellipse2, ellipses.posterior, draws = 100, p.interval = 0.95, n = 100) # a histogram of the overlap hist(bayes95.overlap[,3], 10) # and as above, you can express this a proportion of the non-overlapping area of # the two ellipses, would be bayes.prop.95.over <- (bayes95.overlap[,3] / (bayes95.overlap[,2] + bayes95.overlap[,1] - bayes95.overlap[,3]) ) hist(bayes.prop.95.over, 10)