Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

How to calculate ecd #76

Open
joethorley opened this issue Aug 19, 2024 · 2 comments
Open

How to calculate ecd #76

joethorley opened this issue Aug 19, 2024 · 2 comments

Comments

@joethorley
Copy link
Member

joethorley commented Aug 19, 2024

In small sample size bias vignette it is as follows:

Ag$ecdf <- (rank(Ag$Conc)+.25)/(nrow(Ag)+.5)

which gives

# A tibble: 9 × 6
  Chemical Species                Conc Group        Units  ecdf
  <chr>    <chr>                 <dbl> <fct>        <chr> <dbl>
1 Silver   Oncorhynchus mykiss    0.24 Fish         ug/L  0.132
2 Silver   Lemna gibba            0.63 Plant        ug/L  0.237
3 Silver   Ceriodaphnia dubia     0.78 Invertebrate ug/L  0.342
4 Silver   Pimephales promelas    0.83 Fish         ug/L  0.447
5 Silver   Ictalurus punctatus    1.9  Fish         ug/L  0.553
6 Silver   Daphnia magna          2.12 Invertebrate ug/L  0.658
7 Silver   Hyalella azteca        4    Invertebrate ug/L  0.763
8 Silver   Chironomus tentans    13    Invertebrate ug/L  0.868
9 Silver   Micropterus salmoides 23    Fish         ug/L  0.974
[1] 0.1315789 0.2368421 0.3421053 0.4473684 0.5526316 0.6578947 0.7631579 0.8684211 0.9736842

but in the functions it is

function(x, ties.method = "first") {
  (rank(x, ties.method = ties.method) - 0.5) / length(x)
}

which gives

[1] 0.05555556 0.16666667 0.27777778 0.38888889 0.50000000 0.61111111 0.72222222 0.83333333 0.94444444
@joethorley
Copy link
Member Author

joethorley commented Oct 24, 2024

The following reprex developed with @beckyfisher demonstrates that the vignette code is suboptimal but ppoints may have better properties than the package function.

ecdf_function <- ssdtools::ssd_ecd

print(ssdtools::ssd_ecd)
#> function (x, ties.method = "first") 
#> {
#>     (rank(x, ties.method = ties.method) - 0.5)/length(x)
#> }
#> <bytecode: 0x129c99a88>
#> <environment: namespace:ssdtools>

ecdf_vignette <- function(x, ties.method = "first") {
  (rank(x, ties.method = ties.method) + 0.25) / (length(x) + 0.5)
}

ecdf_function(1)
#> [1] 0.5
ecdf_vignette(1)
#> [1] 0.8333333
stats::ppoints(1)
#> [1] 0.5

ecdf_function(1:2)
#> [1] 0.25 0.75
ecdf_vignette(1:2)
#> [1] 0.5 0.9
stats::ppoints(1:2)
#> [1] 0.2777778 0.7222222

ecdf_function(1:3)
#> [1] 0.1666667 0.5000000 0.8333333
ecdf_vignette(1:3)
#> [1] 0.3571429 0.6428571 0.9285714
stats::ppoints(1:3)
#> [1] 0.1923077 0.5000000 0.8076923

ecdf_function(1:4)
#> [1] 0.125 0.375 0.625 0.875
ecdf_vignette(1:4)
#> [1] 0.2777778 0.5000000 0.7222222 0.9444444
stats::ppoints(1:4)
#> [1] 0.1470588 0.3823529 0.6176471 0.8529412

ecdf_function(c(1,1))
#> [1] 0.25 0.75
ecdf_vignette(c(1,1))
#> [1] 0.5 0.9
stats::ppoints(c(1,1))
#> [1] 0.2777778 0.7222222

Created on 2024-10-23 with reprex v2.1.1

@beckyfisher
Copy link

For completeness, the code for ppoints is:

function (n, a = if (n <= 10) 3/8 else 1/2) 
{
    if (length(n) > 1L) 
        n <- length(n)
    if (n > 0) 
        (1L:n - a)/(n + 1 - 2 * a)
    else numeric()
}

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants