R Under development (unstable) (2025-01-06 r87534 ucrt) -- "Unsuffered Consequences" Copyright (C) 2025 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. > require("fitdistrplus") Loading required package: fitdistrplus Loading required package: MASS Loading required package: survival > > nbboot <- 201 > nbboot <- 10 > ggplotEx <- requireNamespace("ggplot2", quietly = TRUE) > > # (1) Fit of a gamma distribution > # > > set.seed(123) > s1 <- rgamma(20, 3, 2) > f1 <- fitdist(s1, "gamma") > b1 <- bootdist(f1, niter=nbboot, silent=TRUE) > > plot(b1) > quantile(b1) (original) estimated quantiles for each specified probability (non-censored data) p=0.1 p=0.2 p=0.3 p=0.4 p=0.5 p=0.6 p=0.7 estimate 0.4741503 0.6728387 0.848578 1.022138 1.205004 1.408756 1.650833 p=0.8 p=0.9 estimate 1.967001 2.466538 Median of bootstrap estimates p=0.1 p=0.2 p=0.3 p=0.4 p=0.5 p=0.6 p=0.7 estimate 0.5082066 0.6866054 0.8522081 1.018619 1.200685 1.384233 1.577598 p=0.8 p=0.9 estimate 1.7978 2.213969 two-sided 95 % CI of each quantile p=0.1 p=0.2 p=0.3 p=0.4 p=0.5 p=0.6 p=0.7 2.5 % 0.3930161 0.5545033 0.6866142 0.8156814 0.950533 1.099723 1.275833 97.5 % 0.7636052 0.9342519 1.0730947 1.2027269 1.333378 1.479409 1.762129 p=0.8 p=0.9 2.5 % 1.498317 1.827895 97.5 % 2.151868 2.778431 > > par(mfrow=c(1,2)) > CIcdfplot(b1, CI.level=95/100, CI.output = "probability", CI.fill="grey80", CI.col="black") > CIcdfplot(b1, CI.level=95/100, CI.output = "quantile", datacol="blue") > if(ggplotEx) CIcdfplot(b1, CI.level=95/100, CI.output = "probability", CI.fill="grey80", CI.col="black", plotstyle = "ggplot") > if(ggplotEx) CIcdfplot(b1, CI.level=95/100, CI.output = "quantile", datacol="blue", plotstyle = "ggplot") > > par(mfrow=c(1,2)) > CIcdfplot(b1, CI.level=85/100, CI.output = "probability") > CIcdfplot(b1, CI.level=85/100, CI.output = "quantile") > if(ggplotEx) CIcdfplot(b1, CI.level=85/100, CI.output = "probability", plotstyle = "ggplot") > if(ggplotEx) CIcdfplot(b1, CI.level=85/100, CI.output = "quantile", plotstyle = "ggplot") > > par(mfrow=c(1,2)) > CIcdfplot(b1, CI.level=90/100, CI.output = "probability") > CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="black", CI.type = "less", + CI.fill="grey85", verticals=TRUE, datacol="blue", do.points=FALSE) > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "probability", plotstyle = "ggplot") > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="black", CI.type = "less", + CI.fill="grey85", verticals=TRUE, datacol="blue", do.points=FALSE, plotstyle = "ggplot") > > par(mfrow=c(1,2)) > CIcdfplot(b1, CI.level=90/100, CI.output = "probability", CI.type = "greater") > CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="black", CI.type = "greater", + CI.fill="grey90", datacol="blue", datapch=21) > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "probability", CI.type = "greater", plotstyle = "ggplot") > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="black", CI.type = "greater", + CI.fill="grey90", datacol="blue", datapch=21, plotstyle = "ggplot") > > par(mfrow=c(1,1)) > CIcdfplot(b1, CI.level=90/100, CI.output = "probability", CI.col="black", CI.type = "less", CI.fill="grey90") > CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="black", CI.type = "less", CI.fill="grey90", + verticals=TRUE, datacol="blue", do.points=FALSE) > CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="grey90", CI.type = "less", CI.fill="grey90", + verticals=TRUE, datacol="blue", do.points=FALSE, CI.only=TRUE) > CIcdfplot(b1, CI.level=90/100, CI.output = "probability", CI.col="grey85", CI.type = "less", CI.fill="grey90", + CI.only = TRUE) > CIcdfplot(b1, CI.output = "probability", fitlty=3, fitlwd=4) > > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "probability", CI.col="black", CI.type = "less", CI.fill="grey90", plotstyle = "ggplot") > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="black", CI.type = "less", CI.fill="grey90", + verticals=TRUE, datacol="blue", do.points=FALSE, plotstyle = "ggplot") > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="grey90", CI.type = "less", CI.fill="grey90", + verticals=TRUE, datacol="blue", do.points=FALSE, CI.only=TRUE, plotstyle = "ggplot") > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "probability", CI.col="grey85", CI.type = "less", CI.fill="grey90", + CI.only = TRUE, plotstyle = "ggplot") > if(ggplotEx) CIcdfplot(b1, CI.output = "probability", fitlty=3, fitlwd=4, plotstyle = "ggplot") > > # (2) an example from ecotoxicology > # with censored data > # > data(salinity) > log10LC50 <-log10(salinity) > fln <- fitdistcens(log10LC50,"norm") > bln <- bootdistcens(fln, niter=nbboot) > (HC5ln <- quantile(bln,probs = 0.05)) (original) estimated quantiles for each specified probability (censored data) p=0.05 estimate 1.11584 Median of bootstrap estimates p=0.05 estimate 1.110188 two-sided 95 % CI of each quantile p=0.05 2.5 % 1.032713 97.5 % 1.153806 > CIcdfplot(bln, CI.output = "quantile", CI.fill = "lightblue", CI.col = "blue", + xlab = "log10(LC50)",xlim=c(0.5,2),lines01 = TRUE) > if(ggplotEx) CIcdfplot(bln, CI.output = "quantile", CI.fill = "lightblue", CI.col = "blue", + xlab = "log10(LC50)",xlim=c(0.5,2),lines01 = TRUE, plotstyle = "ggplot") > > # zoom around the HC5 with CI on quantiles > CIcdfplot(bln, CI.output = "quantile", CI.fill = "lightblue", CI.col = "blue", + xlab = "log10(LC50)", lines01 = TRUE, xlim = c(0.8, 1.5), ylim = c(0, 0.1)) > abline(h = 0.05, lty = 1) > if(ggplotEx) CIcdfplot(bln, CI.output = "quantile", CI.fill = "lightblue", CI.col = "blue", + xlab = "log10(LC50)", lines01 = TRUE, xlim = c(0.8, 1.5), ylim = c(0, 0.1), plotstyle = "ggplot") + + ggplot2::geom_hline(yintercept = 0.05) > > # zoom around the HC5 with CI on probabilities > CIcdfplot(bln, CI.output = "probability", CI.fill = "lightblue", CI.col = "blue", + xlab = "log10(LC50)", lines01 = TRUE, xlim = c(0.8, 1.5), ylim = c(0, 0.1)) > abline(h = 0.05, lty = 1) > if(ggplotEx) CIcdfplot(bln, CI.output = "probability", CI.fill = "lightblue", CI.col = "blue", + xlab = "log10(LC50)", lines01 = TRUE, xlim = c(0.8, 1.5), ylim = c(0, 0.1), plotstyle = "ggplot") + + ggplot2::geom_hline(yintercept = 0.05) > > > # (3) An example where the difference between "probability" > # and "quantile" is clear on the plot > # > > set.seed(123) > s3 <- rgamma(5, 3, 10) > f3 <- fitdist(s3, "norm") > b3 <- bootdist(f3, niter=nbboot, silent=TRUE) > > par(mfrow=c(1,2)) > CIcdfplot(b3, CI.level=90/100, CI.output = "probability") > CIcdfplot(b3, CI.level=90/100, CI.output = "quantile") > > if(ggplotEx) CIcdfplot(b3, CI.level=90/100, CI.output = "probability", plotstyle = "ggplot") > if(ggplotEx) CIcdfplot(b3, CI.level=90/100, CI.output = "quantile", plotstyle = "ggplot") > > #some ideas from http://edild.github.io/ssd/ > > proc.time() user system elapsed 8.20 0.34 8.48