rm(list=ls())

graphics.off()
library(siar)
# ------------------------------------------------------------------------------

# set the random seed to get the same set of numbers each time 
# this script is run.
set.seed(1)

# create some simulated data in m groups

n <- 10     # number of observations per group
m <- 2      # number of groups

dC <- NULL           # create object dC into which to store data
grp <- NULL          # create object for storing grouping information

mudC <- c(2,6)    # means of the carbon isotopes for the m groups

# loop over the m groups and create the dC values
for (i in 1:m){

  dC <- c(dC,rnorm(n,mudC[i],1))
  grp <- c(grp,rep(i,n))
}          


# bundle the consumer data into a matrix
consumers <- cbind(grp,dC)

# bundle the source data and create a dataframe of the source names,
# means and standard deviations.
sources <-  data.frame(list(item=c("A","B"),meandC=c(-5,10)),sddC=c(1,1),stringsAsFactors=F)


# now run the new model
# NB in this simple example there are no TEF correction factors and we are
# ignoring concentration dependence. See "SIAR for ecologists" for more
# information on these features.
model1 <- siarmcmcdirichletv4(consumers,sources)

# plot the data manually
# (note siarplotdata doesnt work on single isotopes)
dev.new()
plot( consumers[,2], rep(0,nrow(consumers)) , 
        col=consumers[,1], xlim=c(-10,15),
        axes=F,xlab="dC",ylab="",cex=1) # suppress the axes plotting
axis(1,at=seq(-10,15,5))
points(as.numeric(sources[,2]),rep(0,nrow(sources)),col=1:m,pch=15,cex=2)
legend("topright",c("Source 1","Source 2","Consumer 1","Consumer 2"),
         pch=c(15,15,1,1), col=c(1,2,1,2),pt.cex=c(2,2,1,1))

# plot the results for group 1
siarproportionbygroupplot(model1,grp=1)

# plot the results for group 2
siarproportionbygroupplot(model1,grp=2)