From 935f390bc392b0eee83d8ed0a4b49560875c9ffa Mon Sep 17 00:00:00 2001 From: Andrew Parnell Date: Wed, 20 Jan 2021 12:19:24 +0000 Subject: [PATCH] Fixed styler and prepare for CRAN --- DESCRIPTION | 7 +- NEWS.md | 3 +- R/combine_sources.R | 180 ++--- R/compare_groups.R | 274 ++++---- R/compare_sources.R | 248 +++---- R/data.R | 20 +- R/plot.simmr_input.R | 311 +++++---- R/plot.simmr_output.R | 239 ++++--- R/posterior_predictive.simmr_output.R | 175 ++--- R/print.simmr_input.R | 28 +- R/print.simmr_output.R | 14 +- R/prior_viz.simmr_output.R | 164 ++--- R/simmr.R | 55 +- R/simmr_elicit.R | 150 ++-- R/simmr_load.R | 127 ++-- R/simmr_mcmc.R | 487 ++++++------- R/simmr_mcmc_tdf.R | 304 +++++---- R/summary.simmr_output.R | 120 ++-- R/summary.simmr_output_tdf.R | 59 +- cran-comments.md | 4 +- docs/404.html | 8 +- docs/articles/index.html | 8 +- docs/articles/quick_start.html | 124 ++-- .../figure-html/unnamed-chunk-10-1.png | Bin 96414 -> 73626 bytes .../figure-html/unnamed-chunk-11-1.png | Bin 128606 -> 100723 bytes .../figure-html/unnamed-chunk-12-1.png | Bin 73658 -> 46346 bytes .../figure-html/unnamed-chunk-13-1.png | Bin 85657 -> 44821 bytes docs/articles/simmr.html | 640 ++++++++++-------- .../figure-html/unnamed-chunk-10-1.png | Bin 81013 -> 56528 bytes .../figure-html/unnamed-chunk-11-1.png | Bin 127416 -> 101025 bytes .../figure-html/unnamed-chunk-13-1.png | Bin 79475 -> 53471 bytes .../figure-html/unnamed-chunk-14-1.png | Bin 119628 -> 67907 bytes .../figure-html/unnamed-chunk-15-1.png | Bin 58190 -> 34927 bytes .../figure-html/unnamed-chunk-16-1.png | Bin 61125 -> 37675 bytes .../figure-html/unnamed-chunk-22-1.png | Bin 63262 -> 37935 bytes .../figure-html/unnamed-chunk-22-2.png | Bin 81416 -> 53770 bytes .../figure-html/unnamed-chunk-22-3.png | Bin 113213 -> 61304 bytes .../figure-html/unnamed-chunk-23-1.png | Bin 60126 -> 35586 bytes .../figure-html/unnamed-chunk-24-1.png | Bin 60453 -> 35126 bytes .../figure-html/unnamed-chunk-25-1.png | Bin 60655 -> 37392 bytes .../figure-html/unnamed-chunk-26-1.png | Bin 57924 -> 34027 bytes .../figure-html/unnamed-chunk-26-2.png | Bin 86220 -> 43492 bytes .../figure-html/unnamed-chunk-26-3.png | Bin 61018 -> 36138 bytes .../figure-html/unnamed-chunk-31-1.png | Bin 60944 -> 38625 bytes .../figure-html/unnamed-chunk-38-1.png | Bin 126931 -> 100114 bytes .../figure-html/unnamed-chunk-46-1.png | Bin 67015 -> 39605 bytes .../figure-html/unnamed-chunk-47-1.png | Bin 63262 -> 37935 bytes .../figure-html/unnamed-chunk-47-2.png | Bin 65606 -> 41593 bytes .../figure-html/unnamed-chunk-48-1.png | Bin 67142 -> 46183 bytes docs/authors.html | 8 +- docs/index.html | 14 +- docs/news/index.html | 16 +- docs/pkgdown.css | 4 +- docs/pkgdown.yml | 4 +- docs/reference/Rplot001.png | Bin 0 -> 1011 bytes docs/reference/combine_sources.html | 91 +-- docs/reference/compare_groups.html | 83 ++- docs/reference/compare_sources.html | 81 ++- docs/reference/geese_data.html | 14 +- docs/reference/geese_data_day1.html | 14 +- docs/reference/index.html | 8 +- docs/reference/plot.simmr_input.html | 134 ++-- docs/reference/plot.simmr_output.html | 78 ++- docs/reference/posterior_predictive.html | 45 +- docs/reference/print.simmr_input.html | 10 +- docs/reference/print.simmr_output.html | 10 +- docs/reference/prior_viz.html | 63 +- docs/reference/simmr.html | 50 +- docs/reference/simmr_data_1.html | 16 +- docs/reference/simmr_data_2.html | 16 +- docs/reference/simmr_elicit.html | 70 +- docs/reference/simmr_load.html | 59 +- docs/reference/simmr_mcmc.html | 276 ++++---- docs/reference/simmr_mcmc_tdf.html | 123 ++-- docs/reference/square_data.html | 18 +- docs/reference/summary.simmr_output.html | 64 +- docs/reference/summary.simmr_output_tdf.html | 25 +- man/combine_sources.Rd | 48 +- man/compare_groups.Rd | 48 +- man/compare_sources.Rd | 40 +- man/geese_data.Rd | 2 +- man/geese_data_day1.Rd | 2 +- man/plot.simmr_input.Rd | 78 ++- man/plot.simmr_output.Rd | 30 +- man/posterior_predictive.Rd | 24 +- man/prior_viz.Rd | 24 +- man/simmr.Rd | 28 +- man/simmr_data_1.Rd | 4 +- man/simmr_data_2.Rd | 4 +- man/simmr_elicit.Rd | 36 +- man/simmr_load.Rd | 21 +- man/simmr_mcmc.Rd | 197 +++--- man/simmr_mcmc_tdf.Rd | 86 +-- man/square_data.Rd | 6 +- man/summary.simmr_output.Rd | 35 +- tests/testthat/test_simmr_combine_sources.R | 245 ++++--- .../test_simmr_compare_sources_and_groups.R | 225 +++--- tests/testthat/test_simmr_input_plots.R | 123 ++-- tests/testthat/test_simmr_load.R | 78 ++- tests/testthat/test_simmr_mcmc_simple.R | 96 +-- tests/testthat/test_simmr_misc_functions.R | 104 +-- tests/testthat/test_simmr_plot_output.R | 49 +- tests/testthat/test_simmr_print_summary.R | 100 +-- tests/testthat/test_simmr_tdf_estimation.R | 87 +-- 104 files changed, 3790 insertions(+), 3173 deletions(-) create mode 100644 docs/reference/Rplot001.png diff --git a/DESCRIPTION b/DESCRIPTION index 9b26a56..3510dbd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,10 +1,11 @@ Package: simmr Type: Package Title: A Stable Isotope Mixing Model -Version: 0.4.3.9000 -Date: 2020-11-10 +Version: 0.4.4 +Date: 2021-01-21 Author: Andrew Parnell -URL: https://github.com/andrewcparnell/simmr +URL: https://github.com/andrewcparnell/simmr, http://andrewcparnell.github.io/simmr +BugReports: https://github.com/andrewcparnell/simmr/issues Language: en-US Maintainer: Andrew Parnell Description: Fits Stable Isotope Mixing Models (SIMMs) and is meant as a longer term replacement to the previous widely-used package SIAR. SIMMs are used to infer dietary proportions of organisms consuming various food sources from observations on the stable isotope values taken from the organisms' tissue samples. However SIMMs can also be used in other scenarios, such as in sediment mixing or the composition of fatty acids. The main functions are simmr_load and simmr_mcmc. The two vignettes contain a quick start and a full listing of all the features. The methods used are detailed in the papers Parnell et al 2010 , and Parnell et al 2013 . diff --git a/NEWS.md b/NEWS.md index 3b6559d..d806ef8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,7 +1,8 @@ -# simmr 0.4.3.9000 +# simmr 0.4.4 - Updated the simmr_elicit function to provide a more explicit warning for bad input objects - Updated compare_sources and compare_groups to allow exporting of the plot object for editing purposes + - Used styler to correct code style # simmr 0.4.3 diff --git a/R/combine_sources.R b/R/combine_sources.R index 6a2332f..e71ad4d 100644 --- a/R/combine_sources.R +++ b/R/combine_sources.R @@ -31,14 +31,18 @@ #' data(geese_data) #' #' # Load into simmr -#' simmr_1 = with(geese_data_day1, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) +#' simmr_1 <- with( +#' geese_data_day1, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) #' #' # Plot #' plot(simmr_1) @@ -47,117 +51,121 @@ #' simmr_1 #' #' # MCMC run -#' simmr_1_out = simmr_mcmc(simmr_1) +#' simmr_1_out <- simmr_mcmc(simmr_1) #' #' # Print it #' print(simmr_1_out) #' #' # Summary #' summary(simmr_1_out) -#' summary(simmr_1_out,type='diagnostics') -#' summary(simmr_1_out,type='correlations') -#' summary(simmr_1_out,type='statistics') -#' ans = summary(simmr_1_out,type=c('quantiles','statistics')) +#' summary(simmr_1_out, type = "diagnostics") +#' summary(simmr_1_out, type = "correlations") +#' summary(simmr_1_out, type = "statistics") +#' ans <- summary(simmr_1_out, type = c("quantiles", "statistics")) #' #' # Plot #' plot(simmr_1_out) -#' plot(simmr_1_out,type='boxplot') -#' plot(simmr_1_out,type='histogram') -#' plot(simmr_1_out,type='density') -#' plot(simmr_1_out,type='matrix') +#' plot(simmr_1_out, type = "boxplot") +#' plot(simmr_1_out, type = "histogram") +#' plot(simmr_1_out, type = "density") +#' plot(simmr_1_out, type = "matrix") #' -#' simmr_out_combine = combine_sources(simmr_1_out, -#' to_combine=c('Source A','Source D'), -#' new_source_name='Source A+D') +#' simmr_out_combine <- combine_sources(simmr_1_out, +#' to_combine = c("Source A", "Source D"), +#' new_source_name = "Source A+D" +#' ) #' plot(simmr_out_combine$input) -#' plot(simmr_out_combine,type='boxplot',title='simmr output: combined sources') -#' +#' plot(simmr_out_combine, type = "boxplot", title = "simmr output: combined sources") #' } #' -#' @export -combine_sources = function(simmr_out, - to_combine = simmr_out$input$source_names[1:2], - new_source_name = 'combined_source') { - UseMethod('combine_sources') +#' @export +combine_sources <- function(simmr_out, + to_combine = simmr_out$input$source_names[1:2], + new_source_name = "combined_source") { + UseMethod("combine_sources") } -#' @export -combine_sources.simmr_output = function(simmr_out, - to_combine = simmr_out$input$source_names[1:2], - new_source_name = 'combined_source') { +#' @export +combine_sources.simmr_output <- function(simmr_out, + to_combine = simmr_out$input$source_names[1:2], + new_source_name = "combined_source") { # A posteriori combining of sources - + # Check only two sources to be combined - if (length(to_combine) != 2) + if (length(to_combine) != 2) { stop("Currently only two sources can be combined") - + } + # # Check class # if (class(simmr_out) != 'simmr_output') # stop("Only objects of class simmr_output can be run through this function") - + # Find which columns to combine by number - to_combine_cols = sort(match(to_combine, simmr_out$input$source_names)) - if (any(is.na(to_combine_cols))) - stop('1 or more source names not found') - - simmr_new_out = simmr_out - + to_combine_cols <- sort(match(to_combine, simmr_out$input$source_names)) + if (any(is.na(to_combine_cols))) { + stop("1 or more source names not found") + } + + simmr_new_out <- simmr_out + # 1 combine the chosen source means - old_source_means = simmr_out$input$source_means - simmr_new_out$input$source_means = old_source_means[-to_combine_cols[2], , - drop = FALSE] - simmr_new_out$input$source_means[to_combine_cols[1], ] = apply(old_source_means[to_combine_cols, ,drop = FALSE], 2, 'mean') - + old_source_means <- simmr_out$input$source_means + simmr_new_out$input$source_means <- old_source_means[-to_combine_cols[2], , + drop = FALSE + ] + simmr_new_out$input$source_means[to_combine_cols[1], ] <- apply(old_source_means[to_combine_cols, , drop = FALSE], 2, "mean") + # 2 combine the source sds - old_source_sds = simmr_out$input$source_sds - simmr_new_out$input$source_sds = old_source_sds[-to_combine_cols[2], ,drop= FALSE] - simmr_new_out$input$source_sds[to_combine_cols[1], ] = apply(old_source_sds[to_combine_cols, , drop = FALSE], 2, function(x) - sqrt(sum(x ^ 2))) - + old_source_sds <- simmr_out$input$source_sds + simmr_new_out$input$source_sds <- old_source_sds[-to_combine_cols[2], , drop = FALSE] + simmr_new_out$input$source_sds[to_combine_cols[1], ] <- apply(old_source_sds[to_combine_cols, , drop = FALSE], 2, function(x) { + sqrt(sum(x^2)) + }) + # 3 combine the correction means - old_correction_means = simmr_out$input$correction_means - simmr_new_out$input$correction_means = old_correction_means[-to_combine_cols[2], ,drop = FALSE] - simmr_new_out$input$correction_means[to_combine_cols[1], ] = apply(old_correction_means[to_combine_cols, ,drop = FALSE], 2, 'mean') - + old_correction_means <- simmr_out$input$correction_means + simmr_new_out$input$correction_means <- old_correction_means[-to_combine_cols[2], , drop = FALSE] + simmr_new_out$input$correction_means[to_combine_cols[1], ] <- apply(old_correction_means[to_combine_cols, , drop = FALSE], 2, "mean") + # 4 combine the correction sds - old_correction_sds = simmr_out$input$correction_sds - simmr_new_out$input$correction_sds = old_correction_sds[-to_combine_cols[2], ,drop = FALSE] - simmr_new_out$input$correction_sds[to_combine_cols[1], ] = apply(old_correction_sds[to_combine_cols, , drop = FALSE], 2, function(x) - sqrt(sum(x ^ 2))) - + old_correction_sds <- simmr_out$input$correction_sds + simmr_new_out$input$correction_sds <- old_correction_sds[-to_combine_cols[2], , drop = FALSE] + simmr_new_out$input$correction_sds[to_combine_cols[1], ] <- apply(old_correction_sds[to_combine_cols, , drop = FALSE], 2, function(x) { + sqrt(sum(x^2)) + }) + # 5 combine the concentraion means - old_concentration_means = simmr_out$input$concentration_means - simmr_new_out$input$concentration_means = old_concentration_means[-to_combine_cols[2], ,drop = FALSE] - simmr_new_out$input$concentration_means[to_combine_cols[1], ] = apply(old_concentration_means[to_combine_cols, ,drop = FALSE], 2, 'mean') - + old_concentration_means <- simmr_out$input$concentration_means + simmr_new_out$input$concentration_means <- old_concentration_means[-to_combine_cols[2], , drop = FALSE] + simmr_new_out$input$concentration_means[to_combine_cols[1], ] <- apply(old_concentration_means[to_combine_cols, , drop = FALSE], 2, "mean") + # 6 change the source names - old_source_names = simmr_out$input$source_names - simmr_new_out$input$source_names = old_source_names[-to_combine_cols[2]] - simmr_new_out$input$source_names[to_combine_cols[1]] = new_source_name - + old_source_names <- simmr_out$input$source_names + simmr_new_out$input$source_names <- old_source_names[-to_combine_cols[2]] + simmr_new_out$input$source_names[to_combine_cols[1]] <- new_source_name + # 7 Change n_sources - simmr_new_out$input$n_sources = simmr_new_out$input$n_sources - 1 - + simmr_new_out$input$n_sources <- simmr_new_out$input$n_sources - 1 + # 8 Sum across all the output values - n_groups = simmr_out$input$n_groups + n_groups <- simmr_out$input$n_groups for (j in 1:n_groups) { - simmr_new_out$output[[j]] = simmr_out$output[[j]] + simmr_new_out$output[[j]] <- simmr_out$output[[j]] # Change sims.list and sims.matrix # First sims.matrix - sims.matrix = simmr_out$output[[j]]$BUGSoutput$sims.matrix - new_sims.matrix = sims.matrix[,-(to_combine_cols[2]+1)] - new_sims.matrix[,to_combine_cols[1]+1] = sims.matrix[,to_combine_cols[1]+1] + sims.matrix[,to_combine_cols[2]+1] - colnames(new_sims.matrix)[to_combine_cols[1]+1] = new_source_name - simmr_new_out$output[[j]]$BUGSoutput$sims.matrix = new_sims.matrix + sims.matrix <- simmr_out$output[[j]]$BUGSoutput$sims.matrix + new_sims.matrix <- sims.matrix[, -(to_combine_cols[2] + 1)] + new_sims.matrix[, to_combine_cols[1] + 1] <- sims.matrix[, to_combine_cols[1] + 1] + sims.matrix[, to_combine_cols[2] + 1] + colnames(new_sims.matrix)[to_combine_cols[1] + 1] <- new_source_name + simmr_new_out$output[[j]]$BUGSoutput$sims.matrix <- new_sims.matrix # Now sims.list - sims.list = simmr_out$output[[j]]$BUGSoutput$sims.list - new_sims.list = sims.list - new_sims.list$p = sims.list$p[,-to_combine_cols[2]] - new_sims.list$p[,to_combine_cols[1]] = sims.list$p[,to_combine_cols[2]] + sims.list$p[,to_combine_cols[1]] - colnames(new_sims.list$p)[to_combine_cols[1]] = new_source_name - simmr_new_out$output[[j]]$BUGSoutput$sims.list = new_sims.list + sims.list <- simmr_out$output[[j]]$BUGSoutput$sims.list + new_sims.list <- sims.list + new_sims.list$p <- sims.list$p[, -to_combine_cols[2]] + new_sims.list$p[, to_combine_cols[1]] <- sims.list$p[, to_combine_cols[2]] + sims.list$p[, to_combine_cols[1]] + colnames(new_sims.list$p)[to_combine_cols[1]] <- new_source_name + simmr_new_out$output[[j]]$BUGSoutput$sims.list <- new_sims.list } - - class(simmr_new_out) = 'simmr_output' + + class(simmr_new_out) <- "simmr_output" return(simmr_new_out) - } diff --git a/R/compare_groups.R b/R/compare_groups.R index a3d1b4f..edeb4a1 100644 --- a/R/compare_groups.R +++ b/R/compare_groups.R @@ -1,16 +1,16 @@ #' Compare dietary proportions for a single source across different groups -#' +#' #' This function takes in an object of class \code{simmr_output} and creates #' probabilistic comparisons for a given source and a set of at least two #' groups. -#' +#' #' When two groups are specified, the function produces a direct calculation of #' the probability that one group is bigger than the other. When more than two #' groups are given, the function produces a set of most likely probabilistic #' orderings for each combination of groups. The function produces boxplots by #' default and also allows for the storage of the output for further analysis if #' required. -#' +#' #' @param simmr_out An object of class \code{simmr_output} created from #' \code{\link{simmr_mcmc}}. #' @param source_name The name of a source. This should match the names exactly @@ -19,9 +19,9 @@ #' least two groups must be specified. #' @param plot A logical value specifying whether plots should be produced or #' not. -#' +#' #' @import ggplot2 -#' +#' #' @return If there are two groups, a vector containing the differences between #' the two groups proportions for that source. If there are multiple groups, a #' list containing the following fields: \item{Ordering }{The different @@ -33,153 +33,161 @@ #' @examples #' \dontrun{ #' data(geese_data) -#' simmr_in = with(geese_data, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means, -#' group = groups)) -#' +#' simmr_in <- with( +#' geese_data, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means, +#' group = groups +#' ) +#' ) +#' #' # Print #' simmr_in -#' +#' #' # Plot -#' plot(simmr_in,group=1:8,xlab=expression(paste(delta^13, "C (\u2030)",sep="")), -#' ylab=expression(paste(delta^15, "N (\u2030)",sep="")), -#' title='Isospace plot of Inger et al Geese data') -#' +#' plot(simmr_in, +#' group = 1:8, xlab = expression(paste(delta^13, "C (\\u2030)", sep = "")), +#' ylab = expression(paste(delta^15, "N (\\u2030)", sep = "")), +#' title = "Isospace plot of Inger et al Geese data" +#' ) +#' #' # Run MCMC for each group -#' simmr_out = simmr_mcmc(simmr_in) -#' +#' simmr_out <- simmr_mcmc(simmr_in) +#' #' # Print output #' simmr_out -#' +#' #' # Summarise output -#' summary(simmr_out,type='quantiles',group=1) -#' summary(simmr_out,type='quantiles',group=c(1,3)) -#' summary(simmr_out,type=c('quantiles','statistics'),group=c(1,3)) -#' +#' summary(simmr_out, type = "quantiles", group = 1) +#' summary(simmr_out, type = "quantiles", group = c(1, 3)) +#' summary(simmr_out, type = c("quantiles", "statistics"), group = c(1, 3)) +#' #' # Plot - only a single group allowed -#' plot(simmr_out,type='boxplot',group=2,title='simmr output group 2') -#' plot(simmr_out,type=c('density','matrix'),grp=6,title='simmr output group 6') -#' +#' plot(simmr_out, type = "boxplot", group = 2, title = "simmr output group 2") +#' plot(simmr_out, type = c("density", "matrix"), grp = 6, title = "simmr output group 6") +#' #' # Compare groups -#' compare_groups(simmr_out,source='Zostera',groups=1:2) -#' compare_groups(simmr_out,source='Zostera',groups=1:3) -#' compare_groups(simmr_out,source='U.lactuca',groups=c(4:5,7,2)) +#' compare_groups(simmr_out, source = "Zostera", groups = 1:2) +#' compare_groups(simmr_out, source = "Zostera", groups = 1:3) +#' compare_groups(simmr_out, source = "U.lactuca", groups = c(4:5, 7, 2)) #' } -#' +#' #' @export -compare_groups = function(simmr_out, - source_name = simmr_out$input$source_names[1], - groups = 1:2, - plot = TRUE) { - UseMethod('compare_groups') +compare_groups <- function(simmr_out, + source_name = simmr_out$input$source_names[1], + groups = 1:2, + plot = TRUE) { + UseMethod("compare_groups") } #' @export -compare_groups.simmr_output = function(simmr_out, - source_name = simmr_out$input$source_names[1], - groups = 1:2, - plot = TRUE) { +compare_groups.simmr_output <- function(simmr_out, + source_name = simmr_out$input$source_names[1], + groups = 1:2, + plot = TRUE) { -# Function to compare between groups both via textual output and with boxplots -# Things to supply are: -# If two groups are given: -# - provide the probability of one group being bigger than the other -# - give the probability distribution of the difference -# - optional boxplot of two -# If more than two groups are given: -# - provide the top most likely orderings of the groups -# An optional boxplot of the groups - -# Throw an error if only one group is specified -if(length(groups)==1) stop("Please use plot(...) or summary(...) if you just want to look at one group.") - -# Throw an error if the source name given doesn't match the source names -if(!source_name%in%simmr_out$input$source_names) stop("This source name not found in the current source names. Be sure to check case and spelling") - -# Get group names -group_names = levels(simmr_out$input$group)[groups] - -# Start with two groups version -if(length(groups)==2) { - # Get the output for this particular source on these two groups - post_mat1 = simmr_out$output[[groups[1]]]$BUGSoutput$sims.matrix - post_mat2 = simmr_out$output[[groups[2]]]$BUGSoutput$sims.matrix - match_name = match(source_name, colnames(post_mat1)) - out_all_grp_1 = post_mat1[,match_name] - out_all_grp_2 = post_mat2[,match_name] - # Produce the difference between the two - out_diff = out_all_grp_1 - out_all_grp_2 + # Function to compare between groups both via textual output and with boxplots + # Things to supply are: + # If two groups are given: + # - provide the probability of one group being bigger than the other + # - give the probability distribution of the difference + # - optional boxplot of two + # If more than two groups are given: + # - provide the top most likely orderings of the groups + # An optional boxplot of the groups - cat(paste("Prob ( proportion of",source_name,'in group',group_names[1],'> proportion of',source_name,'in group',group_names[2],') =',round(mean(out_diff>0),3))) - - if(plot) { - # Stupid fix for packaging ggplot things - Group = Proportion = NULL - df = data.frame(Proportion=c(out_all_grp_1,out_all_grp_2),Group=c(rep(paste(group_names[1]),length(out_all_grp_1)),rep(paste(group_names[2]),length(out_all_grp_2)))) - p = ggplot(df,aes(x=Group,y=Proportion,fill=Group)) + geom_boxplot(alpha=0.5,outlier.size=0) + theme_bw() + theme(legend.position='none') + ggtitle(paste("Comparison of dietary proportions for groups",group_names[1],'and',group_names[2],'for source',source_name)) - print(p) - } - -} + # Throw an error if only one group is specified + if (length(groups) == 1) stop("Please use plot(...) or summary(...) if you just want to look at one group.") -# Now for more groups -if(length(groups)>2) { - # Get the output for all the groups - post_mat = simmr_out$output[[groups[1]]]$BUGSoutput$sims.matrix - match_name = match(source_name, colnames(post_mat)) - len = length(post_mat[,match_name]) - out_all = matrix(NA,nrow=len,ncol=length(groups)) - for(j in 1:length(groups)) out_all[,j] = simmr_out$output[[groups[j]]]$BUGSoutput$sims.matrix[,match_name] - colnames(out_all) = paste(group_names) - - # Now find the ordering of each one - ordering_num = t(apply(out_all,1,order,decreasing=TRUE)) - Ordering = rep(NA,length=nrow(ordering_num)) - for(i in 1:length(Ordering)) Ordering[i] = paste0(group_names[ordering_num[i,]],collapse=" > ") - cat('Most popular orderings are as follows:\n') - tab = t(t(sort(table(Ordering,dnn=NULL),decreasing=TRUE))) - colnames(tab) = 'Probability' - # Do not print all of it if too long - if(nrow(tab)>30) { - print(round(tab[1:30,]/length(Ordering),4)) - } else { - print(round(tab/length(Ordering),4)) - } - - if(plot) { - # Stupid fix for packaging ggplot things - Group = Proportion = NULL - df = reshape2::melt(out_all)[,2:3] - colnames(df) = c('Group','Proportion') - p = ggplot(df,aes(x=factor(Group),y=Proportion,fill=factor(Group))) + - scale_fill_viridis(discrete=TRUE) + - geom_boxplot(alpha=0.5,outlier.size=0) + - xlab('Group') + - theme_bw() + theme(legend.position='none') + - ggtitle(paste("Comparison of dietary proportions for source",source_name)) - print(p) + # Throw an error if the source name given doesn't match the source names + if (!source_name %in% simmr_out$input$source_names) stop("This source name not found in the current source names. Be sure to check case and spelling") + + # Get group names + group_names <- levels(simmr_out$input$group)[groups] + + # Start with two groups version + if (length(groups) == 2) { + # Get the output for this particular source on these two groups + post_mat1 <- simmr_out$output[[groups[1]]]$BUGSoutput$sims.matrix + post_mat2 <- simmr_out$output[[groups[2]]]$BUGSoutput$sims.matrix + match_name <- match(source_name, colnames(post_mat1)) + out_all_grp_1 <- post_mat1[, match_name] + out_all_grp_2 <- post_mat2[, match_name] + # Produce the difference between the two + out_diff <- out_all_grp_1 - out_all_grp_2 + + cat(paste("Prob ( proportion of", source_name, "in group", group_names[1], "> proportion of", source_name, "in group", group_names[2], ") =", round(mean(out_diff > 0), 3))) + + if (plot) { + # Stupid fix for packaging ggplot things + Group <- Proportion <- NULL + df <- data.frame(Proportion = c(out_all_grp_1, out_all_grp_2), Group = c(rep(paste(group_names[1]), length(out_all_grp_1)), rep(paste(group_names[2]), length(out_all_grp_2)))) + p <- ggplot(df, aes(x = Group, y = Proportion, fill = Group)) + + geom_boxplot(alpha = 0.5, outlier.size = 0) + + theme_bw() + + theme(legend.position = "none") + + ggtitle(paste("Comparison of dietary proportions for groups", group_names[1], "and", group_names[2], "for source", source_name)) + print(p) + } } - -} -# Return output -if(length(groups)==2) { - if(plot) { - invisible(list(out_diff = out_diff, plot = p)) - } else { - invisible(list(out_diff = out_diff)) + # Now for more groups + if (length(groups) > 2) { + # Get the output for all the groups + post_mat <- simmr_out$output[[groups[1]]]$BUGSoutput$sims.matrix + match_name <- match(source_name, colnames(post_mat)) + len <- length(post_mat[, match_name]) + out_all <- matrix(NA, nrow = len, ncol = length(groups)) + for (j in 1:length(groups)) out_all[, j] <- simmr_out$output[[groups[j]]]$BUGSoutput$sims.matrix[, match_name] + colnames(out_all) <- paste(group_names) + + # Now find the ordering of each one + ordering_num <- t(apply(out_all, 1, order, decreasing = TRUE)) + Ordering <- rep(NA, length = nrow(ordering_num)) + for (i in 1:length(Ordering)) Ordering[i] <- paste0(group_names[ordering_num[i, ]], collapse = " > ") + cat("Most popular orderings are as follows:\n") + tab <- t(t(sort(table(Ordering, dnn = NULL), decreasing = TRUE))) + colnames(tab) <- "Probability" + # Do not print all of it if too long + if (nrow(tab) > 30) { + print(round(tab[1:30, ] / length(Ordering), 4)) + } else { + print(round(tab / length(Ordering), 4)) + } + + if (plot) { + # Stupid fix for packaging ggplot things + Group <- Proportion <- NULL + df <- reshape2::melt(out_all)[, 2:3] + colnames(df) <- c("Group", "Proportion") + p <- ggplot(df, aes(x = factor(Group), y = Proportion, fill = factor(Group))) + + scale_fill_viridis(discrete = TRUE) + + geom_boxplot(alpha = 0.5, outlier.size = 0) + + xlab("Group") + + theme_bw() + + theme(legend.position = "none") + + ggtitle(paste("Comparison of dietary proportions for source", source_name)) + print(p) + } } -} else { - if(plot) { - invisible(list(Ordering=Ordering,out_all=out_all, plot = p)) + + # Return output + if (length(groups) == 2) { + if (plot) { + invisible(list(out_diff = out_diff, plot = p)) + } else { + invisible(list(out_diff = out_diff)) + } } else { - invisible(list(Ordering=Ordering,out_all=out_all)) + if (plot) { + invisible(list(Ordering = Ordering, out_all = out_all, plot = p)) + } else { + invisible(list(Ordering = Ordering, out_all = out_all)) + } } -} - } diff --git a/R/compare_sources.R b/R/compare_sources.R index b24e049..4bffd2a 100644 --- a/R/compare_sources.R +++ b/R/compare_sources.R @@ -1,16 +1,16 @@ #' Compare dietary proportions between multiple sources -#' +#' #' This function takes in an object of class \code{simmr_output} and creates #' probabilistic comparisons between the supplied sources. The group number can #' also be specified. -#' +#' #' When two sources are specified, the function produces a direct calculation #' of the probability that the dietary proportion for one source is bigger than #' the other. When more than two sources are given, the function produces a set #' of most likely probabilistic orderings for each combination of sources. The #' function produces boxplots by default and also allows for the storage of the #' output for further analysis if required. -#' +#' #' @param simmr_out An object of class \code{simmr_output} created from #' \code{\link{simmr_mcmc}}. #' @param source_names The names of at least two sources. These should match @@ -19,10 +19,10 @@ #' specified assumes the first or only group #' @param plot A logical value specifying whether plots should be produced or #' not. -#' +#' #' @import ggplot2 #' @importFrom reshape2 "melt" -#' +#' #' @return If there are two sources, a vector containing the differences #' between the two dietary proportion proportions for these two sources. If #' there are multiple sources, a list containing the following fields: @@ -35,148 +35,154 @@ #' \dontrun{ #' # Data set: 10 obs on 2 isos, 4 sources, with tefs and concdep #' data(geese_data_day1) -#' simmr_1 = with(geese_data_day1, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) -#' +#' simmr_1 <- with( +#' geese_data_day1, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) +#' #' # Plot #' plot(simmr_1) -#' +#' #' # Print #' simmr_1 -#' +#' #' # MCMC run -#' simmr_1_out = simmr_mcmc(simmr_1) -#' +#' simmr_1_out <- simmr_mcmc(simmr_1) +#' #' # Print it #' print(simmr_1_out) -#' +#' #' # Summary #' summary(simmr_1_out) -#' summary(simmr_1_out,type='diagnostics') -#' summary(simmr_1_out,type='correlations') -#' summary(simmr_1_out,type='statistics') -#' ans = summary(simmr_1_out,type=c('quantiles','statistics')) -#' +#' summary(simmr_1_out, type = "diagnostics") +#' summary(simmr_1_out, type = "correlations") +#' summary(simmr_1_out, type = "statistics") +#' ans <- summary(simmr_1_out, type = c("quantiles", "statistics")) +#' #' # Plot -#' plot(simmr_1_out,type='boxplot') -#' plot(simmr_1_out,type='histogram') -#' plot(simmr_1_out,type='density') -#' plot(simmr_1_out,type='matrix') -#' +#' plot(simmr_1_out, type = "boxplot") +#' plot(simmr_1_out, type = "histogram") +#' plot(simmr_1_out, type = "density") +#' plot(simmr_1_out, type = "matrix") +#' #' # Compare two sources -#' compare_sources(simmr_1_out,source_names=c('Source B','Source D')) -#' +#' compare_sources(simmr_1_out, source_names = c("Source B", "Source D")) +#' #' # Compare multiple sources #' compare_sources(simmr_1_out) #' } -#' +#' #' @export -compare_sources = function(simmr_out, - source_names = simmr_out$input$source_names, - group = 1, - plot = TRUE) { - UseMethod('compare_sources') -} +compare_sources <- function(simmr_out, + source_names = simmr_out$input$source_names, + group = 1, + plot = TRUE) { + UseMethod("compare_sources") +} #' @export -compare_sources.simmr_output = function(simmr_out, - source_names = simmr_out$input$source_names, - group = 1, - plot = TRUE) { +compare_sources.simmr_output <- function(simmr_out, + source_names = simmr_out$input$source_names, + group = 1, + plot = TRUE) { -# Function to compare between sources within a group both via textual output and with boxplots -# Things to ly are: -# If two sources are given: -# - provide the probability of one group having higher dietary proportion than the other -# - give the probability distribution of the difference -# - optional boxplot of two -# If more than two sources are given: -# - provide the top most likely orderings of the sources -# An optional boxplot of the sources - -# Throw an error if only one group is specified -if(length(source_names)==1) stop("Use compare_between_groups if you want to compare a single source between groups.") - -# Throw an error if the source name given doesn't match the source names -if(!all(source_names%in%simmr_out$input$source_names)) stop("Some source names not found in the current source names. Be sure to check case and spelling") - -# Throw an error if more than one group specified -if(length(group) > 1) stop("Only one group allowed") - -# Start with two groups version -if(length(source_names)==2) { - # Get the output for this particular source on these two groups - match_names = match(source_names, simmr_out$input$source_names) - out_all_src_1 = simmr_out$output[[group]]$BUGSoutput$sims.list$p[,match_names[1]] - out_all_src_2 = simmr_out$output[[group]]$BUGSoutput$sims.list$p[,match_names[2]] - # Produce the difference between the two - out_diff = out_all_src_1 - out_all_src_2 - cat(paste("Prob ( proportion of",source_names[1],'> proportion of',source_names[2],') =',round(mean(out_diff>0),3))) - - if(plot) { - # Stupid fix for packaging ggplot things - Source = Proportion = NULL - df = data.frame(Proportion=c(out_all_src_1,out_all_src_2),Source=c(rep(source_names[1],length(out_all_src_1)),rep(source_names[2],length(out_all_src_2)))) - p = ggplot(df,aes(x=Source,y=Proportion,fill=Source)) + geom_boxplot(alpha=0.5,outlier.size=0) + theme_bw() + theme(legend.position='none') + ggtitle(paste("Comparison of dietary proportions for sources",source_names[1],'and',source_names[2])) - print(p) - } - -} + # Function to compare between sources within a group both via textual output and with boxplots + # Things to ly are: + # If two sources are given: + # - provide the probability of one group having higher dietary proportion than the other + # - give the probability distribution of the difference + # - optional boxplot of two + # If more than two sources are given: + # - provide the top most likely orderings of the sources + # An optional boxplot of the sources -# Now for more sources -if(length(source_names)>2) { - # Get the output for all the sources - match_names = match(source_names, simmr_out$input$source_names) - out_all = simmr_out$output[[group]]$BUGSoutput$sims.list$p[,match_names] - - # Now find the ordering of each one - ordering_num = t(apply(out_all,1,order,decreasing=TRUE)) - Ordering = rep(NA,length=nrow(ordering_num)) - for(i in 1:length(Ordering)) Ordering[i] = paste0(source_names[ordering_num[i,]],collapse=" > ") - if(simmr_out$input$n_groups > 1) cat("Results for group:", group, '\n') - cat('Most popular orderings are as follows:\n') - tab = t(t(sort(table(Ordering,dnn=NULL),decreasing=TRUE))) - colnames(tab) = 'Probability' - # Do not print all of it if too long - if(nrow(tab)>30) { - print(round(tab[1:30,]/length(Ordering),4)) - } else { - print(round(tab/length(Ordering),4)) + # Throw an error if only one group is specified + if (length(source_names) == 1) stop("Use compare_between_groups if you want to compare a single source between groups.") + + # Throw an error if the source name given doesn't match the source names + if (!all(source_names %in% simmr_out$input$source_names)) stop("Some source names not found in the current source names. Be sure to check case and spelling") + + # Throw an error if more than one group specified + if (length(group) > 1) stop("Only one group allowed") + + # Start with two groups version + if (length(source_names) == 2) { + # Get the output for this particular source on these two groups + match_names <- match(source_names, simmr_out$input$source_names) + out_all_src_1 <- simmr_out$output[[group]]$BUGSoutput$sims.list$p[, match_names[1]] + out_all_src_2 <- simmr_out$output[[group]]$BUGSoutput$sims.list$p[, match_names[2]] + # Produce the difference between the two + out_diff <- out_all_src_1 - out_all_src_2 + cat(paste("Prob ( proportion of", source_names[1], "> proportion of", source_names[2], ") =", round(mean(out_diff > 0), 3))) + + if (plot) { + # Stupid fix for packaging ggplot things + Source <- Proportion <- NULL + df <- data.frame(Proportion = c(out_all_src_1, out_all_src_2), Source = c(rep(source_names[1], length(out_all_src_1)), rep(source_names[2], length(out_all_src_2)))) + p <- ggplot(df, aes(x = Source, y = Proportion, fill = Source)) + + geom_boxplot(alpha = 0.5, outlier.size = 0) + + theme_bw() + + theme(legend.position = "none") + + ggtitle(paste("Comparison of dietary proportions for sources", source_names[1], "and", source_names[2])) + print(p) + } } - - if(plot) { - # Stupid fix for packaging ggplot things - Source = Proportion = NULL - df = reshape2::melt(out_all)[,2:3] - colnames(df) = c('Source','Proportion') - p = ggplot(df,aes(x=Source,y=Proportion,fill=Source)) + - scale_fill_viridis(discrete=TRUE) + - geom_boxplot(alpha=0.5,outlier.size=0) + - theme_bw() + - theme(legend.position='none') + - ggtitle(paste("Comparison of dietary proportions between sources")) - print(p) + + # Now for more sources + if (length(source_names) > 2) { + # Get the output for all the sources + match_names <- match(source_names, simmr_out$input$source_names) + out_all <- simmr_out$output[[group]]$BUGSoutput$sims.list$p[, match_names] + + # Now find the ordering of each one + ordering_num <- t(apply(out_all, 1, order, decreasing = TRUE)) + Ordering <- rep(NA, length = nrow(ordering_num)) + for (i in 1:length(Ordering)) Ordering[i] <- paste0(source_names[ordering_num[i, ]], collapse = " > ") + if (simmr_out$input$n_groups > 1) cat("Results for group:", group, "\n") + cat("Most popular orderings are as follows:\n") + tab <- t(t(sort(table(Ordering, dnn = NULL), decreasing = TRUE))) + colnames(tab) <- "Probability" + # Do not print all of it if too long + if (nrow(tab) > 30) { + print(round(tab[1:30, ] / length(Ordering), 4)) + } else { + print(round(tab / length(Ordering), 4)) + } + + if (plot) { + # Stupid fix for packaging ggplot things + Source <- Proportion <- NULL + df <- reshape2::melt(out_all)[, 2:3] + colnames(df) <- c("Source", "Proportion") + p <- ggplot(df, aes(x = Source, y = Proportion, fill = Source)) + + scale_fill_viridis(discrete = TRUE) + + geom_boxplot(alpha = 0.5, outlier.size = 0) + + theme_bw() + + theme(legend.position = "none") + + ggtitle(paste("Comparison of dietary proportions between sources")) + print(p) + } } - -} -# Return output - if(length(source_names)==2) { - if(plot) { + # Return output + if (length(source_names) == 2) { + if (plot) { invisible(list(out_diff = out_diff, plot = p)) } else { invisible(list(out_diff = out_diff)) } } else { - if(plot) { - invisible(list(Ordering=Ordering,out_all=out_all, plot = p)) + if (plot) { + invisible(list(Ordering = Ordering, out_all = out_all, plot = p)) } else { - invisible(list(Ordering=Ordering,out_all=out_all)) + invisible(list(Ordering = Ordering, out_all = out_all)) } - } + } } diff --git a/R/data.R b/R/data.R index e86fa3a..6a6e949 100644 --- a/R/data.R +++ b/R/data.R @@ -13,9 +13,9 @@ #' \item{correction_sds}{A matrix of TEFs sd values for the tracers in the same order as \code{mixtures} above} #' \item{concentration_means}{A matrix of concentration dependence mean values for the tracers in the same order as \code{mixtures} above} #' ... -#' +#' #' @seealso \code{\link{simmr_mcmc}} for an example where it is used -#' +#' #' } "simmr_data_1" @@ -34,9 +34,9 @@ #' \item{correction_sds}{A matrix of TEFs sd values for the tracers in the same order as \code{mixtures} above} #' \item{concentration_means}{A matrix of concentration dependence mean values for the tracers in the same order as \code{mixtures} above} #' ... -#' +#' #' @seealso \code{\link{simmr_mcmc}} for an example where it is used -#' +#' #' } "simmr_data_2" @@ -57,7 +57,7 @@ #' \item{correction_sds}{A matrix of TEFs sd values for the tracers in the same order as \code{mixtures} above} #' \item{concentration_means}{A matrix of concentration dependence mean values for the tracers in the same order as \code{mixtures} above} #' ... -#' +#' #' @seealso \code{\link{simmr_mcmc}} for an example where it is used #' } "geese_data" @@ -80,10 +80,10 @@ #' \item{correction_sds}{A matrix of TEFs sd values for the tracers in the same order as \code{mixtures} above} #' \item{concentration_means}{A matrix of concentration dependence mean values for the tracers in the same order as \code{mixtures} above} #' ... -#' +#' #' @seealso \code{\link{simmr_mcmc}} for an example where it is used -#' -#' +#' +#' #' } "square_data" @@ -104,7 +104,7 @@ #' \item{correction_sds}{A matrix of TEFs sd values for the tracers in the same order as \code{mixtures} above} #' \item{concentration_means}{A matrix of concentration dependence mean values for the tracers in the same order as \code{mixtures} above} #' ... -#' +#' #' @seealso \code{\link{simmr_mcmc}} for an example where it is used #' } -"geese_data_day1" \ No newline at end of file +"geese_data_day1" diff --git a/R/plot.simmr_input.R b/R/plot.simmr_input.R index d24c797..1847a4f 100644 --- a/R/plot.simmr_input.R +++ b/R/plot.simmr_input.R @@ -1,15 +1,15 @@ #' Plot the \code{simmr_input} data created from \code{simmr_load} -#' +#' #' This function creates iso-space (AKA tracer-space or delta-space) plots. #' They are vital in determining whether the data are suitable for running in a #' SIMM. -#' +#' #' It is desirable to have the vast majority of the mixture observations to be #' inside the convex hull defined by the food sources. When there are more than #' two tracers (as in one of the examples below) it is recommended to plot all #' the different pairs of the food sources. See the vignette for further #' details of richer plots. -#' +#' #' @param x An object created via the function \code{\link{simmr_load}} #' @param tracers The choice of tracers to plot. If there are more than two #' tracers, it is recommended to plot every pair of tracers to determine @@ -28,185 +28,198 @@ #' black and white #' @param ggargs Extra arguments to be included in the ggplot (e.g. axis limits) #' @param ... Not used -#' +#' #' @import ggplot2 #' @import viridis -#' +#' #' @author Andrew Parnell #' @seealso See \code{\link{plot.simmr_output}} for plotting the output of a #' simmr run. See \code{\link{simmr_mcmc}} for running a simmr object once the #' iso-space is deemed acceptable. #' @examples -#' +#' #' # A simple example with 10 observations, 4 food sources and 2 tracers #' data(geese_data_day1) -#' simmr_1 = with(geese_data_day1, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) -#' +#' simmr_1 <- with( +#' geese_data_day1, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) +#' #' # Plot #' plot(simmr_1) -#' +#' #' ### A more complicated example with 30 obs, 3 tracers and 4 sources #' data(simmr_data_2) -#' simmr_3 = with(simmr_data_2, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) -#' -#' # Plot 3 times - first default d13C vs d15N +#' simmr_3 <- with( +#' simmr_data_2, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) +#' +#' # Plot 3 times - first default d13C vs d15N #' plot(simmr_3) #' # Now plot d15N vs d34S -#' plot(simmr_3,tracers=c(2,3)) +#' plot(simmr_3, tracers = c(2, 3)) #' # and finally d13C vs d34S -#' plot(simmr_3,tracers=c(1,3)) +#' plot(simmr_3, tracers = c(1, 3)) #' # See vignette('simmr') for fancier x-axis labels -#' +#' #' # An example with multiple groups - the Geese data from Inger et al 2006 #' data(geese_data) -#' simmr_4 = with(geese_data, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means, -#' group = groups)) -#' +#' simmr_4 <- with( +#' geese_data, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means, +#' group = groups +#' ) +#' ) +#' #' # Print #' simmr_4 -#' +#' #' # Plot #' plot(simmr_4, -#' xlab=expression(paste(delta^13, "C (\u2030)",sep="")), -#' ylab=expression(paste(delta^15, "N (\u2030)",sep="")), -#' title='Isospace plot of Inger et al Geese data')#' -#' +#' xlab = expression(paste(delta^13, "C (\u2030)", sep = "")), +#' ylab = expression(paste(delta^15, "N (\u2030)", sep = "")), +#' title = "Isospace plot of Inger et al Geese data" +#' ) #' #' @export plot.simmr_input <- -function(x, - tracers = c(1,2), - title = 'Tracers plot', - xlab = colnames(x$mixtures)[tracers[1]], - ylab = colnames(x$mixtures)[tracers[2]], - sigmas = 1, - group = 1:x$n_groups, - mix_name = 'Mixtures', - ggargs = NULL, - colour = TRUE, - ...) { + function(x, + tracers = c(1, 2), + title = "Tracers plot", + xlab = colnames(x$mixtures)[tracers[1]], + ylab = colnames(x$mixtures)[tracers[2]], + sigmas = 1, + group = 1:x$n_groups, + mix_name = "Mixtures", + ggargs = NULL, + colour = TRUE, + ...) { -# Get mixtures to match current group(s) -curr_rows = which(x$group_int%in%group) -curr_mix = x$mixtures[curr_rows,,drop=FALSE] -curr_n_groups = length(group) + # Get mixtures to match current group(s) + curr_rows <- which(x$group_int %in% group) + curr_mix <- x$mixtures[curr_rows, , drop = FALSE] + curr_n_groups <- length(group) -# Throw error if too many groups (can only handle max 6 before it runs out of shapes) -#if((length(group)+x$n_sources)>6) stop("Too many groups specified. Total number of groups plus number of sources cannot exceed 6") + # Throw error if too many groups (can only handle max 6 before it runs out of shapes) + # if((length(group)+x$n_sources)>6) stop("Too many groups specified. Total number of groups plus number of sources cannot exceed 6") -# Throw error if plotting only one isotope -#if(ncol(curr_mix)==1) stop("This function only works for two or more tracers") + # Throw error if plotting only one isotope + # if(ncol(curr_mix)==1) stop("This function only works for two or more tracers") -# First get the mean corrected sources and the sd corrected sources -source_means_c = x$source_means + x$correction_means -source_sds_c = sqrt(x$source_sds^2 + x$correction_sds^2) + # First get the mean corrected sources and the sd corrected sources + source_means_c <- x$source_means + x$correction_means + source_sds_c <- sqrt(x$source_sds^2 + x$correction_sds^2) -# Set up data frame for ggplot - have to do it this stupid way because of cran -x2=c(source_means_c[,tracers[1]],curr_mix[,tracers[1]]) -x_lower=c(source_means_c[,tracers[1]]-sigmas*source_sds_c[,tracers[1]],curr_mix[,tracers[1]]) -x_upper=c(source_means_c[,tracers[1]]+sigmas*source_sds_c[,tracers[1]],curr_mix[,tracers[1]]) + # Set up data frame for ggplot - have to do it this stupid way because of cran + x2 <- c(source_means_c[, tracers[1]], curr_mix[, tracers[1]]) + x_lower <- c(source_means_c[, tracers[1]] - sigmas * source_sds_c[, tracers[1]], curr_mix[, tracers[1]]) + x_upper <- c(source_means_c[, tracers[1]] + sigmas * source_sds_c[, tracers[1]], curr_mix[, tracers[1]]) -if(ncol(curr_mix)>1) { - y=c(source_means_c[,tracers[2]],curr_mix[,tracers[2]]) - y_lower=c(source_means_c[,tracers[2]]-sigmas*source_sds_c[,tracers[2]],curr_mix[,tracers[2]]) - y_upper=c(source_means_c[,tracers[2]]+sigmas*source_sds_c[,tracers[2]],curr_mix[,tracers[2]]) -} + if (ncol(curr_mix) > 1) { + y <- c(source_means_c[, tracers[2]], curr_mix[, tracers[2]]) + y_lower <- c(source_means_c[, tracers[2]] - sigmas * source_sds_c[, tracers[2]], curr_mix[, tracers[2]]) + y_upper <- c(source_means_c[, tracers[2]] + sigmas * source_sds_c[, tracers[2]], curr_mix[, tracers[2]]) + } -if(x$n_groups==1) { - Source=factor(c(x$source_names,rep(mix_name,nrow(curr_mix))),levels=c(mix_name,x$source_names)) -} else { - Source=factor(c(x$source_names,paste(mix_name,x$group[curr_rows])),levels=c(paste(mix_name,unique(x$group[curr_rows])),x$source_names)) -} -size=c(rep(0.5,x$n_sources),rep(0.5,nrow(curr_mix))) + if (x$n_groups == 1) { + Source <- factor(c(x$source_names, rep(mix_name, nrow(curr_mix))), levels = c(mix_name, x$source_names)) + } else { + Source <- factor(c(x$source_names, paste(mix_name, x$group[curr_rows])), levels = c(paste(mix_name, unique(x$group[curr_rows])), x$source_names)) + } + size <- c(rep(0.5, x$n_sources), rep(0.5, nrow(curr_mix))) -if(ncol(curr_mix) == 1) { - df=data.frame(x=x2,x_lower,x_upper,Source,size, y = Source) -} else { - df=data.frame(x=x2,y=y,x_lower,y_lower,x_upper,y_upper,Source,size) -} + if (ncol(curr_mix) == 1) { + df <- data.frame(x = x2, x_lower, x_upper, Source, size, y = Source) + } else { + df <- data.frame(x = x2, y = y, x_lower, y_lower, x_upper, y_upper, Source, size) + } -# Plot for bivariate mixtures -if(ncol(curr_mix) > 1) { - if(colour) { - g=ggplot(data=df, aes(x = x,y = y,colour=Source)) + - scale_color_viridis(discrete=TRUE) + - theme_bw() + - labs(x=xlab,y=ylab,title=title) + - geom_errorbarh(aes(xmax=x_upper,xmin=x_lower,height=0)) + - #geom_pointrange(aes(x=x,y=y,ymax=y_upper,ymin=y_lower,height=0.2,shape=Source)) + - geom_pointrange(aes(x=x,y=y,ymax=y_upper,ymin=y_lower,shape=Source)) + - scale_shape_manual(values=1:nlevels(df$Source)) + - theme(legend.title=element_blank(),legend.key = element_blank()) + - guides(color=guide_legend(override.aes=list(linetype=c(rep(0,curr_n_groups),rep(1,x$n_sources))))) + - ggargs - } else { - g=ggplot(data=df, aes(x = x,y = y,colour=Source)) + - theme_bw() + - labs(x=xlab,y=ylab,title=title) + - geom_errorbarh(aes(xmax=x_upper,xmin=x_lower,height=0)) + - #geom_pointrange(aes(x=x,y=y,ymax=y_upper,ymin=y_lower,height=0.2,shape=Source)) + - geom_pointrange(aes(x=x,y=y,ymax=y_upper,ymin=y_lower,shape=Source)) + - scale_shape_manual(values=1:nlevels(df$Source)) + - theme(legend.title=element_blank(),legend.key = element_blank()) + - guides(color=guide_legend(override.aes=list(linetype=c(rep(0,curr_n_groups),rep(1,x$n_sources))))) + - scale_colour_grey() + - ggargs - } -} - -# Plot for univariate mixtures -if(ncol(curr_mix) == 1) { - if(colour) { - g = ggplot(data=df, aes(x = x, y = y, colour = Source)) + - scale_color_viridis(discrete=TRUE) + - theme_bw() + - theme(axis.title.y=element_blank()) + - labs(x = xlab, title = title) + - geom_errorbarh(aes(xmax=x_upper,xmin=x_lower,height=0)) + - #geom_pointrange(aes(x=x,y=y,ymax=y_upper,ymin=y_lower,height=0.2,shape=Source)) + - geom_point(aes(shape = Source)) + - scale_shape_manual(values=1:nlevels(df$Source)) + - theme(legend.position = 'None') + - guides(color=guide_legend(override.aes=list(linetype=c(rep(0,curr_n_groups),rep(1,x$n_sources))))) + - ggargs - } else { - g = ggplot(data=df, aes(x = x, y = y, colour = Source)) + - scale_color_grey() + - theme_bw() + - theme(axis.title.y=element_blank(), - axis.text.y=element_blank(), - axis.ticks.y=element_blank()) + - labs(x = xlab, title = title) + - geom_errorbarh(aes(xmax=x_upper,xmin=x_lower,height=0)) + - #geom_pointrange(aes(x=x,y=y,ymax=y_upper,ymin=y_lower,height=0.2,shape=Source)) + - geom_point(aes(shape = Source)) + - scale_shape_manual(values=1:nlevels(df$Source)) + - theme(legend.title=element_blank(),legend.key = element_blank()) + - guides(color=guide_legend(override.aes=list(linetype=c(rep(0,curr_n_groups),rep(1,x$n_sources))))) + - ggargs - } -} + # Plot for bivariate mixtures + if (ncol(curr_mix) > 1) { + if (colour) { + g <- ggplot(data = df, aes(x = x, y = y, colour = Source)) + + scale_color_viridis(discrete = TRUE) + + theme_bw() + + labs(x = xlab, y = ylab, title = title) + + geom_errorbarh(aes(xmax = x_upper, xmin = x_lower, height = 0)) + + # geom_pointrange(aes(x=x,y=y,ymax=y_upper,ymin=y_lower,height=0.2,shape=Source)) + + geom_pointrange(aes(x = x, y = y, ymax = y_upper, ymin = y_lower, shape = Source)) + + scale_shape_manual(values = 1:nlevels(df$Source)) + + theme(legend.title = element_blank(), legend.key = element_blank()) + + guides(color = guide_legend(override.aes = list(linetype = c(rep(0, curr_n_groups), rep(1, x$n_sources))))) + + ggargs + } else { + g <- ggplot(data = df, aes(x = x, y = y, colour = Source)) + + theme_bw() + + labs(x = xlab, y = ylab, title = title) + + geom_errorbarh(aes(xmax = x_upper, xmin = x_lower, height = 0)) + + # geom_pointrange(aes(x=x,y=y,ymax=y_upper,ymin=y_lower,height=0.2,shape=Source)) + + geom_pointrange(aes(x = x, y = y, ymax = y_upper, ymin = y_lower, shape = Source)) + + scale_shape_manual(values = 1:nlevels(df$Source)) + + theme(legend.title = element_blank(), legend.key = element_blank()) + + guides(color = guide_legend(override.aes = list(linetype = c(rep(0, curr_n_groups), rep(1, x$n_sources))))) + + scale_colour_grey() + + ggargs + } + } -invisible(g) + # Plot for univariate mixtures + if (ncol(curr_mix) == 1) { + if (colour) { + g <- ggplot(data = df, aes(x = x, y = y, colour = Source)) + + scale_color_viridis(discrete = TRUE) + + theme_bw() + + theme(axis.title.y = element_blank()) + + labs(x = xlab, title = title) + + geom_errorbarh(aes(xmax = x_upper, xmin = x_lower, height = 0)) + + # geom_pointrange(aes(x=x,y=y,ymax=y_upper,ymin=y_lower,height=0.2,shape=Source)) + + geom_point(aes(shape = Source)) + + scale_shape_manual(values = 1:nlevels(df$Source)) + + theme(legend.position = "None") + + guides(color = guide_legend(override.aes = list(linetype = c(rep(0, curr_n_groups), rep(1, x$n_sources))))) + + ggargs + } else { + g <- ggplot(data = df, aes(x = x, y = y, colour = Source)) + + scale_color_grey() + + theme_bw() + + theme( + axis.title.y = element_blank(), + axis.text.y = element_blank(), + axis.ticks.y = element_blank() + ) + + labs(x = xlab, title = title) + + geom_errorbarh(aes(xmax = x_upper, xmin = x_lower, height = 0)) + + # geom_pointrange(aes(x=x,y=y,ymax=y_upper,ymin=y_lower,height=0.2,shape=Source)) + + geom_point(aes(shape = Source)) + + scale_shape_manual(values = 1:nlevels(df$Source)) + + theme(legend.title = element_blank(), legend.key = element_blank()) + + guides(color = guide_legend(override.aes = list(linetype = c(rep(0, curr_n_groups), rep(1, x$n_sources))))) + + ggargs + } + } -} + invisible(g) + } diff --git a/R/plot.simmr_output.R b/R/plot.simmr_output.R index 611b648..5bd5011 100644 --- a/R/plot.simmr_output.R +++ b/R/plot.simmr_output.R @@ -41,136 +41,153 @@ #' data(geese_data) #' #' # Load into simmr -#' simmr_1 = with(geese_data_day1, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) +#' simmr_1 <- with( +#' geese_data_day1, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) #' # Plot #' plot(simmr_1) #' #' #' # MCMC run -#' simmr_1_out = simmr_mcmc(simmr_1) +#' simmr_1_out <- simmr_mcmc(simmr_1) #' #' # Plot #' plot(simmr_1_out) # Creates all 4 plots -#' plot(simmr_1_out,type='boxplot') -#' plot(simmr_1_out,type='histogram') -#' plot(simmr_1_out,type='density') -#' plot(simmr_1_out,type='matrix') +#' plot(simmr_1_out, type = "boxplot") +#' plot(simmr_1_out, type = "histogram") +#' plot(simmr_1_out, type = "density") +#' plot(simmr_1_out, type = "matrix") #' } #' @export plot.simmr_output <- -function(x, - type = c('isospace', - 'histogram', - 'density', - 'matrix', - 'boxplot'), - group = 1, - binwidth = 0.05, - alpha = 0.5, - title = if(length(group)==1){ 'simmr output plot'} else {paste('simmr output plot: group',group)}, - ggargs = NULL, - ...) { + function(x, + type = c( + "isospace", + "histogram", + "density", + "matrix", + "boxplot" + ), + group = 1, + binwidth = 0.05, + alpha = 0.5, + title = if (length(group) == 1) { + "simmr output plot" + } else { + paste("simmr output plot: group", group) + }, + ggargs = NULL, + ...) { - # Get the specified type - type=match.arg(type,several.ok=TRUE) + # Get the specified type + type <- match.arg(type, several.ok = TRUE) - # Iso-space plot is special as all groups go on one plot - # Add in extra dots here as they can be sent to this plot function - if('isospace' %in% type) graphics::plot(x$input,group=group,title=title,...) + # Iso-space plot is special as all groups go on one plot + # Add in extra dots here as they can be sent to this plot function + if ("isospace" %in% type) graphics::plot(x$input, group = group, title = title, ...) - # Get group names - group_names = levels(x$input$group)[group] - - for(i in 1:length(group)) { + # Get group names + group_names <- levels(x$input$group)[group] - # Prep data - out_all = x$output[[group[i]]]$BUGSoutput$sims.list$p - colnames(out_all) = x$input$source_names - df = reshape2::melt(out_all) - colnames(df) = c('Num','Source','Proportion') - if ('histogram'%in%type) { - g=ggplot(df,aes_string(x="Proportion",y="..density..", - fill="Source")) + - scale_fill_viridis(discrete=TRUE) + - geom_histogram(binwidth=binwidth,alpha=alpha) + - theme_bw() + - ggtitle(title[i]) + - facet_wrap("~ Source") + - theme(legend.position='none') + - ggargs - print(g) - } - - if ('density'%in%type) { - g=ggplot(df,aes_string(x="Proportion",y="..density..", - fill="Source")) + - scale_fill_viridis(discrete=TRUE) + - geom_density(alpha=alpha,linetype=0) + - theme_bw() + - theme(legend.position='none') + - ggtitle(title[i]) + - ylab("Density") + - facet_wrap("~ Source") + - ggargs - print(g) - } + for (i in 1:length(group)) { - if ('boxplot'%in%type) { - g=ggplot(df,aes_string(y="Proportion",x="Source", - fill="Source",alpha="alpha")) + - scale_fill_viridis(discrete=TRUE) + - geom_boxplot(alpha=alpha,notch=TRUE,outlier.size=0) + - theme_bw() + - ggtitle(title[i]) + - theme(legend.position='none') + - coord_flip() + - ggargs - print(g) - } - - # if ('convergence'%in%type) { - # coda::gelman.plot(x$output[[group[i]]],transform=TRUE) - # } + # Prep data + out_all <- x$output[[group[i]]]$BUGSoutput$sims.list$p + colnames(out_all) <- x$input$source_names + df <- reshape2::melt(out_all) + colnames(df) <- c("Num", "Source", "Proportion") + if ("histogram" %in% type) { + g <- ggplot(df, aes_string( + x = "Proportion", y = "..density..", + fill = "Source" + )) + + scale_fill_viridis(discrete = TRUE) + + geom_histogram(binwidth = binwidth, alpha = alpha) + + theme_bw() + + ggtitle(title[i]) + + facet_wrap("~ Source") + + theme(legend.position = "none") + + ggargs + print(g) + } - if ('matrix'%in%type) { - # These taken from the help(pairs) file - panel.hist <- function(x, ...) { - usr <- graphics::par("usr"); on.exit(graphics::par(usr)) - graphics::par(usr = c(usr[1:2], 0, 1.5) ) - h <- graphics::hist(x, plot = FALSE) - breaks <- h$breaks; nB <- length(breaks) - y <- h$counts; y <- y/max(y) - graphics::rect(breaks[-nB], 0, breaks[-1], y, col = "lightblue", ...) + if ("density" %in% type) { + g <- ggplot(df, aes_string( + x = "Proportion", y = "..density..", + fill = "Source" + )) + + scale_fill_viridis(discrete = TRUE) + + geom_density(alpha = alpha, linetype = 0) + + theme_bw() + + theme(legend.position = "none") + + ggtitle(title[i]) + + ylab("Density") + + facet_wrap("~ Source") + + ggargs + print(g) } - panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) - { - usr <- graphics::par("usr"); on.exit(graphics::par(usr)) - graphics::par(usr = c(0, 1, 0, 1)) - r <- stats::cor(x, y) - txt <- format(c(r, 0.123456789), digits = digits)[1] - txt <- paste0(prefix, txt) - if(missing(cex.cor)) cex.cor <- 0.8/graphics::strwidth(txt) - graphics::text(0.5, 0.5, txt, cex = cex.cor * abs(r)) + + if ("boxplot" %in% type) { + g <- ggplot(df, aes_string( + y = "Proportion", x = "Source", + fill = "Source", alpha = "alpha" + )) + + scale_fill_viridis(discrete = TRUE) + + geom_boxplot(alpha = alpha, notch = TRUE, outlier.size = 0) + + theme_bw() + + ggtitle(title[i]) + + theme(legend.position = "none") + + coord_flip() + + ggargs + print(g) } - panel.contour <- function(x, y, ...) - { - usr <- graphics::par("usr"); on.exit(graphics::par(usr)) - graphics::par(usr = c(usr[1:2], 0, 1.5) ) - kd <- MASS::kde2d(x,y) - kdmax <- max(kd$z) - graphics::contour(kd,add=TRUE,drawlabels=FALSE,levels=c(kdmax*0.1,kdmax*0.25,kdmax*0.5,kdmax*0.75,kdmax*0.9)) + + # if ('convergence'%in%type) { + # coda::gelman.plot(x$output[[group[i]]],transform=TRUE) + # } + + if ("matrix" %in% type) { + # These taken from the help(pairs) file + panel.hist <- function(x, ...) { + usr <- graphics::par("usr") + on.exit(graphics::par(usr)) + graphics::par(usr = c(usr[1:2], 0, 1.5)) + h <- graphics::hist(x, plot = FALSE) + breaks <- h$breaks + nB <- length(breaks) + y <- h$counts + y <- y / max(y) + graphics::rect(breaks[-nB], 0, breaks[-1], y, col = "lightblue", ...) + } + panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) { + usr <- graphics::par("usr") + on.exit(graphics::par(usr)) + graphics::par(usr = c(0, 1, 0, 1)) + r <- stats::cor(x, y) + txt <- format(c(r, 0.123456789), digits = digits)[1] + txt <- paste0(prefix, txt) + if (missing(cex.cor)) cex.cor <- 0.8 / graphics::strwidth(txt) + graphics::text(0.5, 0.5, txt, cex = cex.cor * abs(r)) + } + panel.contour <- function(x, y, ...) { + usr <- graphics::par("usr") + on.exit(graphics::par(usr)) + graphics::par(usr = c(usr[1:2], 0, 1.5)) + kd <- MASS::kde2d(x, y) + kdmax <- max(kd$z) + graphics::contour(kd, add = TRUE, drawlabels = FALSE, levels = c(kdmax * 0.1, kdmax * 0.25, kdmax * 0.5, kdmax * 0.75, kdmax * 0.9)) + } + graphics::pairs(out_all, xlim = c(0, 1), ylim = c(0, 1), main = title[i], diag.panel = panel.hist, lower.panel = panel.cor, upper.panel = panel.contour) } - graphics::pairs(out_all,xlim=c(0,1),ylim=c(0,1),main=title[i],diag.panel=panel.hist,lower.panel=panel.cor,upper.panel=panel.contour) } - + if (exists("g")) invisible(g) } - if(exists('g')) invisible(g) - -} diff --git a/R/posterior_predictive.simmr_output.R b/R/posterior_predictive.simmr_output.R index cdb50ee..88379a1 100644 --- a/R/posterior_predictive.simmr_output.R +++ b/R/posterior_predictive.simmr_output.R @@ -5,107 +5,116 @@ #' @param simmr_out A run of the simmr model from \code{\link{simmr_mcmc}} #' @param group Which group to run it for (currently only numeric rather than group names) #' @param prob The probability interval for the posterior predictives. The default is 0.5 (i.e. 50pc intervals) -#' @param plot_ppc Whether to create a bayesplot of the posterior predictive or not. +#' @param plot_ppc Whether to create a bayesplot of the posterior predictive or not. #' #' @importFrom bayesplot ppc_intervals #' #' @export -#' +#' #' @examples #' \dontrun{ #' data(geese_data_day1) -#' simmr_1 = with(geese_data_day1, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) -#' +#' simmr_1 <- with( +#' geese_data_day1, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) +#' #' # Plot #' plot(simmr_1) -#' +#' #' # Print #' simmr_1 -#' +#' #' # MCMC run -#' simmr_1_out = simmr_mcmc(simmr_1) +#' simmr_1_out <- simmr_mcmc(simmr_1) #' #' # Prior predictive -#' post_pred = posterior_predictive(simmr_1_out) +#' post_pred <- posterior_predictive(simmr_1_out) #' } -posterior_predictive = function(simmr_out, - group = 1, - prob = 0.5, - plot_ppc = TRUE) { - UseMethod('posterior_predictive') -} +posterior_predictive <- function(simmr_out, + group = 1, + prob = 0.5, + plot_ppc = TRUE) { + UseMethod("posterior_predictive") +} #' @export -posterior_predictive.simmr_output = function(simmr_out, - group = 1, - prob = 0.5, - plot_ppc = TRUE) { - -# Can't do more than 1 group for now -if(length(group) > 1) stop("Multiple groups not supported") -# Get the original jags script -model_string_old = simmr_out$output[[group]]$model$model() - -# Plug in y_pred -copy_lines = model_string_old[6] -copy_lines = sub("y\\[i","y_pred\\[i",copy_lines) -model_string_new = c(model_string_old[1:6],copy_lines,model_string_old[7:length(model_string_old)]) - -# Re-Run in JAGS -output = R2jags::jags(data=simmr_out$output[[group]]$model$data(), - parameters.to.save = c("y_pred"), - model.file = textConnection(paste0(model_string_new, collapse = '\n')), - n.chains = simmr_out$output[[group]]$BUGSoutput$n.chains, - n.iter = simmr_out$output[[group]]$BUGSoutput$n.iter, - n.burnin = simmr_out$output[[group]]$BUGSoutput$n.burnin, - n.thin = simmr_out$output[[group]]$BUGSoutput$n.thin) +posterior_predictive.simmr_output <- function(simmr_out, + group = 1, + prob = 0.5, + plot_ppc = TRUE) { -y_post_pred = output$BUGSoutput$sims.list$y_pred + # Can't do more than 1 group for now + if (length(group) > 1) stop("Multiple groups not supported") + # Get the original jags script + model_string_old <- simmr_out$output[[group]]$model$model() -# Make is look nicer -low_prob = 0.5 - prob/2 -high_prob = 0.5 + prob/2 -y_post_pred_ci = apply(y_post_pred, - 2:3, - 'quantile', - prob = c(low_prob, high_prob)) -y_post_pred_out = data.frame( - interval = matrix(y_post_pred_ci, - ncol = simmr_out$input$n_tracers, - byrow = TRUE), - data = as.vector(simmr_out$input$mixtures[simmr_out$input$group_int == group,]) -) + # Plug in y_pred + copy_lines <- model_string_old[6] + copy_lines <- sub("y\\[i", "y_pred\\[i", copy_lines) + model_string_new <- c(model_string_old[1:6], copy_lines, model_string_old[7:length(model_string_old)]) -y_post_pred_out$outside = y_post_pred_out[,3] > y_post_pred_out[,2] | - y_post_pred_out[,3] < y_post_pred_out[,1] -prop_outside = mean(y_post_pred_out$outside) + # Re-Run in JAGS + output <- R2jags::jags( + data = simmr_out$output[[group]]$model$data(), + parameters.to.save = c("y_pred"), + model.file = textConnection(paste0(model_string_new, collapse = "\n")), + n.chains = simmr_out$output[[group]]$BUGSoutput$n.chains, + n.iter = simmr_out$output[[group]]$BUGSoutput$n.iter, + n.burnin = simmr_out$output[[group]]$BUGSoutput$n.burnin, + n.thin = simmr_out$output[[group]]$BUGSoutput$n.thin + ) -if(plot_ppc) { - y_rep = y_post_pred - dim(y_rep) = c(dim(y_post_pred)[1], dim(y_post_pred)[2]*dim(y_post_pred)[3]) - curr_rows = which(simmr_out$input$group_int==group) - curr_mix = simmr_out$input$mixtures[curr_rows,,drop=FALSE] - g = ppc_intervals( - y = as.vector(curr_mix), - yrep = y_rep, - x = rep(1:nrow(curr_mix), simmr_out$input$n_tracers), - prob = prob, - fatten = 1 - ) + ggplot2::ylab("Tracer value") + - ggplot2::xlab('Observation') + - ggplot2::ggtitle(paste0(prob*100, '% posterior predictive')) + - ggplot2::theme_bw() + - ggplot2::scale_x_continuous(breaks = 1:simmr_out$input$n_obs) - print(g) -} -# Return the simulations -invisible(list(table = y_post_pred_out, - prop_outside = prop_outside)) + y_post_pred <- output$BUGSoutput$sims.list$y_pred + + # Make is look nicer + low_prob <- 0.5 - prob / 2 + high_prob <- 0.5 + prob / 2 + y_post_pred_ci <- apply(y_post_pred, + 2:3, + "quantile", + prob = c(low_prob, high_prob) + ) + y_post_pred_out <- data.frame( + interval = matrix(y_post_pred_ci, + ncol = simmr_out$input$n_tracers, + byrow = TRUE + ), + data = as.vector(simmr_out$input$mixtures[simmr_out$input$group_int == group, ]) + ) -} \ No newline at end of file + y_post_pred_out$outside <- y_post_pred_out[, 3] > y_post_pred_out[, 2] | + y_post_pred_out[, 3] < y_post_pred_out[, 1] + prop_outside <- mean(y_post_pred_out$outside) + + if (plot_ppc) { + y_rep <- y_post_pred + dim(y_rep) <- c(dim(y_post_pred)[1], dim(y_post_pred)[2] * dim(y_post_pred)[3]) + curr_rows <- which(simmr_out$input$group_int == group) + curr_mix <- simmr_out$input$mixtures[curr_rows, , drop = FALSE] + g <- ppc_intervals( + y = as.vector(curr_mix), + yrep = y_rep, + x = rep(1:nrow(curr_mix), simmr_out$input$n_tracers), + prob = prob, + fatten = 1 + ) + ggplot2::ylab("Tracer value") + + ggplot2::xlab("Observation") + + ggplot2::ggtitle(paste0(prob * 100, "% posterior predictive")) + + ggplot2::theme_bw() + + ggplot2::scale_x_continuous(breaks = 1:simmr_out$input$n_obs) + print(g) + } + # Return the simulations + invisible(list( + table = y_post_pred_out, + prop_outside = prop_outside + )) +} diff --git a/R/print.simmr_input.R b/R/print.simmr_input.R index 4a7b78f..dd35ca5 100644 --- a/R/print.simmr_input.R +++ b/R/print.simmr_input.R @@ -3,19 +3,19 @@ #' @param x An object of class \code{simmr_input} #' @param ... Other arguments (not supported) #' -#' @return A neat presentation of your simmr object. +#' @return A neat presentation of your simmr object. #' @export print.simmr_input <- -function(x,...) { - cat('This is a valid simmr input object with ') - cat(paste(x$n_obs,'observations, ')) - cat(paste(x$n_tracers,'tracers, and ')) - cat(paste(x$n_sources,'sources.\n')) - if(x$n_groups>1) cat(paste('There are',x$n_groups,'groups.\n')) - cat('The source names are: ') - cat(x$source_names,sep=', ') - cat('.\n') - cat('The tracer names are: ') - cat(colnames(x$mixtures), sep = ', ') - cat('.\n\n') -} + function(x, ...) { + cat("This is a valid simmr input object with ") + cat(paste(x$n_obs, "observations, ")) + cat(paste(x$n_tracers, "tracers, and ")) + cat(paste(x$n_sources, "sources.\n")) + if (x$n_groups > 1) cat(paste("There are", x$n_groups, "groups.\n")) + cat("The source names are: ") + cat(x$source_names, sep = ", ") + cat(".\n") + cat("The tracer names are: ") + cat(colnames(x$mixtures), sep = ", ") + cat(".\n\n") + } diff --git a/R/print.simmr_output.R b/R/print.simmr_output.R index 5764c0a..9e30d64 100644 --- a/R/print.simmr_output.R +++ b/R/print.simmr_output.R @@ -4,13 +4,13 @@ #' @param ... Other arguments (not supported) #' #' @return Returns a neat summary of the object -#' +#' #' @seealso \code{\link{simmr_mcmc}} for creating \code{simmr_output} objects #' @export print.simmr_output <- -function(x,...) { - print(x$input) - cat('The input data has been run via simmr_mcmc and has produced ') - cat(nrow(x$output[[1]]$BUGSoutput$sims.matrix),'iterations over',x$output[[1]]$BUGSoutput$n.chains,'MCMC chains.') - cat('\n\n') -} + function(x, ...) { + print(x$input) + cat("The input data has been run via simmr_mcmc and has produced ") + cat(nrow(x$output[[1]]$BUGSoutput$sims.matrix), "iterations over", x$output[[1]]$BUGSoutput$n.chains, "MCMC chains.") + cat("\n\n") + } diff --git a/R/prior_viz.simmr_output.R b/R/prior_viz.simmr_output.R index 3c50c6d..b894749 100644 --- a/R/prior_viz.simmr_output.R +++ b/R/prior_viz.simmr_output.R @@ -1,7 +1,7 @@ #' Plot the prior distribution for a simmr run #' #' This function takes the output from \code{\link{simmr_mcmc}} and plots the prior distribution to enable visual inspection. This can be used by itself or as part of \code{\link{posterior_predictive}} to visually evaluate the influence of the prior on the posterior distribution. -#' +#' #' @param simmr_out A run of the simmr model from \code{\link{simmr_mcmc}} #' @param group Which group to run it for (currently only numeric rather than group names) #' @param plot Whether to create a density plot of the prior or not. The simulated prior values are returned as part of the object @@ -11,112 +11,122 @@ #' #' @export -#' -#' @examples +#' +#' @examples #' \dontrun{ #' data(geese_data_day1) -#' simmr_1 = with(geese_data_day1, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) -#' +#' simmr_1 <- with( +#' geese_data_day1, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) +#' #' # Plot #' plot(simmr_1) -#' +#' #' # Print #' simmr_1 -#' +#' #' # MCMC run -#' simmr_1_out = simmr_mcmc(simmr_1) +#' simmr_1_out <- simmr_mcmc(simmr_1) #' #' # Prior predictive -#' prior = prior_viz(simmr_1_out) +#' prior <- prior_viz(simmr_1_out) #' head(prior) #' summary(prior) #' } -prior_viz = function(simmr_out, - group = 1, - plot = TRUE, - include_posterior = TRUE, - n_sims = 10000, - ggargs = NULL) { - UseMethod('prior_viz') -} +prior_viz <- function(simmr_out, + group = 1, + plot = TRUE, + include_posterior = TRUE, + n_sims = 10000, + ggargs = NULL) { + UseMethod("prior_viz") +} #' @export -prior_viz.simmr_output = function(simmr_out, - group = 1, - plot = TRUE, - include_posterior = TRUE, - n_sims = 10000, - ggargs = NULL) { +prior_viz.simmr_output <- function(simmr_out, + group = 1, + plot = TRUE, + include_posterior = TRUE, + n_sims = 10000, + ggargs = NULL) { # Can't do more than 1 group for now - if(length(group) > 1) stop("Multiple groups not supported") - + if (length(group) > 1) stop("Multiple groups not supported") + # Get group_name - n_groups = simmr_out$input$n_groups - plot_title_1 = "Prior distributions" - plot_title_2 = "Prior and posterior distributions" - if(n_groups > 1) { - group_name = levels(simmr_out$input$group)[group] - plot_title_1 = paste0(plot_title_1, ': ', group_name) - plot_title_2 = paste0(plot_title_2, ': ', group_name) + n_groups <- simmr_out$input$n_groups + plot_title_1 <- "Prior distributions" + plot_title_2 <- "Prior and posterior distributions" + if (n_groups > 1) { + group_name <- levels(simmr_out$input$group)[group] + plot_title_1 <- paste0(plot_title_1, ": ", group_name) + plot_title_2 <- paste0(plot_title_2, ": ", group_name) } - - # Plot and/or output the prior - mu_f_mean = simmr_out$output[[1]]$model$data()$mu_f_mean - sigma_f_sd = simmr_out$output[[1]]$model$data()$sigma_f_sd - n_sources = simmr_out$input$n_sources - + + # Plot and/or output the prior + mu_f_mean <- simmr_out$output[[1]]$model$data()$mu_f_mean + sigma_f_sd <- simmr_out$output[[1]]$model$data()$sigma_f_sd + n_sources <- simmr_out$input$n_sources + # Now simulate some ps - p_prior_sim = matrix(NA, ncol = n_sources, nrow = n_sims) - for(i in 1:n_sims) { - f = rnorm(n_sources,mean = mu_f_mean, sd = sigma_f_sd) - p_prior_sim[i,] = exp(f)/sum(exp(f)) + p_prior_sim <- matrix(NA, ncol = n_sources, nrow = n_sims) + for (i in 1:n_sims) { + f <- rnorm(n_sources, mean = mu_f_mean, sd = sigma_f_sd) + p_prior_sim[i, ] <- exp(f) / sum(exp(f)) } - colnames(p_prior_sim) = simmr_out$input$source_names - if(plot) { - df = reshape2::melt(p_prior_sim) - colnames(df) = c('Num','Source','Proportion') - df$Type = "Prior" - out_all = simmr_out$output[[group]]$BUGSoutput$sims.list$p - df2 = reshape2::melt(out_all) - colnames(df2) = c('Num','Source','Proportion') - df2$Type = "Posterior" - df_all = rbind(df2, df) - - if(include_posterior) { - g=ggplot(df_all, - aes_string(x="Proportion", - y="..density..", - fill="Source", - linetype="Type")) + - scale_fill_viridis(discrete=TRUE) + + colnames(p_prior_sim) <- simmr_out$input$source_names + if (plot) { + df <- reshape2::melt(p_prior_sim) + colnames(df) <- c("Num", "Source", "Proportion") + df$Type <- "Prior" + out_all <- simmr_out$output[[group]]$BUGSoutput$sims.list$p + df2 <- reshape2::melt(out_all) + colnames(df2) <- c("Num", "Source", "Proportion") + df2$Type <- "Posterior" + df_all <- rbind(df2, df) + if (include_posterior) { + g <- ggplot( + df_all, + aes_string( + x = "Proportion", + y = "..density..", + fill = "Source", + linetype = "Type" + ) + ) + + scale_fill_viridis(discrete = TRUE) + geom_density(alpha = 0.5) + theme_bw() + - ggtitle(plot_title_2) + + ggtitle(plot_title_2) + ylab("Density") + facet_wrap("~ Source") } else { - g=ggplot(df, - aes_string(x="Proportion", - y="..density..", - fill="Source")) + - scale_fill_viridis(discrete=TRUE) + + g <- ggplot( + df, + aes_string( + x = "Proportion", + y = "..density..", + fill = "Source" + ) + ) + + scale_fill_viridis(discrete = TRUE) + geom_density(alpha = 0.5, linetype = 0) + theme_bw() + - ggtitle(plot_title_1) + + ggtitle(plot_title_1) + ylab("Density") + facet_wrap("~ Source") } print(g) + ggargs } - + # Return the simulations invisible(p_prior_sim) - -} \ No newline at end of file +} diff --git a/R/simmr.R b/R/simmr.R index 3fe412e..c47b64d 100644 --- a/R/simmr.R +++ b/R/simmr.R @@ -1,5 +1,5 @@ #' simmr: A package for fitting stable isotope mixing models via JAGS in R -#' +#' #' This package runs a #' simple Stable Isotope Mixing Model (SIMM) and is meant as a longer term #' replacement to the previous function SIAR.. These are used to infer dietary @@ -9,56 +9,57 @@ #' composition of fatty acids. The main functions are \code{\link{simmr_load}} #' and \code{\link{simmr_mcmc}}. The help files contain examples of the use of #' this package. See also the vignette for a longer walkthrough. -#' +#' #' An even longer term replacement for properly running SIMMs is MixSIAR, which #' allows for more detailed random effects and the inclusion of covariates. -#' +#' #' @name simmr #' @docType package #' @author Andrew Parnell -#' +#' #' @references Andrew C. Parnell, Donald L. Phillips, Stuart Bearhop, Brice X. #' Semmens, Eric J. Ward, Jonathan W. Moore, Andrew L. Jackson, Jonathan Grey, #' David J. Kelly, and Richard Inger. Bayesian stable isotope mixing models. #' Environmetrics, 24(6):387–399, 2013. -#' +#' #' Andrew C Parnell, Richard Inger, Stuart Bearhop, and Andrew L Jackson. #' Source partitioning using stable isotopes: coping with too much variation. #' PLoS ONE, 5(3):5, 2010. #' @keywords multivariate #' @examples -#' +#' #' \dontrun{ #' # A first example with 2 tracers (isotopes), 10 observations, and 4 food sources #' data(geese_data_day1) -#' simmr_in = with(geese_data_day1, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) -#' +#' simmr_in <- with( +#' geese_data_day1, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) +#' #' # Plot #' plot(simmr_in) -#' +#' #' # MCMC run -#' simmr_out = simmr_mcmc(simmr_in) -#' +#' simmr_out <- simmr_mcmc(simmr_in) +#' #' # Check convergence - values should all be close to 1 -#' summary(simmr_out, type = 'diagnostics') -#' +#' summary(simmr_out, type = "diagnostics") +#' #' # Look at output -#' summary(simmr_out, type = 'statistics') -#' +#' summary(simmr_out, type = "statistics") +#' #' # Look at influence of priors #' prior_viz(simmr_out) -#' +#' #' # Plot output -#' plot(simmr_out, type = 'histogram') +#' plot(simmr_out, type = "histogram") #' } NULL - - - diff --git a/R/simmr_elicit.R b/R/simmr_elicit.R index b52cb36..d7d9892 100644 --- a/R/simmr_elicit.R +++ b/R/simmr_elicit.R @@ -37,7 +37,7 @@ #' \code{control.prior} in \code{\link{simmr_mcmc}}} #' #' @author Andrew Parnell -#' +#' #' @importFrom stats rnorm optim #' @importFrom compositions clrInv #' @importFrom boot logit @@ -47,74 +47,81 @@ #' \dontrun{ #' # Data set: 10 observations, 2 tracers, 4 sources #' data(geese_data_day1) -#' simmr_1 = with(geese_data_day1, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) -#' +#' simmr_1 <- with( +#' geese_data_day1, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) +#' #' # MCMC run -#' simmr_1_out = simmr_mcmc(simmr_1) -#' +#' simmr_1_out <- simmr_mcmc(simmr_1) +#' #' # Look at the prior influence #' prior_viz(simmr_1_out) -#' +#' #' # Summary -#' summary(simmr_1_out,'quantiles') +#' summary(simmr_1_out, "quantiles") #' # A bit vague: #' # 2.5% 25% 50% 75% 97.5% #' # Source A 0.029 0.115 0.203 0.312 0.498 #' # Source B 0.146 0.232 0.284 0.338 0.453 #' # Source C 0.216 0.255 0.275 0.296 0.342 #' # Source D 0.032 0.123 0.205 0.299 0.465 -#' -#' # Now suppose I had prior information that: -#' # proportion means = 0.5,0.2,0.2,0.1 +#' +#' # Now suppose I had prior information that: +#' # proportion means = 0.5,0.2,0.2,0.1 #' # proportion sds = 0.08,0.02,0.01,0.02 -#' prior = simmr_elicit(4, c(0.5,0.2,0.2,0.1), c(0.08,0.02,0.01,0.02)) -#' -#' simmr_1a_out = simmr_mcmc(simmr_1,prior_control=list(means=prior$mean,sd=prior$sd)) -#' +#' prior <- simmr_elicit(4, c(0.5, 0.2, 0.2, 0.1), c(0.08, 0.02, 0.01, 0.02)) +#' +#' simmr_1a_out <- simmr_mcmc(simmr_1, prior_control = list(means = prior$mean, sd = prior$sd)) +#' #' #' # Look at the prior influence now #' prior_viz(simmr_1a_out) -#' -#' summary(simmr_1a_out,'quantiles') +#' +#' summary(simmr_1a_out, "quantiles") #' # Much more precise: #' # 2.5% 25% 50% 75% 97.5% #' # Source A 0.441 0.494 0.523 0.553 0.610 #' # Source B 0.144 0.173 0.188 0.204 0.236 #' # Source C 0.160 0.183 0.196 0.207 0.228 -#' # Source D 0.060 0.079 0.091 0.105 0.135 +#' # Source D 0.060 0.079 0.091 0.105 0.135 #' } #' @export simmr_elicit simmr_elicit <- function(n_sources, - proportion_means = rep(1/n_sources, n_sources), + proportion_means = rep(1 / n_sources, n_sources), proportion_sds = rep(0.1, n_sources), n_sims = 1000) { # proportion_means must be a vector of length n_sources - if (length(proportion_means) != n_sources) - stop('proportion_means must be of same length as the number of sources') + if (length(proportion_means) != n_sources) { + stop("proportion_means must be of same length as the number of sources") + } # proportion_sds must be a vector of length n_sources - if (length(proportion_sds) != n_sources) - stop('proportion_sds must be of same length as the number of sources') - if(any(proportion_sds==0)) + if (length(proportion_sds) != n_sources) { + stop("proportion_sds must be of same length as the number of sources") + } + if (any(proportion_sds == 0)) { stop("No proportion_sds should be 0 as this will mean that food source is not being consumed.") - - low_cis = proportion_means - 2*proportion_sds - high_cis = proportion_means + 2*proportion_sds - if(any(low_cis < 0) | any(high_cis > 1) ) warning("Some proportion sds are large and lie at the edge of the simplex. Check results are reasonable by running prior_viz on the object created by simmr_mcmc.") - - cat('Running elicitation optimisation routine...\n') - + } + + low_cis <- proportion_means - 2 * proportion_sds + high_cis <- proportion_means + 2 * proportion_sds + if (any(low_cis < 0) | any(high_cis > 1)) warning("Some proportion sds are large and lie at the edge of the simplex. Check results are reasonable by running prior_viz on the object created by simmr_mcmc.") + + cat("Running elicitation optimisation routine...\n") + # Perform an optimisation to match the standard deviations to a normal distribution - clr_opt_mean = function(pars) { - means = c(pars[(1:(n_sources - 1))], -sum(pars[(1:(n_sources - 1))])) + clr_opt_mean <- function(pars) { + means <- c(pars[(1:(n_sources - 1))], -sum(pars[(1:(n_sources - 1))])) # Generate n_sims observations from this normal distribution - f = matrix( + f <- matrix( stats::rnorm( n_sims * n_sources, mean = means, @@ -124,14 +131,16 @@ simmr_elicit <- nrow = n_sims, byrow = TRUE ) - p = as.matrix(compositions::clrInv(f)) + p <- as.matrix(compositions::clrInv(f)) return(sum(( - boot::logit(apply(p, 2, 'mean')) - boot::logit(proportion_means) - ) ^ 2)) + boot::logit(apply(p, 2, "mean")) - boot::logit(proportion_means) + )^2)) } - - opt_mean = stats::optim(c(rep(0, (n_sources - 1))), clr_opt_mean, method = - "SANN") + + opt_mean <- stats::optim(c(rep(0, (n_sources - 1))), clr_opt_mean, + method = + "SANN" + ) if (opt_mean$convergence != 0) { warning( "Optimisation for means did not converge properly. Please either increase n_sims or adjust proportion_means and proportion_sds" @@ -139,25 +148,25 @@ simmr_elicit <- } else { cat("Mean optimisation successful.\n") } - best_mean = c(opt_mean$par, -sum(opt_mean$par)) - + best_mean <- c(opt_mean$par, -sum(opt_mean$par)) + options(warn = -1) - clr_opt_sd = function(pars) { - sd = pars[1:(n_sources)] + clr_opt_sd <- function(pars) { + sd <- pars[1:(n_sources)] # Generate n_sims observations from this normal distribution - f = matrix( + f <- matrix( stats::rnorm(n_sims * n_sources, mean = best_mean, sd = sd), ncol = n_sources, nrow = n_sims, byrow = TRUE ) - p = as.matrix(compositions::clrInv(f)) + p <- as.matrix(compositions::clrInv(f)) return(sum(( - boot::logit(apply(p, 2, 'sd')) - boot::logit(proportion_sds) - ) ^ 2)) + boot::logit(apply(p, 2, "sd")) - boot::logit(proportion_sds) + )^2)) } - - opt_sd = stats::optim(rep(1, (n_sources)), clr_opt_sd, method = "SANN") + + opt_sd <- stats::optim(rep(1, (n_sources)), clr_opt_sd, method = "SANN") if (opt_sd$convergence != 0) { warning( "Optimisation for stand deviations did not converge properly. Please either increase n_sims or adjust proportion_means and proportion_sds" @@ -165,25 +174,24 @@ simmr_elicit <- } else { cat("Standard deviation optimisation successful.\n") } - best_sd = opt_sd$par + best_sd <- opt_sd$par options(warn = 0) - - best_f = matrix( + + best_f <- matrix( stats::rnorm(n_sims * n_sources, mean = best_mean, sd = best_sd), ncol = n_sources, nrow = n_sims, byrow = TRUE ) - best_p = as.matrix(compositions::clrInv(best_f)) - - cat('Best fit estimates provide proportion means of:\n') - cat(round(apply(best_p, 2, 'mean'), 3)) - cat('\n') - cat('... and best fit standard deviations of:\n') - cat(round(apply(best_p, 2, 'sd'), 3)) - cat('\n') - cat('Check these match the input values before proceeding with a model run.\n') - + best_p <- as.matrix(compositions::clrInv(best_f)) + + cat("Best fit estimates provide proportion means of:\n") + cat(round(apply(best_p, 2, "mean"), 3)) + cat("\n") + cat("... and best fit standard deviations of:\n") + cat(round(apply(best_p, 2, "sd"), 3)) + cat("\n") + cat("Check these match the input values before proceeding with a model run.\n") + return(list(mean = best_mean, sd = best_sd)) - -} + } diff --git a/R/simmr_load.R b/R/simmr_load.R index ba00b88..d42d59a 100644 --- a/R/simmr_load.R +++ b/R/simmr_load.R @@ -39,10 +39,10 @@ #' columns is the number of tracers. These should be between 0 and 1. If not #' provided these are all set to 1. #' @param group A grouping variable. These can be a character or factor variable -#' +#' #' @import checkmate #' -#' +#' #' @return An object of class \code{simmr_input} with the following elements: #' \item{mixtures }{The mixture data} \item{source_neams }{Source means} #' \item{sources_sds }{Source standard deviations} \item{correction_means @@ -57,84 +57,96 @@ #' #' # A simple example with 10 observations, 2 tracers and 4 sources #' #' data(geese_data_day1) -#' simmr_1 = with(geese_data_day1, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) +#' simmr_1 <- with( +#' geese_data_day1, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) #' #' print(simmr_1) -#' #' @export simmr_load -simmr_load = function(mixtures, - source_names, - source_means, - source_sds, - correction_means = NULL, - correction_sds = NULL, - concentration_means = NULL, - group = NULL) { +simmr_load <- function(mixtures, + source_names, + source_means, + source_sds, + correction_means = NULL, + correction_sds = NULL, + concentration_means = NULL, + group = NULL) { # Function to load in data for simmr and check whether it's appropriate for running through simmr_mcmc - + # Go through each object and check that it matches the requirements - + # Mixtures must be a matrix - the number of rows is the number of observations and the number of columns is the number of tracers assert_matrix(mixtures) - n_obs = nrow(mixtures) - n_tracers = ncol(mixtures) - + n_obs <- nrow(mixtures) + n_tracers <- ncol(mixtures) + # Add column names if they're not there - if (is.null(colnames(mixtures))) - colnames(mixtures) = paste0('tracer', 1:n_tracers) - + if (is.null(colnames(mixtures))) { + colnames(mixtures) <- paste0("tracer", 1:n_tracers) + } + # source_names must be a character vector - the length of it is the number of sources assert_character(source_names) - n_sources = length(source_names) - + n_sources <- length(source_names) + # source_means and source_sds must both be matrices where the number of rows is n_sources (in the same order as source_names) and the number of columns is n_tracers assert_matrix(source_means, nrows = n_sources, ncols = n_tracers) assert_matrix(source_sds, nrows = n_sources, ncols = n_tracers) - assert_matrix(correction_means, nrows = n_sources, - ncols = n_tracers, - null.ok = ifelse(is.null(correction_sds), - TRUE, FALSE)) - assert_matrix(correction_sds, nrows = n_sources, - ncols = n_tracers, - null.ok = ifelse(is.null(correction_means), - TRUE, FALSE)) - assert_matrix(concentration_means, nrows = n_sources, - ncols = n_tracers, null.ok = TRUE) + assert_matrix(correction_means, + nrows = n_sources, + ncols = n_tracers, + null.ok = ifelse(is.null(correction_sds), + TRUE, FALSE + ) + ) + assert_matrix(correction_sds, + nrows = n_sources, + ncols = n_tracers, + null.ok = ifelse(is.null(correction_means), + TRUE, FALSE + ) + ) + assert_matrix(concentration_means, + nrows = n_sources, + ncols = n_tracers, null.ok = TRUE + ) # Fill in correction means if (is.null(correction_means)) { - correction_means = matrix(0, ncol = n_tracers, nrow = n_sources) - correction_sds = matrix(0, ncol = n_tracers, nrow = n_sources) + correction_means <- matrix(0, ncol = n_tracers, nrow = n_sources) + correction_sds <- matrix(0, ncol = n_tracers, nrow = n_sources) } - + # concentration_means must be a matrix where all elements are less than 1 - if(is.null(concentration_means)) { - concentration_means = matrix(1, ncol = n_tracers, nrow = n_sources) + if (is.null(concentration_means)) { + concentration_means <- matrix(1, ncol = n_tracers, nrow = n_sources) } else { assert_true(all(concentration_means < 1) & all(concentration_means > 0)) } - + # Check the groups are the right length and structure if given if (!is.null(group)) { # Allow for the group to be a factor variable - group = as.factor(group) - group_int = as.integer(group) + group <- as.factor(group) + group_int <- as.integer(group) assert_vector(group, len = n_obs) } else { - group = factor(rep(1, n_obs)) - group_int = as.integer(group) + group <- factor(rep(1, n_obs)) + group_int <- as.integer(group) } - n_groups = length(unique(group)) - + n_groups <- length(unique(group)) + # Prepare output and give class - out = list( + out <- list( mixtures = mixtures, source_names = source_names, source_means = source_means, @@ -149,14 +161,13 @@ simmr_load = function(mixtures, n_sources = n_sources, n_groups = n_groups ) - + # Look through to see whether there are any missing values in anything - if (any(unlist(lapply(out, 'is.na')))) { + if (any(unlist(lapply(out, "is.na")))) { warning("Missing values provided for some values. Check your inputs") } - - class(out) = 'simmr_input' - + + class(out) <- "simmr_input" + return(out) - } diff --git a/R/simmr_mcmc.R b/R/simmr_mcmc.R index bd2a632..a84f3a5 100644 --- a/R/simmr_mcmc.R +++ b/R/simmr_mcmc.R @@ -1,17 +1,17 @@ #' Run a \code{simmr_input} object through the main simmr Markov chain Monte #' Carlo (MCMC) function -#' +#' #' This is the main function of simmr. It takes a \code{simmr_input} object #' created via \code{\link{simmr_load}}, runs an MCMC to determine the dietary #' proportions, and then outputs a \code{simmr_output} object for further #' analysis and plotting via \code{\link{summary.simmr_output}} and #' \code{\link{plot.simmr_output}}. -#' +#' #' If, after running \code{\link{simmr_mcmc}} the convergence diagnostics in #' \code{\link{summary.simmr_output}} are not satisfactory, the values of #' \code{iter}, \code{burn} and \code{thin} in \code{mcmc_control} should be #' increased by a factor of 10. -#' +#' #' @param simmr_in An object created via the function \code{\link{simmr_load}} #' @param prior_control A list of values including arguments named \code{means} #' and \code{sd} which represent the prior means and standard deviations of the @@ -28,260 +28,291 @@ #' \code{mcmc.list} from the coda package. These can be analysed using the #' \code{\link{summary.simmr_output}} and \code{\link{plot.simmr_output}} #' functions.} -#' +#' #' @author Andrew Parnell -#' +#' #' @seealso \code{\link{simmr_load}} for creating objects suitable for this #' function, \code{\link{plot.simmr_input}} for creating isospace plots, #' \code{\link{summary.simmr_output}} for summarising output, and #' \code{\link{plot.simmr_output}} for plotting output. -#' +#' #' @references Andrew C. Parnell, Donald L. Phillips, Stuart Bearhop, Brice X. #' Semmens, Eric J. Ward, Jonathan W. Moore, Andrew L. Jackson, Jonathan Grey, #' David J. Kelly, and Richard Inger. Bayesian stable isotope mixing models. #' Environmetrics, 24(6):387–399, 2013. -#' +#' #' Andrew C Parnell, Richard Inger, Stuart Bearhop, and Andrew L Jackson. #' Source partitioning using stable isotopes: coping with too much variation. #' PLoS ONE, 5(3):5, 2010. -#' +#' #' @importFrom R2jags jags -#' +#' #' @examples #' \dontrun{ #' ## See the package vignette for a detailed run through of these 4 examples -#' +#' #' # Data set 1: 10 obs on 2 isos, 4 sources, with tefs and concdep #' data(geese_data_day1) -#' simmr_1 = with(geese_data_day1, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) -#' +#' simmr_1 <- with( +#' geese_data_day1, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) +#' #' # Plot #' plot(simmr_1) -#' +#' #' # Print #' simmr_1 -#' +#' #' # MCMC run -#' simmr_1_out = simmr_mcmc(simmr_1) -#' +#' simmr_1_out <- simmr_mcmc(simmr_1) +#' #' # Print it #' print(simmr_1_out) -#' +#' #' # Summary -#' summary(simmr_1_out,type='diagnostics') -#' summary(simmr_1_out,type='correlations') -#' summary(simmr_1_out,type='statistics') -#' ans = summary(simmr_1_out,type=c('quantiles','statistics')) -#' +#' summary(simmr_1_out, type = "diagnostics") +#' summary(simmr_1_out, type = "correlations") +#' summary(simmr_1_out, type = "statistics") +#' ans <- summary(simmr_1_out, type = c("quantiles", "statistics")) +#' #' # Plot -#' plot(simmr_1_out,type='boxplot') -#' plot(simmr_1_out,type='histogram') -#' plot(simmr_1_out,type='density') -#' plot(simmr_1_out,type='matrix') -#' +#' plot(simmr_1_out, type = "boxplot") +#' plot(simmr_1_out, type = "histogram") +#' plot(simmr_1_out, type = "density") +#' plot(simmr_1_out, type = "matrix") +#' #' # Compare two sources -#' compare_sources(simmr_1_out,source_names=c('Source A','Source D')) -#' +#' compare_sources(simmr_1_out, source_names = c("Source A", "Source D")) +#' #' # Compare multiple sources #' compare_sources(simmr_1_out) -#' +#' #' ##################################################################################### -#' +#' #' # A version with just one observation #' data(geese_data_day1) -#' simmr_2 = with(geese_data_day1, -#' simmr_load(mixtures=mixtures[1,,drop = FALSE], -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) -#' +#' simmr_2 <- with( +#' geese_data_day1, +#' simmr_load( +#' mixtures = mixtures[1, , drop = FALSE], +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) +#' #' # Plot #' plot(simmr_2) -#' +#' #' # MCMC run - automatically detects the single observation -#' simmr_2_out = simmr_mcmc(simmr_2) -#' +#' simmr_2_out <- simmr_mcmc(simmr_2) +#' #' # Print it #' print(simmr_2_out) -#' +#' #' # Summary #' summary(simmr_2_out) -#' summary(simmr_2_out,type='diagnostics') -#' ans = summary(simmr_2_out,type=c('quantiles')) -#' +#' summary(simmr_2_out, type = "diagnostics") +#' ans <- summary(simmr_2_out, type = c("quantiles")) +#' #' # Plot #' plot(simmr_2_out) -#' plot(simmr_2_out,type='boxplot') -#' plot(simmr_2_out,type='histogram') -#' plot(simmr_2_out,type='density') -#' plot(simmr_2_out,type='matrix') -#' +#' plot(simmr_2_out, type = "boxplot") +#' plot(simmr_2_out, type = "histogram") +#' plot(simmr_2_out, type = "density") +#' plot(simmr_2_out, type = "matrix") +#' #' ##################################################################################### -#' +#' #' # Data set 2: 3 isotopes (d13C, d15N and d34S), 30 observations, 4 sources #' data(simmr_data_2) -#' simmr_3 = with(simmr_data_2, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) -#' +#' simmr_3 <- with( +#' simmr_data_2, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) +#' #' # Get summary #' print(simmr_3) -#' +#' #' # Plot 3 times #' plot(simmr_3) -#' plot(simmr_3,tracers=c(2,3)) -#' plot(simmr_3,tracers=c(1,3)) +#' plot(simmr_3, tracers = c(2, 3)) +#' plot(simmr_3, tracers = c(1, 3)) #' # See vignette('simmr') for fancier axis labels -#' +#' #' # MCMC run -#' simmr_3_out = simmr_mcmc(simmr_3) -#' +#' simmr_3_out <- simmr_mcmc(simmr_3) +#' #' # Print it #' print(simmr_3_out) -#' +#' #' # Summary #' summary(simmr_3_out) -#' summary(simmr_3_out,type='diagnostics') -#' summary(simmr_3_out,type='quantiles') -#' summary(simmr_3_out,type='correlations') -#' +#' summary(simmr_3_out, type = "diagnostics") +#' summary(simmr_3_out, type = "quantiles") +#' summary(simmr_3_out, type = "correlations") +#' #' # Plot #' plot(simmr_3_out) -#' plot(simmr_3_out,type='boxplot') -#' plot(simmr_3_out,type='histogram') -#' plot(simmr_3_out,type='density') -#' plot(simmr_3_out,type='matrix') -#' +#' plot(simmr_3_out, type = "boxplot") +#' plot(simmr_3_out, type = "histogram") +#' plot(simmr_3_out, type = "density") +#' plot(simmr_3_out, type = "matrix") +#' #' ##################################################################################### -#' +#' #' # Data set 4 - identified by Fry (2014) as a failing of SIMMs #' # See the vignette for more interpreation of these data and the output -#' +#' #' # The data #' data(square_data) -#' simmr_4 = with(square_data, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds)) -#' +#' simmr_4 <- with( +#' square_data, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds +#' ) +#' ) +#' #' # Get summary #' print(simmr_4) -#' -#' # Plot +#' +#' # Plot #' plot(simmr_4) -#' +#' #' # MCMC run - needs slightly longer -#' simmr_4_out = simmr_mcmc(simmr_4) -#' +#' simmr_4_out <- simmr_mcmc(simmr_4) +#' #' # Print it #' print(simmr_4_out) -#' +#' #' # Summary #' summary(simmr_4_out) -#' summary(simmr_4_out,type='diagnostics') -#' ans = summary(simmr_4_out,type=c('quantiles','statistics')) -#' +#' summary(simmr_4_out, type = "diagnostics") +#' ans <- summary(simmr_4_out, type = c("quantiles", "statistics")) +#' #' # Plot #' plot(simmr_4_out) -#' plot(simmr_4_out,type='boxplot') -#' plot(simmr_4_out,type='histogram') -#' plot(simmr_4_out,type='density') -#' plot(simmr_4_out,type='matrix') # Look at the massive correlations here -#' +#' plot(simmr_4_out, type = "boxplot") +#' plot(simmr_4_out, type = "histogram") +#' plot(simmr_4_out, type = "density") +#' plot(simmr_4_out, type = "matrix") # Look at the massive correlations here +#' #' ##################################################################################### -#' +#' #' # Data set 5 - Multiple groups Geese data from Inger et al 2006 -#' +#' #' # Do this in raw data format - Note that there's quite a few mixtures! #' data(geese_data) -#' simmr_5 = with(geese_data, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means, -#' group = groups)) -#' +#' simmr_5 <- with( +#' geese_data, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means, +#' group = groups +#' ) +#' ) +#' #' # Plot #' plot(simmr_5, -#' xlab=expression(paste(delta^13, "C (\u2030)",sep="")), -#' ylab=expression(paste(delta^15, "N (\u2030)",sep="")), -#' title='Isospace plot of Inger et al Geese data') -#' +#' xlab = expression(paste(delta^13, "C (\\u2030)", sep = "")), +#' ylab = expression(paste(delta^15, "N (\\u2030)", sep = "")), +#' title = "Isospace plot of Inger et al Geese data" +#' ) +#' #' # Run MCMC for each group -#' simmr_5_out = simmr_mcmc(simmr_5) -#' +#' simmr_5_out <- simmr_mcmc(simmr_5) +#' #' # Summarise output -#' summary(simmr_5_out,type='quantiles',group=1) -#' summary(simmr_5_out,type='quantiles',group=c(1,3)) -#' summary(simmr_5_out,type=c('quantiles','statistics'),group=c(1,3)) -#' +#' summary(simmr_5_out, type = "quantiles", group = 1) +#' summary(simmr_5_out, type = "quantiles", group = c(1, 3)) +#' summary(simmr_5_out, type = c("quantiles", "statistics"), group = c(1, 3)) +#' #' # Plot - only a single group allowed -#' plot(simmr_5_out,type='boxplot',group=2,title='simmr output group 2') -#' plot(simmr_5_out,type=c('density','matrix'),grp=6,title='simmr output group 6') -#' +#' plot(simmr_5_out, type = "boxplot", group = 2, title = "simmr output group 2") +#' plot(simmr_5_out, type = c("density", "matrix"), grp = 6, title = "simmr output group 6") +#' #' # Compare sources within a group -#' compare_sources(simmr_5_out,source_names=c('Zostera','U.lactuca'),group=2) -#' compare_sources(simmr_5_out,group=2) -#' +#' compare_sources(simmr_5_out, source_names = c("Zostera", "U.lactuca"), group = 2) +#' compare_sources(simmr_5_out, group = 2) +#' #' # Compare between groups -#' compare_groups(simmr_5_out,source='Zostera',groups=1:2) -#' compare_groups(simmr_5_out,source='Zostera',groups=1:3) -#' compare_groups(simmr_5_out,source='U.lactuca',groups=c(4:5,7,2)) -#' -#' +#' compare_groups(simmr_5_out, source = "Zostera", groups = 1:2) +#' compare_groups(simmr_5_out, source = "Zostera", groups = 1:3) +#' compare_groups(simmr_5_out, source = "U.lactuca", groups = c(4:5, 7, 2)) #' } -#' +#' #' @export -simmr_mcmc = function(simmr_in, - prior_control=list(means=rep(0, - simmr_in$n_sources), - sd=rep(1, - simmr_in$n_sources)), - mcmc_control=list(iter=10000, - burn=1000, - thin=10, - n.chain=4)) { - UseMethod('simmr_mcmc') -} +simmr_mcmc <- function(simmr_in, + prior_control = list( + means = rep( + 0, + simmr_in$n_sources + ), + sd = rep( + 1, + simmr_in$n_sources + ) + ), + mcmc_control = list( + iter = 10000, + burn = 1000, + thin = 10, + n.chain = 4 + )) { + UseMethod("simmr_mcmc") +} #' @export -simmr_mcmc.simmr_input = function(simmr_in, - prior_control=list(means=rep(0,simmr_in$n_sources), - sd=rep(1,simmr_in$n_sources)), - mcmc_control=list(iter=10000, - burn=1000, - thin=10, - n.chain=4)) { +simmr_mcmc.simmr_input <- function(simmr_in, + prior_control = list( + means = rep(0, simmr_in$n_sources), + sd = rep(1, simmr_in$n_sources) + ), + mcmc_control = list( + iter = 10000, + burn = 1000, + thin = 10, + n.chain = 4 + )) { -# Main function to run simmr through JAGS -# if(class(simmr_in)!='simmr_input') stop("Input argument simmr_in must have come from simmr_load") + # Main function to run simmr through JAGS + # if(class(simmr_in)!='simmr_input') stop("Input argument simmr_in must have come from simmr_load") -# Throw warning if n.chain =1 -if(mcmc_control$n.chain==1) warning("Running only 1 MCMC chain will cause an error in the convergence diagnostics") + # Throw warning if n.chain =1 + if (mcmc_control$n.chain == 1) warning("Running only 1 MCMC chain will cause an error in the convergence diagnostics") -# Throw a warning if less than 4 observations in a group - 1 is ok as it wil do a solo run -if(min(table(simmr_in$group))>1 & min(table(simmr_in$group))<4) warning("At least 1 group has less than 4 observations - either put each observation in an individual group or use informative prior information") + # Throw a warning if less than 4 observations in a group - 1 is ok as it wil do a solo run + if (min(table(simmr_in$group)) > 1 & min(table(simmr_in$group)) < 4) warning("At least 1 group has less than 4 observations - either put each observation in an individual group or use informative prior information") -# Set up the model string -model_string = " + # Set up the model string + model_string <- " model{ # Likelihood for (j in 1:J) { @@ -304,68 +335,70 @@ model{ } " -output = vector('list',length=simmr_in$n_groups) -names(output) = levels(simmr_in$group) + output <- vector("list", length = simmr_in$n_groups) + names(output) <- levels(simmr_in$group) -# Loop through all the groups -for(i in 1:simmr_in$n_groups) { - if(simmr_in$n_groups>1) cat(paste("\nRunning for group",levels(simmr_in$group)[i],'\n\n')) - - curr_rows = which(simmr_in$group_int==i) - curr_mix = simmr_in$mixtures[curr_rows,,drop=FALSE] - - # Determine if a single observation or not - if(nrow(curr_mix)==1) { - cat('Only 1 mixture value, performing a simmr solo run...\n') - solo=TRUE - } else { - solo=FALSE - } - - # Create data object - data = with(simmr_in,list( - y=curr_mix, - s_mean=source_means, - s_sd=source_sds, - N=nrow(curr_mix), - J=n_tracers, - c_mean=correction_means, - c_sd = correction_sds, - q=concentration_means, - K=n_sources, - mu_f_mean=prior_control$means, - sigma_f_sd=prior_control$sd, - sig_upp=ifelse(solo,0.001,1000))) - - # Run in JAGS - output[[i]] = R2jags::jags(data=data, - parameters.to.save = c("p", "sigma"), - model.file = textConnection(model_string), - n.chains = mcmc_control$n.chain, - n.iter = mcmc_control$iter, - n.burnin = mcmc_control$burn, - n.thin = mcmc_control$thin) - - # Get the names right and interpretable everywhere - # Set the posterior names right - n_tracers = simmr_in$n_tracers - n_sources = simmr_in$n_sources - s_names = simmr_in$source_names - colnames(output[[i]]$BUGSoutput$sims.matrix)[2:(2 + n_sources - 1)] = s_names - colnames(output[[i]]$BUGSoutput$sims.list$p) = s_names - # Also do it in the summary - rownames(output[[i]]$BUGSoutput$summary)[2:(2 + n_sources - 1)] = s_names - sd_names = paste0('sd[',colnames(simmr_in$mixtures),']') - colnames(output[[i]]$BUGSoutput$sims.matrix)[(n_sources + 2):(n_sources + n_tracers + 1)] = sd_names - colnames(output[[i]]$BUGSoutput$sims.list$sigma) = sd_names - rownames(output[[i]]$BUGSoutput$summary)[(n_sources + 2):(n_sources + n_tracers + 1)] = sd_names -} + # Loop through all the groups + for (i in 1:simmr_in$n_groups) { + if (simmr_in$n_groups > 1) cat(paste("\nRunning for group", levels(simmr_in$group)[i], "\n\n")) + + curr_rows <- which(simmr_in$group_int == i) + curr_mix <- simmr_in$mixtures[curr_rows, , drop = FALSE] -output_all = vector('list') -output_all$input = simmr_in -output_all$output = output -class(output_all) = 'simmr_output' + # Determine if a single observation or not + if (nrow(curr_mix) == 1) { + cat("Only 1 mixture value, performing a simmr solo run...\n") + solo <- TRUE + } else { + solo <- FALSE + } + + # Create data object + data <- with(simmr_in, list( + y = curr_mix, + s_mean = source_means, + s_sd = source_sds, + N = nrow(curr_mix), + J = n_tracers, + c_mean = correction_means, + c_sd = correction_sds, + q = concentration_means, + K = n_sources, + mu_f_mean = prior_control$means, + sigma_f_sd = prior_control$sd, + sig_upp = ifelse(solo, 0.001, 1000) + )) + + # Run in JAGS + output[[i]] <- R2jags::jags( + data = data, + parameters.to.save = c("p", "sigma"), + model.file = textConnection(model_string), + n.chains = mcmc_control$n.chain, + n.iter = mcmc_control$iter, + n.burnin = mcmc_control$burn, + n.thin = mcmc_control$thin + ) + + # Get the names right and interpretable everywhere + # Set the posterior names right + n_tracers <- simmr_in$n_tracers + n_sources <- simmr_in$n_sources + s_names <- simmr_in$source_names + colnames(output[[i]]$BUGSoutput$sims.matrix)[2:(2 + n_sources - 1)] <- s_names + colnames(output[[i]]$BUGSoutput$sims.list$p) <- s_names + # Also do it in the summary + rownames(output[[i]]$BUGSoutput$summary)[2:(2 + n_sources - 1)] <- s_names + sd_names <- paste0("sd[", colnames(simmr_in$mixtures), "]") + colnames(output[[i]]$BUGSoutput$sims.matrix)[(n_sources + 2):(n_sources + n_tracers + 1)] <- sd_names + colnames(output[[i]]$BUGSoutput$sims.list$sigma) <- sd_names + rownames(output[[i]]$BUGSoutput$summary)[(n_sources + 2):(n_sources + n_tracers + 1)] <- sd_names + } -return(output_all) + output_all <- vector("list") + output_all$input <- simmr_in + output_all$output <- output + class(output_all) <- "simmr_output" + return(output_all) } diff --git a/R/simmr_mcmc_tdf.R b/R/simmr_mcmc_tdf.R index c15d8de..a802cf3 100644 --- a/R/simmr_mcmc_tdf.R +++ b/R/simmr_mcmc_tdf.R @@ -1,158 +1,190 @@ #' Estimate correction factors from stable isotope data with known dietary #' proportions -#' -#' This function runs a slightly different version of the main +#' +#' This function runs a slightly different version of the main #' \code{\link{simmr_mcmc}} function with the key difference that it estimates #' the correction factors (sometimes called trophic enrichment or trophic #' discrimination factors; TEFs/TDFs) for a given set of dietary proportions. -#' -#' The idea is that this code can be used for feeding studies where an +#' +#' The idea is that this code can be used for feeding studies where an #' organism is fed a known proportional diet with a view to estimating -#' the correction factors to be used in a later stable isotope mixing -#' model when the organisms are observed in the field. -#' -#' The main argument of the function is an object created from +#' the correction factors to be used in a later stable isotope mixing +#' model when the organisms are observed in the field. +#' +#' The main argument of the function is an object created from #' \code{\link{simmr_load}} which contains mixture data on a number of tracers -#' and food source means and standard deviations. Any correction factors +#' and food source means and standard deviations. Any correction factors #' included in this object will be ignored. The known dietary proportions should be provided for each individual (i.e. should be a matrix with the same number of rows as \code{mix}). It is advisable to have multiple different dietary proportion values as part of the feeding experimental design -#' +#' #' The output of the function is a posterior distribution on the correction #' factors for each food source. Just like the output from #' \code{\link{simmr_mcmc}}, this should be checked for convergence. Examples #' are included below to help assist with this check and further plots -#' +#' #' If, after running \code{\link{simmr_mcmc_tdf}} the convergence diagnostics in #' \code{\link{summary.simmr_output_tdf}} are not satisfactory, the values of #' \code{iter}, \code{burn} and \code{thin} in \code{mcmc_control} should be #' increased by e.g. a factor of 10. -#' +#' #' @param simmr_in An object created via the function \code{\link{simmr_load}} #' @param p The known dietary proportions for the feeding study. Dietary proportions should be given per individual (even if they are all identical) #' @param prior_control A list of values including arguments named \code{means} #' and \code{sd} which represent the prior means and standard deviations of the -#' correction factors. These can usually be left at their default values unless +#' correction factors. These can usually be left at their default values unless #' you wish to include to include prior information on them. #' @param mcmc_control A list of values including arguments named \code{iter} #' (number of iterations), \code{burn} (size of burn-in), \code{thin} (amount #' of thinning), and \code{n.chain} (number of MCMC chains). #' @return An object of class \code{simmr_tdf} with two named top-level -#' components: +#' components: #' \item{input}{The \code{simmr_input} object given to the -#' \code{simmr_mcmc} function} +#' \code{simmr_mcmc} function} #' \item{output}{A set of MCMC chains of class -#' \code{mcmc.list} from the coda package. These can be analysed using +#' \code{mcmc.list} from the coda package. These can be analysed using #' \code{\link{summary.simmr_output_tdf}}}. -#' +#' #' @author Andrew Parnell -#' +#' #' @seealso \code{\link{simmr_load}} for creating objects suitable for this -#' function, \code{\link{simmr_mcmc}} for estimating dietary proportions, +#' function, \code{\link{simmr_mcmc}} for estimating dietary proportions, #' \code{\link{plot.simmr_input}} for creating isospace plots, #' \code{\link{summary.simmr_output_tdf}} for summarising output -#' +#' #' @references Andrew C. Parnell, Donald L. Phillips, Stuart Bearhop, Brice X. #' Semmens, Eric J. Ward, Jonathan W. Moore, Andrew L. Jackson, Jonathan Grey, #' David J. Kelly, and Richard Inger. Bayesian stable isotope mixing models. #' Environmetrics, 24(6):387–399, 2013. -#' +#' #' Andrew C Parnell, Richard Inger, Stuart Bearhop, and Andrew L Jackson. #' Source partitioning using stable isotopes: coping with too much variation. #' PLoS ONE, 5(3):5, 2010. -#' -#' +#' +#' #' @examples #' \dontrun{ #' ## Example of estimating TDFs for a simple system with known dietary proportions -#' +#' #' # Data set 1: 10 obs on 2 isos, 4 sources, with tefs and concdep #' # Assume p = c(0.25, 0.25, 0.25, 0.25) -#' +#' #' # The data #' data(simmr_data_1) #' # Load into simmr -#' simmr_1 = with(simmr_data_1, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) -#' +#' simmr_1 <- with( +#' simmr_data_1, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) +#' #' # Plot #' plot(simmr_tdf) -#' +#' #' # MCMC run -#' simmr_tdf_out = simmr_mcmc_tdf(simmr_tdf, -#' p = matrix(rep(1/simmr_tdf$n_sources, -#' simmr_tdf$n_sources), -#' ncol = simmr_tdf$n_sources, -#' nrow = simmr_tdf$n_obs, -#' byrow = TRUE)) -#' +#' simmr_tdf_out <- simmr_mcmc_tdf(simmr_tdf, +#' p = matrix(rep( +#' 1 / simmr_tdf$n_sources, +#' simmr_tdf$n_sources +#' ), +#' ncol = simmr_tdf$n_sources, +#' nrow = simmr_tdf$n_obs, +#' byrow = TRUE +#' ) +#' ) +#' #' # Summary -#' summary(simmr_tdf_out,type='diagnostics') -#' summary(simmr_tdf_out,type='quantiles') -#' -#' # Now put these corrections back into the model and check the +#' summary(simmr_tdf_out, type = "diagnostics") +#' summary(simmr_tdf_out, type = "quantiles") +#' +#' # Now put these corrections back into the model and check the #' # iso-space plots and dietary output -#' simmr_tdf_2 = simmr_load(mixtures=mix, -#' source_names=s_names, -#' source_means=s_means, -#' source_sds=s_sds, -#' correction_means = simmr_tdf_out$c_mean_est, -#' correction_sds = simmr_tdf_out$c_sd_est, -#' concentration_means = conc) -#' +#' simmr_tdf_2 <- simmr_load( +#' mixtures = mix, +#' source_names = s_names, +#' source_means = s_means, +#' source_sds = s_sds, +#' correction_means = simmr_tdf_out$c_mean_est, +#' correction_sds = simmr_tdf_out$c_sd_est, +#' concentration_means = conc +#' ) +#' #' # Plot with corrections now #' plot(simmr_tdf_2) -#' -#' simmr_tdf_2_out = simmr_mcmc(simmr_tdf_2) -#' summary(simmr_tdf_2_out, type = 'diagnostics') -#' plot(simmr_tdf_2_out, type = 'boxplot') +#' +#' simmr_tdf_2_out <- simmr_mcmc(simmr_tdf_2) +#' summary(simmr_tdf_2_out, type = "diagnostics") +#' plot(simmr_tdf_2_out, type = "boxplot") #' } -#' +#' #' @export -simmr_mcmc_tdf = function(simmr_in, - p = matrix(rep(1/simmr_in$n_sources, - simmr_in$n_sources), - ncol = simmr_in$n_sources, - nrow = simmr_in$n_obs, - byrow = TRUE), - prior_control=list(c_mean_est=rep(2, - simmr_in$n_tracers), - c_sd_est=rep(2, - simmr_in$n_tracers)), - mcmc_control=list(iter=10000, - burn=1000, - thin=10, - n.chain=4)) { - UseMethod('simmr_mcmc_tdf') -} +simmr_mcmc_tdf <- function(simmr_in, + p = matrix(rep( + 1 / simmr_in$n_sources, + simmr_in$n_sources + ), + ncol = simmr_in$n_sources, + nrow = simmr_in$n_obs, + byrow = TRUE + ), + prior_control = list( + c_mean_est = rep( + 2, + simmr_in$n_tracers + ), + c_sd_est = rep( + 2, + simmr_in$n_tracers + ) + ), + mcmc_control = list( + iter = 10000, + burn = 1000, + thin = 10, + n.chain = 4 + )) { + UseMethod("simmr_mcmc_tdf") +} #' @export -simmr_mcmc_tdf.simmr_input = function(simmr_in, - p = matrix(rep(1/simmr_in$n_sources, - simmr_in$n_sources), - ncol = simmr_in$n_sources, - nrow = simmr_in$n_obs, - byrow = TRUE), - prior_control=list(c_mean_est=rep(2, - simmr_in$n_tracers), - c_sd_est=rep(2, - simmr_in$n_tracers)), - mcmc_control=list(iter=10000, - burn=1000, - thin=10, - n.chain=4)) { +simmr_mcmc_tdf.simmr_input <- function(simmr_in, + p = matrix(rep( + 1 / simmr_in$n_sources, + simmr_in$n_sources + ), + ncol = simmr_in$n_sources, + nrow = simmr_in$n_obs, + byrow = TRUE + ), + prior_control = list( + c_mean_est = rep( + 2, + simmr_in$n_tracers + ), + c_sd_est = rep( + 2, + simmr_in$n_tracers + ) + ), + mcmc_control = list( + iter = 10000, + burn = 1000, + thin = 10, + n.chain = 4 + )) { -# Throw warning if n.chain =1 -if(mcmc_control$n.chain==1) warning("Running only 1 MCMC chain will cause an error in the convergence diagnostics") + # Throw warning if n.chain =1 + if (mcmc_control$n.chain == 1) warning("Running only 1 MCMC chain will cause an error in the convergence diagnostics") -# Throw a warning if less than 4 observations in a group - 1 is ok as it wil do a solo run -if(min(table(simmr_in$group))>1 & min(table(simmr_in$group))<4) warning("At least 1 group has less than 4 observations - either put each observation in an individual group or use informative prior information") + # Throw a warning if less than 4 observations in a group - 1 is ok as it wil do a solo run + if (min(table(simmr_in$group)) > 1 & min(table(simmr_in$group)) < 4) warning("At least 1 group has less than 4 observations - either put each observation in an individual group or use informative prior information") -# Set up the model string -model_string = " + # Set up the model string + model_string <- " model { # Likelihood for (j in 1:J) { @@ -180,50 +212,52 @@ model { } " -if(simmr_in$n_groups > 1) stop("TDF calculation currently only works for single group data") - -# Loop through all the groups -curr_rows = which(simmr_in$group_int==1) -curr_mix = simmr_in$mixtures[curr_rows,,drop=FALSE] + if (simmr_in$n_groups > 1) stop("TDF calculation currently only works for single group data") -# Determine if a single observation or not -if(nrow(curr_mix)==1) { - cat('Only 1 mixture value, performing a simmr solo run...\n') - solo=TRUE -} else { - solo=FALSE -} + # Loop through all the groups + curr_rows <- which(simmr_in$group_int == 1) + curr_mix <- simmr_in$mixtures[curr_rows, , drop = FALSE] -# Create data object -data = with(simmr_in,list( - y=curr_mix, - p = p, - s_mean=source_means, - s_sd=source_sds, - N=nrow(curr_mix), - J=n_tracers, - q=concentration_means, - K=n_sources, - c_sd_est = prior_control$c_sd_est, - c_mean_est = prior_control$c_mean_est, - sig_upp=ifelse(solo,0.001,1000))) + # Determine if a single observation or not + if (nrow(curr_mix) == 1) { + cat("Only 1 mixture value, performing a simmr solo run...\n") + solo <- TRUE + } else { + solo <- FALSE + } -# Run in JAGS -output = R2jags::jags(data=data, - parameters.to.save = c("c_mean", "c_sd"), - model.file = textConnection(model_string), - n.chains = mcmc_control$n.chain, - n.iter = mcmc_control$iter, - n.burnin = mcmc_control$burn, - n.thin = mcmc_control$thin) + # Create data object + data <- with(simmr_in, list( + y = curr_mix, + p = p, + s_mean = source_means, + s_sd = source_sds, + N = nrow(curr_mix), + J = n_tracers, + q = concentration_means, + K = n_sources, + c_sd_est = prior_control$c_sd_est, + c_mean_est = prior_control$c_mean_est, + sig_upp = ifelse(solo, 0.001, 1000) + )) -output_all = vector('list') -output_all$input = simmr_in -output_all$output = output -output_all$c_mean_est = output$BUGSoutput$median$c_mean -output_all$c_sd_est = output$BUGSoutput$median$c_sd -class(output_all) = c('simmr_output_tdf') + # Run in JAGS + output <- R2jags::jags( + data = data, + parameters.to.save = c("c_mean", "c_sd"), + model.file = textConnection(model_string), + n.chains = mcmc_control$n.chain, + n.iter = mcmc_control$iter, + n.burnin = mcmc_control$burn, + n.thin = mcmc_control$thin + ) -return(output_all) + output_all <- vector("list") + output_all$input <- simmr_in + output_all$output <- output + output_all$c_mean_est <- output$BUGSoutput$median$c_mean + output_all$c_sd_est <- output$BUGSoutput$median$c_sd + class(output_all) <- c("simmr_output_tdf") + return(output_all) } diff --git a/R/summary.simmr_output.R b/R/summary.simmr_output.R index 49332f2..67fefd3 100644 --- a/R/summary.simmr_output.R +++ b/R/summary.simmr_output.R @@ -1,28 +1,28 @@ #' Summarises the output created with \code{\link{simmr_mcmc}} -#' +#' #' Produces textual summaries and convergence diagnostics for an object created #' with \code{\link{simmr_mcmc}}. The different options are: 'diagnostics' #' which produces Brooks-Gelman-Rubin diagnostics to assess MCMC convergence, #' 'quantiles' which produces credible intervals for the parameters, #' 'statistics' which produces means and standard deviations, and #' 'correlations' which produces correlations between the parameters. -#' +#' #' The quantile output allows easy calculation of 95 per cent credible #' intervals of the posterior dietary proportions. The correlations, along with #' the matrix plot in \code{\link{plot.simmr_output}} allow the user to judge #' which sources are non-identifiable. The Gelman diagnostic values should be #' close to 1 to ensure satisfactory convergence. -#' +#' #' When multiple groups are included, the output automatically includes the #' results for all groups. -#' +#' #' @param object An object of class \code{simmr_output} produced by the #' function \code{\link{simmr_mcmc}} #' @param type The type of output required. At least none of 'diagnostics', #' 'quantiles', 'statistics', or 'correlations'. #' @param group Which group or groups the output is required for. #' @param ... Not used -#' @return An list containing the following components: \item{gelman }{The +#' @return A list containing the following components: \item{gelman }{The #' convergence diagnostics} \item{quantiles }{The quantiles of each parameter #' from the posterior distribution} \item{statistics }{The means and standard #' deviations of each parameter} \item{correlations }{The posterior @@ -34,95 +34,99 @@ #' function, and many more examples. See also \code{\link{simmr_load}} for #' creating simmr objects, \code{\link{plot.simmr_input}} for creating isospace #' plots, \code{\link{plot.simmr_output}} for plotting output. -#' +#' #' @importFrom stats sd cor -#' +#' #' @examples #' \dontrun{ #' # A simple example with 10 observations, 2 tracers and 4 sources -#' +#' #' # The data #' data(geese_data_day1) -#' simmr_1 = with(geese_data_day1, -#' simmr_load(mixtures=mixtures, -#' source_names=source_names, -#' source_means=source_means, -#' source_sds=source_sds, -#' correction_means=correction_means, -#' correction_sds=correction_sds, -#' concentration_means = concentration_means)) -#' +#' simmr_1 <- with( +#' geese_data_day1, +#' simmr_load( +#' mixtures = mixtures, +#' source_names = source_names, +#' source_means = source_means, +#' source_sds = source_sds, +#' correction_means = correction_means, +#' correction_sds = correction_sds, +#' concentration_means = concentration_means +#' ) +#' ) +#' #' # Plot #' plot(simmr_1) -#' -#' +#' +#' #' # MCMC run -#' simmr_1_out = simmr_mcmc(simmr_1) -#' -#' # Summarise -#' summary(simmr_1_out) # This outputs all the summaries -#' summary(simmr_1_out,type='diagnostics') # Just the diagnostics +#' simmr_1_out <- simmr_mcmc(simmr_1) +#' +#' # Summarise +#' summary(simmr_1_out) # This outputs all the summaries +#' summary(simmr_1_out, type = "diagnostics") # Just the diagnostics #' # Store the output in an -#' ans = summary(simmr_1_out, -#' type=c('quantiles','statistics')) object +#' ans <- summary(simmr_1_out, +#' type = c("quantiles", "statistics") +#' ) #' } #' @export -summary.simmr_output = - function(object,type=c('diagnostics','quantiles','statistics','correlations'),group=1,...) { +summary.simmr_output <- + function(object, type = c("diagnostics", "quantiles", "statistics", "correlations"), group = 1, ...) { # Get the specified type - type=match.arg(type,several.ok=TRUE) + type <- match.arg(type, several.ok = TRUE) # Set up containers - out_bgr = out_quantiles = out_statistics = out_cor = vector('list',length=length(group)) - group_names = levels(object$input$group) - names(out_bgr) = paste0('group_',group) - names(out_quantiles) = paste0('group_',group) - names(out_statistics) = paste0('group_',group) - names(out_cor) = paste0('group_',group) + out_bgr <- out_quantiles <- out_statistics <- out_cor <- vector("list", length = length(group)) + group_names <- levels(object$input$group) + names(out_bgr) <- paste0("group_", group) + names(out_quantiles) <- paste0("group_", group) + names(out_statistics) <- paste0("group_", group) + names(out_cor) <- paste0("group_", group) # Loop through groups - for(i in 1:length(group)) { - - cat(paste("\nSummary for",group_names[group[i]],'\n')) - out_all = object$output[[group[i]]]$BUGSoutput$sims.matrix + for (i in 1:length(group)) { + cat(paste("\nSummary for", group_names[group[i]], "\n")) + out_all <- object$output[[group[i]]]$BUGSoutput$sims.matrix # Get objects - out_bgr[[i]] = object$output[[i]]$BUGSoutput$summary[,'Rhat'] - out_quantiles[[i]] = t(apply(out_all,2,'quantile',probs=c(0.025,0.25,0.5,0.75,0.975))) + out_bgr[[i]] <- object$output[[i]]$BUGSoutput$summary[, "Rhat"] + out_quantiles[[i]] <- t(apply(out_all, 2, "quantile", probs = c(0.025, 0.25, 0.5, 0.75, 0.975))) # coda:::summary.mcmc.list(object$output)$quantiles - out_statistics[[i]] = t(apply(out_all,2,function(x) {return(c(mean=mean(x),sd=stats::sd(x)))})) + out_statistics[[i]] <- t(apply(out_all, 2, function(x) { + return(c(mean = mean(x), sd = stats::sd(x))) + })) # coda:::summary.mcmc.list(object$output)$statistics[,1:2] - out_cor[[i]] = stats::cor(out_all) + out_cor[[i]] <- stats::cor(out_all) - if ('diagnostics'%in%type) { + if ("diagnostics" %in% type) { # Print out gelman diagnostics of the output - cat('Gelman diagnostics - these values should all be close to 1.\n') - cat('If not, try a longer run of simmr_mcmc.\n') - print(round(out_bgr[[i]],2)) + cat("Gelman diagnostics - these values should all be close to 1.\n") + cat("If not, try a longer run of simmr_mcmc.\n") + print(round(out_bgr[[i]], 2)) } - if ('quantiles'%in%type) { + if ("quantiles" %in% type) { # Print out quantiles argument - print(round(out_quantiles[[i]],3)) + print(round(out_quantiles[[i]], 3)) } - if ('statistics'%in%type) { + if ("statistics" %in% type) { # Print out quantiles argument - print(round(out_statistics[[i]],3)) + print(round(out_statistics[[i]], 3)) } - if ('correlations'%in%type) { + if ("correlations" %in% type) { # Print out quantiles argument - print(round(out_cor[[i]],3)) + print(round(out_cor[[i]], 3)) } - } - if(object$input$n_groups==1) { - invisible(list(gelman=out_bgr[[1]],quantiles=out_quantiles[[1]],statistics=out_statistics[[1]],correlations=out_cor[[1]])) + if (object$input$n_groups == 1) { + invisible(list(gelman = out_bgr[[1]], quantiles = out_quantiles[[1]], statistics = out_statistics[[1]], correlations = out_cor[[1]])) } else { - invisible(list(gelman=out_bgr,quantiles=out_quantiles,statistics=out_statistics,correlations=out_cor)) + invisible(list(gelman = out_bgr, quantiles = out_quantiles, statistics = out_statistics, correlations = out_cor)) } - } diff --git a/R/summary.simmr_output_tdf.R b/R/summary.simmr_output_tdf.R index a1c68da..c900a5c 100644 --- a/R/summary.simmr_output_tdf.R +++ b/R/summary.simmr_output_tdf.R @@ -1,17 +1,17 @@ #' Summarises the output created with \code{\link{simmr_mcmc_tdf}} -#' +#' #' Produces textual summaries and convergence diagnostics for an object created #' with \code{\link{simmr_mcmc_tdf}}. The different options are: 'diagnostics' #' which produces Brooks-Gelman-Rubin diagnostics to assess MCMC convergence, #' 'quantiles' which produces credible intervals for the parameters, #' 'statistics' which produces means and standard deviations, and #' 'correlations' which produces correlations between the parameters. -#' +#' #' The quantile output allows easy calculation of 95 per cent credible #' intervals of the posterior dietary proportions. The Gelman diagnostic values should be close to 1 to ensure satisfactory convergence. -#' +#' #' Multiple groups are not currently supported for estimating TDFs. -#' +#' #' @param object An object of class \code{simmr_output_df} produced by the #' function \code{\link{simmr_mcmc_tdf}} #' @param type The type of output required. At least none of 'diagnostics', @@ -24,54 +24,55 @@ #' correlations between the parameters} Note that this object is reported #' silently so will be discarded unless the function is called with an object #' as in the example below. -#' +#' #' @author Andrew Parnell -#' +#' #' @seealso See \code{\link{simmr_mcmc_tdf}} for creating objects suitable for this #' function, and many more examples. See also \code{\link{simmr_load}} for #' creating simmr objects, \code{\link{plot.simmr_input}} for creating isospace #' plots, \code{\link{plot.simmr_output}} for plotting output. -#' +#' #' @importFrom stats sd cor -#' +#' #' @export -summary.simmr_output_tdf = - function(object,type=c('diagnostics','quantiles','statistics','correlations'),...) { +summary.simmr_output_tdf <- + function(object, type = c("diagnostics", "quantiles", "statistics", "correlations"), ...) { # Get the specified type - type=match.arg(type,several.ok=TRUE) + type <- match.arg(type, several.ok = TRUE) - out_all = object$output$BUGSoutput$sims.matrix + out_all <- object$output$BUGSoutput$sims.matrix # Get objects - out_bgr = object$output$BUGSoutput$summary[,'Rhat'] - out_quantiles = t(apply(out_all,2,'quantile',probs=c(0.025,0.25,0.5,0.75,0.975))) + out_bgr <- object$output$BUGSoutput$summary[, "Rhat"] + out_quantiles <- t(apply(out_all, 2, "quantile", probs = c(0.025, 0.25, 0.5, 0.75, 0.975))) # coda:::summary.mcmc.list(object$output)$quantiles - out_statistics = t(apply(out_all,2,function(x) {return(c(mean=mean(x),sd=stats::sd(x)))})) + out_statistics <- t(apply(out_all, 2, function(x) { + return(c(mean = mean(x), sd = stats::sd(x))) + })) # coda:::summary.mcmc.list(object$output)$statistics[,1:2] - out_cor = stats::cor(out_all) + out_cor <- stats::cor(out_all) - if ('diagnostics'%in%type) { + if ("diagnostics" %in% type) { # Print out gelman diagnostics of the output - cat('Gelman diagnostics - these values should all be close to 1.\n') - cat('If not, try a longer run of simmr_mcmc_tdf.\n') - print(round(out_bgr,2)) + cat("Gelman diagnostics - these values should all be close to 1.\n") + cat("If not, try a longer run of simmr_mcmc_tdf.\n") + print(round(out_bgr, 2)) } - if ('quantiles'%in%type) { + if ("quantiles" %in% type) { # Print out quantiles argument - print(round(out_quantiles,3)) + print(round(out_quantiles, 3)) } - if ('statistics'%in%type) { + if ("statistics" %in% type) { # Print out quantiles argument - print(round(out_statistics,3)) + print(round(out_statistics, 3)) } - if ('correlations'%in%type) { + if ("correlations" %in% type) { # Print out quantiles argument - print(round(out_cor,3)) + print(round(out_cor, 3)) } - invisible(list(gelman=out_bgr,quantiles=out_quantiles,statistics=out_statistics,correlations=out_cor)) - -} + invisible(list(gelman = out_bgr, quantiles = out_quantiles, statistics = out_statistics, correlations = out_cor)) + } diff --git a/cran-comments.md b/cran-comments.md index 6fd51cd..fa4800a 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,6 +1,6 @@ ## Test environments -* local OS X install, R 4.0.0 -* travis +* local OS X install, R 4.0.3 +* GitHub actions * rhub * win-builder (devel and release) diff --git a/docs/404.html b/docs/404.html index 9306fe4..ccf8094 100644 --- a/docs/404.html +++ b/docs/404.html @@ -71,7 +71,7 @@ simmr - 0.4.3.9000 + 0.4.4 @@ -79,7 +79,7 @@