--- title: "IsoMCv21_beta" author: "RS" date: "11 de abril de 2017" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ## R Markdown # This is an R Markdown document. Markdown is a simple formatting syntax for # authoring HTML, PDF, and MS Word documents. For more details on using R # Markdown see # http:rmarkdown.rstudio.com. # When you click the **Knit** button a document will be generated that includes # both content as well as the output of any embedded R code chunks within the # document. You can embed an R code chunk like this: ``` ### Monte Carlo Analysis of Isotope Mixing Models. ### Outputs: 95% confidence intervals of mangrove contributions for each species # (R. Schwamborn, December 2009) ### A. Single, One-Time Operations # 1. Import Data for Input Iso <- read.table("C:/Dokumente und Einstellungen/Ralf/Eigene Dateien/R Software/OutputdataforIsotopeMixingModelsBA.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE, stringsAsFactors = F) # 2. Define Mixing equations (Source A = mangrove) mix.funA <- function(mix, sA, sB) { contrib.A <- (mix-sB) / (sA-sB) contrib.A } mix.funB <- function(mix, sA, sB) { contrib.B <- (mix-sA) / (sB-sA) contrib.B } # 3. load fUtlities (for skewness and kurtosis) library(fUtilities) ### B. LOOP operations for(species.no in 4:57) { # 3. Read Data from Table ### INPUT : define species number species.no species.name <- names(Iso[(species.no)]) Iso.means <- c(Iso[1,(species.no)], Iso[2,(species.no)], Iso[3,(species.no)]) Iso.sd <- c(Iso[7,(species.no)], Iso[8,(species.no)], Iso[9,(species.no)]) Iso.mix.vect <- rnorm( 4000 , mean = Iso.means[1], sd = Iso.sd[1]) Iso.sA.vect <- rnorm( 4000 , mean = Iso.means[2], sd = Iso.sd[2]) Iso.sB.vect <- rnorm( 4000 , mean = Iso.means[3], sd = Iso.sd[3]) # Contribution of Source A = Mangrove Contribution ContribA.outp.vect <- mix.funA(Iso.mix.vect, Iso.sA.vect, Iso.sB.vect) mean(ContribA.outp.vect) median(ContribA.outp.vect) quantile(ContribA.outp.vect, probs=c(0.5, .025, 0.975)) # 95% percentiles quantile(ContribA.outp.vect, probs=c(0.5, 0.159, .841)) # 68.2% percentiles (approx. equiv. to SD) quantsA <- quantile(ContribA.outp.vect, probs=c(0.5, .025,0.975, 0.159, .841)) 0.5 +(0.682/2) # Upper 68.2% percentile 0.5 -(0.682/2) # Lower 68.2% percentile hist(ContribA.outp.vect, breaks=200, freq = F, main = species.name, sub = "red: median, blue: lower and upper 95% quantiles", xlab = "Mangrove Contribution (zero to 1)") abline(v= quantsA[2], col = "blue") abline(v= quantsA[3], col = "blue") abline(v= median(ContribA.outp.vect), col = "red") text(quantile(ContribA.outp.vect, probs=c(.02)), 400, "Mangrove Contribution (zero to 1)", col = "blue") text(quantile(ContribA.outp.vect, probs=c(.02)), 350, "95% confidence intervals:") text(quantile(ContribA.outp.vect, probs=c(.02)), 300, quantsA[1]) text(quantile(ContribA.outp.vect, probs=c(.02)), 250, "to") text(quantile(ContribA.outp.vect, probs=c(.02)), 200, quantsA[2]) text(quantile(ContribA.outp.vect, probs=c(.02)), 150, "median:") text(quantile(ContribA.outp.vect, probs=c(.02)), 100, median(ContribA.outp.vect)) x.axis <- seq(-2, 2, length=200) points(x.axis, 40 * (dnorm(x.axis, mean= mean(ContribA.outp.vect), sd= sd(ContribA.outp.vect))), col = "red") abline(h=0, col="gray") # 4. Attach Quantiles and Median to output matrix # 4.1 generate output,1 column, with 95 Conf. int. results out.mat1 <- quantsA (name.mat <- data.frame(cbind(species.name))) (out.mat2 <- data.frame(name.mat, out.mat1)) out.mat2[2] out.mat3 <- out.mat2[2] out.mat3[1,1] out.mat3[2,1] out.mat3[3,1] out.mat3[4,1] out.mat3[5,1] # 4.2 generate output,1 column, with skewness and kurtosis, output of the Shapiro-Wilks-test for normality ContribA.skew <- skewness(ContribA.outp.vect) ContribA.kurt <- kurtosis(ContribA.outp.vect) S.W.test.output <- shapiro.test(ContribA.outp.vect) S.W.test.p.value <- S.W.test.output[2] out.mat4 <- data.frame(ContribA.skew[1], ContribA.kurt[1], S.W.test.p.value[1]) out.mat4 # 4.3 concatenate output to full Results table, x columns, with 95 Conf. int. results, skewness and kurtosis Iso[10, species.no] Iso[11, species.no] Iso[12, species.no] Iso[10, species.no] <- (out.mat3[1,1]) Iso[11, species.no] <- (out.mat3[2,1]) Iso[12, species.no] <- (out.mat3[3,1]) #Iso[10, 1] <- "median contribution of Source A" #Iso[11, 1] <- c("lower 95% percentiles for contribution of Source A") #Iso[12, 1] <- c("upper 95% percentiles for contribution of Source A") Iso[13, species.no] Iso[14, species.no] Iso[15, species.no] Iso[13, species.no] <- (out.mat4[1,1]) Iso[14, species.no] <- (out.mat4[1,2]) Iso[15, species.no] <- (out.mat4[1,3]) #Iso[13, 1] <- c("skewness of distribution for contribution of Source A") #Iso[14, 1] <- c("skewness of distribution for contribution of Source A") #Iso[15, 1] <- c("p-value of th Shapiro-Wilks-test") Iso[16, species.no] <- (out.mat3[4,1]) Iso[17, species.no] <- (out.mat3[5,1]) #Iso[16, 1] <- c("lower 68.2% percentiles for contribution of Source A") #Iso[17, 1] <- c("upper 68.2% percentiles for contribution of Source A") # 4.3b. New Row Names (2nd Column of Data Matrix): Iso[2] <- factor(c( "mean isotopic signature for mixture (e.g., d13C)", #row 1# "mean isotopic signature for source A", #row 2# "mean isotopic signature for source B", #row 3# "number of mixture samples analyzed", #row 4# "number of source A samples analyzed", #row 5# "number of source B samples analyzed", #row 6# "standard deviation of signatures from samples of M", #row 7# "standard deviation of signatures from samples of A", #row 8# "standard deviation of signatures from samples of B", #row 9# "median contribution of Source A", #row 10 (new)# "lower 95% percentiles for contribution of Source A", #row 11 (new)# "upper 95% percentiles for contribution of Source A", #row 12 (new)# "skewness of distribution for contribution of Source A", #row 13 (new)# "skewness of distribution for contribution of Source A", #row 14 (new)# "p-value of th Shapiro-Wilks-test", #row 15 (new)# "lower 68.2% percentiles for contribution of Source A", #row 16 (new)# "upper 68.2% percentiles for contribution of Source A" #row 17 (new)# )) Iso[species.no] quants } Iso # 5. Write output file write.table(Iso, file = "C:/Dokumente und Einstellungen/Ralf/Eigene Dateien/R Software/OUTPUT-IsoMC21.csv", sep=",", row.names = T) ### tests for normality: shapiro.test(ContribA.outp.vect) shapiro.test(Iso.sA.vect) skewness(Iso.sA.vect) kurtosis(Iso.sA.vect) shapiro.test(ContribA.outp.vect) skewness(ContribA.outp.vect) kurtosis(ContribA.outp.vect) ``` ## Including Plots You can also embed plots, for example: ```{r pressure, echo=FALSE} plot(pressure) ``` Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.