Merge pull request #372 from poissonconsulting/dev

Switched to `ssd_min_pmix()`

.gitignore

@@ -12,6 +12,7 @@
 
 docs/
 inst/docs/
-scripts/
 inst/doc
 paper/paper.pdf
+scripts/boron_samples1
+scripts/checks.md

DESCRIPTION

@@ -73,7 +73,6 @@ Imports:
 Suggests:
     actuar,
     covr,
-    car,
     doFuture,
     dplyr,
     EnvStats,
@@ -83,7 +82,6 @@ Suggests:
     future,
     glue,
     grDevices,
-    ggpubr,
     knitr,
     magrittr,
     mle.tools,
@@ -110,4 +108,4 @@ Encoding: UTF-8
 Language: en-US
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2

NAMESPACE

@@ -55,7 +55,6 @@ export(geom_xribbon)
 export(glance)
 export(is.fitdists)
 export(is_censored)
-export(licensing_md)
 export(logLik)
 export(nobs)
 export(npars)
@@ -97,7 +96,9 @@ export(ssd_hc_burrlioz)
 export(ssd_hp)
 export(ssd_hp_bcanz)
 export(ssd_is_censored)
+export(ssd_licensing_md)
 export(ssd_match_moments)
+export(ssd_min_pmix)
 export(ssd_pal)
 export(ssd_pburrIII3)
 export(ssd_pgamma)

R/fit.R

@@ -168,7 +168,7 @@ ssd_fit_dists <- function(
     computable = TRUE,
     at_boundary_ok = FALSE,
     all_dists = FALSE,
-    min_pmix = 0,
+    min_pmix = ssd_min_pmix(nrow(data)),
     range_shape1 = c(0.05, 20),
     range_shape2 = range_shape1,
     control = list(),

R/licensing.R

@@ -4,8 +4,8 @@
 #' and documentation
 #' @export
 #' @examples
-#' licensing_md()
-licensing_md <- function() {
+#' ssd_licensing_md()
+ssd_licensing_md <- function() {
 "## Licensing
   Copyright 2018-2024 Province of British Columbia\\
   Copyright 2021 Environment and Climate Change Canada\\

R/min-pmix.R

@@ -0,0 +1,16 @@
+#' Calculate Minimum Proportion in Mixture Models
+#'
+#' @inheritParams params
+#'
+#' @return A number between 0 and 0.5 of the minimum proportion in mixture models.
+#' @seealso [ssd_fit_dists()]
+#' @export
+#'
+#' @examples
+#' ssd_min_pmix(6)
+#' ssd_min_pmix(50)
+ssd_min_pmix <- function(n) {
+  chk_whole_number(n)
+  chk_gt(n)
+  max(min(3/n, 0.5), 0.1)
+}

R/multi.R

@@ -198,6 +198,9 @@ ssd_qmulti <- function(
 #' # multi 
 #' set.seed(50)
 #' hist(ssd_rmulti(1000), breaks = 100)
+#' 
+#' fits <- ssd_fit_dists(ssddata::ccme_boron)
+#' do.call("ssd_rmulti", c(n = 10, estimates(fits)))
 ssd_rmulti <- function(
     n,     
     burrIII3.weight = 0, 

R/params.R

@@ -60,6 +60,7 @@
 #' proportion of bootstrap samples that must successfully fit (return a likelihood) 
 #' to report the confidence intervals.
 #' @param min_pmix A number between 0 and 0.5 specifying the minimum proportion in mixture models.
+#' @param n A whole number of the effective number of rows of data.
 #' @param npars A whole numeric vector specifying which distributions to include based on the number of parameters.
 #' @param all_estimates A flag specifying whether to calculate estimates for all implemented distributions.
 #' @param ci_method A string specifying which method to use for estimating the bootstrap values. 

R/predict.R

@@ -65,6 +65,7 @@ predict.fitdists <- function(
     level = level,
     nboot = nboot, 
     min_pboot = min_pboot,
+    multi_est = multi_est,
     average = average, 
     delta = delta, 
     parametric = parametric,

README.Rmd

@@ -144,5 +144,5 @@ If you would like to contribute to the package, please see our
 Please note that the ssdtools project is released with a [Contributor Code of Conduct](https://contributor-covenant.org/version/2/1/CODE_OF_CONDUCT.html). By contributing to this project, you agree to abide by its terms.
 
 ```{r, results = "asis", echo = FALSE}
-cat(licensing_md())
+cat(ssd_licensing_md())
 ```

_pkgdown.yml

@@ -83,7 +83,7 @@ reference:
   - '`comma_signif`'
   - '`dist_data`'
   - '`is.fitdists`'
-  - '`licensing_md`'
+  - '`ssd_licensing_md`'
   - '`npars`'
   - '`ssd_data`'
   - '`ssd_ecd`'

scripts/sims.R

@@ -0,0 +1,56 @@
+library(ssddata)
+library(ssdtools)
+library(purrr)
+library(stringr)
+library(future.apply)
+
+rm(list = ls())
+
+plan("multisession")
+
+nsims <- 200
+ndata <- 1024
+dists <- c("lnorm", "lgumbel")
+
+fit <- ssd_fit_dists(ssddata::ccme_glyphosate, dists = dists)
+ssd_plot_cdf(fit)
+true <- ssd_hc(fit)$est
+true_ave <- ssd_hc(fit, multi_est = FALSE)$est
+ssd_hc(fit, average = FALSE)
+
+estimates <- estimates(fit, all_estimates = TRUE)
+
+set.seed(42)
+data <- do.call("ssd_rmulti", c(n = nsims * ndata, estimates))
+data <- matrix(data, ncol = nsims)
+
+e5 <- future_apply(data, MARGIN = 2, function(x) quantile(x, 0.05))
+hist(e5)
+e1 <- future_apply(data, MARGIN = 2, function(x) quantile(x, 0.01))
+hist(e1)
+e20 <- future_apply(data, MARGIN = 2, function(x) quantile(x, 0.2))
+hist(e20)
+emax <- future_apply(data, MARGIN = 2, function(x) max(x))
+hist(log(emax))
+
+log(max(ssd_rlgumbel(round(204800 * 0.719), locationlog = 8.357433, scalelog = 1.223079)))
+
+ssd_plot_data(data.frame(Conc = data[,10]))
+
+fits <- future_apply(data, MARGIN = 2, function(x) ssd_fit_dists(data.frame(Conc = x), dists = dists))
+
+hc5s <- vapply(fits, function(x) { ssd_hc(x, delta = Inf)$est }, 1)
+
+props <- (hc5s - true) / true
+hist(props)
+
+which_biased <- which.min(props)
+fit_biased <- fits[[which_biased]]
+ssd_hc(fit_biased, average = FALSE, delta = Inf)
+ssd_plot_cdf(fit_biased, average = FALSE, delta = Inf)
+
+which_not_biased <- which.min(abs(props-0))
+fit_not_biased <- fits[[which_not_biased]]
+ssd_hc(fit_not_biased, average = FALSE, delta = Inf)
+ssd_plot_cdf(fit_not_biased, average = FALSE, delta = Inf)
+

tests/testthat/_snaps/data/boron_unstable.csv

@@ -3,8 +3,8 @@ gompertz,location,0.0394073,0.0119791
 gompertz,shape,0.00260432,0.00997275
 invpareto,scale,75.2608,0
 invpareto,shape,0.568403,0.107418
-llogis_llogis,locationlog1,0.896677,0.379787
-llogis_llogis,locationlog2,3.14911,0.288357
-llogis_llogis,pmix,0.255182,0.144381
-llogis_llogis,scalelog1,0.317596,0.168008
-llogis_llogis,scalelog2,0.496754,0.13515
+llogis_llogis,locationlog1,0.896787,0.379717
+llogis_llogis,locationlog2,3.14918,0.288289
+llogis_llogis,pmix,0.255219,0.14435
+llogis_llogis,scalelog1,0.317647,0.167982
+llogis_llogis,scalelog2,0.49671,0.135114

tests/testthat/_snaps/fit/tidy_stable_anon_e.csv

@@ -7,10 +7,10 @@ llogis,locationlog,6.53111,0.964136
 llogis,scalelog,2.18789,0.423632
 lnorm,meanlog,6.2366,0.858092
 lnorm,sdlog,3.538,0.606761
-lnorm_lnorm,meanlog1,2.22857,0.541818
-lnorm_lnorm,meanlog2,9.04838,0.266599
-lnorm_lnorm,pmix,0.412305,0.119536
-lnorm_lnorm,sdlog1,1.42155,0.397676
-lnorm_lnorm,sdlog2,0.840691,0.188728
+lnorm_lnorm,meanlog1,2.22869,0.541909
+lnorm_lnorm,meanlog2,9.04841,0.266623
+lnorm_lnorm,pmix,0.412379,0.119543
+lnorm_lnorm,sdlog1,1.42176,0.39782
+lnorm_lnorm,sdlog2,0.84076,0.188761
 weibull,scale,2701.04,1846.65
 weibull,shape,0.372736,0.0768177

tests/testthat/_snaps/fit/tidy_stable_computable.csv

@@ -5,8 +5,8 @@ llogis,locationlog,0.012,0.278
 llogis,scalelog,0.0109,0.253
 lnorm,meanlog,0.0233,1.48
 lnorm,sdlog,0.0345,2.19
-lnorm_lnorm,meanlog1,0.0236,NA
-lnorm_lnorm,meanlog2,0.879,NA
-lnorm_lnorm,pmix,1,NA
-lnorm_lnorm,sdlog1,0.035,NA
-lnorm_lnorm,sdlog2,0.623,NA
+lnorm_lnorm,meanlog1,0.0199,NA
+lnorm_lnorm,meanlog2,0.329,NA
+lnorm_lnorm,pmix,0.885,NA
+lnorm_lnorm,sdlog1,0.0346,NA
+lnorm_lnorm,sdlog2,1.57e-4,NA

tests/testthat/_snaps/fit/tidy_stable_rescale.csv

@@ -7,10 +7,10 @@ llogis,locationlog,0.497053,0.248257
 llogis,scalelog,0.740426,0.114375
 lnorm,meanlog,0.432508,0.234627
 lnorm,sdlog,1.24153,0.165905
-lnorm_lnorm,meanlog1,-1.17973,0.318005
-lnorm_lnorm,meanlog2,1.07186,0.252917
-lnorm_lnorm,pmix,0.283997,0.123401
-lnorm_lnorm,sdlog1,0.554526,0.212462
-lnorm_lnorm,sdlog2,0.768811,0.193728
+lnorm_lnorm,meanlog1,-1.17973,0.317999
+lnorm_lnorm,meanlog2,1.07186,0.252923
+lnorm_lnorm,pmix,0.283988,0.123401
+lnorm_lnorm,sdlog1,0.554502,0.212446
+lnorm_lnorm,sdlog2,0.768838,0.19374
 weibull,scale,2.79651,0.577423
 weibull,shape,0.9661,0.14542

tests/testthat/_snaps/fit/tidy_unstable_anon_e.csv

@@ -1,8 +1,8 @@
 dist,term,est,se
 invpareto,scale,38106.2,0
 invpareto,shape,0.231937,0.0562529
-llogis_llogis,locationlog1,2.26675,0.557708
-llogis_llogis,locationlog2,9.14324,0.257326
-llogis_llogis,pmix,0.41598,0.121011
-llogis_llogis,scalelog1,0.839584,0.296992
-llogis_llogis,scalelog2,0.461779,0.129323
+llogis_llogis,locationlog1,2.26716,0.557743
+llogis_llogis,locationlog2,9.14327,0.257326
+llogis_llogis,pmix,0.415978,0.121012
+llogis_llogis,scalelog1,0.839608,0.297029
+llogis_llogis,scalelog2,0.461775,0.129324

tests/testthat/_snaps/hc-root.md

@@ -38,7 +38,7 @@
         dist    proportion   est    se   lcl   ucl    wt method    nboot pboot samples
         <chr>        <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>     <dbl> <dbl> <I<lis>
       1 average       0.05  1.26 0.622 0.492  2.13     1 parametr~    10     1 <dbl>  
-      2 average       0.1   2.38 0.928 1.17   3.59     1 parametr~    10     1 <dbl>  
+      2 average       0.1   2.38 0.926 1.17   3.58     1 parametr~    10     1 <dbl>  
 
 # hc multi_ci lnorm ci
 

tests/testthat/_snaps/hc/censored_5ll.csv

@@ -1,2 +1,2 @@
 dist,proportion,est,se,lcl,ucl,wt,method,nboot,pboot,samples
-average,0.05,1.36598,NA,NA,NA,1,parametric,0,NA,numeric(0)
+average,0.05,1.36601,NA,NA,NA,1,parametric,0,NA,numeric(0)

tests/testthat/_snaps/hc/hc.csv

@@ -3,5 +3,5 @@ gamma,0.05,1.07428,0.983875,0.709111,3.19472,0.356574,parametric,10,1,"c(`000000
 lgumbel,0.05,1.76939,0.401763,1.39032,2.63032,0.0134466,parametric,10,1,"c(`000000001_lgumbel` = 2.2833, `000000002_lgumbel` = 1.48094, `000000003_lgumbel` = 2.12922, `000000004_lgumbel` = 2.36856, `000000005_lgumbel` = 1.81154, `000000006_lgumbel` = 1.96535, `000000007_lgumbel` = 1.36401, `000000008_lgumbel` = 1.93797, `000000009_lgumbel` = 2.70632, `000000010_lgumbel` = 2.09232)"
 llogis,0.05,1.56226,0.834776,0.851639,3.37274,0.0656452,parametric,10,1,"c(`000000001_llogis` = 0.751505, `000000002_llogis` = 3.04268, `000000003_llogis` = 2.10953, `000000004_llogis` = 2.22634, `000000005_llogis` = 1.30249, `000000006_llogis` = 2.52802, `000000007_llogis` = 3.46857, `000000008_llogis` = 2.04533, `000000009_llogis` = 1.85618, `000000010_llogis` = 1.19654)"
 lnorm,0.05,1.68118,0.823602,0.889403,3.02109,0.177236,parametric,10,1,"c(`000000001_lnorm` = 3.09184, `000000002_lnorm` = 2.42899, `000000003_lnorm` = 1.325, `000000004_lnorm` = 1.61081, `000000005_lnorm` = 2.6033, `000000006_lnorm` = 0.865976, `000000007_lnorm` = 2.77742, `000000008_lnorm` = 1.19715, `000000009_lnorm` = 2.45546, `000000010_lnorm` = 0.970096)"
-lnorm_lnorm,0.05,1.54136,0.32909,0.952197,1.82189,0.0296268,parametric,10,1,"c(`000000001_lnorm_lnorm` = 1.64162, `000000002_lnorm_lnorm` = 1.67905, `000000003_lnorm_lnorm` = 1.80875, `000000004_lnorm_lnorm` = 0.921876, `000000005_lnorm_lnorm` = 1.68355, `000000006_lnorm_lnorm` = 1.28521, `000000007_lnorm_lnorm` = 1.82571, `000000008_lnorm_lnorm` = 1.05664, `000000009_lnorm_lnorm` = 1.20995, `000000010_lnorm_lnorm` = 1.67582)"
+lnorm_lnorm,0.05,1.54136,0.329201,0.952193,1.82182,0.0296268,parametric,10,1,"c(`000000001_lnorm_lnorm` = 1.64164, `000000002_lnorm_lnorm` = 1.67905, `000000003_lnorm_lnorm` = 1.80859, `000000004_lnorm_lnorm` = 0.921876, `000000005_lnorm_lnorm` = 1.68554, `000000006_lnorm_lnorm` = 1.28523, `000000007_lnorm_lnorm` = 1.82566, `000000008_lnorm_lnorm` = 1.05662, `000000009_lnorm_lnorm` = 1.20995, `000000010_lnorm_lnorm` = 1.67571)"
 weibull,0.05,1.08673,0.895134,0.819219,3.40858,0.357472,parametric,10,1,"c(`000000001_weibull` = 1.67077, `000000002_weibull` = 0.93999, `000000003_weibull` = 1.45323, `000000004_weibull` = 3.60435, `000000005_weibull` = 1.0464, `000000006_weibull` = 1.48364, `000000007_weibull` = 2.08463, `000000008_weibull` = 1.05416, `000000009_weibull` = 2.73428, `000000010_weibull` = 0.784157)"

tests/testthat/_snaps/hc/hc_save_to1_not_multi_default.csv

@@ -1,2 +1,2 @@
 dist,proportion,est,se,lcl,ucl,wt,method,nboot,pboot,samples
-average,0.05,1.24152,NA,1.62613,1.62613,1,parametric,1,1,"c(`000000001_gamma` = 0.987785, `000000001_lgumbel` = 2.2833, `000000001_llogis` = 0.751505, `000000001_lnorm` = 3.09184, `000000001_lnorm_lnorm` = 1.64162, `000000001_weibull` = 1.67077)"
+average,0.05,1.24152,NA,1.62613,1.62613,1,parametric,1,1,"c(`000000001_gamma` = 0.987785, `000000001_lgumbel` = 2.2833, `000000001_llogis` = 0.751505, `000000001_lnorm` = 3.09184, `000000001_lnorm_lnorm` = 1.64164, `000000001_weibull` = 1.67077)"

tests/testthat/_snaps/hc/hc_save_to_not_multi_default.csv

@@ -1,2 +1,2 @@
 dist,proportion,est,se,lcl,ucl,wt,method,nboot,pboot,samples
-average,0.05,1.24152,NA,1.62613,1.62613,1,parametric,1,1,"c(`000000001_gamma` = 0.987785, `000000001_lgumbel` = 2.2833, `000000001_llogis` = 0.751505, `000000001_lnorm` = 3.09184, `000000001_lnorm_lnorm` = 1.64162, `000000001_weibull` = 1.67077)"
+average,0.05,1.24152,NA,1.62613,1.62613,1,parametric,1,1,"c(`000000001_gamma` = 0.987785, `000000001_lgumbel` = 2.2833, `000000001_llogis` = 0.751505, `000000001_lnorm` = 3.09184, `000000001_lnorm_lnorm` = 1.64164, `000000001_weibull` = 1.67077)"

tests/testthat/_snaps/hc/hc_unweighted2.csv

@@ -1,5 +1,5 @@
 dist,proportion,est,se,lcl,ucl,wt,method,nboot,samples
-average,0.05,1.24152,0.886732,0.806144,3.20382,1,parametric,10,"c(`000000001_gamma` = 0.987785, `000000002_gamma` = 3.15112, `000000003_gamma` = 0.84656, `000000004_gamma` = 0.70604, `000000005_gamma` = 1.98672, `000000006_gamma` = 0.873452, `000000007_gamma` = 1.18798, `000000008_gamma` = 0.836688, `000000009_gamma` = 0.719689, `000000010_gamma` = 3.20738, `000000001_lgumbel` = 2.2833, `000000002_lgumbel` = 1.48094, `000000003_lgumbel` = 2.12922, `000000004_lgumbel` = 2.36856, `000000005_lgumbel` = 1.81154, `000000006_lgumbel` = 1.96535, `000000007_lgumbel` = 1.36401, 
+average,0.05,1.24152,0.886735,0.806144,3.20382,1,parametric,10,"c(`000000001_gamma` = 0.987785, `000000002_gamma` = 3.15112, `000000003_gamma` = 0.84656, `000000004_gamma` = 0.70604, `000000005_gamma` = 1.98672, `000000006_gamma` = 0.873452, `000000007_gamma` = 1.18798, `000000008_gamma` = 0.836688, `000000009_gamma` = 0.719689, `000000010_gamma` = 3.20738, `000000001_lgumbel` = 2.2833, `000000002_lgumbel` = 1.48094, `000000003_lgumbel` = 2.12922, `000000004_lgumbel` = 2.36856, `000000005_lgumbel` = 1.81154, `000000006_lgumbel` = 1.96535, `000000007_lgumbel` = 1.36401, 
 `000000008_lgumbel` = 1.93797, `000000009_lgumbel` = 2.70632, `000000010_lgumbel` = 2.09232, `000000001_llogis` = 0.751505, `000000002_llogis` = 3.04268, `000000003_llogis` = 2.10953, `000000004_llogis` = 2.22634, `000000005_llogis` = 1.30249, `000000006_llogis` = 2.52802, `000000007_llogis` = 3.46857, `000000008_llogis` = 2.04533, `000000009_llogis` = 1.85618, `000000010_llogis` = 1.19654, `000000001_lnorm` = 3.09184, `000000002_lnorm` = 2.42899, `000000003_lnorm` = 1.325, `000000004_lnorm` = 1.61081, 
-`000000005_lnorm` = 2.6033, `000000006_lnorm` = 0.865976, `000000007_lnorm` = 2.77742, `000000008_lnorm` = 1.19715, `000000009_lnorm` = 2.45546, `000000010_lnorm` = 0.970096, `000000001_lnorm_lnorm` = 1.64162, `000000002_lnorm_lnorm` = 1.67905, `000000003_lnorm_lnorm` = 1.80875, `000000004_lnorm_lnorm` = 0.921876, `000000005_lnorm_lnorm` = 1.68355, `000000006_lnorm_lnorm` = 1.28521, `000000007_lnorm_lnorm` = 1.82571, `000000008_lnorm_lnorm` = 1.05664, `000000009_lnorm_lnorm` = 1.20995, `000000010_lnorm_lnorm` = 1.67582, 
+`000000005_lnorm` = 2.6033, `000000006_lnorm` = 0.865976, `000000007_lnorm` = 2.77742, `000000008_lnorm` = 1.19715, `000000009_lnorm` = 2.45546, `000000010_lnorm` = 0.970096, `000000001_lnorm_lnorm` = 1.64164, `000000002_lnorm_lnorm` = 1.67905, `000000003_lnorm_lnorm` = 1.80859, `000000004_lnorm_lnorm` = 0.921876, `000000005_lnorm_lnorm` = 1.68554, `000000006_lnorm_lnorm` = 1.28523, `000000007_lnorm_lnorm` = 1.82566, `000000008_lnorm_lnorm` = 1.05662, `000000009_lnorm_lnorm` = 1.20995, `000000010_lnorm_lnorm` = 1.67571, 
 `000000001_weibull` = 1.67077, `000000002_weibull` = 0.93999, `000000003_weibull` = 1.45323, `000000004_weibull` = 3.60435, `000000005_weibull` = 1.0464, `000000006_weibull` = 1.48364, `000000007_weibull` = 2.08463, `000000008_weibull` = 1.05416, `000000009_weibull` = 2.73428, `000000010_weibull` = 0.784157)"

tests/testthat/_snaps/hp/hp.csv

@@ -3,5 +3,5 @@ gamma,1,4.67759,2.14332,0.993989,6.55631,0.356574,parametric,10,1,numeric(0)
 lgumbel,1,0.856452,0.732699,0.0636536,2.21195,0.0134466,parametric,10,1,numeric(0)
 llogis,1,2.80047,1.89766,0.460372,6.09523,0.0656452,parametric,10,1,numeric(0)
 lnorm,1,1.95429,2.09645,0.409858,5.86466,0.177236,parametric,10,1,numeric(0)
-lnorm_lnorm,1,1.23428,1.92465,0.251344,5.72089,0.0296268,parametric,10,1,numeric(0)
+lnorm_lnorm,1,1.23428,1.92474,0.251385,5.72094,0.0296268,parametric,10,1,numeric(0)
 weibull,1,4.62301,1.81696,0.79621,6.01478,0.357472,parametric,10,1,numeric(0)