A second contribution is a discussion of model averaging applied to SSDs. The literature cited is useful. It would have been helpful to include a few detailed examples showing the utility or dangers of such an approach. It would also have been helpful to include more details on the software packages cited, at least for the two they say they focused on.

8 Responses

  • RESOLUTION

    Adopted suggestions to incorporate statement about Schwarz and Tillmanns (2019) :

    
    Schwarz and Tillmans (2019) used data sets extracted from the CCME Guidelines for the Protection of Aquatic Life for boron (CCME 2009) and silver (CCME 2015) to assess and compare the results from model averaged and single SSDs. Among other things, they concluded that model averaging can reduce the uncertainty associated with fitting distributions to small data sets as well as providing some immunity to perturbations in HCx values due to the influence of a single sensitive data point.

  • An example (e.g. using the ERISS ammonia dataset) could be put together and included as Supplementary material, to avoid word limit issues. However, it would still need someone’s time to draft it.

  • I found an email thread from when we were exploring the issue of example datasets. I think we had decided on total ammonia as the best option before David decided it was getting too much for the paper.

    My code for exploring this was:
    library(ssdtools)
    library(ggplot2)
    source(“lnormlnorm.R”)
    require(tidyverse)

    # Total ammonia —-
    # this is the output of dput, which is used to create a data.frame from data entered in interactive spreadsheet
    data <- data.frame(Concentration = c(16.12, 3.27, 0.47, 4.61, 4.23, 1.13,
                                            0.5, 0.6))#,
                                            #495, 431.5, 87, 27.9, 25.4, 172.6, 37.6))
    # fix unacceptable column names
    colnames(data) <- make.names(colnames(data))
    data$Species <- NA

    #data <- rbind(data,data)

    # fit distributions
    dist <- ssd_fit_dists(data, left = ‘Concentration’, dists = c(‘llogis’, ‘gamma’, ‘lnorm’, ‘lgumbel’, ‘weibull’,
                                                                 ‘lnormlnorm’))
    # plot distributions
    ssdtools::ssd_plot_cdf(dist)
    # goodness of fit table
    gof.1 <- ssd_gof(dist)

    # width and height are in inches, dpi (dots per inch) sets resolution
    ggsave(‘fit_dist_plot_total_amonia.png’, width = 8 , height = 6 , dpi = 300)

    # plot model average
    # to add confidence intervals set ci = TRUE in predict and ssd_plot
    # we recommend using nboot = 10000 in predict, although this may take several minutes to run
    pred <- predict(dist, nboot = 10L)
    ssd_plot(data, pred, left = ‘Concentration’, label = ‘Species’, color = NULL, shape = NULL, hc = 5L, ci = FALSE,
            shift_x = 1.3, xlab = ‘Concentration’, ylab = ‘Percent of Species Affected’) +
     ggtitle(”) +
     theme(panel.border = element_blank(),
           panel.grid.major = element_blank(),
           panel.grid.minor = element_blank(),
           panel.background = element_rect(fill = NA, colour=’black’),
           axis.text = element_text(color = ‘black’),
           legend.key = element_rect(fill = NA, colour = NA)) +
     expand_limits(x = 16.12) +
     scale_color_brewer(palette = ‘Dark2’, name = ‘-none-‘) +
     scale_shape(name = NULL)

    # save plot
    # width and height are in inches, dpi (dots per inch) sets resolution
    ggsave(‘model_average_plot.png’, width = 8 , height = 6 , dpi = 600)

    # get confidence limits
    # use the nboot argument in ssd_hc to set the number of bootstrap samples
    dist <- ssd_fit_dists(data, left = ‘Concentration’, dists = c(‘llogis’, ‘gamma’, ‘lnorm’, ‘lgumbel’, ‘weibull’))
    hc.vals <- ssd_hc(dist, percent = 5L, ci = TRUE, nboot = 500L)
    hc5.all <- data.frame(do.call(“rbind”,lapply(dist,FUN=ssd_hc, percent = 5L, ci = TRUE, nboot = 500L)))
    hc1.all <- data.frame(do.call(“rbind”,lapply(dist,FUN=ssd_hc, percent = 1L, ci = TRUE, nboot = 500L)))

    hc.out <- data.frame(‘dist’=hc5.all$dist, ‘HC5’=hc5.all$est, ‘hc1’=hc1.all$est)

    stats.out <- data.frame(gof.1) %>%
     left_join(hc.out)

    stats.out
    write.csv(stats.out, “Total ammonia – gof with lnormlnorm.csv”)

  • Given the word limits in the journal and the number of topics covered in the MS adding detailed examples would be very difficult. We could add a model averaged estimate for the Fig1 data – although that example doesn’t actually show much difference in the HC5 estimates across the four distributions:

      percent   est    se   lcl   ucl dist   
*   <dbl> <dbl> <dbl> <dbl> <dbl> <chr>  
1       5  4.23  2.22  1.83 10.3  average
2       5  4.05  1.85  2.10  9.12 lnorm  
3       5  4.32  2.43  1.49 10.2  weibull
4       5  4.53  2.01  1.99  9.70 llogis 
5       5  4.00  2.27  1.99 10.5  gamma

    The weights for that example are:

    # A tibble: 4 x 9
  dist       ad    ks   cvm   aic  aicc   bic delta weight
  <chr>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
1 lnorm      NA    NA    NA  44.5  48.5  44.1 2.07   0.152
2 weibull    NA    NA    NA  42.5  46.5  42.0 0      0.427
3 llogis     NA    NA    NA  44.5  48.5  44.1 2.04   0.154
4 gamma      NA    NA    NA  43.4  47.4  43.0 0.944  0.267

    Another relatively simple option would be the boron dataset because it is already available on the web, but again – not a lot of difference across the distributions:

      percent   est    se   lcl   ucl dist   
*   <dbl> <dbl> <dbl> <dbl> <dbl> <chr>  
1       5  1.22 0.817 0.450  3.52 average
2       5  1.68 0.709 0.917  3.71 lnorm  
3       5  1.09 0.744 0.405  3.37 weibull
4       5  1.56 0.714 0.665  3.42 llogis

    We do plan a further MS based on specific examples and simulations to explore these issues.
    Rather than adding an example in this MS, perhaps we could just expand the detail around the ssdtools package and refer to the testing that was done in Schwarz and Tilmlmanns (2019). A couple of key lines in their executive summary that stand out are “The model averaging approach showed greater stability than single distributions when data points from the middle or right side of the distribution were added” and that model averaging can “also reduce the uncertainty associated to fitting distributions to small data sets”.

    • Perhaps citing Schwarz and Tillmanns and providing a few sentences that describe how model averaging was evaluated and the results.

Comments are closed.