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.777572 2.103596 2.43865 two-sided 95 % CI of each quantile p=0.05 p=0.1 p=0.2 2.5 % 1.472471 1.818170 2.233205 97.5 % 2.370334 2.639452 2.977245 > 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.777572 2.103596 2.43865 left bound of one-sided 95 % CI of each quantile p=0.05 p=0.1 p=0.2 5 % 1.524112 1.86079 2.261322 > 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.777572 2.103596 2.43865 two-sided 90 % CI of each quantile p=0.05 p=0.1 p=0.2 5 % 1.524112 1.860790 2.261322 95 % 2.264801 2.525721 2.865501 > > # (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.139796 1.214227 1.298036 two-sided 95 % CI of each quantile p=0.05 p=0.1 p=0.2 2.5 % 1.057824 1.135910 1.230466 97.5 % 1.179856 1.254932 1.345843 > 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.139796 1.214227 1.298036 left bound of one-sided 95 % CI of each quantile p=0.05 p=0.1 p=0.2 5 % 1.0582 1.136305 1.230884 > 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.139796 1.214227 1.298036 two-sided 90 % CI of each quantile p=0.05 p=0.1 p=0.2 5 % 1.058200 1.136305 1.230884 95 % 1.170991 1.246258 1.337401 > > > # (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.792414 2.133804 2.539207 two-sided 95 % CI of each quantile p=0.05 p=0.1 p=0.2 2.5 % 0.9760929 1.349996 1.797993 97.5 % 2.6927550 2.900403 3.151848 > > # (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.961062 two-sided 95 % CI of each quantile p=0.05 2.5 % 1.195361 97.5 % 2.683243 > # in ATV > 10^(HC5ll$basequant) numeric(0) > 10^(HC5ll$quantCI) p=0.05 2.5 % 15.68055 97.5 % 482.21783 > > proc.time() user system elapsed 1.75 0.25 1.96