Problems with CIs for ssd_plot with censored data

[beckyfisher]

library(ssdtools)
library(tidyverse)
#> Warning: package 'ggplot2' was built under R version 4.3.3
example_dat <- ssddata::ccme_boron |>
dplyr::mutate(left=Conc, right=Conc)

left_censored_example <- example_dat
#left_censored_example$left[c(3,6,8)] <- NA

left_censored_dists <- ssd_fit_dists(left_censored_example,
dists = ssd_dists_bcanz(n = 2),
left = "left", right = "right")
ssd_hc(left_censored_dists, average = FALSE)
#> # A tibble: 5 × 11
#> dist proportion est se lcl ucl wt method nboot pboot samples
#> <I
#> 1 gamma 0.05 1.07 NA NA NA 0.367 paramet… 0 NA
#> 2 lgumbel 0.05 1.77 NA NA NA 0.0139 paramet… 0 NA
#> 3 llogis 0.05 1.56 NA NA NA 0.0676 paramet… 0 NA
#> 4 lnorm 0.05 1.68 NA NA NA 0.183 paramet… 0 NA
#> 5 weibull 0.05 1.09 NA NA NA 0.368 paramet… 0 NA
ssd_hc(left_censored_dists)
#> # A tibble: 1 × 11
#> dist proportion est se lcl ucl wt method nboot pboot samples
#> <I
#> 1 average 0.05 1.24 NA NA NA 1 parametr… 0 NA
ssd_gof(left_censored_dists)
#> # A tibble: 5 × 9
#> dist ad ks cvm aic aicc bic delta weight
#>
#> 1 gamma 0.440 0.117 0.0554 238. 238. 240. 0.005 0.367
#> 2 lgumbel 0.829 0.158 0.134 244. 245. 247. 6.56 0.014
#> 3 llogis 0.487 0.0994 0.0595 241. 241. 244. 3.39 0.068
#> 4 lnorm 0.507 0.107 0.0703 239. 240. 242. 1.40 0.183
#> 5 weibull 0.434 0.117 0.0542 238. 238. 240. 0 0.368

set.seed(99)
left_censored_pred <- predict(left_censored_dists, ci = TRUE, nboot = 100)

ssd_plot(left_censored_example, left_censored_pred,
left = "left", right = "right",
xlab = "Concentration (mg/L)"
)

Created on 2024-04-11 with reprex v2.1.0

[joethorley]

Thanks for reporting.

With resultant plot without censoring.

library(ssdtools)
library(tidyverse)
example_dat <- ssddata::ccme_boron |>
  dplyr::mutate(left=Conc, right=Conc)

left_censored_example <- example_dat
#left_censored_example$left[c(3,6,8)] <- NA

left_censored_dists <- ssd_fit_dists(left_censored_example,
                                     dists = ssd_dists_bcanz(n = 2),
                                     left = "left", right = "right")
ssd_hc(left_censored_dists, average = FALSE)
#> # A tibble: 5 × 11
#>   dist    proportion   est    se   lcl   ucl     wt method   nboot pboot samples
#>   <chr>        <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <chr>    <int> <dbl> <I<lis>
#> 1 gamma         0.05  1.07    NA    NA    NA 0.367  paramet…     0    NA <dbl>  
#> 2 lgumbel       0.05  1.77    NA    NA    NA 0.0139 paramet…     0    NA <dbl>  
#> 3 llogis        0.05  1.56    NA    NA    NA 0.0676 paramet…     0    NA <dbl>  
#> 4 lnorm         0.05  1.68    NA    NA    NA 0.183  paramet…     0    NA <dbl>  
#> 5 weibull       0.05  1.09    NA    NA    NA 0.368  paramet…     0    NA <dbl>

ssd_hc(left_censored_dists)
#> # A tibble: 1 × 11
#>   dist    proportion   est    se   lcl   ucl    wt method    nboot pboot samples
#>   <chr>        <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>     <int> <dbl> <I<lis>
#> 1 average       0.05  1.24    NA    NA    NA     1 parametr…     0   NaN <dbl>

ssd_gof(left_censored_dists)
#> # A tibble: 5 × 9
#>   dist       ad     ks    cvm   aic  aicc   bic delta weight
#>   <chr>   <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
#> 1 gamma   0.440 0.117  0.0554  238.  238.  240. 0.005  0.367
#> 2 lgumbel 0.829 0.158  0.134   244.  245.  247. 6.56   0.014
#> 3 llogis  0.487 0.0994 0.0595  241.  241.  244. 3.39   0.068
#> 4 lnorm   0.507 0.107  0.0703  239.  240.  242. 1.40   0.183
#> 5 weibull 0.434 0.117  0.0542  238.  238.  240. 0      0.368

set.seed(99)
left_censored_pred <- predict(left_censored_dists, ci = TRUE, nboot = 100)

ssd_plot(left_censored_example, left_censored_pred,
         left = "left", right = "right",
         xlab = "Concentration (mg/L)"
)

^{Created on 2024-04-11 with reprex v2.1.0}

[joethorley]

with censoring

library(ssdtools)
library(tidyverse)
example_dat <- ssddata::ccme_boron |>
  dplyr::mutate(left=Conc, right=Conc)

left_censored_example <- example_dat
left_censored_example$left[c(3,6,8)] <- NA

left_censored_dists <- ssd_fit_dists(left_censored_example,
                                     dists = ssd_dists_bcanz(n = 2),
                                     left = "left", right = "right")
ssd_hc(left_censored_dists, average = FALSE)
#> # A tibble: 5 × 11
#>   dist    proportion   est    se   lcl   ucl     wt method   nboot pboot samples
#>   <chr>        <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <chr>    <int> <dbl> <I<lis>
#> 1 gamma         0.05 0.674    NA    NA    NA 0.376  paramet…     0    NA <dbl>  
#> 2 lgumbel       0.05 1.51     NA    NA    NA 0.0221 paramet…     0    NA <dbl>  
#> 3 llogis        0.05 1.15     NA    NA    NA 0.0590 paramet…     0    NA <dbl>  
#> 4 lnorm         0.05 1.32     NA    NA    NA 0.176  paramet…     0    NA <dbl>  
#> 5 weibull       0.05 0.752    NA    NA    NA 0.367  paramet…     0    NA <dbl>

ssd_hc(left_censored_dists)
#> # A tibble: 1 × 11
#>   dist    proportion   est    se   lcl   ucl    wt method    nboot pboot samples
#>   <chr>        <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>     <int> <dbl> <I<lis>
#> 1 average       0.05 0.859    NA    NA    NA     1 parametr…     0   NaN <dbl>

ssd_gof(left_censored_dists)
#> # A tibble: 5 × 9
#>   dist       ad    ks   cvm   aic  aicc   bic delta weight
#>   <chr>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
#> 1 gamma      NA    NA    NA  222.    NA    NA 0      0.376
#> 2 lgumbel    NA    NA    NA  228.    NA    NA 5.67   0.022
#> 3 llogis     NA    NA    NA  226.    NA    NA 3.70   0.059
#> 4 lnorm      NA    NA    NA  224.    NA    NA 1.52   0.176
#> 5 weibull    NA    NA    NA  222.    NA    NA 0.046  0.367

set.seed(99)
left_censored_pred <- predict(left_censored_dists, ci = TRUE, nboot = 100)

ssd_plot(left_censored_example, left_censored_pred,
         left = "left", right = "right",
         xlab = "Concentration (mg/L)"
)

^{Created on 2024-04-11 with reprex v2.1.0}