R/iso.demo.R

In rogerclarkgc/estableiso0.99: Stable isotope analysis which using euclideanmetric model

iso.demo <- function(){
  cat("*********************DEMO*********************\n\n")
  cat("This is a demo using data from siar package\n")
  cat("The model used in this R script is simpler than siar, but also effective\n")
  cat("Press <Enter> to continue\n ")
  readline()
  invisible()
  geese <- matrix(c(10.22, 10.37, 10.44, 10.52, 10.19, 
                    10.45, 9.91, 11.27, 9.34, -11.36, -11.88, -10.6, -11.25, 
                    -11.66, -10.41, -10.88, -14.73, -11.52), ncol = 2, nrow = 9)
  colnames(geese) <- c("d15NPl", "d13CPl")
  geese <- as.data.frame(geese)
  sources <- data.frame(d15NPl = c(6.50, 4.42, 11.19, 9.82), d13CPl = c(-11.17, -30.88, -10.19, -15.01), sd.N = c(1.4594632, 2.2680709, 1.1124385, 0.8271039), 
                        mean.N = c(6.488984, 4.432160, 11.192613, 9.816280), sd.C = c(1.2149562, 0.6413182, 1.9593306, 1.1724677), mean.C = c(-11.17023, -30.87984, -11.17090, -14.05701))
  rownames(sources) <- c("Zostera", "Grass", "U.lactuca", "Enteromorpha")
  cat("The target mixture isotope data is from siar package, which is called geese1demo\n")
  cat("and must has following format:\n")
  print(geese)
  cat("Press <Enter> to see the data\n ")
  readline()
  invisible()
  cat("The source isotope data is from siar package, which is called sourcesdemo\n")
  cat("and must has following format:\n")
  print(sources)
  cat("Press <Enter> to continue\n")
  readline()
  invisible()
  cat("These data can be loaded by typing iso.menu() at the Console \n")
  cat("Your also can use function iso.loaddata() to load your data\n")
  cat("Press <Enter> to see the source and target plot\n")
  readline()
  invisible()
  cat("Using funtion iso.plot() to draw the plot.\n")
  cat("Producing plot..... done\n")
  iso.plot(datalist = list(sources = sources, mixture = geese))
  cat("Press <Enter> to run the model and see the result>")
  readline()
  invisible()
  choose <- menu(c("Using Montel Carlo sampling method only(faster).", "Using Montel Carlo Markov Chain sampling method and probability threshold(better robust, but slowly)"), title = "Choose a method")
  if(choose == 1){
    datalist <- list(sources = sources, mixture = apply(geese, 2, mean))
    start <- Sys.time()
    cat("---step.1 Creating random sample list---\n")
    sample.list <- iso.euclidean2(source.matrix = datalist$sources, mixture.matrix = datalist$mixture)
    cat("---step.1 finished!---\n")
    cat("---step.2 Creating resample list---\n")
    mcm.table <- iso.mcm(source.list = sample.list, numbers = length(sample.list))
    cat("---step.2 finished---\n")
    cat("---step.3 Drawing histograms...---\n")
    iso.histograms(mcm.table, numbers = length(sample.list))
    end <- Sys.time()
    cat("All steps done!\n")
    cat("Time use: \n")
    print(end - start)
    cat("\nThis demo is finished\n")
    cat("Thanks for using\n")
    cat("Press <Enter> to continue\n")
    readline()
    invisible()
  }
  if(choose == 2){
    start <- Sys.time()
    cat("---step.1 Creating random sample list and set the probability threshold---\n")
    threshold.list <- iso.threshold(source.matrix = sources, mixture.matrix = apply(geese, 2, mean), simplemode = FALSE, correctiso.matrix = NULL)
    cat("---step.1 finished!---\n")
    cat("---step.2 Creating resample list use MCMC method and probability threshold\n")
    mcmc.list <- iso.mcmc(threshold.result = threshold.list, numbers = length(threshold.list$sample_list), iter = 1000)
    apportion.table <- iso.apportion(mcmc.result = mcmc.list, run = 1000)
    cat("---step.2 finished---\n")
    cat("---step.3 Drawing histograms...---\n")
    iso.histograms(apportion.table, numbers = length(threshold.list$sample_list))
    end <- Sys.time()
    cat("All steps done!\n")
    cat("Time use: ")
    print(end - start)
    cat("Press <Enter> to continue\n")
    readline()
    invisible()
    
  }
}