# IsoBoot version 3.2 BETA
### Bootstrap resampling analysis of stable isotope mixing models with 95% confidence intervals.
## Designed for Food Web Analysis. Inputs include trophic level and isotopic fractionation of consumers.
# (Copyright: R. Schwamborn, 2009)
# INSTRUCTIONS:
# 1. Download and install "R" (www.r-project.org), then open this textfile in any editor.
# 2. Write your data into "INPUT VALUES", copy the whole text and paste it into the "R" Console.
### INPUT:
# M: mean isotopic signature for mixture (e.g. 13C of a consumer),
# A: mean isotopic signature for source A (e.g. 13C of mangroves),
# B: mean isotopic signature for source B (e.g. 13C of algae),
# sM: standard deviation of signatures from samples of M,
# sA: standard deviation of signatures from samples of A,
# sB: standard deviation of signatures from samples of B,
# nM: n of samples of mixture M,
# nA: n of samples of source A,
# nB: n of samples of source B,
# nF: n of samples used to calculate F,
### Note: unequal n allowed.
# TL: Trophic level of consumer M (TL = 2 for herbivores, 3 for carnivores, zero if unknown),
# F: Isotopic fractionation per trophic level (zero if unknown),
# sF: Standard deviation of F (zero if unknown),
# E: Measurement precision (zero or 0.3 if unknown),
# nMC: Length of the original Monte Carlo simulation matrix (default: 20000 random numbers per variable),
# NBoot: Number of bootstrap resampling runs (default: 999 runs).
### OUTPUT:
# Percent contribution of source A consumer M.
# Upper and lower limits of the 95% confidence interval for this source contribution, based bootstrap resampling.
# Graph with the frequency distribution of the resampled contributions for source A.
### INPUT VALUES:
M <- (-25.54)
A <- (-28.42)
B <- (-20.63)
sM <- (1.19)
sA <- (1.03)
sB <- (1.19)
nM <- (3) # n of samples for M
nA <- (3) # n of samples for source A
nB <- (3) # n of samples for source B
nF <- (2) # n of samples used for F (default: 2)
TL <- (2) # trophic level of M (default: 2)
F <- (0) # Isotopic fractionation per TL (default: zero or 1.6)
sF <- (0) # s.d. of F (default: 0)
E <- (0.3) # measurement precision (default: 0.3)
nMC <- (20000) # Length of the original Monte Carlo simulation matrix.
NBoot <- (999) # Number of bootstrap resampling runs.
### NOTES ###
# Note (nF): For nF, default is 2. Suggestion: use nF = 2 when F was obtained from a pair of single data, or n-1, when F was obtained from a regression of 13C vs TL.
# Note (E): It is not necessary to include a value for measurement precision E. Is is not used in the basic calculations, since it is implicit in the variability of the measured variables. E is used to check for the relative effects of sources of variability only.
### Calculations ###
# 1. General Mixing Function: Calculates Exact Contribution of Source A to Mixture M ("ContribA", zero to one):
# 1.a Defines the mixing function.
mix.funA <- function(M, A, B)
{ contrib.A <- (M-B) / (A-B)
contrib.A }
# Contribution of source A to mixture M, no correction for F:
mix.funA(M, A, B)
## [1] 0.6302953
# 1.b Correction for F and trophic level
# Equation for source A, considering F:
mix.funFA <- function(M, A, B, F, TL)
{ contrib.A <- ((M +((-F) *(TL-1)))-B) / (A-B)
contrib.A }
# Contribution of source A to mixture M, with correction for F:
mix.funFA(M, A, B, F, TL)
## [1] 0.6302953
# 2. Generates four Monte Carlo vectors with random numbers from original data,
# using Mean and Standard Devs of samples.
# length = 20,000 random numbers
M.vec <- rnorm(nMC, M, (sM))
A.vec <- rnorm(nMC, A, (sA))
B.vec <- rnorm(nMC, B, (sB))
F.vec <- rnorm(nMC, mean = F, sd = sF)
# hist(F.vec)
data.input <- data.frame(F.vec, M.vec, A.vec, B.vec)
#############
# 3. Resampling using the sample() function
## 3.a Define the sampling function
resample.fun <- function(d) {
M.sample <- sample(d$M.vec, nM)
Mean.M.sample <- mean(M.sample)
A.sample <- sample(d$A.vec, nA)
Mean.A.sample <- mean(A.sample)
B.sample <- sample(d$B.vec, nB)
Mean.B.sample <- mean(B.sample)
F.sample <- sample(d$F.vec, nB)
Mean.F.sample <- mean(F.sample)
outp.sample.mat <- data.frame(Mean.M.sample, Mean.A.sample, Mean.B.sample, Mean.F.sample)
}
(subsample1 <- resample.fun(data.input))
## Mean.M.sample Mean.A.sample Mean.B.sample Mean.F.sample
## 1 -26.52027 -28.74843 -20.00138 0
#############
## 3.b Defines the function that generates a single value for ContribA from single averages for M, A, B, F
mix.funFA((subsample1$Mean.M.sample),(subsample1$Mean.A.sample),(subsample1$Mean.B.sample), (subsample1$Mean.F.sample), TL)
## [1] 0.7452678
##############
## 3.c Construct a single function out of both, that takes n samples for each variable and generates a single value for ContribA from three means.
sample.and.calculate.fun <- function(d)
{
resample.and.mean.fun <- function(d) {
M.sample <- sample(d$M.vec, nM)
Mean.M.sample <- mean(M.sample)
A.sample <- sample(d$A.vec, nA)
Mean.A.sample <- mean(A.sample)
B.sample <- sample(d$B.vec, nB)
Mean.B.sample <- mean(B.sample)
F.sample <- sample(d$F.vec, nB)
Mean.F.sample <- mean(F.sample)
outp.sample.mat <- data.frame(Mean.M.sample, Mean.A.sample, Mean.B.sample, Mean.F.sample)
}
## take samples, calculate means
subsample1 <- resample.and.mean.fun(data.input)
## calculate a contribA value with this sample
mix.funFA((subsample1$Mean.M.sample),(subsample1$Mean.A.sample),(subsample1$Mean.B.sample), (subsample1$Mean.F.sample), TL)
}
sample.and.calculate.fun(data.input)
## [1] 0.4807765
# 4. Resampling: Apply this function N times, taking n samples out of each column, using the initial normally distributed variables (n = 20000)
N.rep <- NBoot # N of repaeted resamplings
outp.mat1 <- data.frame(c(1:(N.rep)),rep.int(0,(N.rep)))
names(outp.mat1) <- c("sample.no", "contribA")
summary(outp.mat1)
## sample.no contribA
## Min. : 1.0 Min. :0
## 1st Qu.:250.5 1st Qu.:0
## Median :500.0 Median :0
## Mean :500.0 Mean :0
## 3rd Qu.:749.5 3rd Qu.:0
## Max. :999.0 Max. :0
for(line.no in 1:(N.rep)) {
outp.mat1$contribA[line.no] <- sample.and.calculate.fun(data.input)
}
# 5. Outputs
## 5.a: Mean and Median of the resampling distribution
summary(outp.mat1)
## sample.no contribA
## Min. : 1.0 Min. :0.2880
## 1st Qu.:250.5 1st Qu.:0.5517
## Median :500.0 Median :0.6292
## Mean :500.0 Mean :0.6296
## 3rd Qu.:749.5 3rd Qu.:0.7077
## Max. :999.0 Max. :1.0109
mean(outp.mat1$contribA)
## [1] 0.6296151
median(outp.mat1$contribA)
## [1] 0.6291652
## 5.b: Makes a histogram of the Bootstrap Resampling Distribution
hist(outp.mat1$contribA, breaks = 50, main = "Bootstrap Resampling Distribution",
sub = "lines: median and 95% quantiles",
xlab = "Mangrove Contribution (zero to 1)",
xlim = c(0,1),
axes = F )
axis(1, at = c(seq(0, 1, by = 0.1)))
abline(v = median(outp.mat1$contribA), lty = 3)
abline(v = quantile(outp.mat1$contribA, probs=c(.025)))
abline(v = quantile(outp.mat1$contribA, probs=c(.975)))
## 5.c: 95% Percentiles of resampling distribution - Approx. 95% Conf Int.
round(quantile(outp.mat1$contribA, probs=c(.025,0.975)), 2)
## 2.5% 97.5%
## 0.42 0.85
#########################