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 <- 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.4742872 0.67299 0.8487352 1.022297 1.20516 1.408907 1.650974 p=0.8 p=0.9 estimate 1.967124 2.466624 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.508392 0.6868063 0.8523895 1.0188 1.200886 1.384462 1.577764 p=0.8 p=0.9 estimate 1.797957 2.214205 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.393174 0.5546783 0.6867869 0.8158468 0.9506868 1.099861 1.275947 97.5 % 0.763781 0.9344348 1.0732803 1.2029130 1.3335633 1.479578 1.762352 p=0.8 p=0.9 2.5 % 1.498368 1.827862 97.5 % 2.152084 2.778622 > > 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 12.64 0.43 13.07