R Under development (unstable) (2024-07-10 r86888 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(fitdistrplus) Loading required package: MASS Loading required package: survival > nbboot <- 1001 > nbboot <- 10 > > > # (1) Fit of a normal distribution on acute toxicity log-transformed values of endosulfan for > # nonarthropod invertebrates, using maximum likelihood estimation > # to estimate what is called a species sensitivity distribution > # (SSD) in ecotoxicology, followed by estimation of the 5, 10 and 20 percent quantile values of > # the fitted distribution, which are called the 5, 10, 20 percent hazardous concentrations (HC5, HC10, HC20) > # in ecotoxicology, followed with calculations of their confidence intervals with various definitions. > # > data(endosulfan) > ATV <- subset(endosulfan, group == "NonArthroInvert")$ATV > log10ATV <- log10(subset(endosulfan, group == "NonArthroInvert")$ATV) > fln <- fitdist(log10ATV, "norm") > quantile(fln, probs = c(0.05, 0.1, 0.2)) Estimated quantiles for each specified probability (non-censored data) p=0.05 p=0.1 p=0.2 estimate 1.744227 2.080093 2.4868 > bln <- bootdist(fln, niter=nbboot, bootmethod="param") > quantile(bln, probs = c(0.05, 0.1, 0.2)) (original) estimated quantiles for each specified probability (non-censored data) p=0.05 p=0.1 p=0.2 estimate 1.744227 2.080093 2.4868 Median of bootstrap estimates p=0.05 p=0.1 p=0.2 estimate 1.821602 2.136697 2.513664 two-sided 95 % CI of each quantile p=0.05 p=0.1 p=0.2 2.5 % 1.180667 1.591642 2.075168 97.5 % 2.398608 2.653510 2.975204 > quantile(bln, probs = c(0.05, 0.1, 0.2), CI.type = "greater") (original) estimated quantiles for each specified probability (non-censored data) p=0.05 p=0.1 p=0.2 estimate 1.744227 2.080093 2.4868 Median of bootstrap estimates p=0.05 p=0.1 p=0.2 estimate 1.821602 2.136697 2.513664 left bound of one-sided 95 % CI of each quantile p=0.05 p=0.1 p=0.2 5 % 1.218419 1.609704 2.080949 > quantile(bln, probs = c(0.05, 0.1, 0.2), CI.level = 0.9) (original) estimated quantiles for each specified probability (non-censored data) p=0.05 p=0.1 p=0.2 estimate 1.744227 2.080093 2.4868 Median of bootstrap estimates p=0.05 p=0.1 p=0.2 estimate 1.821602 2.136697 2.513664 two-sided 90 % CI of each quantile p=0.05 p=0.1 p=0.2 5 % 1.218419 1.609704 2.080949 95 % 2.363704 2.636803 2.969874 > > # (2) Fit of a distribution on acute salinity log-transformed tolerance > # for riverine macro-invertebrates, using maximum likelihood estimation > # to estimate what is called a species sensitivity distribution > # (SSD) in ecotoxicology, followed by estimation of the 5, 10 and 20 percent quantile values of > # the fitted distribution, which are called the 5, 10, 20 percent hazardous concentrations (HC5, HC10, HC20) > # in ecotoxicology, followed with calculations of their confidence intervals with various definitions. > # > data(salinity) > log10LC50 <-log10(salinity) > flncens <- fitdistcens(log10LC50,"norm") > quantile(flncens, probs = c(0.05, 0.1, 0.2)) Estimated quantiles for each specified probability (censored data) p=0.05 p=0.1 p=0.2 estimate 1.11584 1.194121 1.288913 > blncens <- bootdistcens(flncens,niter=nbboot) > quantile(blncens, probs = c(0.05, 0.1, 0.2)) (original) estimated quantiles for each specified probability (censored data) p=0.05 p=0.1 p=0.2 estimate 1.11584 1.194121 1.288913 Median of bootstrap estimates p=0.05 p=0.1 p=0.2 estimate 1.120514 1.197321 1.289855 two-sided 95 % CI of each quantile p=0.05 p=0.1 p=0.2 2.5 % 1.086488 1.167199 1.261777 97.5 % 1.135722 1.214900 1.316927 > quantile(blncens, probs = c(0.05, 0.1, 0.2), CI.type = "greater") (original) estimated quantiles for each specified probability (censored data) p=0.05 p=0.1 p=0.2 estimate 1.11584 1.194121 1.288913 Median of bootstrap estimates p=0.05 p=0.1 p=0.2 estimate 1.120514 1.197321 1.289855 left bound of one-sided 95 % CI of each quantile p=0.05 p=0.1 p=0.2 5 % 1.087555 1.167574 1.263897 > quantile(blncens, probs = c(0.05, 0.1, 0.2), CI.level = 0.9) (original) estimated quantiles for each specified probability (censored data) p=0.05 p=0.1 p=0.2 estimate 1.11584 1.194121 1.288913 Median of bootstrap estimates p=0.05 p=0.1 p=0.2 estimate 1.120514 1.197321 1.289855 two-sided 90 % CI of each quantile p=0.05 p=0.1 p=0.2 5 % 1.087555 1.167574 1.263897 95 % 1.134726 1.213400 1.313030 > > > # (3) Estimation of quantiles of the fitted distribution (fln) > # and two-sided 95 percent confidence intervals for various > # probabilities using non-parametric bootstrap with 101 iterations > # > bln.np <- bootdist(fln, bootmethod = "nonparam", niter = nbboot) > quantile(bln.np, probs = c(0.05, 0.1, 0.2)) (original) estimated quantiles for each specified probability (non-censored data) p=0.05 p=0.1 p=0.2 estimate 1.744227 2.080093 2.4868 Median of bootstrap estimates p=0.05 p=0.1 p=0.2 estimate 1.734788 2.079622 2.497188 two-sided 95 % CI of each quantile p=0.05 p=0.1 p=0.2 2.5 % 1.189562 1.524928 1.931030 97.5 % 2.464358 2.720569 3.046818 > > # (4) Fit of a loglogistic distribution on the same acute toxicity values and > # estimation of the 5 percent quantile (HC5) of the fitted distribution > # and associated two-sided 95 percent confidence interval > # > fll <- fitdist(log10ATV, "logis") > bll <- bootdist(fll, bootmethod = "param", niter = nbboot) > # in log10(ATV) > HC5ll <- quantile(bll, probs = 0.05) > HC5ll (original) estimated quantiles for each specified probability (non-censored data) p=0.05 estimate 1.780291 Median of bootstrap estimates p=0.05 estimate 1.87385 two-sided 95 % CI of each quantile p=0.05 2.5 % 1.431954 97.5 % 2.223686 > # in ATV > 10^(HC5ll$basequant) numeric(0) > 10^(HC5ll$quantCI) p=0.05 2.5 % 27.0367 97.5 % 167.3732 > > proc.time() user system elapsed 1.81 0.25 2.01