Keep getting this error when trying to plot:

"Error in is.finite(x) :
default method not implemented for type 'language'"

This happens when using the script straight from the Shiny App using Boron as the example as well as when I use my own data.
Any suggestions?

I tried to do again all the command to obtain HC5 and I realized I can’t have same number that showed by Joe and team.

> ssd_hc5(boron_dists)
    prop   est    se   lcl   ucl
   <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.0500  1.25 0.726 0.623  3.19

Hi @joethorley. Interesting package.
Please, change line 42 of ssdtools.R to ggplot2::autoplot or add library(ggplot2)

line for each distribution

to reduce runtime

as tbl data frame

hi there, when i try to fit my data set using the R code, it works but i get errors when i attempt to plot, or make predictions. I've attached the data set. Cu_lowdata.xlsx

This is the code and errors:

Cu_low_dists <- ssd_fit_dists(LowerCuMarData)
autoplot(Cu_low_dists)
Error in pgompertz(q = c(0.7, 1.3, 1.94), shape = 0.391386784692936, scale = 0.313425202466778) :
unused argument (q = c(0.7, 1.3, 1.94))
ssd_gof(Cu_low_dists)
Error in pgompertz(q = c(0.7, 1.3, 1.94), shape = 0.391386784692936, scale = 0.313425202466778) :
unused argument (q = c(0.7, 1.3, 1.94))

It works when i do it in the R shiny app.

Is there a way to exclude some fits when using the R code, as you can in the Shiny app?
Thanks,
Jenni

QQ, PP, residuals etc

Gamma distribution with shape (a) and scale (s) has
"mean and variance are E(X) = as and Var(X) = as^2" according to the help file for dgamma.

Solving for $a$ and $s$ gives
$ a = V(x) / E[x]$
so using method of moments, use a starting value for $a = var(x) / mean(x)$.

Solving for $a$ then gives
$ a = E[X] / s$
and a moment estimator is then $a = mean(x) / s = mean(x) / (V(x)/Mean(x) = mean^2/V(x)$

We have
if (dist$dist == "gamma") {
dist$start <- list(
scale = stats::var(x) / mean(x),
shape = mean(x)^2 / stats::var(x)^2
)
so it appear we have and extra square in the last line. This should read

if (dist$dist == "gamma") {
dist$start <- list(
scale = stats::var(x) / mean(x),
shape = mean(x)^2 / stats::var(x)
)

i.e. remove the ^2 from the var(x) in the initial value for shape.

Fortunately, the fitdist() function appears to be robust to making the shape parameter too small so no obvious user impact.

in tbl data frame

Namely the HC5s appear to be incorrectly assigned to distributions in Table 1

For example

> ssd_hc5(boron_lnorm)
# A tibble: 1 x 5
   prob      est        se       lcl      ucl
  <dbl>    <dbl>     <dbl>     <dbl>    <dbl>
1  0.05 1.681175 0.6848542 0.8618846 3.574394

the ssdca package with lnorm produces an hc5 estimate of 1.68

and below

> ssd_hc5(boron_dists, average = F)
# A tibble: 6 x 7
      dist  prop      est        se       lcl      ucl weight
     <chr> <dbl>    <dbl>     <dbl>     <dbl>    <dbl>  <dbl>
1    lnorm  0.05 1.681175 0.6874052 0.8851202 3.611764  0.133
2     llog  0.05 1.562052 0.7760264 0.6636144 3.652784  0.049
3 gompertz  0.05 1.299302 0.4419386 0.9581149 2.788982  0.271
4  lgumbel  0.05 1.768918 0.5369794 1.1299498 3.231247  0.010
5    gamma  0.05 1.075854 0.8355977 0.2652210 3.476163  0.268
6  weibull  0.05 1.087170 0.7784698 0.3654033 3.292916  0.269

but in Table 1 1.68 is the hc5 value for log-logistic

note aicc values in Table 1 appear to be correctly assigned

> ssd_gof(boron_dists)
# A tibble: 6 x 7
      dist        ad         ks        cvm      aic     aicc      bic
     <chr>     <dbl>      <dbl>      <dbl>    <dbl>    <dbl>    <dbl>
1    lnorm 0.5070335 0.10651430 0.07033164 239.0284 239.5084 241.6928
2     llog 0.4870694 0.09934088 0.05950233 241.0149 241.4949 243.6793
3 gompertz 0.6019563 0.12018745 0.08221977 237.6112 238.0912 240.2756
4  lgumbel 0.8286390 0.15826649 0.13401800 244.1860 244.6660 246.8504
5    gamma 0.4409319 0.11691481 0.05550569 237.6303 238.1103 240.2947
6  weibull 0.4346273 0.11697580 0.05426727 237.6253 238.1053 240.2897

For some extreme sets of data, the gamma distribution reports NA for the standard error. Eg.]

library(ssdtools)
test.csv <- textConnection(
"dummy , Conc , Units
all , 0.1 , ug/L
all , 0.12 , ug/L
all , 0.24 , ug/L
all , 0.42 , ug/L
all , 0.67 , ug/L
all , 0.78 , ug/L
all , 120 , ug/L
all , 2030 , ug/L
all , 9033 , ug/L
all , 15000 , ug/L
all , 15779 , ug/L
all , 20000 , ug/L
all , 31000 , ug/L
all , 40000 , ug/L
all , 105650 , ug/L")
test <- read.csv(test.csv, header=TRUE, strip.white=TRUE, as.is=TRUE)

ssd_fit_dists(test, dists="gamma")

reports

Fitting of the distribution ' gamma ' by maximum likelihood
Parameters:
estimate Std. Error
scale 9.695693e+04 NaN
shape 1.642224e-01 0.04142257

Estimates of scale and shape are ok, but it was unable to compute a standard error for the scale parameter. This appears to be a scaling problem in inverting the hessian. The standard error for the scale parameter is about 67000 and the se for the shape is .04. So the vcov will have terms on the order of 67000^2 and .04^2 which vary by 2.8 x 10^12 (!) so I suspect that it is a numerical issue

If you rescale the data

factor=100
test2 <- test
test2$Conc <- test2$Conc/ factor
ssd_fit_dists(test2, dists="gamma")

all is well
Fitting of the distribution ' gamma ' by maximum likelihood
Parameters:
estimate Std. Error
scale 969.5359604 674.2705097
shape 0.1641717 0.0454135

The pdf for the gamma has a term x/s so if you multiple x by a factor k and multiply s by a factor of k, the pdf is unchanged except for the change in the scale factor.

I would advise user if any estimates have NA for the standard errors and advise them that a rescaling may be appropriate. We could make a specific message for each distribution, i.e. typically advise gamma users to divide by a factor k, or a generic message and the user needs to try both multiplying or dividing the data by a factor until message "disappears".

Estimates are ok, so estimates of HCx are ok. Bootstrapping used to compute SE for HCx so this should also work fine.

Carl

log-normal and gamma

for narrative around different distributions

Check Apache 2.0 license for all code and creative commons license for non-code pages

https://github.com/bcgov/BC-Policy-Framework-For-GitHub/blob/master/BC-Open-Source-Development-Employee-Guide/Licenses.md

From CRAN:

FIr the future: Please elaborate in the Description field what this is about. Is there some reference about the method you can add in the Description field in the form Authors (year) doi:.....?

Can you please remove the old background paper and update it with this document:

http://a100.gov.bc.ca/appsdata/acat/documents/r57400/2_1568399094009_8398900200.pdf

Thanks!

as tbl data frame

upgrade from exploration

consider

  data <- ssdtools::boron_data
  data$Conc <- data$Conc / 100
  ssdtools:::ssd_fit_dist(data, dist = "llog")

AIC with weighted likelihoods??

i.e. replace 0. mykiss with Onchorhynchus mykiss etc

based on AIC

add to readme