IsoMC v.2.3 beta

### IsoMC (v.2.3) 
## Isotope mixing model with 95% interquantile intervals, based on Monte Carlo simulations (20.000 runs)
# (R. Schwamborn, 2009)

# Inputs:
# 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. 
# TL: Trophic level (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)


# Outputs: 
# MF, isotopic signature for mixture, corrected for trophic fractionation MF = M + (F*(TL-1))
# Contrib.A =  contribution of source A to the mixture (ratio, zero to 1) 
# ConfUP, ConfLOW, Approximate upper and lower limits of the 95% confidence interval.

### INPUTS

M <- -25.54
A <- -28.42
B <- -20.63
sM <- 1.19
sA <- 1.03
sB <- 1.19
TL <- 2
F <- 1.6
sF <- 0.4


# Step 1. Mixing functions:

# Basic equations for sources A and B , without F

mix.funA <- function(M, A, B)
{ contrib.A <- (M-B) / (A-B)
  contrib.A }

mix.funB <- function(M, A, B)
{ contrib.B <- (M-A) / (B-A)
  contrib.B }

mix.funA(M, A, B)

## [1] 0.6302953

mix.funB(M, A, B)

## [1] 0.3697047

# 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 }


# Step 2. Monte Carlo simulation with 20000 runs:
 
iso.means <- c(M, A, B )
iso.sd <- c(sM, sA, sB)

iso.mix.vect <- rnorm( 20000 , mean = iso.means[1], sd = iso.sd[1])
iso.sA.vect <- rnorm( 20000 , mean = iso.means[2], sd = iso.sd[2])
iso.sB.vect <- rnorm( 20000 , mean = iso.means[3], sd = iso.sd[3])
iso.sF.vect <- rnorm( 20000 , mean = F, sd = sF)


# Contribution of Source A

iso.contribA.outp.vect <- mix.funFA(iso.mix.vect, iso.sA.vect, iso.sB.vect, iso.sF.vect, TL)

mean(iso.contribA.outp.vect)

## [1] 0.8497003

median(iso.contribA.outp.vect)

## [1] 0.8380125

quantile(iso.contribA.outp.vect, probs=c(.025,0.975))

##      2.5%     97.5% 
## 0.4638475 1.3103126

quantsA <- quantile(iso.contribA.outp.vect, probs=c(.025,0.975))
hist(iso.contribA.outp.vect, breaks=200, 
   main = "Contribution of Source A", sub= "red: median, blue: 95% quantiles", 
   xlab = "Source Contribution (zero to 1)")
abline(v= quantsA[1], col = "blue") 
abline(v= quantsA[2], col = "blue")
abline(v= median(iso.contribA.outp.vect), col = "red")

### OUTPUTS

# 1 graph in R: histograms of Monte Carlo output for source A with 95% quantiles.
# contribution of Source A, without considering trophic fractionation.
mix.funA(M, A, B) 

## [1] 0.6302953

# contribution of Source B, without considering trophic fractionation.
mix.funB(M, A, B)

## [1] 0.3697047

# contribution of Source A, considering F (trophic fractionation).
median(iso.contribA.outp.vect) 

## [1] 0.8380125

# 95% interquantile intervals for the contribution of Source A, considering F
quantile(iso.contribA.outp.vect, probs=c(.025,0.975))

##      2.5%     97.5% 
## 0.4638475 1.3103126

# 95% interquantile range for the contribution of Source A, considering F
(i.q.range <- quantile(iso.contribA.outp.vect, probs=c(.975))-quantile(iso.contribA.outp.vect, probs=c(0.025)))

##     97.5% 
## 0.8464652